MNRAS 527, 8886–8906 (2024) [Link]

1093/mnras/stad3782
Advance Access publication 2023 December 8

α-enhanced astrochemistry: the carbon cycle in extreme galactic

conditions
Thomas G. Bisbas ,1 ‹ Zhi-Yu Zhang ,2,3 ‹ Eda Gjergo ,2,3 Ying-He Zhao ,4,5 Gan Luo ,6,2
Donghui Quan ,1,7 Xue-Jian Jiang ,1 Yichen Sun ,2,3 Theodoros Topkaras,8 Di Li 9,1,10 and Ziyi Guo2,3
1 Research Center for Intelligent Computing Platforms, Zhejiang Lab, Hangzhou 311100, China
2 School of Astronomy and Space Science, Nanjing University, Nanjing 210093, China
3 Key Laboratory of Modern Astronomy and Astrophysics, Nanjing University, Ministry of Education, Nanjing 210093, China
4 Yunnan Observatories, Chinese Academy of Sciences, Kunming 650011, China
5 Key Laboratory of Radio Astronomy and Technology, Chinese Academy of Sciences, A20 Datun Road, Chaoyang District, Beijing 100101, P. R. China
6 Institut de Radioastronomie Millimetrique, 300 rue de la Piscine, Domaine Universitaire de Grenoble, F-38406 Saint-Martin d’Héres, France

Downloaded from [Link] by guest on 21 May 2026

7 Xinjiang Astronomical Observatory, Chinese Academy of Sciences, No. 150 Science 1-Street, Urumqi 830011, P. R. China
8 I. Physikalisches Institut, Universität zu Köln, Zülpicher Straße 77, D-50937 Köln, Germany
9 CAS Key Laboratory of FAST, National Astronomical Observatories, Chinese Academy of Sciences, Beijing 100101, China
10 NAOC-UKZN Computational Astrophysics Centre, University of KwaZulu-Natal, Durban 4000, South Africa

Accepted 2023 December 5. Received 2023 December 5; in original form 2023 September 21

ABSTRACT
Astrochemistry has been widely developed as a power tool to probe the physical properties of the interstellar medium (ISM)
in various conditions of the Milky Way (MW) Galaxy, and in near and distant galaxies. Most current studies conventionally
apply linear scaling to all elemental abundances based on the gas-phase metallicity. However, these elements, including carbon
and oxygen, are enriched differentially by stellar nucleosynthesis and the overall galactic chemical evolution, evident from
α-enhancement in multiple galactic observations such as starbursts, high-redshift star-forming galaxies, and low-metallicity
dwarfs. We perform astrochemical modelling to simulate the impact of an α-enhanced ISM gas cloud on the abundances of
the three phases of carbon (C+ , C, CO) dubbed as ‘the carbon cycle’. The ISM environmental parameters considered include
two cosmic-ray ionization rates (ζ CR = 10−17 and 10−15 s−1 ), two isotropic FUV radiation field strengths (χ /χ 0 = 1 and 102 ),
and (sub-)linear dust-to-gas relations against metallicity, mimicking the ISM conditions of different galaxy types. In galaxies
with [C/O] < 0, CO, C, and C+ , all decrease in both abundances and emission, though with differential biases. The low-J CO
emission is found to be the most stable tracer for the molecular gas, while C and C+ trace H2 gas only under limited conditions,
in line with recent discoveries of [C I]-dark galaxies. We call for caution when using [C II] 158 μm and [C I](1–0) as alternative
H2 -gas tracers for both diffuse and dense gas with non-zero [C/O] ratios.
Key words: radiative transfer – methods: numerical – ISM: abundances – photodissociation region (PDR) – galaxies: ISM.

Seifried et al. 2020; Hu, Sternberg & van Dishoeck 2021). These
1 I N T RO D U C T I O N
help to understand the physical properties and the line emission of
To study the evolution of the multiphase interstellar medium (ISM) the ISM both in Galactic local clouds and distant galaxies, especially
and the star-formation process in galaxies through cosmic time, nowadays since the JWST in combination with the Atacama Large
it is required to understand the astrochemistry that occurs under Millimeter/sub-millimeter Array (ALMA) offer an unprecedented
the different environmental conditions. These conditions act as an view of our Universe.
energy source that powers the ISM, which in turn compensates with Metallicity is one of the most important environmental parameters
a number of cooling functions, leading to its thermal balance. Various of the ISM, as it plays a significant role in controlling its chemistry
computational works model these processes, either in static but (see Maiolino & Mannucci 2019, for a review). Metals are the key
chemically evolved clouds (e.g. Le Petit et al. 2006; Röllig et al. 2007; coolants in the cold phase (Tkin < 104 K) of the ISM, such as the
Bisbas et al. 2012; Ferland et al. 2017; Röllig & Ossenkopf-Okada [C II] 158 μm and [O I] 63 μm fine-structure lines, molecular lines
2022) or by evolving them together with (magneto-)hydrodynamics such as the CO ladder, and dust emission (Tielens 2005; Draine
(e.g. Glover et al. 2010; Walch et al. 2015; Girichidis et al. 2016; 2011; Goldsmith et al. 2012). Denoted as Z, metallicity refers to the
Richings & Schaye 2016; Gong, Ostriker & Wolfire 2017; Bate 2019; mass fraction of elements heavier than hydrogen (X) and helium (Y)
present in a medium, so that X + Y + Z = 1. The galactic metallicity is
often presented by the abundance measurement of a single element;
 E-mail: tbisbas@[Link] (TGB); zzhang@[Link] (ZYZ) either by Oxygen ([O/H], frequently observed for gas), or by Iron

© 2023 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative
Commons Attribution License ([Link] which permits unrestricted reuse, distribution, and reproduction in any medium,
provided the original work is properly cited.
 α-enhanced ISM 8887
([Fe/H], frequently observed for stars). In the MW and also in other Function (IMF), and the evolutionary history of the host galaxy,
galaxies, Z may change as a function of the galactocentric radius with the relative evolutionary patterns of C and O started to be known
the general trend to decrease towards the outer galactic regions (e.g. through observational and theoretical studies. In the MW and in
Wolfire et al. 2003; Wuyts et al. 2016; Kreckel et al. 2019; Esteban nearby galaxies, both C and O abundances as well as their relative
et al. 2022; Matsunaga et al. 2023). ratio C/O have been systematically measured from the surface of
Throughout this work, we will refer to the oxygen and carbon stars of various ages and metallicities, which reflects the conditions
abundances using the square bracket notation (Aller & Greenstein of the ISM throughout history (e.g. Suárez-Andrés et al. 2018).
1960; Pagel 2009; Maiolino & Mannucci 2019), as was also used In particular, throughout the galactic evolution and starting from
above. The notation [A/B] represents the logarithmic ratio of the metallicities as low as [O/H] ≈ −3 (or [Fe/H] ≈ −2) and until
number of atoms of element A to element B, normalized to a reference [O/H] ≈ −1.5, several observations (e.g. Akerman et al. 2004;
abundance, typically the Sun: Fabbian et al. 2009; Berg et al. 2016, 2019) show that the [C/O]
    ratio drops from solar values, reaching a minimum of [C/O] ≈ −0.7
NA NA
[A/B] = log10 − log10 . (1) to −1 (Trainor et al. 2016; Cooke, Pettini & Steidel 2017; Maiolino &
NB NB  Mannucci 2019). This trend (which decreases also as a function of

Downloaded from [Link] by guest on 21 May 2026

For the purposes of the above equation, we adopt solar abun- [Fe/H]) is measured from Carbon-Enhanced Metal-Poor (CEMP)
dances (Asplund et al. 2009; Asplund, Amarsi & Grevesse 2021) stars. CEMP stars are possibly enriched by the first generation
in which 12 + log10 (O/H) = 8.691 and 12 + log10 (C/H) = 8.43. (Population III) stars or by binary stars whose one member is an AGB
Consequently, for [O/H] = 0 and [C/O] = 0, it is understood that resulting in mass transfer (Bonifacio et al. 2015). At [O/H] ≈ −1.5
solar values are assumed. However, throughout the paper, gas-phase (or [Fe/H] > −2), a turnover occurs bringing the [C/O] ratio back
elemental abundances will be considered, which are lower than the to solar values ([C/O] ≈ 0) (Spitoni et al. 2019; Romano et al. 2020;
solar values to account for effects related to depletion on grains (see Delgado Mena et al. 2021). This trend is observed in both MW stars
Section 2.1). and metal-poor gas of local dwarf galaxies (Berg et al. 2019; Romano
2022).
For [C/O] ratios less than the expected ones from solar scaling,
1.1 Chemical evolution of elements in galaxies it means that the relative abundance of oxygen (the most typical α-
Oxygen and carbon are the most abundant metals in the Universe element) is enhanced relative to the expected value. Therefore, the
and they, therefore, constitute the core of chemical evolution of aforementioned conditions can be considered as ‘α-enhanced’.
galaxies. Similar to all other elements with atomic mass A > 12, i.e.
from carbon to iron, C and O are continuously enriched via stellar 1.2 Treatment of elemental abundances and dust in
nucleosynthesis throughout the galactic evolution (see Romano 2022, astrochemical models
for a review), while elements heavier than 56 Fe are mainly produced
via neutron capture (Burbidge et al. 1957) by Asymptotic Giant In astrochemical models concerning small-scale objects, such as
Branch stars (AGB) (Gallino et al. 1998; Cescutti & Matteucci 2022) the circumstellar discs of AGB stars (e.g. Li et al. 2016) and
and core collapse supernovae (Limongi & Chieffi 2018), or by more protoplanetary discs (e.g. Öberg, Murray-Clay & Bergin 2011;
catastrophic events such as binary-neutron star mergers or collapsars van Dishoeck et al. 2023), non-Solar values of the C/O ratio are
(Kasen et al. 2017; Kajino et al. 2019). However, carbon and oxygen considered due to their particular importance in icy chemistry. For
originate from different processes, making their abundances vary example, the Gas phase Elemental abundances in Molecular cloudS
both from galaxy-to-galaxy and also within a galaxy. (‘GEMS’; Fuente et al. 2019) is a dedicated IRAM program to study
The 16 O (O-16) isotope is mostly synthesized from 12 C and α the depletion of key elements, such as S, C, N, and O, in star-forming
particles (i.e. the Helium nucleus), through the reaction of 12 C(α, filaments. However, such systematic studies in molecular clouds and
γ )16 O, which has a low reaction rate in low-temperature hydrostatic on galactic scales, especially considering extragalactic conditions,
He-burning. The major producer of 16 O is the short-lived massive are still lacking.
stars with mass  8 M . Through type-II supernovae explosions, In most astrochemical models of the MW and external galaxies,
oxygen is almost always released back to the ISM within in a short it is frequently assumed that the CNO abundances have a linear cor-
period of time (Romano 2022). relation with metallicity. This assumption means that, for example,
The 12 C (C-12) isotope, on the other hand, is synthesized through the abundances of C and N will decrease by the same amount with
the triple-α capture process during the He-burning stage, in which [O/H], leading to constant [C/O] and [N/O] as a function of [O/H]
three α particles are synthesized to the Hoyle state of 12 C (Hoyle (e.g. Kaufman et al. 1999; Glover & Clark 2016; Gong, Ostriker &
1954; Freer & Fynbo 2014). Apart from massive stars, carbon is also Wolfire 2017; Bisbas, Schruba & van Dishoeck 2019; Gong et al.
and mainly produced in low- and intermediate-mass stars (1 < M < 2020; Bisbas, Tan & Tanaka 2021). It was only recently since the
8 M ), which live much longer than the first ones. Carbon is released impact of negative [C/O] and [N/O] ratios on the larger-scale ISM
into the ISM through type Ia supernovae and through the Asymptotic and star-formation process have been considered in numerical and
Giant Branch (AGB) stars (e.g. Berg et al. 2019). As discussed below, analytical works (see Section 1.3).
this delay in the release compared to oxygen, makes the C/O ratio Similar to the abundance scaling, the dust-to-gas ratio is frequently
to decrease throughout the chemical evolution period, allowing also assumed to decrease linearly with metallicity (see e.g. Bate 2014;
astronomers to use it as a clock. Tanaka & Omukai 2014; Tanaka et al. 2018; Bate 2019; Bisbas, Tan &
Although the abundance evolution of C and O elements involves Tanaka 2021). However, observations (Galametz et al. 2011; Herrera-
various processes, including stellar physics, the stellar Initial Mass Camus et al. 2012; Rémy-Ruyer et al. 2014) show that for low values
of [O/H] (e.g. for the [O/H] = −1.21 modelled here), the dust-
to-gas ratio decreases even further, favouring the photodissociation
1 Themost recent measurements by Pietrow et al. (2023) set the solar oxygen processes of both H2 and CO molecules. This behaviour is important
abundance to 12 + log10 (O/H) = 8.73 ± 0.03. for the C/O ratio all the more since the carbon abundance and the

MNRAS 527, 8886–8906 (2024)

8888 T. G. Bisbas et al.
dust-to-gas relation are related (Mathis 1990; Dwek 1998). The Other environmental parameters, which will not be considered here,
assumption of a linear correlation leads to an overestimation of include X-rays (Maloney, Hollenbach & Tielens 1996; Meijerink,
the dust opacity having consequences in the formation of H2 on Spaans & Israel 2006; Mackey et al. 2019), turbulence (Xie, Allen &
grains. Apart from the fact that the grain abundance controls the Langer 1995), and shocks (Meijerink et al. 2011; Kelly et al. 2017;
FUV attenuation along a column, it is also directly proportional to Cosentino et al. 2019; James et al. 2020).
the rate of photoelectric heating and, thus, important for the thermal Ionized carbon (‘C II’ when referring to its emission or ‘C+ ’ when
balance (Wolfire et al. 2008). referring to its abundance) is an ion with significant importance to
the study of the ISM. Its fine-structure transition, at 157.7 μm, is
one of the brightest cooling lines and is commonly used to trace
1.3 Impact of a non-linear [C/O] with [O/H] in the ISM of
warm neutral gas and ionized gas. [C II] is excited due to inelastic
galaxies
collisions mainly with e− , H I, and H2 (Goldsmith et al. 2012; Lique
The impact of a negative [C/O] in the ISM has been recently et al. 2013), and can be enhanced by interstellar shocks (Draine &
considered in observations and numerical simulations of high- McKee 1993; Appleton et al. 2013). It is emitted by both diffuse and
redshift galaxies. Using ALMA, Harikane et al. (2020) studied dense gas in star-forming regions (Stacey et al. 1991; Brauher, Dale &

Downloaded from [Link] by guest on 21 May 2026

galaxies at a redshift of z > 6 and found that low [C/O] partially Helou 2008; Accurso et al. 2017; Franeck et al. 2018; Cormier et al.
explain the enhanced L[O III] /L[C II] and possibly the L[C II] /SFR (Star 2019). The [C II] 158 μm line is also often adopted as a tracer of Star
Formation Rate) ratios observed in these systems. Cosmological Formation Rate (SFR) of galaxies (De Looze et al. 2011; Pineda,
radiation hydrodynamics simulations of Arata et al. (2020) find that Langer & Goldsmith 2014; Herrera-Camus et al. 2015; Sutter et al.
the L[O III] /L[C II] luminosity ratio decreases with increasing [C/O], 2019; Bisbas et al. 2022; Liang et al. 2023) and as a tracer of the
which is eventually connected with the metal enrichment during molecular gas mass (Accurso et al. 2017; Combes 2018; Zanella
the galaxy evolution. Similar conclusions were addressed in the et al. 2018; Madden et al. 2020), including the dynamical evolution
models of Katz et al. (2022) who showed that low C/O ratios are of colliding clouds (Bisbas et al. 2018; Schneider et al. 2023).
needed to reproduce the observed relation of [CII]158 μ m –SFR and Atomic carbon (‘C I’ when referring to its emission or ‘C’ when
[OIII]88 μ m –SFR in z > 6 galaxies, as well as the observed higher referring to its abundance) is another important tracer of the ISM.
[OIII]88 μ m /[CII]158 μ m line ratio in the Epoch of Reionization. Katz It is emitted at frequencies of 492.2 and 809.3 GHz through the
et al. (2022) further suggest that this results in a top-heavy IMF in [C I] 3 P1 → 3 P0 (hereafter ‘[C I] (1–0)’) and 3 P2 → 3 P1 (hereafter
early-Universe star-forming galaxies. This finding is in agreement ‘[C I] (2–1)’) fine-structure lines, respectively. While the C layer is
with the Zhang et al. (2018) 13 C/18 O isotope ratio observations in believed to be thin and located between C+ and carbon monoxide
high-z galaxies, implying an increased population of massive stars. (CO) in classical one-dimensional PDRs (Draine 2011), it has been
In studying the IMF, the analytical work of Sharda et al. (2023b) found to accurately trace the molecular mass content of ISM gas
showed that negative [C/O] values in metal-poor ISM environments (Papadopoulos, Thi & Viti 2004; Bell, Viti & Williams 2007; Lo
impact the transition point from top-heavy to bottom-heavy IMF, et al. 2014; Offner et al. 2014; Zhang et al. 2014; Glover et al.
shifting it upwards by ∼0.5–1.0 dex in metallicity and that the 2015; Jiao et al. 2017; Papadopoulos, Bisbas & Zhang 2018; Jiao
characteristic mass (which sets the peak of the stellar IMF depending et al. 2019; Gaches, Offner & Bisbas 2019b; Bisbas, Tan & Tanaka
on the ISM environmental conditions; Sharda & Krumholz 2022) 2021; Dunne et al. 2022). This suggests that C emission originates
also increases by a factor of ∼7. This effect is a consequence of from a wider range of densities than those where the C abundance
the cooling processes in the ISM, since both carbon and oxygen peaks. The ratio of the two [C I] lines is often used to investigate the
are major cooling components in the total cooling function. It is properties of the observed gas (Bothwell et al. 2017; Valentino et al.
therefore expected that such negative [C/O] values, regardless of the 2020), including the environmental parameters of the ISM (Bisbas,
metallicity, will have a direct impact on the observables, especially Tan & Tanaka 2021; Bisbas et al. 2023).
12
concerning the carbon cycle as described below. CO (hereafter referred to as CO) is the most abundant molecule
(except for H2 ) in the ISM, and it is the most commonly used tracer of
the molecular gas in the MW and extragalactic objects (Tacconi et al.
1.4 Carbon cycle: background
2008; Genzel et al. 2012; Narayanan et al. 2012; Bolatto, Wolfire &
Photodissociation regions (PDR) are essential for understanding the Leroy 2013; Gong et al. 2020; Luo et al. 2020; Frias Castillo
astrochemistry of molecular clouds and the ISM at large. They are et al. 2023; Montoya Arroyave et al. 2023). It emits at different
the sites where the atomic-to-molecular transition occurs, as well wavelengths depending on the J → J − 1 transition. The collection
as the transition between the three carbon phases (ionized, atomic, of these transitions consists the Spectral Line Energy Distribution
and molecular in the form of carbon monoxide), known as the (SLED), which is used as an important diagnostics for studying the
‘carbon cycle’ (hereafter ‘C-cycle’; see reviews by Hollenbach & chemical and dynamical states of the ISM (Papadopoulos et al. 2010;
Tielens 1999; Wolfire, Vallini & Chevance 2022). However, the Narayanan & Krumholz 2014; Mashian et al. 2015; Rosenberg et al.
exact location of these phase transitions depends on a variety of 2015; Vallini et al. 2018; Klitsch et al. 2022; Stanley et al. 2023)
environmental parameters (Jura 1974; Black & Dalgarno 1977; van and a discriminant between the X-ray and FUV heating (Meijerink,
Dishoeck & Black 1986; Offner et al. 2013; Sternberg et al. 2014; Spaans & Israel 2007; Vallini et al. 2019; Esposito et al. 2022).
Bialy et al. 2015; Bisbas et al. 2023). The most important of those
include the radiation due to far-UV photons emitted by massive stars,
1.5 This work: C-cycle in extreme galactic environments
which photodissociate the molecules of H2 and CO (van Dishoeck &
Black 1988), the ionization rate due to the interaction of ISM gas with The CO molecule has a very high binding energy of ∼7.7 eV, making
charged particles carrying high energies at high column densities it stable against various chemical conditions. Intuitively speaking,
known as ‘cosmic-rays’ (see reviews by Strong, Moskalenko & molecular clouds with [C/O] < 0 would first bond most carbon in
Ptuskin 2007; Grenier, Black & Strong 2015), and the metallicity, the CO molecule. Only dissociation processes, such as UV radiation
which is described above and consists of the main focus of this work. fields, cosmic rays, and shocks, would release carbon into free atoms

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8889
and produce C and C+ . Therefore, to quantitatively evaluate such using the treatment of Cazaux & Tielens (2002, 2004, 2010). The H2
impacts, we perform detailed chemical modelling to simulate various self-shielding is calculated along each HEALPIX ray emanated from
galactic conditions. each computational cell and using the results of Lee et al. (1996).
This work aims to evaluate how an ISM component with a [C/O] < For the three-dimensional models, the ‘Dense’ cloud of Bisbas,
0 ratio affects the abundances and the line emission of the gas-phase Tan & Tanaka (2021) is used, which is a subregion of the binary
C-cycle as well as the location of the H I-to-H2 transition. To our collision simulation of giant molecular clouds presented in Wu
knowledge, the impact of a negative [C/O] on the C-cycle observables et al. (2017). The selected subregion has a resolution of 1123
and the H I-to-H2 transition has not been studied in the past. Particular uniform cells and contains a dense filamentary structure that remains
attention is drawn on how an [C/O] < 0 environment affects the C molecular under most ISM conditions explored and that undergoes
layer and the emission of both its fine-structure lines, and its impact star-formation. It has a size of L = 13.88 pc and a total mass of
on extragalactic observations. M = 5.9 × 104 M with a mean density of nH  = 640 cm−3 . In
A total of nineteen three-dimensional PDR models have been a full hydro-chemical simulation, it is expected that different ISM
performed in this regard using the density distribution of a molecular environmental parameters would affect hydrodynamical pressures
cloud that has been studied in previous works (Bisbas et al. 2017b; resulting in different density and velocity distributions. These would,

Downloaded from [Link] by guest on 21 May 2026

Bisbas, Tan & Tanaka 2021; Gaches, Bisbas & Bialy 2022a; Gaches in turn, impact the abundances and the corresponding line emission.
et al. 2022b) and considering various combinations of cosmic- In the present work, however, the astrochemical models have been
ray ionization rates, dust-to-gas ratios and FUV intensities to best performed in a static cloud; this provides a different and perhaps more
represent different galactic environments. We examine the response educational view of how negative [C/O] ratios affect the observables.
in the abundances of C+ , C, and CO, and we also perform radiative Numerical instabilities may arise in PDR models as a result of
transfer to produce velocity integrated emission maps of [C II] the very non-linear nature of the ordinary differential equations that
158 μm, [C I](1–0) and (2–1), and the first ten CO transitions (J = are solved until thermal balance is reached. In order to minimize
1–0 to 10–9). bi-stabilities (Le Bourlot et al. 1993; Viti et al. 2001; Roueff & Le
This paper is organized as follows. In Section 2, we present the Bourlot 2020; Dufour & Charnley 2021), which result in a significant
numerical approach followed. Section 3 shows the results of our noise in the calculations, a mixture of carbon phase abundances is
calculations including velocity integrated emission maps, column preferred in which the total abundance of carbon is always divided
densities, and line ratios. Section 4 discusses the impact of our results in three parts for C+ , C, and CO, respectively.
on observations. We conclude in Section 5.

2.1 Initial abundances

2 NUMERICAL METHOD
The fiducial (MW) gas-phase carbon abundance3 is taken to be
The PDR calculations in this study have been performed using the xC, MW = 1.4 × 10−4 and refers to the interstellar carbon abundance
3D-PDR2 code (Bisbas et al. 2012). 3D-PDR is a publicly available within the 600 pc solar neighbourhood (Cardelli et al. 1996; Röllig
astrochemical code that treats one- and three-dimensional pho- et al. 2007; Draine 2011). The fiducial gas-phase oxygen abundance,
todissociation regions. By taking into account several heating and is taken to be xO, MW = 3 × 10−4 (Cartledge et al. 2004; Draine 2011).
cooling processes, the code performs thermal balance calculations Both aforementioned abundances are lower than the reported Solar
and terminates when the total heating is approximately equal to the ones (Asplund et al. 2009) due to depletion on grains and have been
total cooling. Once the code is terminated, it outputs the abundances measured by optical/UV absorption lines in diffuse clouds (Cardelli
of species, gas and dust temperatures, heating and cooling functions et al. 1996; Cartledge et al. 2004). In addition, we assume a helium
as well as the level populations of various coolants, all versus the abundance of xHe = 1.00 × 10−1 . The fiducial model has, thus, an
column depth. [O/H] = −0.21 and a [C/O] = −0.07.
A subset of the UMIST2012 chemical network (McElroy et al. In both panels of Fig. 1, the larger shapes of circles, boxes,
2013) consisting of 33 species and 330 reactions is used. For the and triangles show the combinations of [O/H] and [C/O] (left-
purposes of this work, we only consider the elements of carbon, hand panel), and of [O/H] and D (right-hand panel; see 2.2 for
oxygen, helium, and hydrogen in the gas-phase. All abundances are the definition of D) considered here. In all our models, the notation
normalized to the total hydrogen. The intensity of the incident FUV ‘Oi Cj Dk -x-y’ is followed, where i is the identifier of the oxygen
radiation field, χ /χ 0 , is normalized according to the spectral shape of abundance, j of the carbon abundance, and k of the adopted D. The
Draine (1978). Two cosmic-ray ionization rates per H2 molecule are terms ‘x’ and ‘y’ correspond to the logarithm of the cosmic-ray
adopted; a low one with ζCR = 10−17 s−1 and a high one with ζCR = ionization rate (‘17’ for ζ CR = 10−17 or ‘15’ for 10−15 s−1 ) and to
10−15 s−1 . Such higher values of ζ CR are found near supernovae the logarithm of the FUV intensity field (no number for χ /χ 0 = 1 or
(Indriolo 2023) and are expected in galaxies with high star-formation ‘2’ for χ /χ 0 = 102 ), respectively. Table 1 shows a summary of all
rates (Papadopoulos 2010; Indriolo et al. 2018; Gaches, Offner & presented models.
Bisbas 2019a; Yang et al. 2023). In all simulations, the isotropic FUV In the left-hand panel of Fig. 1, individual observations of metal-
radiation field impinges radially and incident to each HEALPIX ray poor dwarf galaxies (Berg et al. 2016, 2019), star-forming galaxies
used in the ray-tracing scheme (Bisbas et al. 2012). In most of the of low-mass and low-metallicity at a redshift of z ∼ 0.1–0.3 known
models, the FUV intensity is taken to be χ /χ 0 = 1; however we also as ‘green pea galaxies’ (Ravindranath et al. 2020) as well as Lyman-
consider an additional and stronger intensity of χ /χ 0 = 102 for the continuum leaking galaxies at a redshift of z ∼ 0.3–0.4 (Izotov
ζCR = 10−15 s−1 models. At all times, the microturbulence velocity is et al. 2023) are shown with small diamonds. The purple diamond
set to vturb = 2.0 km s−1 . The formation of H2 molecule is modelled refers to the recent JWST observations of GN-z11 (Isobe et al.

2 [Link] 3 Hereafter, xsp is to denote the relative abundance of species ‘sp’.

MNRAS 527, 8886–8906 (2024)

8890 T. G. Bisbas et al.

Downloaded from [Link] by guest on 21 May 2026

Figure 1. Left-hand panel: Pairs of [O/H] and [C/O] used in models. Blue circles correspond to models O0 C0 and O0 C1 . Orange squares and green triangles
correspond to models O1 C0 , O1 C1 , O1 C2 according to the assumed D, which is the dust-to-gas ratio normalized to the solar value (see right-hand panel). The
horizontal dotted black line (‘Linear’) corresponds to the assumption of a constant [C/O] (equal to −0.07) and independent on [O/H]. The Sharda et al. (2023b,
c) (equation 2) cubic-fit, the Nicholls et al. (2017) (equation 3) best-fit, and the Garnett et al. (1995) (equation 4) best-fit are shown in dot–dashed, dashed,
and solid lines, respectively. Small brown diamonds refer to Lyman-continuum leaking galaxies at z ∼ 0.3–0.4 of Izotov et al. (2023), yellow and dark blue
to metal-poor dwarf galaxies of Berg et al. (2016) and Berg et al. (2019), respectively and red to green pea galaxies of Ravindranath et al. (2020). The purple
diamond refers to the Isobe et al. (2023) JWST observations of GN-z11. Right-hand panel: Relation of [O/H] versus D. The big blue circle, orange square, and
green triangle correspond to the values used in our models while their colours are in tandem with those of the left-hand panel. Blue and orange models follow
an approximately linear decrease of D as a function of [O/H] and are labelled as ‘D0 ’. The green triangle has a sublinear D and is labelled as ‘D1 ’. Black dotted
line corresponds to a linear connection between [O/H] and D to guide the eye. Grey diamonds correspond to the KINGFISH observations of Rémy-Ruyer et al.
(2014), while magenta diamonds to the Dwarf Galaxy Survey data of Rémy-Ruyer et al. (2013).

Table 1. Initial conditions and PDR parameters of the presented simulations. (2023b) expression.4 This expression is a cubic-fit of the Amarsi,
The first column gives the identifier of each simulation. To ease the reader, Nissen & Skúladóttir (2019) 3D non-LTE calculations in a sample
we represent the key models with the colours and shapes shown in Fig. 1 of 187 stars that exist in the thin and thick discs and the metal-poor
and elsewhere. The ‘x’ notation refers to the negative logarithm of the
halo of the MW. The dashed line shows the Nicholls et al. (2017)
cosmic-ray ionization rate and it is either 17 (for ζCR = 10−17 s−1 ) or 15 (for
ζCR = 10−15 s−1 ). Models including the ‘y’ notation refer to the additional
expression.5 This expression is a best-fit of various observations
calculations with a higher FUV intensity (2 for χ /χ 0 = 102 , otherwise no including disc and halo MW stars and damped Lyα systems (see
number noted). The second and third columns show the [O/H] and [C/O] Berg et al. 2016, 2019, and references therein). The thin solid line
ratios, respectively. The fourth and fifth column refers to the fractional corresponds to the Garnett et al. (1995) best-fit of least squares.6
gas-phase abundances with respect to hydrogen nuclei used in the 3D-PDR
simulations. Simulations O0 C0 D0 , O1 C0 D0 , and O1 C0 D1 follow a linear
decrease of the C/O ratio (they have all C0 ). The sixth column shows D, 2.2 Dust-to-gas mass ratio
which is the dust-to-gas ratio normalized to the solar one of 10−2 . D0 refers The reference value of the dust-to-gas ratio is taken to be 10−2
to an approximately linear decrease of D, whereas D1 to a sublinear. In all
(Sandstrom et al. 2013). We further define with D the dust-to-gas ratio
cases, the bracketed numbers indicate the order of magnitude.

ID [O/H] [C/O] xC xO D 4 The Sharda et al. (2023b) expression is given by:

O0 C0 D0 -x −0.21 −0.07 1.40(−4) 3.00(−4) 1 [C/O]MW = aS [O/H]MW 3 + bS [O/H]MW 2 + cS [O/H]MW + dS , (2)
O0 C1 D0 -x-y −0.21 −0.57 4.43(−5) 3.00(−4) 1
O1 C0 D0 -x −1.21 −0.07 1.40(−5) 3.00(−5) 3(−2) where aS = −0.02, bS = 0.14, cS = 0.6, and dS = −0.09. Note that the
O1 C1 D0 -x −1.21 −0.57 4.43(−6) 3.00(−5) 3(−2) above presented in Sharda et al. (2023b) uses the MW abundances as a
O1 C2 D0 -x-y −1.21 −0.77 2.80(−6) 3.00(−5) 3(−2) reference value (see Sharda et al. 2023c, for correction). The corresponding
O1 C0 D1 -x −1.21 −0.07 1.40(−5) 3.00(−5) 3(−3) curve of Fig. 1 is plotted after some further calculations considering the Solar
O1 C1 D1 -x −1.21 −0.57 4.43(−6) 3.00(−5) 3(−3) abundances adopted here as the reference value.
5 The Nicholls et al. (2017) expression is given by:
O1 C2 D1 -x-y −1.21 −0.77 2.80(−6) 3.00(−5) 3(−3)
 
x
   
[O/H]+log10 xO +bN xC
[C/O] = log10 10aN + 10 H  − log10 (3)
xO 

2023). The horizontal dotted black line refers to the constant [C/O] where aN = −0.8 and bN = 2.72.
6 The Garnett et al. (1995) expression is given by
(equal to −0.07), which corresponds to the commonly assumed linear
     
relationship between carbon and oxygen. In addition, three best-fit xO xC
[C/O] = 0.44 [O/H] + log10 + 1.14 − log10 . (4)
relations are illustrated. The dot–dashed line shows the Sharda et al. xH  xO 

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8891
normalized to the above value, thus D = 1 is the fiducial one. This [C I](1–0), and CO J = 1–0. The latter quantity is defined as
ratio affects both the opacity and the grain photoelectric heating. The 
W dS
right-hand panel of Fig. 1 shows the observed relationship between W  =  , (5)
D and [O/H]. The dotted black line shows the linear connection dS
between [O/H] and D that is frequently adopted. As described in where W is the velocity integrated emission of the line (in units of
the Introduction this linear correlation of the dust-to-gas ratio with K km s−1 ) and S is the area of the whole map, without adopting a
metallicity, may not be valid for very metal-poor galaxies (Rémy- lower observational limit. As can be seen from columns 3 and 4, the
Ruyer et al. 2014). This is further shown with the individual small cloud remains molecular for most of the cases (xH I  < xH2 ), except
triangles, which correspond to observations from the Dwarf Galaxy for the extreme conditions of O1 C2 D0 -15-2 and for all O1 Cj D1 -15-y
Survey (Rémy-Ruyer et al. 2013) and the KINGFISH survey of models (last five rows of Table 2).
nearby galaxies (Rémy-Ruyer et al. 2014). To accommodate this ob- The above results are visualized in Figs 2 and 3 where the
served behaviour and explore its imprint on the C-cycle, we consider aforementioned average abundance ratios and emissions of the C-
two D values: D = 3 × 10−2 (which is a 30 × decrease as opposed cycle are shown, respectively. In particular, Fig. 2 shows the relative
to the 10 × decrease in abundances) to represent an approximately behaviour of the above species for ζCR = 10−17 s−1 (top panel) and

Downloaded from [Link] by guest on 21 May 2026

linear/sublinear scaling, as well as a much lower D of 3 × 10−3 . Our ζCR = 10−15 s−1 (bottom panel). This relative behaviour refers to
calculations assume a standard MRN distribution (Mathis, Rumpl & the normalization of the [C/O] < −0.07 models (C1 and C2 ) with
Nordsieck 1977) as described in Bisbas et al. (2012). the corresponding ones for [C/O] = −0.07 at a fixed [O/H] (e.g.
Models O0 C0 D0 refer to the standard MW carbon and oxygen O0 C1 D0 -17/O0 C0 D0 -17 etc.). For the higher-FUV runs (O0 C1 D0 -15-
abundances. In model O0 C1 D0 , the carbon abundance is reduced 2 and O1 C2 D1 -15-2), these are normalized with their counterparts of
by half an order of magnitude while keeping D constant. Such a χ /χ 0 = 1 (e.g. O0 C1 D0 -15 etc.).
combination between [O/H] and [C/O] has been observed in the It is found that the atomic-to-molecular mass content remains
GN-z11 galaxy using JWST (Isobe et al. 2023), which is located at a unaffected as a function of [C/O] and [O/H] for fixed ζ CR (the
redshift of z  10.6 (Oesch et al. 2016). In particular, the O0 C1 D0 -15 models with higher FUV intensity are described below). In all cases,
model is to represent conditions that are met in metal-rich starburst the C+ , C, and CO abundances decrease by a factor between 3–
galaxies while O0 C1 D0 -15-2 is to explore also the effect of an FUV 5 when compared to the reference [C/O] = −0.07 value. Notably,
enhancement. in the higher ζ CR models the C abundance can be decreased up to
Models O1 C0 D0 use carbon and oxygen abundances that are both approximately 15 × (simulation O0 C1 D0 -15).
reduced by one order of magnitude (linear decrease). This model The stronger FUV intensity in the selected ζCR = 10−15 s−1 sim-
is used as a benchmark for the O1 C1 D0 and O1 C2 D0 , which use ulations affect all aforementioned species. In particular, the stronger
the same (reduced) oxygen abundance but a further suppressed χ /χ 0 slightly increases xHI  by approximately 2.5 × in the O0 C1 D0 -
carbon abundance by ∼3.15 (C1 ) and 5 (C2 ) times corresponding to 15-2 simulation while the cloud remains molecular. However, in the
[C/O] = −0.5 and [C/O] = −0.7, respectively. In all these cases, lower metallicity model of O1 C2 D0 -15-2, the stronger FUV radiation
the normalized dust-to-gas ratio is D = 3 × 10−2 (D0 ). The effect pushes the H I-to-H2 transition towards higher column densities
of a lower, sublinear D is also explored and indicated with ‘D1 ’ in increasing its xHI . This effect is more prominent in the O1 C2 D1 -15-2
the ID of the corresponding model. Recently, Arellano-Córdova et al. model, where the additionally lower dust-to-gas ratio allows the FUV
(2022) reported JWST observations (not shown in the left-hand panel photons to propagate even further, thereby reducing the abundance of
of Fig. 1) of the z > 7 galaxy cluster SMACS J0723.3-7327 having H2 even more. A similar behaviour is also reflected in the C-cycle in
[O/H]  −1.57 and [C/O]  −0.57 (see also fig. 4 of Jones et al. the aforementioned models, where an increase in C+ and a decrease
2023 for a more complete picture). in CO are observed due to the photodissociation of the latter.
Fig. 3 shows the relative emission lines from the C-cycle. These
include the CO SLED for the first 10 transitions, both [C I] and the
2.3 Synthetic observations and radiative transfer
[C II] 158 μm fine-structure lines. For ζCR = 10−17 s−1 (top panel),
To calculate the line emission of the C-cycle and thus obtain the the high-J CO lines become weaker for lower [C/O] ratios. For J =
corresponding synthetic observations (see Haworth et al. 2018, for 6–5 and above, the emission of CO is found to be approximately 2 ×
a review), the equation of radiative transfer is solved along the line- weaker. The [C I] (1–0) line becomes on average 2.5 × weaker as
of-sight. For each C-cycle coolant, the methodology described in [C/O] decreases at a fixed [O/H]. The emission of both [C I] (2–1) and
Bisbas et al. (2017b) is followed, including the updates described in [C II] is also decreasing, albeit at a smaller factor. For ζCR = 10−15 s−1
Bisbas, Tan & Tanaka (2021) to account for the dust contribution. and low FUV intensity (bottom panel), it is found that the relative
We redirect the reader to these works for further details. CO SLEDs of the models O0 C1 D0 -15 and O0 C0 D0 -15 are mostly
affected, with the (low-to-mid)-J lines to become twice as bright.
This effect is largely connected with the efficient destruction of CO
3 R E S U LT S
in the more diffuse gas and also with the increase of the local heating
due to cosmic-rays (Bisbas, Papadopoulos & Viti 2015; Bisbas et al.
3.1 General behaviour
2017a; Bisbas, Tan & Tanaka 2021). This makes the emission of the
Table 2 shows a summary of all simulation results. Column 2 line to be closely associated with a more dense and slightly warmer
refers to the most representative type of objects based on the ISM gas, hence its increase in the (low-to-mid)-J transitions (Gaches et al.
environmental parameters of each simulation. These types include 2022b). All models with [C/O] < −0.07 and χ /χ 0 = 1 appear to be
local clouds in the solar neighbourhood, post-starburst, starburst, dimmer for the CO J = 7–6 transition and above. Both [C I] emission
and dwarf galaxies with different C/O and dust-to-gas ratios. From lines are very suppressed (≈3–10 ×) and [C II] 158 μm becomes also
column 3 onward, the density-weighted average abundances of H I, weaker by a factor of 2–3.
H2 , C+ , C, and CO are shown, as well as the average gas temperature In the O0 C1 D0 -15-2 model, the higher FUV radiation field does not
and the average velocity integrated emission of [C II] 158 μm, affect much the CO SLED. However, the CO SLED is significantly

MNRAS 527, 8886–8906 (2024)

8892 T. G. Bisbas et al.
Table 2. Average abundances and velocity integrated emissions of all simulations. Column 1 shows the IDs of all models. Column 2 refers to the object type
each simulation best represents. Columns 3–7 refer to the average (density-weighted) abundances of H I, H2 , C+ , C, and CO, respectively. The numbers in the
parenthesis denote the order of magnitude. Column 8 refers to the average (density-weighted) gas temperature (K). Columns 9–11 refer to the area average
emissions of [C II] 158 μm, [C I](1–0), and CO J = 1–0, respectively. The average emission is given in units of K km s−1 .

ID Representative type xH I  xH2  xC+  xC  xCO  Tgas  [C II] [C I](1–0) CO(1–0)

O0 C0 D0 −17 Solar Neighbourhood 0.27(−2) 0.49(0) 1.07(−5) 5.42(−6) 1.21(−4) 11.8 0.21 2.53 45.32
O0 C1 D0 −17 Local post-starbursts 0.25(−2) 0.49(0) 4.03(−6) 1.12(−6) 3.87(−5) 14.7 0.14 1.01 51.11
O1 C0 D0 −17 Dwarfs (C/OL-S ) 0.88(−1) 0.45(0) 4.47(−6) 3.55(−6) 5.84(−6) 17.6 0.23 4.59 12.64
O1 C1 D0 −17 Dwarfs (C/OObs ) 0.82(−1) 0.45(0) 1.46(−6) 1.11(−6) 1.80(−6) 24.7 0.20 2.23 10.82
O1 C2 D0 −17 Dwarfs (C/OExt ) 0.80(−1) 0.46(0) 9.32(−7) 7.27(−7) 1.13(−6) 28.2 0.18 1.62 9.49
O1 C0 D1 −17 Dwarfs (C/OL-S , oD) 0.27(0) 0.36(0) 5.61(−6) 5.01(−6) 3.26(−6) 15.2 0.17 4.72 5.36
O1 C1 D1 −17 Dwarfs (C/OExt , oD) 0.25(0) 0.37(0) 1.85(−6) 1.54(−6) 9.90(−7) 20.8 0.16 2.73 4.89
O1 C2 D1 −17 Dwarfs (C/OL-S , oD) 0.25(0) 0.37(0) 1.18(−6) 9.83(−7) 6.20(−7) 23.1 0.15 1.97 4.31
O0 C0 D0 −15 Solar Neighbourhood (hCR) 0.20(−1) 0.49(0) 2.27(−5) 3.29(−5) 8.22(−5) 28.3 2.77 41.42 67.63

Downloaded from [Link] by guest on 21 May 2026

O0 C1 D0 −15 Local ULIRGs (hCR) 0.20(−1) 0.49(0) 6.26(−6) 2.30(−6) 3.53(−5) 41.2 1.26 4.57 124.63
O0 C1 D0 -15-2 Local ULIRGs (hCR, hUV) 0.50(−1) 0.47(0) 9.68(−6) 2.87(−6) 3.13(−5) 86.8 3.93 5.93 91.01
O1 C0 D0 −15 S-B Dwarfs (C/OL-S ,hCR) 0.16(0) 0.41(0) 4.39(−6) 5.09(−6) 4.39(−6) 46.1 2.65 10.24 24.75
O1 C1 D0 −15 S-B Dwarfs (C/OObs ,hCR) 0.17(0) 0.41(0) 1.40(−6) 1.06(−6) 1.91(−6) 61.3 1.36 2.42 17.79
O1 C2 D0 −15 S-B Dwarfs (C/OExt ,hCR) 0.17(0) 0.41(0) 8.86(−7) 6.24(−7) 1.28(−6) 66.0 0.93 1.48 13.80
O1 C2 D0 -15-2 S-B Dwarfs (C/OObs ,hCR,hUV) 0.38(0) 0.30(0) 1.70(−6) 6.95(−7) 3.83(−7) 93.6 3.19 1.56 3.00
O1 C0 D1 −15 Dwarfs (C/OL-S , oD,hCR) 0.40(0) 0.29(0) 5.19(−6) 6.45(−6) 2.24(−6) 34.4 1.50 11.31 10.08
O1 C1 D1 −15 Dwarfs (C/OObs , oD,hCR) 0.40(0) 0.29(0) 1.72(−6) 1.74(−6) 9.22(−7) 41.8 0.99 3.67 7.88
O1 C2 D1 −15 Dwarfs (C/OExt , oD,hCR) 0.41(0) 0.29(0) 1.10(−6) 1.06(−6) 6.21(−7) 44.9 0.77 2.33 6.09
O1 C2 D1 -15-2 Dwarfs (C/OExt , oD,hCR, hUV) 0.73(0) 0.13(0) 2.13(−6) 5.59(−7) 1.00(−7) 46.7 2.01 1.13 0.64
Note.  ULIRGs: Ultra Luminous Infrared galaxies, S-B: Starburst, L-S: Linearly-scaled, Obs: Observed, Ext: Extreme, oD: observed dust-to-gas ratio, hCR: high
cosmic-ray ionization rate, hUV: high FUV intensity

affected in the models of lower [O/H] = −1.21 values, especially As expected, the column densities of the C-cycle species (C+ , C,
for all transitions with J > 2–1 (models O1 C2 D0 -15-2 and O1 C2 D1 - and CO) are reduced with decreasing the [C/O] ratio at fixed [O/H]
15-2). This is because the abundances of CO and H2 molecules (e.g. O0 C0 D0 -17 with O0 C1 D0 -17 etc.). This is also evident from the
decrease in those models resulting in a weaker excitation of CO average abundances shown in Table 2 and also in the top panel of
for mid/high-J transitions. The emission of both [C I] lines is slightly Fig. 2. Interestingly, all these three species are reduced by the same
enhanced for the higher FUV simulations, except for the O1 C2 D1 -15- amount and by observing their column density maps, it can be seen
2 one, since the much lower dust-to-gas ratio contributes (through that no significant changes in the structure of those column density
the propagation of radiation) to the ionization of atomic carbon. maps can be found.
Contrary, the emission of [C II] 158 μm line is enhanced by a factor However, considering fixed [C/O] = −0.07 and a decreasing
of ≈3 at all cases of higher FUV intensity. [O/H] ratio (models O0 C0 D0 -17, O1 C0 D0 -17, and O0 C0 D1 -17), and
In Appendix A, the behaviour of the most commonly used line focusing on the maps of N(C+ ), N(C), and N(H2 ), it can be seen
ratios between CO and [C I] is discussed. that the column densities of C+ and C become more closely related
with the filamentary structure (see also Bisbas, Schruba & van
Dishoeck 2019, for a relevant discussion). This connection is more
3.2 Low cosmic-ray ionization rate evident in the O1 C0 D1 -17 model in which N(CI)  1018.5 cm−2 (the
Fig. 4 shows (from top to bottom) the column densities of H2 , highest C column density in all ζCR = 10−17 s−1 models) implying
H I, C+ , C, and CO of the most representative models. The bottom that in such conditions, atomic carbon may provide more accurate
row shows the density-weighted gas temperature. The first pair of measurements of the molecular mass content in the ISM. On the
columns correspond to models O0 C0, 1 D0 -17, the second pair to other hand, the photodissociation of CO at lower [O/H] (due to
models O1 C0, 2 D0 -17 (all O1 C1 models are omitted in these plots) the further penetration of FUV photons at higher column densities)
and the third pair to models O1 C0, 2 D1 -17, respectively. reduces even more the abundance and therefore the column density
The gas within the dense, star-forming filamentary structure of this molecule. It is, thus, found that CO is abundant only at the
remains always molecular as can be seen from the first row of very dense parts of the filamentary structure where star-formation is
Fig. 4. The H2 column density of that structure does not change, likely to occur.
remaining almost completely unaffected in all cases. However, the Finally, the density-weighted gas temperature seen in the bottom
N(H2 ) of the gas in the more diffuse medium surrounding the filament panel of Fig. 4 shows an interesting behaviour. The average gas
is reduced particularly in the O1 C0 D1 -17 and O1 C2 D1 -17 models temperature increases with decreasing [C/O] at fixed [O/H] (see
due to the lower dust-to-gas ratio, D, as the latter favours the H2 also column 8 in Table 2) since carbon is a major PDR coolant
photodissociation process. On the other hand, the H I column density and therefore its suppression in abundance leads to a reduction in
remains low in the O0 C0, 1 D0 -17 models, increases in the O1 C0, 2 D0 - its cooling efficacy. However, the densest part of the filamentary
17 models, and eventually increases further in the O1 C0, 2 D1 -17 structure remains always within Tgas  ∼ 10–20 K. Perhaps the most
models. Although the column densities of H2 and H I change with striking feature is the gas temperature comparison between the pair
[O/H] for fixed [C/O] (e.g. moving ‘across’ the left-hand panel of of models O1 C0 D0, 1 -17, and the pair of O1 C2 D0, 1 -17. Although
Fig. 1) and also with D, they do not change as a function of [C/O] at Tgas in both these models is higher than the corresponding one of
fixed [O/H] (e.g. moving ‘top-to-bottom’ of that latter panel). the O0 C1 D0 -17 model, O1 C2 D0 -17 appears to have on average the

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8893

Downloaded from [Link] by guest on 21 May 2026

Figure 2. Comparative predictions of abundances for H I, H2 , C+ , C,
CO, and the average gas temperature. The top panel presents models
for ζCR = 10−17 s−1 , while the bottom panel depicts the ζCR = 10−15 s−1
models. Consequently, ‘x’ in the legend identifies the ‘17’ and ‘15’ models for
the top and bottom panels, respectively. Abundances are normalized according
to their designated reference models. In particular, the top panel reference
Figure 3. Relative emission lines of the CO SLED, both [C I] lines and the
models are characterized by [C/O] = −0.07 for a constant [O/H], while
[C II] 158 μm line. Dashed lines indicate models with D = 3 × 10−3 . The
the bottom panel reference models align with the corresponding lower FUV
rest of the notation follows the one of Fig. 2.
intensity in the ISM. Models marked with open circles adopt a sublinear
D = 3 × 10−3 .
Focusing now on the changes observed in the line emission due
to the decrease of [C/O] at fixed [O/H], it is found that [C II]
warmest ISM with a Tgas   28.3 K. The lower D ratio, although 158 μm although it is originated from different parts of the cloud,
allows the FUV photons for reaching higher column densities, it is not strongly influenced by the carbon abundance decrease and its
increases xC+ and xC as described above, which in turn contribute brightness remains relatively unchanged as the abundance decrease
positively to the total cooling. A lower D also reduces the efficiency compensates with the gas temperature increase (see top panel of
of photoelectric heating in places. Notably, the spread of gas tem- Fig. 3). However, [C I] (1–0) decreases with the decreasing [C/O] by
peratures in both O1 C0 D1 -17 and O1 C2 D1 -17 models is smaller than an average factor of ∼2–3 as can be seen from Table 2. Finally, CO
in the other simulations, making overall a cloud of approximately J = 1–0 shows a different picture; although the abundance of CO
uniform Tgas . It is noted that our models show that the photoelectric molecule is lower, the emission of this line becomes stronger. This is
heating and the cosmic-ray heating (the two most dominant heating due to the increase of gas temperature by a few K everywhere in the
mechanisms) are not affected by the decrease of [C/O] ratio but rather cloud (࣠1.5 × on average) and particularly in the dense gas where
by the lower D ratio. this line originates from.
Fig. 5 shows (from top to bottom) the velocity integrated emission
maps of [C II] 158 μm, [C I] (1–0) and CO J = 1–0 for the above
3.3 High cosmic-ray ionization rate
models. As with the behaviour of the C-cycle column densities
described earlier, the [C II] 158 μm and [C I] (1–0) lines are better Fig. 6 shows the column densities of species as described in Fig. 4 but
connected with the N(H2 ) map for lower [O/H] values (O1 models) for the case of higher cosmic-ray ionization rate of ζCR = 10−15 s−1 .
therefore originating from denser gas than the one in the O0 C0, 1 D0 - The majority of the gas remains molecular in the O0 C0, 1 D0 -15
17 models. In the particular O1 C0 D0, 1 -17 models, the additionally and O1 C0, 1, 2 D0 -15 models (top row of Fig. 6), since the H I-to-H2
higher gas temperature results in the brightness increase of the [C II] transition is not heavily affected by the cosmic-ray ionization rate
158 μm and [C I] (1–0) emission lines. (Bialy et al. 2015; Bisbas, Papadopoulos & Viti 2015). However, in

MNRAS 527, 8886–8906 (2024)

8894 T. G. Bisbas et al.

Downloaded from [Link] by guest on 21 May 2026

Figure 4. Column densities of species for ζCR = 10−17 s−1 . The first pair of columns correspond to simulations O0 C0, 1 D0 -17, the second pair to O1 C0, 2 D0 -17
and the third pair to O1 C0, 2 D1 -17 in which D = 3 × 10−3 . From top-to-bottom, the column densities of H2 , H I, C+ , C, and CO are shown, respectively. In
general, the column densities of all these species decrease with decreasing [C/O], except for the H I and H2 , which depend more strongly on the [O/H] ratio.
The bottom row shows the density-weighted gas temperature, which increases with decreasing [C/O].

the models of lower D (O1 C0, 1, 2 D1 -15) the gas remains molecular estimates a ζCR ∼ 5 × 10−16 s−1 , which is slightly lower than the one
only in regions of higher densities, where star-formation is likely to used.
occur. As can be seen in Table 2 for these models, χ H I  > χ H2 , The response of the C-cycle in an ζ CR -enhanced environment has
thus the cloud is mainly atomic. been extensively studied in previous works (Bisbas, Papadopoulos &
The column densities of H I (second row of Fig. 6) are enhanced by Viti 2015; Bisbas et al. 2017a; Bisbas, Tan & Tanaka 2021). Similarly
a factor between ∼2 × (in models O1 C0, 2 D1 -15) and ∼10 × (in mod- to those studies, it is found here that CO is converted to C+ and
els O0 C0, 1 D0 -15). An interesting feature occurs, however, for both mainly to C through the interaction of He+ resulting from cosmic-
O0 C0, 1 D0 -15. In these models, N(H I) appears to be approximately rays. However, the increase of χ C is not as effective for lower [C/O]
constant (N(HI)  8.1 × 1020 cm−2 with a standard deviation of values. In fact, we find that in the O0 C1 D0 -15 model, C+ transitions
σN(HI)  9.6 × 1019 cm−2 ). Notably, this is in good agreement with to CO almost directly without having a C-dominated layer. This
H I measurements in the Perseus cloud (Burkhart et al. 2015) and in turn lowers N(C), as evident from the corresponding panel in
the IC 348 cloud (Luo et al. 2023). Interestingly, the latter work Fig. 6 (fourth row, second column). In particular, it is estimated that

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8895

Downloaded from [Link] by guest on 21 May 2026

Figure 5. Emission maps for ζCR = 10−17 s−1 . The first pair of columns correspond to simulations O0 C0, 1 D0 -17, the second pair to O1 C0, 2 D0 -17 and the
third pair to O1 C0, 2 D1 -17. From top-to-bottom, the velocity integrated emissions of [C II] 158 μm, [C I](1–0) and CO J = 1–0 are shown, respectively. As [C/O]
decreases with respect to the reference value, the emission of [C I](1–0) decreases substantially while CO J = 1–0 remains bright. Note, however, that as [O/H]
decreases, [C I](1–0) increases. The effect is stronger for sublinear D due to the additional CO photodissociation.

in the O0 C0 D0 -15 model, ∼57.7 per cent of the total carbon is in due to the enhanced ζ CR value, the higher cosmic-ray heating is
CO, ∼23.1 per cent in C and ∼19.1 per cent in C+ , whereas in the consequently increasing the average temperatures 2–3 ×, affecting,
O0 C1 D0 -15 model these percentages are ∼80.5, ∼5.2, and ∼14.3 in turn, all corresponding line emissions described below.
per cent, respectively. Similarly, although the cloud in O1 C0 D0 - Fig. 7 shows the velocity integrated emission maps of [C II]
15 is C-dominated, it becomes CO-dominated in O1 C2 D0 -15. The 158 μm, [C I] (1–0), and CO J = 1–0. In general, all these emission
additionally lower dust-to-gas ratio in the O1 C0 D1 -15 and O1 C2 D1 - lines are much stronger in the O0 C0 D0 -15 model when compared
15 models favour both C+ and C carbon phases. Indeed, the cloud to the O0 C0 D0 -17 one, as was also discussed in Bisbas, Tan &
remains poor in CO, which is found only in the very high density Tanaka (2021); Bisbas et al. (2023). The most striking feature found
areas (nH  3 × 104 cm−3 ). in this figure, is that both [C II] 158 μm and [C I] (1–0) become
By performing further analysis, it is found that for ζCR = significantly dimmer when [C/O] decreases for fixed [O/H]. For
10−17 s−1 , C peaks for a narrow range of densities. The width of instance, recent NOEMA observations of [C II] in the GN-z11 galaxy
this range does not change for lower [C/O] values, so the column (Fudamoto et al. 2023) find a weak emission of this line, in agreement
density of C decreases approximately linearly with the decreasing C with our models for negative values of [C/O]. In regards to the
abundance. However, for ζCR = 10−15 s−1 and for [C/O] = −0.07, [C I] (1–0) line and even though the cosmic-ray ionization rate
a wider range of densities are C-rich (Bisbas, Papadopoulos & Viti is high, it is weak and remains much weaker than CO J = 1–0.
2015). For lower [C/O], the width of this range becomes significantly The latter emission line on the other hand becomes approximately
narrower. Chemical analysis show that this occurs because lower twice brighter in the O0 C1 D0 -15 model than in O0 C0 D0 -15 and
[C/O] values at high ζ CR , influence the formation pathway of C. In remains approximately unaffected by the lower [C/O] in the low-
particular, it is found that C is mainly formed via the reactions: metallicity models (O1 ) and even for sublinear values of D. This in
turn favours the CO J = 1–0 line as a molecular gas tracer in such
H + CH → C + H2 , (R1) an environment. The fact that under these conditions the molecular
gas can remain bright in CO J = 1–0 but dark in [C I](1–0) may have
CH+ −
3 + e → C + H2 + H, (R2) important consequences in observations of high-redshift star-forming
+ galaxies which can potentially have lower [C/O] than expected (see
instead from the C recombination with free electron which has
a formation rate of 90 per cent for [C/O] ∼ 0. For [C/O] < 0, Section 4).
reactions R1 and R2 can dominate the C formation pathway by
∼50–70 per cent, leading to a C-poor gas, since the abundances of
3.4 High cosmic-ray ionization rate and FUV intensity
CH and CH+ 3 are low.
The bottom row of Fig. 6 shows the density-weighted gas temper- We will now explore the case, where both ζ CR and FUV intensity
ature. The pattern discussed for the corresponding panels in Fig. 4 are enhanced. Fig. 8 shows how the column density maps presented
is also repeated here; the average gas temperature increases with in Section 3.3 change if an FUV field with intensity χ /χ 0 = 102
decreasing [O/H] and [C/O] (see also Tgas  in Fig. 2, where the is applied (for ζCR = 10−15 s−1 ). The models calculated consider the
relative ratios are essentially the same for both ζ CR ). However, lower [C/O] values for the explored [O/H]. The combination of high-

MNRAS 527, 8886–8906 (2024)

8896 T. G. Bisbas et al.

Downloaded from [Link] by guest on 21 May 2026

Figure 6. As in Fig. 4 but for ζCR = 10−15 s−1 . The particular abundance of C decreases severely in model O0 C1 D0 -15 and substantially in models O1 C2 D0 -15
and O1 C2 D1 -15, even though a higher cosmic-ray ionization rate is present. It is therefore found that knowledge of the [C/O] ratio is highly important for a
better understanding of the C-cycle emission under extreme environmental conditions. As expected, the column densities of CO are significantly reduced in
comparison to the ζCR = 10−17 s−1 ionization rate (Bisbas, Papadopoulos & Viti 2015; Bisbas et al. 2017a).

ζ CR and high-FUV represents better the conditions of star-forming structure. As discussed in Section 3.3, the column density of H I in the
metal-rich and metal-poor galaxies. In general, it is expected that O0 C1 D0 -15-2 model is approximately uniform with N(HI)  1.4 ×
a galaxy of higher star-formation rate, will also contain on average 1021 cm−2 and with a standard deviation of σN(HI)  3.5 × 1020 cm−2
higher FUV radiation fields and higher cosmic-ray ionization rates. (see also O0 C0 D0 -15 and O0 C1 D0 -15 in Fig. 6). Considering all
While the exact relationship between these quantities is not yet simulations presented in this work, the H I column density obtains its
observationally determined, they are frequently assumed to scale highest average value of N(HI)  2.1 × 1022 cm−2 in the O1 C2 D1 -
in a linear fashion (e.g. Papadopoulos 2010). 15-2 model.
It is found that lower [O/H] values lead to an increase in atomic In regards to the C-cycle, it is found that the cloud is CO-rich in
fraction. More specifically, when the dust-to-gas ratio D is lower, the the O0 C1 D0 -15-2 model but becomes C+ -rich in the O1 C2 D0 -15-2
cloud becomes H I-dominated. In the particular case of O1 C2 D1 -15- and O1 C2 D1 -15-2 models. However, the column density maps of all
2, the higher χ /χ 0 and the higher ζ CR , dissociates the H2 molecule carbon phases agree well with the structure of the N(H2 ) maps. It
almost everywhere in the cloud, except for the dense filamentary is interesting to note that in both O1 C2 D0 -15-2 and O1 C2 D1 -15-2

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8897

Downloaded from [Link] by guest on 21 May 2026

Figure 7. As in Fig. 5 but for ζCR = 10−15 s−1 . The behaviour of the C-cycle emission remains similar albeit all lines are much brighter due to the higher
cosmic-ray heating.

simulations, the molecular gas contains a small amount of CO, thus 2020) and [C/O] < 0. Remarkably, its brightness remains approxi-
making the H2 -rich cloud to be severely CO-dark. mately constant even for sublinear D (O1 C2 D0, 1 -15-2 simulations).
The gas temperatures are shown in the bottom row. The variations For the above simulations of high ζ CR and high χ /χ 0 , it is also found
in Tgas in the different parts of the cloud and between the aforemen- that [O I] 63 μm increases with decreasing [C/O] and decreasing D,
tioned models appear to be somehow counter-intuitive. In particular, contributing significantly to the total cooling.
the diffuse part of the cloud surrounding the dense filamentary struc-
ture obtains the highest Tgas value for the O0 C1 D0 -15-2 simulation.
Although in O1 C2 D0 -15-2 and O1 C2 D1 -15-2 models both metallicity 3.5 Relation of line intensities with H2 column density
and D are lower, the gas temperature in the diffuse part is not higher
as expected. This is better seen in the O1 C2 D1 -15-2 simulation. This Fig. 10 shows the mean values of the velocity integrated emission
effect is due to the increase of C+ abundance, and consequently its of [C II] 158 μm (top row), [C I] (1–0) (middle row), and CO J =
line emission, which cools down the gas (see also Bisbas et al. 2017a). 1–0 (bottom row) versus the H2 column density for ζCR = 10−17 s−1 .
On the other hand, the dense gas along the filament is warmer in the [O/H] and D decrease from left to right. The different line styles
lower [O/H] models as a result of photoelectric heating due to the correspond to different [C/O] values. With the exception of [C II]
lower D compared to the O0 C1 D0 -15-2 model. A double behaviour 158 μm for [O/H] = −0.21, it is found that the emission of all
is therefore seen as [O/H] decreases for negative [C/O]; a decrease lines increases with N(H2 ) and eventually saturates. The different
in Tgas for densities corresponding to diffuse gas and an increase in behaviour of the [C II] emission in the top left-hand panel is due to
Tgas for the higher density gas. its different origin from the ISM, other than the molecular phase,
Fig. 9 shows the corresponding velocity integrated emission maps. and thus no correlation with N(H2 ) can be found. Looking at all
The emission lines follow closely the pattern of the corresponding panels, perhaps the most significant feature is the decrease of the
column density maps of the C-cycle. CO J = 1–0 becomes very W(C I 1–0) emission with [C/O] for fixed [O/H], as has been also
bright and extended in the O0 C1 D0 -15-2 simulation and is much discussed earlier. W(C II) remains weak for all [O/H] and W(C I 1–0),
brighter than [C I] (1–0) and [C II] 158 μm. The most interesting although increases its brightness for [O/H] = −1.21 and for an H2
feature, however, can be seen in the O1 C2 D0 -15-2 and particularly column density of a few ×1022 cm−2 , it becomes weaker with the
in the O1 C2 D1 -15-2 models. In both, these simulations, CO J = 1–0 [C/O] decrease.
originates only from the very dense clumps of the cloud while it is Interestingly, the correlation of W(CO 1–0) with N(H2 ) remains
very dark otherwise. This makes the H2 -rich part of the cloud to be nearly unchanged with [C/O]. However, some dependency on [O/H]
essentially CO-dark. Although [C I] (1–0) shows a good correlation can be seen. Even more remarkably, comparing the results between
with the N(H2 ) map in O1 C2 D0 -15-2, it also becomes very weak in the the second and third column (bottom row) of Fig. 10, the W(CO 1–0)
O1 C2 D1 -15-2 simulation, making the molecular structure essentially versus N(H2 ) relation does not depend on the normalized dust-to-gas
[C I]-dark. ratio, D.
In both these models, [C II] 158 μm appears to give the most rea- Fig. 11 shows the aforementioned correlation of emission lines
sonable and detailed correlation with the molecular cloud structure, with the H2 column densities for ζCR = 10−15 s−1 . The models with
thus favouring the observation of this line in low-metallicity galaxies χ /χ 0 = 102 are also plotted with dotted lines. With the exception
with high star-formation rates (Cormier et al. 2019; Madden et al. of these higher-FUV simulations, the general behaviour of [C II]

MNRAS 527, 8886–8906 (2024)

8898 T. G. Bisbas et al.

Downloaded from [Link] by guest on 21 May 2026

Figure 9. As in Fig. 7 but with a higher FUV intensity of χ /χ 0 = 102 . In the
O0 C1 D0 -15-2 model, CO (1–0) is very bright and well associated with N(H2 )
(see Fig. 8) whereas in the O1 C2 D0 -15-2 and O1 C2 D1 -15-2, the filament is
bright in [C II] 158 μm. In the O1 C2 D1 -15-2 simulation, the molecular cloud
is dark in both [C I](1–0) and CO(1–0). Note the change in the extent of the
colour-bar for the [C I](1–0) line between this figure and Fig. 7.

it is affected mostly when a higher radiation field is applied. It is

interesting to note that for [O/H] = −1.21 and regardless of the
dust-to-gas ratio, W(CO 1–0) builds in an approximately linear
fashion. However, as the corresponding panels in Fig. 11 show, it
originates from the very dense and clumpy medium leading to a
CO-dark molecular gas (see also Figs 7 and 9).
Given that the W(CO 1–0)–N(H2 ) correlation can be considered
independent on the [C/O] ratio for fixed [O/H], it is expected that the
corresponding CO-to-H2 conversion factor will behave in a similar
fashion. This finding also favours the alternative methodology of
Seifried et al. (2020) in calculating the H2 gas mass content for solar
metallicity as their formulas (see their equations 12–14) may hold
also for [C/O] < 0.

Figure 8. As in Fig. 6 but with a higher FUV intensity of χ /χ 0 = 102 . In the

O1 C2 D0 -15 and O1 C2 D1 -15-2 models, N(CO) is severely suppressed even
3.6 Distribution of C+ -, C-, and CO-rich gas masses across
in the densest part of the molecular filament. various galactic environments
How do the gas masses of C+ , C, and CO distribute in each chemical
158 μm, [C I] (1–0), and CO J = 1–0 is similar to the one described phase of the cloud? To answer this question, a definition of ‘pure
above for ζCR = 10−17 s−1 . The models with χ /χ 0 = 102 show two atomic’, ‘PDR’, and ‘molecular’ regions needs to be introduced,
important differences in the W(C II) and W(CO 1–0) emission when based on the relative abundances of H I and H2 . It is defined, here,
compared to the fiducial models of lower FUV intensity. Since the as ‘pure atomic’ (or simply ‘atomic’), the region that satisfies the
higher radiation photodissociates CO and increases the abundance of relation xH I ≥ 0.9, as ‘PDR’ the region that satisfies the relation xH I
C+ , an increase in W(C II) at lower N(H2 ) is observed. At the same < 0.9 and xH2 < 0.49, and as ‘molecular’ the region with xH2 ≥ 0.49.
time, W(CO 1–0) builds for much higher N(H2 ) in the O1 C2 D0 -15- The choice of the latter criterion is taken after consideration of the
2 and O1 C2 D1 -15-2 models. A similar behaviour, albeit at a lower H I narrow-line self-absorption observations of Galactic infrared dark
scale, can be also observed for W(C I 1–0). clouds (Zuo et al. 2018), as the early stage of high-column density,
Like above with ζCR = 10−17 s−1 , it appears that CO J = 1–0 does H2 -dominated clouds. However, a more relaxed criterion satisfying
not depend on the [C/O] value for any [O/H] and D considered. This the relation of xH2 ≥ 0.45 is also explored (see Appendix B).
means that the W(CO 1–0) versus N(H2 ) relation is not a function of
[C/O] neither for the low ζ CR nor for the high ζ CR rate considered.
3.6.1 Distribution of C+
We find that this is because the optical depth of CO (1–0) decreases
with decreasing [C/O] resulting in a weak change in the W(CO 1– The upper panel of Fig. 12 shows the percentage fraction of MC+ .
0)–N(H2 ) relation as a function of [C/O]. However, it seems that Blue colour corresponds to the molecular region, orange to the

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8899

Downloaded from [Link] by guest on 21 May 2026

Figure 10. Mean values of the velocity integrated emission (in units of K km s−1 ) versus N(H2 ) for ζCR = 10−17 s−1 . From top-to-bottom the emission of
[C II], [C I](1–0), and CO(1–0) are shown for all simulations considered, in units of K km s−1 . Solid lines show the [C/O] = −0.07 case whereas dashed and
dot–dashed the [C/O] = −0.57 and −0.77, respectively. The first column is for [O/H] = −0.21, the second is for [O/H] = −1.21 and D = 3 × 10−2 and the
third for [O/H] = −1.21 and D = 3 × 10−3 .

Figure 11. As in Fig. 10 for ζCR = 10−15 s−1 . The dotted lines are for the models with χ /χ 0 = 102 . The W(CO 1–0) versus N(H2 ) relationship does not depend
on the [C/O] ratio.

MNRAS 527, 8886–8906 (2024)

8900 T. G. Bisbas et al.
In all other cases, a trace amount of C+ mass is located in the
molecular gas. Instead, in such models, MC+ is located in PDRs at a
fraction of ∼50–60 per cent with the remaining to be associated
with the atomic phase. It is noted, here, that our simulations
do not model the ionized phase; should that have occurred, the
fraction corresponding to the ‘pure atomic’ would also include MC+
associated with that phase. For instance, the SHINING collaboration
(Herrera-Camus et al. 2018) find that 60–90 per cent of [C II] may
arise from neutral gas including gas of low-ionization.

3.6.2 Distribution of C
The middle panel of Fig. 12 shows the percentage fraction of MC
using the above criterion. It is found that MC is associated either with

Downloaded from [Link] by guest on 21 May 2026

the molecular phase or the PDR depending on the environmental
condition. The atomic phase contains trace amount (1 per cent)
of C-rich gas. More specifically, for ζCR = 10−17 s−1 and a linear
decrease of D (models O0 C0, 1 D0 -17, and O1 C0, 1, 2 D0 -17), a fraction
of MC exceeding 90 per cent of the total C is associated with the
molecular phase. The corresponding fraction for the PDR phase is
10 per cent. Once a sublinear D is adopted (models O1 C0, 1, 2 D1 -
17), the photodissociation of CO at higher optical depths decreases
the fraction of MC associated with the molecular phase, giving an
approximate fifty-fifty per cent between this phase and the PDR
phase.
For the high cosmic-ray ionization rate of ζCR = 10−15 s−1 , it is
found that for metallicities close to solar (model O0 C0 D0 -15), the
majority of MC originates from the molecular phase with a fraction
of 75 per cent. For [C/O] < 0 at this metallicity, this fraction can
be reduced down to approximately 60 per cent (model 1A 15). For
the low-metallicity models (O1 C0, 1, 2 D0 -15), including the linear and
sublinear decrease of the dust-to-gas ratio D, C is abundant on the
surface of PDRs whereas the molecular gas is very C-poor. It can
be therefore argued that in such systems the [C I] emission may be
associated with the PDR surface rather than directly with the H2 rich
gas.
Higher FUV intensities appear to increase the abundance of C
inside the molecular gas. This is evident from the O0 C1 D0 -15-
2, O1 C2 D0 -15-2, and O1 C2 D1 -15-2 models when compared to the
corresponding ones of lower FUV intensity. In particular, it is found
that the two orders of magnitude higher FUV intensity considered,
results in an increase of MC that is associated with the molecular
phase by ∼25 per cent for solar metallicities and by 10–50 per cent
for low-metallicities, the latter depending on the adopted value of D.

Figure 12. Distribution of C+ mass (top), C mass (middle), and CO mass

(bottom) in each model. Green colour refers to the atomic gas (χ H I ≥ 0.9), 3.6.3 Distribution of CO
orange colour to PDRs (χ H I < 0.9 and χ H2 < 0.49), and blue colour to
molecular gas (χ H2 ≥ 0.49). The numbers shown represent the percentage of The bottom panel of Fig. 12 shows the corresponding mass distribu-
the carbon mass associated with each different phase. Fractions less than 5 tion of CO. Overall, CO is associated with molecular gas to a very
per cent are not labelled. See Section 3.6 for discussion. high percentage, as expected from PDR theory. However, for the
models of low [O/H], low D and high ζ CR , it is found that part of
PDR and green to the pure atomic. The general picture is that only the total CO gas mass may be considerably connected with PDRs
two models (O0 C0, 1 D0 -17) appear to contain a significant fraction to 40 per cent. This connection increases from ∼40 to ∼60 per
of C+ mass in the molecular phase of approximately 60 per cent. cent with decreasing [C/O] (models O1 C0, 1, 2 D1 -15). By relaxing the
For low metallicities at low ζ CR (O1 C0, 1, 2 D0 -17), it is found that criterion for defining PDRs and molecular regions as introduced in
approximately 20 per cent of the total C+ mass is associated with the the beginning of this section (i.e. from xH2 ≥ 0.49 to xH2 ≥ 0.45,
molecular gas phase. A similar-to-higher percentage (20–30 per cent) see Appendix B), we find that the aforementioned percentages are
is also found when the higher ζCR = 10−15 s−1 is applied in the two reduced to 15 per cent. This abrupt decrease implies that CO-rich
models of solar metallicities (so for O0 C0, 1 D0 -15) and depending on gas may exist at the very edge of PDRs and that the transition from
the [C/O] ratio. Such a fraction is found also when the FUV intensity C+ and/or C to CO is extremely sharp. For such extreme cases,
is increased (O0 C1 D0 -15-2). fraction of the CO emission may originate from such regions.

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8901
3.7 Model limitations evolving giant molecular cloud for typical conditions met in the
MW. In their work, they find that carbon is subthermally excited,
The present models are to investigate the trends in the emission
implying that the observed [C I] emission seen in real clouds should
lines and the abundances of the carbon cycle as well as the location
underestimate the derived column density of that species. This
of the atomic-to-molecular transition resulting from [C/O] < 0 and
underestimation should arise from the fact that subthermal excitation
sublinear dust-to-gas relations as a function of metallicity ([O/H]). It
reduces the brightness of the line. Observations of [C I] lines in high-
is expected that the dynamical evolution will differ under the various
redshift strongly lensed galaxies by Harrington et al. (2021) indicate
ISM environmental parameters explored (e.g. Gong et al. 2020;
that they are generally subthermally excited. Similar, results were
Hu, Sternberg & van Dishoeck 2021), affecting the aforementioned
concluded also by Papadopoulos, Dunne & Maddox (2022) who
abundances and line emissions. However, it is not expected that the
studied a sample of 106 galaxies. It should be therefore interesting to
presented trends will be altered if dynamical evolution is included
see how the [C I](1–0) emission line behaves in the presented work
(see Offner et al. 2013; Bisbas et al. 2023), and thus an understanding
and whether there is a relationship with the [C/O] ratio.
of how negative [C/O] ratios impact the observables can be reached.
To calculate the amount of subthermally excited [C I](1–0) emis-
Numerical codes incorporating (magneto-)hydrodynamics coupled
sion in our models, an approach similar to Glover et al. (2015) is
with chemistry (such as Glover et al. 2010; Walch et al. 2015; Gong,

Downloaded from [Link] by guest on 21 May 2026

adopted. Through the level populations outputted by 3D-PDR, the
Ostriker & Wolfire 2017; Smith et al. 2020; Hu, Sternberg & van
excitation temperature, Tex , is calculated locally in each cell using
Dishoeck 2021) will be needed to provide a more accurate insight
the expression
into the presented outcomes.
  −1
hν10 s1
Tex = − ln , (6)
kB 3s0
4 DISCUSSION
where ν 10 is the [C I](1–0) frequency, s0 and s1 , the level populations
The observed chemical evolution of galaxies shows that carbon in the ground and first levels, respectively, and the factor three results
and oxygen have been originated from different stellar processes, from the statistical weight ratio between these levels. The column of
resulting in negative [C/O] ratios. While negative ratios have been the excitation temperature (Tex, w ) is then calculated along each line-
considered in studies related to star-formation and the IMF, their of-sight, weighted by the local carbon number density. Similarly, the
impact on the observables using the carbon cycle emission lines TC, w gas temperature weighted by nC is calculated alongside (see
and the location of the H I-to-H2 transition has not been addressed Glover et al. 2015, for more details).
in the past. Through three-dimensional PDR and radiative transfer The map of the TC, w /Tex, w ratio can be then constructed (see
modelling, it has been demonstrated how such [C/O] values change top panels of Fig. 13). For values of TC, w /Tex, w  1, the excitation
the emission of [C II] 158 μm, [C I] (1–0) and (2–1), and the first temperature of [C I](1–0) is approximately equal to the local gas tem-
ten CO transitions in a star-forming cloud, which is embedded in perature. In this situation, local thermodynamical equilibrium (LTE)
different cosmic-ray ionization rates and FUV intensities. is reached and the line becomes thermalized. This occurs in the parts
Concerning the H I-to-H2 transition, for a given [O/H] metallicity of the cloud containing densities that are higher than the critical one of
negative [C/O] ratios do not shift its location, and as a result, the [C I](1–0) (approximately 103 cm−3 for cold gas, see Papadopoulos,
molecular mass content does not change. Instead, this transition Thi & Viti 2004; Glover et al. 2015). On the other hand, for values
is more sensitive to the ISM environmental parameters concerning TC, w /Tex, w > 1, the column of the excitation temperature is lower
the value of [O/H], the cosmic-ray ionization rate and the FUV than the one for the gas temperature implying that radiative de-
intensity as was previously studied (e.g. Bisbas, Tan & Tanaka excitation is more rapid than collisional de-excitation LTE is then
2021). Sublinear dust-to-gas ratios tend to increase the H I abundance not satisfied and the [C I](1–0) line is subthermally excited. The top
therefore shifting the H I-to-H2 transition at higher H-nucleus column panels of Fig. 13 show the aforementioned behaviour in two models;
densities. the O0 C0 D0 -17 (left-hand panel) and O1 C2 D1 -15 (right-hand panel).
Since carbon is a major coolant in PDRs in ISM environments, In the left one, TC, w /Tex, w  1.5 in most of the cloud (except for
where [C/O] < 0, the average (density-weighted) gas temperature two high-density regions), implying that the majority of [C I](1–
is found to increase. The effect is better seen in the diffuse gas 0) emission is subthermally excited (albeit it is weak for the ISM
which can increase Tgas as much as two times. However, for ISM conditions of O0 C0 D0 -17 as discussed in the previous section). In
environments with sublinear dust-to-gas ratios, FUV photons create the right one, corresponding to low-metallicity gas ([O/H] = −1.21)
more C+ , which is an efficient coolant. Although the photoelectric with [C/O] = −0.77, sublinear dust-to-gas ratio and high ζ CR value,
effect operates for a wider range of column densities, the additional the majority of [C I](1–0) emission is very thermalized.
C+ decreases the overall gas temperature (e.g. models O1 C2 D0, 1 - While there is no precise condition defining the transition from
15-2). As discussed earlier, the photoelectric heating also decreases thermal-to-subthermal emission, we adopt here the ratio of 1.3 to
(especially in lower density regions) as a result of the lower D values. define it, meaning a 30 per cent difference between the gas and the
In turn, this has consequences in the line emission of all C-cycle excitation temperature. With this criterion, the percentage emission
coolants. of subthermally excited [C I](1–0) is calculated for all models and
Our discussion below is mainly focused on the impact of [C/O] < presented in the bottom panel of Fig. 13. In general, it is found
0 on observations of the [C I](1–0) emission. The connection of that subthermal emission of [C I] occurs only for metal-rich clouds
[C/O] < 0 with ‘[C II]-deficit’ ISM is also examined. (models with [O/H] = −0.21) regardless to the ISM conditions (and
except for O0 C0 D0 -15). For low-metallicity conditions, [C I](1–0)
emission is generally thermalized. There is an interesting pattern
4.1 Subthermal excitation of [C I](1–0)
seen, however. For fixed [O/H], the percentage of subthermal
Numerical models studying the atomic carbon emission line have emission of the carbon line is increasing with decreasing [C/O], albeit
been the main focus of various groups in the last few years. Glover not substantially for low-metallicity gas. This trend is marked with
et al. (2015) examined the [C I] fine structure lines of a dynamically blue arrows, separated for each group of models. Considering that

MNRAS 527, 8886–8906 (2024)

8902 T. G. Bisbas et al.
NGC 6052 (Mrk 297, Arp 209) is a low-redshift (z = 0.01581)
system of colliding galaxies with an SFR  18 M yr−1 (Michiyama
et al. 2020). The latter group reported nondetection of [C I](1–0) and
a bright emission of the CO J = 4–3, resulting in a luminosity ratio
of LCI 1−0 /LCO 4−3 < 0.08. In their analysis, Michiyama et al. (2020)
used PDR models by assuming MW abundances and they found that
such a ratio can be reproduced if the observed number densities are
nH > 105 cm−3 , which may be reasonably explained from the strong
shocks produced due to the merging process. However, our models
can reproduce the observed ratio if a low C/O ratio is considered.
More specifically, NGC 6052 is observed to be metal rich, with 12
+ log (O/H) = 8.22–8.85 (Sage, Loose & Salzer 1993; James et al.
2002; Shi et al. 2005; Rupke, Veilleux & Baker 2008) corresponding
to [O/H] = −0.47 to +0.16, which matches, on average, with O0

Downloaded from [Link] by guest on 21 May 2026

models. Models O0 C1 D0 -15 and O0 C1 D0 -15-2 contain a strong
emission of CO J = 4–3, and stronger than the O0 C0 D0 -15 of higher
C/O by a factor of 2, as indicated in the bottom panel of Fig. 3 (see
red solid line) and a weak [C I](1–0) emission (see red filled circle),
which is 10 × fainter than O0 C0 D0 -15. Similarly, the abundance of
C is very suppressed for O0 C1 D0 -15 (Fig. 2 bottom panel). The ratio
of these two lines for O0 C2 D0 -15 and O0 C1 D0 -15-2 (see Fig. A1)
is close to the observed luminosity ratio.7 The remaining question
is whether this galaxy has a [C/O] ratio that is much lower than
that of the MW. An insight on that can be taken by Sage, Loose &
Salzer (1993) who report a high ratio of 12 CO/13 CO > 22 and Sage,
Mauersberger & Henkel (1991) who report a high abundance of
18
O in local starburst galaxies. Both of these works claim that such
a ratio may indicate a higher population of massive stars (such as
Figure 13. Top: Maps of the TC, w /Tex, w ratio for the O0 C0 D0 -17 (left)
discussed in Zhang et al. 2018) than in the MW. If 12 C is partly
and O1 C2 D1 -15 (right) models. For ratios ∼1, [C I](1–0) is thermalized,
whereas for >1 it is subthermally excited. Bottom: Percentage emission
produced in intermediate and low-mass stars (see Section 1.1), it
of subthermally excited [C I](1–0) line in each model. For fixed [O/H], the may in turn produce low C/O ratios. Due to the ongoing collision, an
percentage is increasing with decreasing [C/O]. The arrows point these trends increase in the star-formation rate leads to an increase in the cosmic-
in each group of models. ray ionization rate and FUV intensities (e.g. Papadopoulos 2010;
Kashiyama & Mészáros 2014; Lahén et al. 2020), consistent with
the presented models.
the atomic carbon line ratio does not strongly depend on the [C/O]
NGC 7679 (Mrk 534) is a low-redshift (z = 0.00177) Seyfert 2
ratio (see Fig. A1), it is expected that the [C I](2–1) line emission
galaxy whose metallicity has not yet been reasonably estimated. Its
will also follow the behaviour of [C I](1–0) in regards to whether it
star-formation rate is SFR ∼ 10–21 M yr−1 (De Looze et al. 2014;
is subthermally excited or not.
Davies et al. 2016). Recently, Michiyama et al. (2021) reported a
Therefore, it can be demonstrated that for solar metallicities
low [C I](1–0)/CO(4–3) luminosity ratio of LCI /LCO < 0.08, similar
the [C I] emission is generally subthermally excited, in agreement
to NGC 6052. They further reported young starbursts, a finding that
with Papadopoulos, Dunne & Maddox (2022). However, in low-
supports the hypothesis that a larger population of massive stars may
metallicity star-forming clouds and galaxies, [C I](1–0) emission is
exist. NGC 7679 is undergoing interaction, possibly with NGC 7682
generally found to be thermalized. Its intensity is weak and it is a
(Michiyama et al. 2021). Through SNe type II, these conditions
natural consequence of negative [C/O] ratios in the observed ISM
may create α-enhanced environments resulting in low C/O ratios.
rather than being associated with subthermal emission. In turn, this
Like with NGC 6052 discussed above, our models show that an
behaviour is directly related with the galactic chemical evolution.
environment reminiscent to that of O0 C1 D0 -15 and O0 C1 D0 -15-2
models may explain the very weak detection of [C I] (1–0) line.
In studying the kinematics of star-forming galaxies at cosmic-
4.2 Three cases of [C I]-dark galaxies
noon, Lelli et al. (2023) presented ALMA observations of the zC-
Observations of galaxies that are bright in CO but dark in [C I] 400569 galaxy, located at z = 2.24, in mid-J and [C I](1–0) lines.
lines have been reported in the past (Michiyama et al. 2020, 2021; This solar-metallicity galaxy has an elevated SFR = 81 M yr−1 and
Harrington et al. 2021; Dunne et al. 2022; Lelli et al. 2023). it is, thus, expected it will contain higher cosmic-ray energy densities
Explanations for the low LCI /LCO luminosity ratios include the and possibly also high average FUV radiation field strengths when
existence of extended high-density regions that are not bright in compared to the MW. While these conditions favor the [C I](1–0)
[C I] lines, CI-poor environments, missing fluxes, and even specific emission (Bisbas, Tan & Tanaka 2021), Lelli et al. (2023) reported a
processes followed during data analysis. However, the possibility of weak emission of this line which did not allow to perform kinematic
low C/O ratios has not been examined as a potential cause for [C I]- studies. It remains to see whether zC-400569 has [C/O] < 0, as this
dark galaxies. Our attention is drawn for three candidates that may
contain [C/O] < 0 gas: the low-redshift NGC 6052 and NGC 7679
galaxies, which are discussed in Michiyama et al. (2020, 2021), and 7 See Solomon & Vanden Bout (2005) for the equality between the luminosity

the high-redshift zC-400569 discussed in Lelli et al. (2023). ratio and the velocity integrated emission ratio.

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8903
would be in accordance to the results of the present work. Should this intensities, are predicted to appear bright in [C II] 158 μm and CO(1–
be true and considering also the relevance with NGC 6052 and NGC 0) but dark in [C I](1–0). Similarly, low-metallicity dwarf galaxies
7679, it can be demonstrated that galaxies undergoing starburst may are predicted to be in principle [C II] 158 μm bright with very low
become [C I]-dark as a consequence of negative [C/O] ratios resulting emission in CO(1–0) arising from dense clumpy structures only.
from larger populations of massive stars (see, however, Sharda et al. In these galaxies, [C I](1–0) is also expected to be thermalized but
2023a, for potential -but localized- differential metal removal through very weak. Low C/O abundance ratios contribute insignificantly in
galactic winds). making the ISM gas [C II] deficit. It is also expected that O-bearing
molecules will be affected both in terms of their relative abundances
and line emission.
4.3 Contribution of [C/O] < 0 on the [C II]-deficit gas
The impact of negative [C/O] ratios on the aforementioned
Observations examining the [C II]/FIR (far-IR) ratio have shown that abundances and line emission appears to be a complicated and
it has a declining trend with increasing the FIR luminosity (Malhotra non-linear problem. For a better understanding, especially on the
et al. 2001; Luhman et al. 2003; Casey, Narayanan & Cooray 2014; large-scale ISM studies, full hydro-chemical models incorporating
Zhao, Yan & Tsai 2016). This trend is known in the literature as dynamical and chemical evolution simultaneously are needed. It is

Downloaded from [Link] by guest on 21 May 2026

‘[C II] deficit’. Many studies have focused on examining its origin also recommended for observational researchers to report, whenever
˜ & Oh 2016; Narayanan & Krumholz 2017; Smith et al.
(Munoz possible, the C/O abundance ratio, as this is of immense importance
2017; Rybak et al. 2019; Bisbas et al. 2022) and have suggested for performing robust PDR models.
various different physical and chemical processes. The role played
by a [C/O] < 0 gas on the [C II] deficit has not been extensively
studied in the past, however some first insight for its contribution AC K N OW L E D G E M E N T S
was discussed in Harikane et al. (2020), who found that low carbon- The authors thank the anonymous referee for their comments which
to-oxygen ratios are not sufficient to explain the observed trends. It improved the clarity of the work. The authors thank Piyush Sharda
should be therefore interesting to explore whether this claim holds for the discussion and suggestions. ZYZ, EG, GL, YS and ZG
for our set of models. acknowledge support from the National Natural Science Foundation
Assuming optically thin emission for both [C II] and dust (FIR) of China (NSFC) under grants Nos. 12173016, 12041305, the
and following Bisbas et al. (2022) for their calculation, we find in Program for Innovative Talents, Entrepreneur in Jiangsu, and the
our models that low C/O ratios can contribute ࣠3 times to the [C II] science research grants from the China Manned Space Project with
deficit (i.e. to decrease the [C II]/FIR ratio by up to approximately Nos. CMS-CSST-2021-A08 and CMS-CSST-2021-A07. YHZ is
three times). This is a small factor compared to the 1–2 orders grateful for support from the NSFC (Grant No. 12173079). DQ and
of magnitude necessary to consider it as a major factor. It is thus XJJ acknowledge support from the NSFC (Grant No. 12373026).
demonstrated that low C/O ratios play a minor and hence insignificant
role in contributing to the [C II] deficit, in agreement with the findings
of Harikane et al. (2020). Although we do not model the ionized DATA AVA I L A B I L I T Y
phase as mentioned also earlier, it is expected that the above result The data underlying this article will be shared on reasonable request
will also hold even when H II regions are present, since cosmological to the corresponding author.
hydrodynamical simulations of a dwarf galaxy merger (Bisbas et al.
2022) found that that they may account on average ∼10 per cent of
the total C+ luminosity. REFERENCES
Accurso G., Saintonge A., Bisbas T. G., Viti S., 2017, MNRAS, 464, 3315
5 CONCLUSIONS Akerman C. J., Carigi L., Nissen P. E., Pettini M., Asplund M., 2004, A&A,
414, 931
This work explored the impact of negative [C/O] ratios, implying Aller L. H., Greenstein J. L., 1960, ApJS, 5, 139
a subsolar scaling between carbon and oxygen for low-metallicity Amarsi A. M., Nissen P. E., Skúladóttir Á., 2019, A&A, 630, A104
gas, on the location of the atomic-to-molecular transition and on Appleton P. N. et al., 2013, ApJ, 777, 66
the abundances and line emission of the carbon cycle in a star- Arata S., Yajima H., Nagamine K., Abe M., Khochfar S., 2020, MNRAS,
forming cloud. The 3D-PDR astrochemical simulations revealed that 498, 5541
Arellano-Córdova K. Z. et al., 2022, ApJ, 940, L23
consideration of the relative carbon-to-oxygen abundances in both
Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, ARA&A, 47, 481
metal-rich and especially metal-poor ISM conditions is of highly im- Asplund M., Amarsi A. M., Grevesse N., 2021, A&A, 653, A141
portance. The conclusion is that the [C/O] ratio cannot be neglected Bate M. R., 2014, MNRAS, 442, 285
from astrochemical calculations concerning the diffuse and dense Bate M. R., 2019, MNRAS, 484, 2341
ISM. For studies related to (extra-)galactic clouds, the [C/O] ratio Bell T. A., Viti S., Williams D. A., 2007, MNRAS, 378, 983
constitutes an additional and equally important ISM environmental Berg D. A., Skillman E. D., Henry R. B. C., Erb D. K., Carigi L., 2016, ApJ,
parameter on top of the three basic and most commonly used (the 827, 126
cosmic-ray ionization rate, the FUV intensity, and the metallicity in Berg D. A., Erb D. K., Henry R. B. C., Skillman E. D., McQuinn K. B. W.,
the form of the [O/H] ratio) in astrochemical modelling. It should be 2019, ApJ, 874, 93
Bialy S., Sternberg A., Lee M.-Y., Le Petit F., Roueff E., 2015, ApJ, 809, 122
noted that for low-metallicity clouds and galaxies, special attention
Bisbas T. G., Bell T. A., Viti S., Yates J., Barlow M. J., 2012, MNRAS, 427,
must be paid also to the dust-to-gas mass ratio.
2100
The use of [C I](1–0) emission as H2 gas mass tracer in galaxies Bisbas T. G., Papadopoulos P. P., Viti S., 2015, ApJ, 803, 37
with [C/O] < 0 is difficult. For such a medium, it appears that Bisbas T. G., van Dishoeck E. F., Papadopoulos P. P., Szűcs L., Bialy S.,
the [C I] lines may become very weak even when high cosmic-ray Zhang Z.-Y., 2017a, ApJ, 839, 90
energy densities are present. Metal-rich star-forming galaxies of low Bisbas T. G., Tanaka K. E. I., Tan J. C., Wu B., Nakamura F., 2017b, ApJ,
[C/O] ratios, which are expected to contain enhanced ζ CR and FUV 850, 23

MNRAS 527, 8886–8906 (2024)

8904 T. G. Bisbas et al.
Bisbas T. G. et al., 2018, MNRAS, 478, L54 Glover S. C. O., Clark P. C., Micic M., Molina F., 2015, MNRAS, 448, 1607
Bisbas T. G., Schruba A., van Dishoeck E. F., 2019, MNRAS, 485, 3097 Goldsmith P. F., Langer W. D., Pineda J. L., Velusamy T., 2012, ApJS, 203,
Bisbas T. G., Tan J. C., Tanaka K. E. I., 2021, MNRAS, 502, 2701 13
Bisbas T. G. et al., 2022, ApJ, 934, 115 Gong M., Ostriker E. C., Wolfire M. G., 2017, ApJ, 843, 38
Bisbas T. G., van Dishoeck E. F., Hu C.-Y., Schruba A., 2023, MNRAS, 519, Gong M., Ostriker E. C., Kim C.-G., Kim J.-G., 2020, ApJ, 903, 142
729 Grenier I. A., Black J. H., Strong A. W., 2015, ARA&A, 53, 199
Black J. H., Dalgarno A., 1977, ApJS, 34, 405 Gullberg B. et al., 2016, A&A, 591, A73
Bolatto A. D., Wolfire M., Leroy A. K., 2013, ARA&A, 51, 207 Harikane Y. et al., 2020, ApJ, 896, 93
Bonifacio P. et al., 2015, A&A, 579, A28 Harrington K. C. et al., 2021, ApJ, 908, 95
Bothwell M. S. et al., 2017, MNRAS, 466, 2825 Haworth T. J., Glover S. C. O., Koepferl C. M., Bisbas T. G., Dale J. E., 2018,
Brauher J. R., Dale D. A., Helou G., 2008, ApJS, 178, 280 New A Rev., 82, 1
Burbidge E. M., Burbidge G. R., Fowler W. A., Hoyle F., 1957, Rev. Mod. Herrera-Camus R. et al., 2012, ApJ, 752, 112
Phys., 29, 547 Herrera-Camus R. et al., 2015, ApJ, 800, 1
Burkhart B., Lee M.-Y., Murray C. E., Stanimirović S., 2015, ApJ, 811, L28 Herrera-Camus R. et al., 2018, ApJ, 861, 94
Cardelli J. A., Meyer D. M., Jura M., Savage B. D., 1996, ApJ, 467, 334 Hollenbach D. J., Tielens A. G. G. M., 1999, Rev. Mod. Phys., 71, 173

Downloaded from [Link] by guest on 21 May 2026

Cartledge S. I. B., Lauroesch J. T., Meyer D. M., Sofia U. J., 2004, ApJ, 613, Hoyle F., 1954, ApJS, 1, 121
1037 Hu C.-Y., Sternberg A., van Dishoeck E. F., 2021, ApJ, 920, 44
Casey C. M., Narayanan D., Cooray A., 2014, Phys. Rep., 541, 45 Indriolo N., 2023, ApJ, 950, 64
Cazaux S., Tielens A. G. G. M., 2002, ApJ, 575, L29 Indriolo N., Bergin E. A., Falgarone E., Godard B., Zwaan M. A., Neufeld
Cazaux S., Tielens A. G. G. M., 2004, ApJ, 604, 222 D. A., Wolfire M. G., 2018, ApJ, 865, 127
Cazaux S., Tielens A. G. G. M., 2010, ApJ, 715, 698 Isobe Y. et al., 2023, preprint (arXiv:2307.00710)
Cescutti G., Matteucci F., 2022, Universe, 8, 173 Israel F. P., Baas F., 2002, A&A, 383, 82
Combes F., 2018, A&A Rev., 26, 5 Izotov Y. I., Schaerer D., Worseck G., Berg D., Chisholm J., Ravindranath
Cooke R. J., Pettini M., Steidel C. C., 2017, MNRAS, 467, 802 S., Thuan T. X., 2023, MNRAS
Cormier D. et al., 2019, A&A, 626, A23 James A., Dunne L., Eales S., Edmunds M. G., 2002, MNRAS, 335, 753
Cosentino G. et al., 2019, ApJ, 881, L42 James T. A., Viti S., Holdship J., Jiménez-Serra I., 2020, A&A, 634, A17
Davies R. L. et al., 2016, MNRAS, 462, 1616 Jiao Q., Zhao Y., Zhu M., Lu N., Gao Y., Zhang Z.-Y., 2017, ApJ, 840, L18
De Looze I., Baes M., Bendo G. J., Cortese L., Fritz J., 2011, MNRAS, 416, Jiao Q. et al., 2019, ApJ, 880, 133
2712 Jones T. et al., 2023, ApJ, 951, L17
De Looze I. et al., 2014, A&A, 568, A62 Jura M., 1974, ApJ, 191, 375
Delgado Mena E., Adibekyan V., Santos N. C., Tsantaki M., González Kajino T., Aoki W., Balantekin A. B., Diehl R., Famiano M. A., Mathews G.
Hernández J. I., Sousa S. G., Bertrán de Lis S., 2021, A&A, 655, A99 J., 2019, Part. Nucl. Phys., 107, 109
Draine B. T., 1978, ApJS, 36, 595 Kasen D., Metzger B., Barnes J., Quataert E., Ramirez-Ruiz E., 2017, Nature,
Draine B. T., 2011, Physics of the Interstellar and Intergalactic Medium. 551, 80
Princeton Univ. Press, Princeton Kashiyama K., Mészáros P., 2014, ApJ, 790, L14
Draine B. T., McKee C. F., 1993, ARA&A, 31, 373 Katz H. et al., 2022, MNRAS, 510, 5603
Dufour G., Charnley S. B., 2021, ApJ, 909, 171 Kaufman M. J., Wolfire M. G., Hollenbach D. J., Luhman M. L., 1999, ApJ,
Dunne L., Maddox S. J., Papadopoulos P. P., Ivison R. J., Gomez H. L., 2022, 527, 795
MNRAS, 517, 962 Kelly G., Viti S., Garcı́a-Burillo S., Fuente A., Usero A., Krips M., Neri R.,
Dwek E., 1998, ApJ, 501, 643 2017, A&A, 597, A11
Esposito F., Vallini L., Pozzi F., Casasola V., Mingozzi M., Vignali C., Klitsch A. et al., 2022, MNRAS, 514, 2346
Gruppioni C., Salvestrini F., 2022, MNRAS, 512, 686 Kreckel K. et al., 2019, ApJ, 887, 80
Esteban C., Méndez-Delgado J. E., Garcı́a-Rojas J., Arellano-Córdova K. Z., Krips M. et al., 2016, A&A, 592, L3
2022, ApJ, 931, 92 Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S.,
Fabbian D., Nissen P. E., Asplund M., Pettini M., Akerman C., 2009, A&A, Steinwandel U. P., Moster B. P., 2020, ApJ, 891, 2
500, 1143 Le Bourlot J., Pineau des Forets G., Roueff E., Schilke P., 1993, ApJ, 416,
Ferland G. J. et al., 2017, RMxAA, 53, 385 L87
Franeck A. et al., 2018, MNRAS, 481, 4277 Le Petit F., Nehmé C., Le Bourlot J., Roueff E., 2006, ApJS, 164, 506
Freer M., Fynbo H. O. U., 2014, Prog. Part. Nucl. Phys., 78, 1 Lee H. H., Herbst E., Pineau des Forets G., Roueff E., Le Bourlot J., 1996,
Frias Castillo M. et al., 2023, ApJ, 945, 128 A&A, 311, 690
Fudamoto Y. et al., 2023, preprint (arXiv:2309.02493) Lelli F. et al., 2023, A&A, 672, A106
Fuente A. et al., 2019, A&A, 624, A105 Li X., Millar T. J., Heays A. N., Walsh C., van Dishoeck E. F., Cherchneff I.,
Gaches B. A. L., Offner S. S. R., Bisbas T. G., 2019a, ApJ, 878, 105 2016, A&A, 588, A4
Gaches B. A. L., Offner S. S. R., Bisbas T. G., 2019b, ApJ, 883, 190 Liang L. et al., 2023, preprint (arXiv:2301.04149)
Gaches B. A. L., Bisbas T. G., Bialy S., 2022a, A&A, 658, A151 Limongi M., Chieffi A., 2018, ApJS, 237, 13
Gaches B. A. L., Bialy S., Bisbas T. G., Padovani M., Seifried D., Walch S., Lique F., Werfelli G., Halvick P., Stoecklin T., Faure A., Wiesenfeld L.,
2022b, A&A, 664, A150 Dagdigian P. J., 2013, J. Chem. Phys., 138, 204314
Galametz M., Madden S. C., Galliano F., Hony S., Bendo G. J., Sauvage M., Lo N. et al., 2014, ApJ, 797, L17
2011, A&A, 532, A56 Luhman M. L., Satyapal S., Fischer J., Wolfire M. G., Sturm E., Dudley C.
Gallino R., Arlandini C., Busso M., Lugaro M., Travaglio C., Straniero O., C., Lutz D., Genzel R., 2003, ApJ, 594, 758
Chieffi A., Limongi M., 1998, ApJ, 497, 388 Luo G. et al., 2020, ApJ, 889, L4
Garnett D. R., Skillman E. D., Dufour R. J., Peimbert M., Torres-Peimbert Luo G. et al., 2023, ApJ, 942, 101
S., Terlevich R., Terlevich E., Shields G. A., 1995, ApJ, 443, 64 Mackey J., Walch S., Seifried D., Glover S. C. O., Wünsch R., Aharonian F.,
Genzel R. et al., 2012, ApJ, 746, 69 2019, MNRAS, 486, 1094
Girichidis P. et al., 2016, MNRAS, 456, 3432 Madden S. C. et al., 2020, A&A, 643, A141
Glover S. C. O., Clark P. C., 2016, MNRAS, 456, 3596 Maiolino R., Mannucci F., 2019, A&A Rev., 27, 3
Glover S. C. O., Federrath C., Mac Low M. M., Klessen R. S., 2010, MNRAS, Malhotra S. et al., 2001, ApJ, 561, 766
404, 2 Maloney P. R., Hollenbach D. J., Tielens A. G. G. M., 1996, ApJ, 466, 561

MNRAS 527, 8886–8906 (2024)

 α-enhanced ISM 8905
Mashian N. et al., 2015, ApJ, 802, 81 Stacey G. J., Geis N., Genzel R., Lugten J. B., Poglitsch A., Sternberg A.,
Mathis J. S., 1990, ARA&A, 28, 37 Townes C. H., 1991, ApJ, 373, 423
Mathis J. S., Rumpl W., Nordsieck K. H., 1977, ApJ, 217, 425 Stanley F. et al., 2023, ApJ, 945, 24
Matsunaga N. et al., 2023, ApJ, 954, 198 Sternberg A., Le Petit F., Roueff E., Le Bourlot J., 2014, ApJ, 790, 10
McElroy D., Walsh C., Markwick A. J., Cordiner M. A., Smith K., Millar T. Strong A. W., Moskalenko I. V., Ptuskin V. S., 2007, Annu. Rev. Nucl. Part.
J., 2013, A&A, 550, A36 Sci., 57, 285
Meijerink R., Spaans M., Israel F. P., 2006, ApJ, 650, L103 Suárez-Andrés L., Israelian G., González Hernández J. I., Adibekyan V. Z.,
Meijerink R., Spaans M., Israel F. P., 2007, A&A, 461, 793 Delgado Mena E., Santos N. C., Sousa S. G., 2018, A&A, 614, A84
Meijerink R., Spaans M., Loenen A. F., van der Werf P. P., 2011, A&A, 525, Sutter J. et al., 2019, ApJ, 886, 60
A119 Tacconi L. J. et al., 2008, ApJ, 680, 246
Michiyama T. et al., 2020, ApJ, 897, L19 Tanaka K. E. I., Omukai K., 2014, MNRAS, 439, 1884
Michiyama T. et al., 2021, ApJS, 257, 28 Tanaka K. E. I., Tan J. C., Zhang Y., Hosokawa T., 2018, ApJ, 861, 68
Montoya Arroyave I. et al., 2023, A&A, 673, A13 Tielens A. G. G. M., 2005, The Physics and Chemistry of the Interstellar
˜ J. A., Oh S. P., 2016, MNRAS, 463, 2085
Munoz Medium. Cambridge Univ. Press, Cambridge
Narayanan D., Krumholz M. R., 2014, MNRAS, 442, 1411 Trainor R. F., Strom A. L., Steidel C. C., Rudie G. C., 2016, ApJ, 832, 171

Downloaded from [Link] by guest on 21 May 2026

Narayanan D., Krumholz M. R., 2017, MNRAS, 467, 50 Valentino F. et al., 2018, ApJ, 869, 27
Narayanan D., Krumholz M. R., Ostriker E. C., Hernquist L., 2012, MNRAS, Valentino F. et al., 2020, ApJ, 890, 24
421, 3127 Vallini L., Pallottini A., Ferrara A., Gallerani S., Sobacchi E., Behrens C.,
Nicholls D. C., Sutherland R. S., Dopita M. A., Kewley L. J., Groves B. A., 2018, MNRAS, 473, 271
2017, MNRAS, 466, 4403 Vallini L., Tielens A. G. G. M., Pallottini A., Gallerani S., Gruppioni C.,
Öberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16 Carniani S., Pozzi F., Talia M., 2019, MNRAS, 490, 4502
Oesch P. A. et al., 2016, ApJ, 819, 129 van Dishoeck E. F., Black J. H., 1986, ApJS, 62, 109
Offner S. S. R., Bisbas T. G., Viti S., Bell T. A., 2013, ApJ, 770, 49 van Dishoeck E. F., Black J. H., 1988, ApJ, 334, 771
Offner S. S. R., Bisbas T. G., Bell T. A., Viti S., 2014, MNRAS, 440, van Dishoeck E. F. et al., 2023, Faraday Discussions, 245, 52
L81 Viti S., Roueff E., Hartquist T. W., Pineau des Foretsˆ G., Williams D. A.,
Pagel B. E. J., 2009, Nucleosynthesis and Chemical Evolution of Galaxies. 2001, A&A, 370, 557
Cambridge Univ. Press, Cambridge Walch S. et al., 2015, MNRAS, 454, 238
Papadopoulos P. P., 2010, ApJ, 720, 226 Wolfire M. G., McKee C. F., Hollenbach D., Tielens A. G. G. M., 2003, ApJ,
Papadopoulos P. P., Thi W. F., Viti S., 2004, MNRAS, 351, 147 587, 278
Papadopoulos P. P., van der Werf P., Isaak K., Xilouris E. M., 2010, ApJ, 715, Wolfire M. G., Tielens A. G. G. M., Hollenbach D., Kaufman M. J., 2008,
775 ApJ, 680, 384
Papadopoulos P. P., Bisbas T. G., Zhang Z.-Y., 2018, MNRAS, 478, 1716 Wolfire M. G., Vallini L., Chevance M., 2022, ARA&A, 60, 247
Papadopoulos P., Dunne L., Maddox S., 2022, MNRAS, 510, 725 Wu B., Tan J. C., Nakamura F., Van Loo S., Christie D., Collins D., 2017,
Pietrow A. G. M., Hoppe R., Bergemann M., Calvo F., A&A, 2023, 672, L6 ApJ, 835, 137
Pineda J. L., Langer W. D., Goldsmith P. F., 2014, A&A, 570, A121 Wuyts E. et al., 2016, ApJ, 827, 74
Ravindranath S., Monroe T., Jaskot A., Ferguson H. C., Tumlinson J., 2020, Xie T., Allen M., Langer W. D., 1995, ApJ, 440, 674
ApJ, 896, 170 Yang C. et al., 2023, A&A, 680, A95
Rémy-Ruyer A. et al., 2013, A&A, 557, A95 Zanella A. et al., 2018, MNRAS, 481, 1976
Rémy-Ruyer A. et al., 2014, A&A, 563, A31 Zhang Z.-Y. et al., 2014, A&A, 568, A122
Richings A. J., Schaye J., 2016, MNRAS, 458, 270 Zhang Z.-Y., Romano D., Ivison R. J., Papadopoulos P. P., Matteucci F., 2018,
Röllig M., Ossenkopf-Okada V., 2022, A&A, 664, A67 Nature, 558, 260
Röllig M. et al., 2007, A&A, 467, 187 Zhao Y., Yan L., Tsai C.-W., 2016, ApJ, 824, 146
Romano D., 2022, A&A Rev., 30, 7 Zuo P., Li D., Peek J. E. G., Chang Q., Zhang X., Chapman N., Goldsmith P.
Romano D., Franchini M., Grisoni V., Spitoni E., Matteucci F., Morossi C., F., Zhang Z.-Y., 2018, ApJ, 867, 13
2020, A&A, 639, A37
Rosenberg M. J. F. et al., 2015, ApJ, 801, 72
Roueff E., Le Bourlot J., 2020, A&A, 643, A121 A P P E N D I X A : L I N E R AT I O S
Rupke D. S. N., Veilleux S., Baker A. J., 2008, ApJ, 674, 172 Line ratios are frequently used in observations to study the ISM
Rybak M. et al., 2019, ApJ, 876, 112
conditions and constrain the environmental parameters in MW
Sage L. J., Mauersberger R., Henkel C., 1991, A&A, 249, 31
clouds and distant extragalactic objects (e.g. Israel & Baas 2002;
Sage L. J., Loose H. H., Salzer J. J., 1993, A&A, 273, 6
Sandstrom K. M. et al., 2013, ApJ, 777, 5 Gullberg et al. 2016; Krips et al. 2016; Valentino et al. 2018, 2020;
Schneider N. et al., 2023, Nat. Astron., 7, 546 Papadopoulos, Dunne & Maddox 2022). Here, we explore the ratios
Seifried D., Haid S., Walch S., Borchert E. M. A., Bisbas T. G., 2020, of CO(1–0)/[C I](1–0) the lines of which are both important coolants
MNRAS, 492, 1465 of the dense gas, CO(4-3)/[C I](1–0) and CO(7-6)/[C I](2-1) as the
Sharda P., Krumholz M. R., 2022, MNRAS, 509, 1959 frequency separation between these lines is small and are often
Sharda P., Ginzburg O., Krumholz M. R., Forbes J. C., Wisnioski E., Mingozzi observed together especially with ALMA, and the atomic carbon
M., Zovaro H. R. M., Dekel A., 2023a, preprint (arXiv:2303.15853) line ratio, [C I](1–0)/[C I](2-1). Fig. A1 shows the emission line
Sharda P., Amarsi A. M., Grasha K., Krumholz M. R., Yong D., Chiaki G., ratios of the aforementioned combinations of CO and [C I] for the
Roy A., Nordlander T., 2023b, MNRAS, 518, 3985
environmental parameters explored.
Sharda P., Amarsi A. M., Grasha K., Krumholz M. R., Yong D., Chiaki G.,
It is found that for ζCR = 10−17 s−1 , all above line ratios remain
Roy A., Nordlander T., 2023c, MNRAS, 525, 3316
Shi F., Kong X., Li C., Cheng F. Z., 2005, A&A, 437, 849 approximately constant as a function of [C/O], except for the CO(1–
Smith J. D. T. et al., 2017, ApJ, 834, 5 0)/[C I](1–0) and CO(4-3)/[C I](1–0) for the O0 C1 D0 -17 model in
Smith R. J. et al., 2020, MNRAS, 492, 1594 which they show a slight increase when compared to the O0 C0 D0 -17.
Solomon P. M., Vanden Bout P. A., 2005, ARA&A, 43, 677 The most interesting features, however, occur for ζCR = 10−15 s−1 .
Spitoni E., Silva Aguirre V., Matteucci F., Calura F., Grisoni V., 2019, A&A, Here, it is found that all line ratios, except for the atomic carbon line
623, A60 ratio of [C I](1–0)/[C I](2–1), are increased with decreasing [C/O]

MNRAS 527, 8886–8906 (2024)

8906 T. G. Bisbas et al.

Downloaded from [Link] by guest on 21 May 2026

Figure A1. Emission line ratios between CO and [C I] for all explored
simulation cases and that are commonly used in observations.

for fixed [O/H]. The increase in FUV radiation dissociates CO and

ionizes C, therefore affecting their emission. The imprint of that can
be seen in the second column of Fig. A1 (dashed lines), where a
decreasing trend is observed between the high and low FUV models.
The atomic carbon line ratio remains overall remarkably constant
as a function of [C/O] in all cases, even when a higher FUV radiation
field is applied. However, this ratio increases with the increasing ζ CR
as can be seen when comparing it between the left and the right
columns of the above figure. This is in accordance with the findings
of Bisbas, Tan & Tanaka (2021); Bisbas et al. (2023), who suggest
that the atomic carbon line ratio can be used to infer the cosmic-ray
ionization rate, particularly in extragalactic star-forming systems, as
it is not strongly affected by other changes in the ISM environmental Figure B1. As in Fig. 12 but with the molecular gas defined as the gas with
parameters. xH I ≥ 0.45, and the PDR with xH I < 0.9 and xH2 < 0.45, accordingly.

more prominent, however, for the C species, as in all environmental

A P P E N D I X B : M A S S D I S T R I B U T I O N O F C +, C ,
parameters explored and except for the lower D at ζCR = 10−15 s−1 ,
AND CO WITH MOLECULAR GAS DEFINED AS
it is associated with >60 or even >90 per cent with the H2 -rich gas.
χH2 ≥ 0 . 4 5
For CO, the bottom panel of Fig. B1 indicates that CO is always
As discussed, the χ H2 ≥ 0.49 condition for defining the molecular associated with the molecular region as defined with xH2 ≥ 0.45. It
phase was based on observations of Galactic infrared dark clouds can be, therefore, argued that the transition layer between C and CO
(Zuo et al. 2018). It would be interesting, however, to explore how is sharp and very close to the border of H2 -rich gas. A similar picture
the above results change if this condition is relaxed to χ H2 ≥ 0.45. occurs for the C+ to C transition, albeit at lower column densities.
Fig. B1 shows how MC+ , MC , and MCO are now associated with the
PDR and molecular gas phases. For C+ and C, it is found that the
contribution from the molecular gas phase is increased. The effect is This paper has been typeset from a TEX/LATEX file prepared by the author.

© 2023 The Author(s).

Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License
([Link] which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
MNRAS 527, 8886–8906 (2024)