Hydrodynamic simulations of the Disc of Gas Around Supermassive black holes (HDGAS) – II. The transition from neutral atomic to molecular gas phases

ABSTRACT

We use HDGAS (Hydrodynamic simulations of the Disc of Gas Around Supermassive black holes) hydrodynamic simulations to study the impact of active galactic nucleus (AGN) feedback on the conversion of atomic gas to molecular gas within the circumnuclear disc of a typical AGN-dominated galaxy. The comparison of C i, C ii, and CO line intensities and their ratios in the HDGAS post-processing radiative transfer analysis reveals the complex interplay between AGN activity, cold molecular gas properties, and the physical processes governing the evolution of star formation in galaxies. Our results demonstrate that the C i/CO intensity ratio serves as a reliable indicator of the atomic-to-molecular gas transition. We present the probability distribution function and abundance trends of various metal species related to molecular H₂ gas, highlighting differences in clumpiness and intensity maps between AGN-feedback and NoAGN models. The profile of the integrated intensity (moment-0) maps shows that the AGN-feedback model exhibits a lower C i/CO intensity ratio in the vicinity of the supermassive black hole (<50 pc), indicating a smaller atomic gas abundance and the presence of positive AGN feedback. Our simulations have successfully predicted the presence of faint-CO emissions extending to larger radii from the galactic centre. We also explore the relationships between C ii/CO and C i/C ii intensity ratios, as well as the ratios versus CO intensity, which provides insights into the ‘CO-dark’ issues. One notable feature in the later time-scale of the AGN model is the presence of a ‘CO-dark’ region, where the intensity of CO emission (⁠$I_{\mathrm{CO}}$⁠) is depleted relative to the H${}_2$ column density (⁠$N_{{\mathrm{H}}_2}$⁠) compared to the NoAGN model.

1 INTRODUCTION

Molecular gas is both a product and a driver of star formation (e.g. Kennicutt & Evans 2012). It forms from atomic hydrogen in cold and shielded regions of the interstellar medium (ISM), where it can collapse and fragment to form stars. It impacts the star formation process by determining the physical and chemical conditions of the ISM, which can be examined through the observation of different atomic and molecular gas tracers. These tracers are employed to investigate the galactic disc and circumnuclear disc (CND) of galaxies (Kamenetzky et al. 2011; Combes et al. 2014; Takano et al. 2014; Viti et al. 2014; Salvestrini et al. 2022).

Atoms can serve as tracers of relatively dense molecular gas regions in different galaxies. Examples of such atoms include neutral atomic carbon (C i) and singly ionized carbon (C ii; Papadopoulos, Thi & Viti 2004; Bell, Viti & Williams 2007; Glover et al. 2015; Bisbas et al. 2017). The impact of cosmic rays (CRs) on the abundance and distribution of CO, C i, and C ii within molecular clouds of star-forming galaxies has been the subject of recent investigations. Glover et al. (2015) explored the utility of C i emission in turbulent molecular clouds for investigating their structure and dynamics. Through numerical simulations of chemical, thermal, and dynamical evolution, the authors modelled C i emission under various cloud properties and observational conditions. The findings indicated widespread C i emission within the cloud due to turbulence-induced density variations, enabling radiation to penetrate deeply. The C i emission was demonstrated to accurately trace the column density of the cloud across a broad range of visual extinctions, surpassing the reliability of CO for low-extinction regions. The study further revealed that C i excitation temperatures generally remained lower than the kinetic temperatures, suggesting that the carbon atoms were not in thermal equilibrium with the gas. The authors discussed techniques for estimating C i excitation temperatures and atomic carbon column density based on C i line observations, while considering potential estimation errors. Other studies revealed that elevated CR ionization rates can lead to a substantial reduction in CO abundance, thereby complicating the tracing of ${\mathrm{H}}_2$ mass and cloud structure. Additionally, C i and C ii were found to become more prevalent and dominant at higher CR ionization rates (Bisbas et al. 2017; Izumi et al. 2018).

Saito et al. (2022) compare the C i/CO intensity ratios in different regions of NGC 1068, a nearby starburst/active galactic nucleus (AGN) galaxy, and find that they vary from 0.1 to 1.0, with higher values in the outflow regions and lower values in the disc regions. The study interprets these variations as the result of different excitation and optical depth effects, as well as different chemical abundances and destruction rates of C i and CO. They suggest that C i is more abundant and less affected by photodissociation than CO in the outflow regions, where the radiation field and the gas density are higher. They also suggest that C i is more optically thin than CO in the outflow regions, where the gas column density is lower. A kinematic model indicated that the C i enhancement was predominantly driven by the interaction between the outflow and the disc, which compressed and heated the gas. This implies that the AGN feedback has a negative effect on star formation in the CND and several recent studies indicate that the CND of NGC 1068 is a region dominated by X-ray radiation (XDR; Usero et al. 2004; Krips et al. 2008; Pérez-Beaupuits et al. 2009; García-Burillo et al. 2010).

Another example is the nearby Seyfert-2 type galaxy, Circinus, which hosts one of the closest AGNs. Izumi et al. (2018) revealed that the C i(1–0)/CO(3–2) ratio at the AGN position was unusually high, suggesting an XDR-type chemistry (Izumi et al. 2018; Miyamoto et al. 2018). This means that the AGN radiation can destroy molecules and create atomic gas, which can affect the star formation efficiency and rate (Park et al. 2022). The kinematics of both lines showed the multiphase nature in the torus region; i.e. the diffuse atomic gas was more spatially extended along the vertical direction of the disc than the dense molecular gas. This was explained by a radiation-driven fountain model, which proposed that atomic outflows were driven by radiation pressure from the AGN and created a geometrically thick atomic disc. This supported the validity of the radiation-driven fountain scheme in explaining the physical origin of the AGN torus (Uzuo et al. 2021).

CO-dark molecular gas represents a significant component of the ISM in galaxies, particularly in regions with intense radiation fields such as those near supermassive black holes (SMBHs) in the central nuclear region of galaxies (Montoya Arroyave et al. 2023). In Seyfert-2 type AGN-dominated galaxies, the presence of an SMBH can influence the surrounding gas, potentially increasing the proportion of CO-dark gas. As shown by Madden et al. (2020), this gas can be effectively traced by C ii and C i emission lines, which highlight the reservoirs of molecular hydrogen (H${}_2$⁠) that CO fails to detect. The C ii 158 micron line, in particular, is bright and can trace the total H${}_2$ mass, revealing up to 70–100 per cent of the H${}_2$ mass that is not traced by CO in dwarf galaxies. Conversely, in galaxies without AGN, the star formation histories and the properties of the molecular gas can differ.

AGN feedback encompasses the various processes through which the active nucleus of a galaxy interacts with its surrounding gas (Croton et al. 2006; Booth & Schaye 2009; Dubois et al. 2013; Bower et al. 2017; Raouf et al. 2017, 2019; Raouf, Purabbas & Fazel Hesar 2024), influencing its physical and chemical properties. In this study, our focus is on investigating the influence of AGN feedback on the transition from atomic to molecular gas within a CND surrounding an SMBH. Our investigation builds upon recent hydrodynamics simulation studies, HDGAS (Hydrodynamic simulations of the Disc of Gas Around Supermassive black holes; Raouf et al. 2023), to further explore the effects of AGN feedback on the CND scale. The transition from neutral atomic to a molecular gas phase is especially interesting, as it reflects the balance between the various processes that form and destroy molecules, such as heating, cooling, shielding, and dissociation. The transition also affects the observability and detectability of the CND, as different gas phases emit or absorb different types of radiation, such as radio, infrared, or X-rays (Izumi et al. 2018).

Our analysis concentrates on the atomic and molecular species encompassed within the chimes chemistry network, including C i, C ii, and CO. The structure of the paper is as follows: the simulation and methods are described in Section 2. In Section 3, we introduce the observational galaxy templates used for comparison. Section 4 presents our results and analysis, while Section 5 offers a summary and discussion.

2 THE hdgas SIMULATIONS

In our study, we utilize hydrodynamic simulations called HDGAS to explore the ISM within the CND of AGN-dominated galaxies, considering the influence of mechanical feedback from the AGN. These simulations incorporate the chimes non-equilibrium chemistry network to accurately represent radiative cooling and AGN heating. Focusing on Seyfert-2 type galaxies, such as NGC 1068, we develop a model of a gas disc surrounding a black hole (BH). By comparing this model to one without AGN feedback, we investigate the effects of AGN feedback on the central region of the galactic disc. For a comprehensive understanding of the HDGAS hydrodynamic simulation model, we refer interested readers to Raouf et al. (2023) (referred to as Paper-I).

The simulations were run using gizmo (Hopkins 2015), a hydrodynamics and gravity code that implements several different hydrodynamics solvers. The initial conditions consist of gas-rich nuclear disc containing a BH in its centre with $n={10}^6$ gas particles, where each particle initially has a mass of $\sim 100{\mathrm{M}}_{\odot }$⁠. This simulation is coupled with stellar feedback and ISM physics from Feedback In Realistic Environments -2 (FIRE-2) as described in Hopkins et al. (2018). The formation of stars is only possible in cold, molecular, and locally self-gravitating regions that have a number density above $n_{\mathrm{H}}={10}^4{\mathrm{c}\mathrm{m}}^{-3}$ (as used in the most recent FIRE simulation studies; e.g. Torrey et al. 2020). With the chimes non-equilibrium chemistry and cooling model (Richings, Schaye & Oppenheimer 2014a, b), chemical abundances for 157 species are calculated over time, including all ionization states of 11 elements that are important for cooling, as well as for 20 molecules, including CO, HCO${}^{+}$⁠, and ${\mathrm{H}}_2\mathrm{O}$⁠. By integrating the temperature in time with the 157 rate equations, the chimes module calculates cooling and heating rates.

We use the BH accretion model implemented by Hopkins & Quataert (2011). In this approach, BH accretion is assumed to be determined by the gravitational torques (see Hopkins et al. 2016, for more details of accretion methods). In this study, AGNs exhibit blueshifted Broad Absorption Line (BAL) or troughs when their line of sight is intercepted by a high-speed outflow, likely originated from their accretion discs (which is unresolved in the simulation; see Hopkins & Quataert 2011). Through the mechanical feedback process we implement in this simulation, wind mass and kinetic luminosity are continuously ‘injected’ into the gas surrounding the SMBH where the outflow is isotropic (Hopkins et al. 2016). It is assumed that some fraction of the photon momentum drives a wind at the resolution scale around the BH (Murray & Chiang 1995). Hence, the accreted gas is blown out as a wind with velocity $v_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}$⁠. The wind is defined by two parameters, the mass loading of $\beta \equiv {\dot{M}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}/{\dot{M}}_{\mathrm{B}\mathrm{H}}$ and the velocity $v_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}$ that equivalently relate to momentum loading (⁠${\dot{p}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}={\eta }_{p}L/c$⁠) and energy loading (⁠${\dot{E}}_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}={\eta }_{E}L$⁠) of the wind. In this study, our simulation has been carried out using fixed energy and momentum loading factor for AGN model (with value for $\beta =6$ and $v_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}=5000\mathrm{km}{\mathrm{s}}^{-1}$ based on observations and theoretical models; see e.g. Moe et al. 2009; Dunn et al. 2010; Hamann et al. 2011; Borguet et al. 2013). In the NoAGN model, we have set both $\beta$ and $v_{\mathrm{w}\mathrm{i}\mathrm{n}\mathrm{d}}$ to zero.

The maps of emission lines are generated through the post-processing of simulation snapshots using version 0.40 of the publicly available Monte Carlo radiative transfer code radmc-3d (Dullemond et al. 2012), utilizing the abundances of ions and molecules calculated during the simulations with the chimes chemistry module. All line images are produced with a passband ranging from −100 to +100 km s⁻¹ and 400 wavelength points evenly distributed around the line’s centre. In this study, we generated emission maps for C i, C ii, and CO. Specifically, we utilized the first J-line of atomic carbon (C i) at a frequency of 492 GHz, which is crucial for probing the warm, dense gas regions typical of AGN environments. We also examined ionized carbon (C ii) at 1900 GHz, a key fine-structure line in the far-infrared spectrum. Additionally, we analysed carbon monoxide [CO(2–1)] at 230 GHz, which is commonly used to trace molecular gas, essential for star formation. Using NGC 1068 as an example, we use the distance of 14 Mpc (so that 70 pc equals 1 arcsec in each moment map) and inclination of 41°.

3 GALAXY TEMPLATES

We have used submillimetre Atacama Large Millimeter/submillimeter Array (ALMA) observations of a selection of galaxies to compare our simulated maps. We specifically focused on four galaxies¹: NGC 1808, NGC 7469, NGC 3627, and NGC 4321. Note that all the selected galaxies are nearby and comparable in distance to NGC 1068.

NGC 1808 is a galaxy exhibiting intense starburst activity. Observations of C i(1–0) and CO(2–1) have revealed a high column density of atomic carbon in the CND (Salak et al. 2019). Furthermore, the C i/CO intensity ratios suggest that the luminosity of C i can serve as a reliable tracer of molecular gas mass, particularly in resolved starburst nuclei (Salak et al. 2019). In the case of NGC 7469, a type 1 Seyfert galaxy, the circumnuclear gas disc (CND) exhibits a ring-like structure and a two-arm/bi-symmetric spiral pattern. Surrounding the CND is a starbursting ring. By studying the emissions of C i(1–0) and CO(1–0), Nguyen et al. (2021) measure the mass of the SMBH via the C i emission, which seems to be a better mass probe than the CO. NGC 3627, also known as M66, is a galaxy that forms part of the Leo Triplet. ALMA observations (Liu et al. 2023) have shown nearly uniform C i/CO line ratios across the majority of its star-forming disc. There is also an excellent spatial correspondence between the C i and CO emissions. The abundance ratio of [C i/CO] in NGC 3627 has been estimated to be approximately 0.1, consistent with previous large-scale studies conducted on the Milky Way (Liu et al. 2023). Similarly, observations of NGC 4321, also referred to as M100, within its inner discs have revealed nearly uniform line ratios of C i and CO emissions across its star-forming discs. Liu et al. (2023) have observed a mild decrease in the C i/CO ratio with increasing metallicity. The abundance ratio of [C i/CO] in NGC 4321 has also been estimated to be around 0.1, indicating a consistent chemistry with NGC 3627.

Note that these observations generally have a lower resolution (with an average spatial resolution greater than 100 pc) than our simulations, which may affect our interpretations. However, our aim is to highlight the overall trend of our results in relation to the observations, rather than to focus on the specific details of comparison.

4 ANALYSIS AND RESULTS

4.1 Atomic and molecular gas abundances

In this study, we analyse the distributions and probability density functions (PDFs) of C i, C ii, and CO for different regions of the CND with radii lower than 50 and 100 pc. The results are shown in Fig. 1 for $t=7$ Myr, which is the time-scale during which we can confirm the effects of the AGN. This snapshot is optimal since the gas depletion time-scale is under 1 Myr (García-Burillo et al. 2014), making a simulation duration of 10 Myr ideal. By normalizing the PDFs to the total hydrogen number density (⁠$n_{\mathrm{H}}$⁠), we can directly compare the variations in molecular and atomic gas densities, providing an understanding of their respective contributions and the overall gas distribution within the system. Although the density peaks in the PDFs do not show significant differences between the AGN and NoAGN models, the skewness and tails in the high- and low-density regions are crucial indicators of AGN influences. In the region $r<50\text{pc}$⁠, the peak of the logarithm of CO abundance relative to hydrogen number density $\left(n_{\text{CO}}/n_{\text{H}}\right)$ is $-3.7$ for both the AGN and NoAGN models. This similarity suggests that the conditions leading to CO formation are comparable in both scenarios, despite the presence of an AGN in one model. However, as we look at the carbon species, we see a divergence: for C i, the logarithmic abundance peaks at $-5.5$ for NoAGN, compared to $-5$ for AGN. This indicates that the presence of an AGN potentially enhances the C i abundance, possibly due to increased heating or photodissociation processes that favour C i formation. For C ii, the peaks are even lower, with values of $-8.3$ for NoAGN and $-8$ for AGN. This further highlights the influence of the AGN, which likely provides additional energy that facilitates the ionization processes necessary for C ii formation. The differences in C i and C ii suggest that the AGN environment promotes a higher degree of ionization, affecting the chemical pathways of carbon species. In the region $r<100\text{pc}$⁠, the CO PDF peaks at $-4$ for NoAGN, while the PDF peak is at $-3.5$ for the AGN model. This shift indicates that the AGN environment may contribute to a higher CO abundance, perhaps through shock heating or enhanced star formation activity associated with the AGN. Furthermore, the peaks of the PDF for C i in both models are $-5.3$ for NoAGN and $-5$ for AGN, reinforcing the idea that the AGN fosters conditions that are more conducive to C i formation. The differences in PDF peaks reflect the varying chemical environments influenced by the AGN presence. Finally, the peaks of the PDF in C ii for both the AGN and NoAGN models are $-7.9$⁠, suggesting that while the AGN impacts the other carbon species, the C ii formation might reach a saturation point in the conditions present, resulting in similar C ii levels across both models.

Furthermore, we examine the skewness² of the distributions in Fig. 2 to assess the asymmetry between the AGN and NoAGN models for time-scales between 3 and 8 Myr. A positive skewness value (⁠${\delta }_s>0$⁠) indicates a longer tail on the right side of the distribution, suggesting higher density regions. Conversely, a negative skewness value (⁠${\delta }_s<0$⁠) corresponds to a longer tail on the left side, representing lower density regions.

Fig. 1 illustrates a bimodal distribution of C i density (purple lines) in the PDFs for both the AGN and NoAGN models across various radii. The density exhibits fluctuations in low-density areas (early time; t < 5 Myr) and high-density areas (late time; t > 5 Myr), suggesting an increase in C i density over time for the AGN model (see Fig. 2; blue lines). In the CO skewness [${\delta }_s(\text{AGN-NoAGN})$], the temporal trend indicates higher density tails for the AGN model at both 50 and 100 pc scales of the disc. Specifically, at t = 7 Myr (shown in Fig. 1), the AGN model is indeed located in lower density regions compared to the NoAGN model. However, we also emphasize that the overall trend from t = 3 to t = 8 Myr reveals distinct behaviours for the AGN and NoAGN models, especially in regions greater than 50 pc. This distinction highlights the influence of AGN feedback on the gas distribution. For C ii, the deviations in the PDFs are notably smaller compared to those for C i. As illustrated in Fig. 2, C ii predominantly exists in the lower density mode. The AGN model shows a significant change in skewness, transitioning from ${\delta }_s(\text{AGN},t=5\text{ Myr})=-1.54$ to ${\delta }_s(\text{AGN},t=6\text{ Myr})=0.84$⁠. In contrast, the NoAGN model exhibits a smaller change, with skewness values of ${\delta }_s(\text{NoAGN},t=5\text{ Myr})=-1.10$ and ${\delta }_s(\text{NoAGN},t=6\text{ Myr})=-1.03$⁠. This shift in skewness for the AGN model indicates a transition to higher density regions at 6 Myr, particularly in the inner regions (⁠$<$50 pc).

Furthermore, we investigate the correlations between C i, C ii, and CO abundances normalized to ${\mathrm{H}}_2$ within the CND. The results are visualized in Fig. 3, shedding light on the distinct effects of AGN feedback on the molecular composition of CO at different radii. In the top panels of Fig. 3, we show the C i/H${}_2$ ratio against the CO/H${}_2$ ratio. The data points exhibit a correlation at low CO abundances (CO/H${}_2<{10}^{-6}$⁠). However, an intriguing anticorrelation trend (between C i and CO) emerges at high CO/H${}_2$ ratios (⁠$>{10}^{-6}$⁠) for both the AGN and NoAGN models. Further, the AGN model displays significantly higher C i/H${}_2$ ratios compared to the NoAGN model, particularly within radii less than 50 pc. This suggests that AGN feedback is particularly effective in enhancing the abundance of C i in the inner regions of the disc, likely due to increased radiation fields, heating of the ISM, gas mixing from outflows, and the formation of extensive photodissociation regions (PDRs). In the bottom panel of Fig. 3, we examine the C ii/H${}_2$ ratio as a function of the CO/H${}_2$ ratio. We see a steeper anticorrelation for CO/H${}_2>{10}^{-6}$ compared to C i/H${}_2$⁠. Also, we see that at lower CO abundance levels, the abundance of C ii remains relatively constant and the AGN model exhibits higher levels of ionized C ii gas abundances at different time-scales (most of the times). The intense photoionization or dissociation processes induced by the AGN likely contribute to this effect.

Moreover, the higher CO/H${}_2$ ratios observed in the NoAGN model suggest that the AGN is influencing these ratios, particularly in both regions beyond 100 and 50 pc from the disc, where we do not observe $\log ({\text{CO/H}}_2)$ values exceeding approximately −3 in the AGN model. This indicates that the presence of the AGN leads to a suppression of the CO/H${}_2$ ratio in these areas, causing the solid lines to terminate at lower CO/H${}_2$ values compared to the dashed lines.

4.2 Integrated intensity maps and profiles

Fig. 4 shows the spatial distribution of C i, C ii, and CO integrated intensities, as well as their ratios (C i/CO, C ii/CO, and C i/C ii) in the CND for both the AGN and NoAGN models. The integrated intensity maps can reveal structures such as molecular spiral arms, rings, or clumps, which may be associated with ongoing star formation or other phenomena related to the interaction between the AGN and the surrounding gas. The C i/CO ratio maps (t = 7 Myr) generally illustrate the extent of CO photodissociation and/or C ii recombination, highlighting their effects on the molecular gas abundance in the PDRs. Higher C i/CO ratios indicate decreased CO molecule levels due to photodissociation by AGN radiation.

In the NoAGN model, the high-intensity regions within 50 pc correspond to areas with a greater concentration of molecular gas, likely linked to active star-forming regions or dense molecular clouds. These regions exhibit stronger C i emission due to the higher abundance of carbon-bearing molecules undergoing CO photodissociation. Conversely, the AGN model reveals low-intensity regions within the same radius, indicating lower molecular gas densities or reduced C i and CO abundance. These areas may be influenced by AGN radiation, shocks, or other physical processes that affect the excitation or destruction of atomic C i and CO molecules, highlighting that the factors responsible for CO destruction differ from those that ionize C i (see C ii in the second panels of Fig. 4). Fig. 5 illustrates the radial profiles of C i and CO intensity for both the AGN and NoAGN models. The NoAGN model displays higher intensity values for radii below 50 pc compared to the AGN model. In the lower panels of Fig. 4, the C i/CO(2–1) ratio intensity maps reveal characteristics of the photon-dominated region (PDR). In Fig. 5, the NoAGN model displays a higher C i/CO ratio at radii below 50 pc for t  $>$3 Myr, suggesting a larger interface between neutral atomic gas and molecular gas in the absence of AGN influence. Furthermore, while the overall trend does indicate an increase in the C i/CO ratio as a function of radius, the AGN model exhibits a higher C i/CO(2–1) ratio at t = 3 Myr across all radii compared to the NoAGN model and other time points.

At later times (⁠$t=8$ Myr), the AGN model shows a higher C i/CO ratio for radii between 30 and 100 pc, indicating that the presence of an AGN influences the gas dynamics in this region. A high C i/CO ratio reflects the balance between different carbon phases in the ISM. Specifically, a higher ratio suggests a greater presence of C i compared to CO, which may arise from processes such as increased star formation or heating associated with AGN activity. The colour-shaded lines illustrate the observed trends of CO, C i, and C i/CO profiles for samples without AGN, specifically NGC 3627 and NGC 4321, as well as for AGN samples such as NGC 1808 and NGC 7469, within the central 200 pc regions (see Section 3). The AGN samples show higher C i and CO intensities across all regions (limited by resolution) and exhibit a consistent trend in C i/CO ratios similar to simulations at different time-scales beyond 50 pc, although they are generally higher compared to the without AGN samples.

The radial profile of C ii and C ii/CO and C i/C ii integrated intensity maps illustrate in the bottom of Fig. 4. The C ii profile shows distinct differences between the AGN and NoAGN models. For the NoAGN model, the C ii emission is higher in the central region across all time-scales. However, in the outer regions beyond $\sim$50 pc, the C ii profiles for both models drop off inversely at later time-steps, indicating a decrease in C ii emission in the outer galaxy over time. The C ii/CO(2–1) ratio profile also exhibits interesting trends. At earlier time-steps, the AGN model has a higher C ii/CO(2–1) ratio, especially in the central regions. However, at later time-steps, this ratio increases in the outer parts (beyond $\sim$50 pc) for both the AGN and NoAGN models (t = 7, 8 Myr), suggesting a decrease in the relative C ii to CO(2–1) emission in the outer galaxy over time, regardless of the presence of an AGN.

Lastly, the C ii/C i ratio profile reveals that in the earlier time-steps, the AGN model has a lower C ii/C i ratio compared to the NoAGN model, particularly in the central regions. This trend reverses at larger radii (beyond $\sim 30$ pc), where the AGN model exhibits a higher C ii/C i ratio than the NoAGN model at later time-steps (⁠$t=7,8$ Myr). This indicates that the relative amount of C ii to C i emission is initially lower in the AGN model’s central region, but becomes higher in the outer parts of the galaxy at later times.

4.3 The brightnesses in star-forming and starburst disc environments

In Fig. 6, we present scatter points representing the observed brightnesses of C i and CO, along with the C i/CO line ratios colour-coded by distance to the centre of disc (r). A correlation is revealed for different time-scales (t = 3–8 Myr) between C i brightness as a function of CO, aligning with the range of observations from starburst (SB) regions (Salak et al. 2019) and the CNDs (Liu et al. 2023). Notably, C i brightness is higher in AGN models, which is consistent with trends observed in galaxies hosting AGN. In general, the observation trend is confined within the range of $1<\log (I_{\text{CO}})<3$⁠. However, the simulation shows an extension to lower CO values in the range $-1<\log (I_{\text{CO}})<3$ for both the AGN and NoAGN models. The scatter of C i intensity brightness is approximately $-2\text{dex}$⁠, while the C i/CO ratio exhibits a scatter of around $-2.5\text{dex}$⁠. The lower C i/CO ratios are found at greater distances from the centre (⁠$r>200\text{kpc}$⁠). Time evolution leads to a lower C i/CO ratio at larger distances from the disc centre (⁠$r>200\text{kpc}$⁠) for both the AGN and NoAGN models.

As shown in the bottom panels of Fig. 6, the simulations predict higher C i/CO ratios at fainter CO pixels, especially in the outer regions. However, it is important to note that the outer regions of the simulations actually display lower ratios than the observations. The upward tails observed at the CO-faint end in the C i versus CO plot are primarily a result of the observed resolution limits of each galaxy. In contrast, the simulations provide a more accurate representation of this behaviour, reflecting the underlying gas dynamics without the same detection constraints. Additionally, even at earlier time-steps, the simulations indicate that both C i and C i/CO ratios are higher in the AGN model compared to the NoAGN model. Overall, the simulations show both higher and lower C i/CO ratios at fainter CO pixels, contributing to the scatter observed in the data.

Additionally, systematic differences in C i/CO ratios are observed among these galaxies and within their galactic environments. The majority of molecular gas in the 3–7 kpc star-forming discs of NGC 3627 and NGC 4321 exhibits a distribution of C i/CO ratios centred around 0.1 (indicated by the dotted line in the bottom panels of Fig. 6) (Liu et al. 2023). It is important to note that this apparent tightness in the distribution, as shown in Fig. 6, results from the binning process applied to the original data points, which helps to emphasize the underlying trends while smoothing out individual variations.

This distribution remains relatively flat over an order of magnitude in CO(2–1), which is referred to as the star-forming discs’ regime (as described in Fig. 5; Liu et al. 2023). In contrast, the C i/CO ratio increases to about 0.2 within a galactocentric radius of a few hundred parsecs to approximately 1 kpc in the starburst regions of NGC 1808 and NGC 7469, marking the starburst disc regime, as indicated by the dotted line at the identity value of $I_{\mathrm{C}I}$ and $I_{\mathrm{CO}}$⁠.

4.4 The median integrated intensity

The top panel of Fig. 7 displays the median integrated intensity measurements from radii less than 50 and 100 pc from the centre of the disc. These measurements, which include species such as CO, C i, and C ii, were collected over a time span of 3–8 Myr. Overall, the data indicate a declining trend in intensity over time, with this decline being more pronounced at larger radii (<100 pc). This suggests that the intensity of these species decreases as we move further from the centre of the disc in both the AGN and NoAGN models. In both the inner and outer regions, the median intensities of C i, C ii, and CO are lower in the AGN model for time-scales below 7 Myr compared to the NoAGN model. However, in the AGN model, after 7 Myr, the inner region shows an increase in the median intensity of CO, followed by similar increases for the inner and outer regions of C i and C ii at 8 Myr. This suggests that the presence of an AGN influences the distribution and abundance of these atomic species within the disc.

The bottom panel of Fig. 7 illustrates the median integrated intensity ratios obtained from radii smaller than 50 and 100 pc from the centre of the disc for C i/CO, C ii/CO, and C i/C ii in both the AGN and NoAGN models. Overall, the trends reveal a gradual decrease over time for the C i/CO and C i/C ii ratios, indicating a diminishing abundance of C i relative to C ii as the system evolves. Notably, the C i/C ii ratio shows an increase after 6 Myr for radii smaller than 50 pc. Additionally, the C ii/CO ratios exhibit relatively increasing patterns over time. In the inner regions, the AGN models display lower median values for C i/CO, C ii/CO, and C i/C ii for time-scales below 8 Myr compared to the NoAGN model. However, there is an increase in the C i/CO ratio at 8 Myr for the AGN model. As discussed in Section 4.2, the C i/CO ratio is higher in the NoAGN model compared to the AGN model before 7 Myr, suggesting a more prominent presence of C i relative to CO in the absence of AGN activity. This trend reverses at t = 8 Myr in the inner regions (<50 pc). Furthermore, the C ii/CO ratios are predominantly higher for the NoAGN model, indicating a relatively greater abundance of C ii compared to CO when an AGN is not present. Interestingly, the AGN model shows an increasing trend for C i after 6 Myr.

Table 1 summarizes the trends of increases and decreases in the inner and outer regions for various species and ratios in both the AGN and NoAGN models. It categorizes trends in median intensities of C i, C ii, CO, and their ratios by inner regions (⁠$<$50 pc), outer regions (⁠$>$50 pc), and overall trends over time. Changes are indicated as increases (⁠$\uparrow$⁠), decreases (⁠$\downarrow$⁠), or no change (⁠$\to$⁠) for the C i/CO, C ii/CO, and C i/C ii ratios. The table highlights that the C i/CO and C i/C ii ratios generally decline over time, while the C ii/CO ratios tend to increase. The AGN model shows lower median values for C i/CO, C ii/CO, and C i/C ii before 8 Myr, with a reversal in the C i/CO ratio at that time, suggesting that AGN activity influences changes in abundances.

4.5 CO-dark regions

In Fig. 8, the intensity ratios of CO(2–1), CO(1–0), C i, and C ii to the ${\mathrm{H}}_2$ column density (⁠$N_{{\mathrm{H}}_2}$⁠) reveal the relative abundance and behaviour of these carbon-bearing species within the ISM and molecular clouds. The above figure illustrates that in the NoAGN model, during the early time-scale (t < 6 Myr), the column density remains above ${10}^{21}{\text{cm}}^{-3}$⁠, with CO, C i, and C ii decreasing as $N_{{\mathrm{H}}_2}$ increases. Starting at t = 6 Myr, the intensities of CO, C i, and C ii rise with $N_{{\mathrm{H}}_2}$ for column densities below ${10}^{21}{\text{cm}}^{-3}$⁠, with CO and C i showing a steeper increase at later times (t = 7, 8 Myr) before declining for higher column densities. This behaviour suggests that the NoAGN environment initially limits carbon species formation but becomes more favourable as conditions evolve, particularly at lower densities. In contrast, the AGN model exhibits different trends. Here, CO consistently decreases across all time-scales, with a more pronounced decline after t > 7 Myr, indicating that the AGN environment may suppress CO formation more effectively over time. C i also shows a gradual decrease, with a sharp decline at $N_{{\mathrm{H}}_2}$ below ${10}^{21}{\text{cm}}^{-3}$ before stabilizing. Similarly, C ii decreases across all time-scales with increasing $N_{{\mathrm{H}}_2}$⁠, reflecting the AGN’s impact on ionization processes.

In the context of the CND around an AGN, the shape and trends of these ratios can be affected by the AGN’s mechanical feedback. A notable feature in the later stages (t = 7, 8 Myr) of the AGN model is the emergence of a ‘CO-dark’ region, where the intensity of CO emission (⁠$I_{\mathrm{CO}}$⁠) is diminished relative to the H${}_2$ column density compared to the NoAGN model. This CO-dark region likely arises from the AGN’s effects, which can disrupt and dissociate the molecular gas, resulting in a decreased CO abundance relative to other carbon-bearing species, particularly evident at $t=8$ Myr. This disruption mainly impacts lower excitation CO transitions, while higher J CO transitions [e.g. CO(2–1)] may show a less pronounced depletion or even an increase in their intensity ratios, as the mechanical feedback can also heat and excite the remaining molecular gas.

In contrast, as shown in the top two panels of Fig. 8, the $I_{\mathrm{C}I}/N_{{\mathrm{H}}_2}$ and $I_{\mathrm{C}II}/N_{{\mathrm{H}}_2}$ ratios exhibit distinct behaviours. The scatter in the $I_{\mathrm{C}I}/N_{{\mathrm{H}}_2}$ ratios is approximately 2 dex larger than that observed for the $I_{\mathrm{C}II}/N_{{\mathrm{H}}_2}$ ratios across different time-scales, indicating that the C i emission responds more variably to environmental changes. Furthermore, AGN activity becomes increasingly efficient at reducing the C i emission at later time-scales in the $I_{\mathrm{C}I}/N_{{\mathrm{H}}_2}$ ratio compared to the NoAGN model. In contrast, the $I_{\mathrm{C}II}/N_{{\mathrm{H}}_2}$ ratio in the AGN scenario exhibits lower scatter, influenced by local conditions, when compared to the NoAGN model. This could indicate the enhanced abundance of atomic and ionized carbon (CO gets dissociated) compared to molecular hydrogen, as the mechanical feedback from the AGN alters the chemical and physical conditions in the CND (the increase in atomic carbon in CO-dark regions has been shown by various studies; Bisbas, Papadopoulos & Viti 2015; Franeck et al. 2018).

Fig. 9 also shows that the relationships between these CO J-line intensities and the H${}_2$ column density (⁠$N_{{\mathrm{H}}_2}$⁠) can reveal important information about the CO excitation and the impact of AGN feedback. In a typical molecular cloud environment, one would expect to see a certain pattern in the CO J-line intensities as a function of $N_{{\mathrm{H}}_2}$⁠. The lower J CO transitions, such as CO(1–0) and CO(2–1), would typically exhibit a linear or near-linear relationship with $N_{{\mathrm{H}}_2}$⁠, reflecting the relatively uniform CO excitation conditions. However, in the CND influenced by the mechanical feedback from the AGN, the CO J-line intensity versus $N_{{\mathrm{H}}_2}$ relationships deviate from this expected pattern, particularly within the CO-dark region. The CO(1–0) and CO(2–1) line intensities show a depletion relative to the H${}_2$ column density, indicating a disruption and dissociation of the molecular gas by the AGN’s mechanical feedback, which preferentially affects the lower excitation CO transitions. As illustrated in Fig. 9, the CO lines in AGN-affected regions compared to NoAGN regions demonstrate a depletion. This depletion suggests a reduction in the abundance of CO molecules in areas influenced by AGN activity. We will conduct further analysis on the CO intensity ratios in our future study. In contrast, the higher J CO transitions, such as CO(3–2) and CO(4–3), exhibit a less pronounced depletion or even an increase in their intensity versus $N_{{\mathrm{H}}_2}$ relationships. This can be attributed to the AGN feedback selectively heating and exciting the remaining molecular gas, leading to a relative enhancement of the higher J CO lines compared to the lower J transitions.

5 DISCUSSION

In this study, we utilized the HDGAS hydrodynamic simulation (Raouf et al. 2023) to investigate the impact of AGN feedback on the ISM by analysing various atomic to molecular line ratios, including C i/CO, C ii/CO, and C i/C ii, as well as the abundances of different species, to explore the transition from neutral atomic gas to molecular gas in models with and without AGN influence. Comparing the AGN and NoAGN models, we show significant differences in the intensity maps and radial profiles of C i, C ii, and CO and their ratios. The results show that the C i/CO ratio serves as a useful indicator of the transition from atomic to molecular gas.

The probability distribution functions (PDFs) of these species demonstrate the influence of an AGN on the relative abundances within various scales of the CND (50 and 100 pc). The C i PDF indicates a time-dependent increase in C i density specifically for the AGN model. Meanwhile, the CO PDFs reveal that the AGN model exhibits higher density tails at both the 50 and 100 pc scales of the disc, highlighting the AGN’s impact on the distribution of CO as well (as shown in Raouf et al. 2023, regarging positive feedback from AGN). The AGN model exhibits elevated levels of C i and C ii gas abundances across various time-scales. We observe notable differences in carbon species abundance between the AGN and NoAGN models, particularly in the region $r<50\text{pc}$⁠. Both models exhibit a peak logarithm of CO abundance relative to hydrogen density $\left(n_{\text{CO}}/n_{\text{H}}\right)$ at $-3.7$⁠, indicating similar conditions for CO formation. However, the AGN model shows a higher peak for C i at $-5$⁠, compared to $-5.5$ for NoAGN, suggesting that the AGN enhances C i abundance through increased heating and photodissociation. C ii peaks at $-8$ in the AGN model and $-8.3$ in the NoAGN, reflecting greater ionization in the AGN environment. At $r<100\text{pc}$⁠, the CO PDF peaks at $-4$ for NoAGN and $-3.5$ for AGN, indicating higher CO abundance in the AGN model, likely due to shock heating and enhanced star formation. These results underscore the significant influence of an AGN on the chemical pathways and abundances of carbon species.

In the NoAGN model, the intensity maps show higher values of C i and CO at radii below 50 pc and time-scales under 7 Myr, suggesting a more pronounced interface between ionized and molecular gas, likely due to local star formation processes. In contrast, the AGN model reveals increased C i and CO intensities in the outer regions (greater than 50 pc) at time-scales exceeding 7 Myr, indicating that AGN feedback contributes to additional heating and excitation of the remaining molecular gas. The AGN model also exhibits a higher C i/CO ratio at larger radii at later time points, reflecting the AGN’s influence on photodissociation and the chemical composition of the gas. Meanwhile, the NoAGN model shows a dominant C i/CO ratio at smaller radii, highlighting a different balance of physical processes. This aligns with observations of galaxies such as NGC 1808 and NGC 7469, which show higher C i/CO intensity ratios in their outer regions. While our simulations predict CO emissions extending to larger radii, the observed upward tails in C i versus CO ratios are mainly due to detection limits. Overall, our models indicate that CO is higher in the AGN model, particularly in outer regions (⁠$>$50 pc) and inner region at later times (⁠$>$7 Myr). Further, the observed trends are limited to the range $1<\log (I_{\text{CO}})<3$⁠, while simulations extend to lower CO values, $-1<\log (I_{\text{CO}})<3$⁠, suggesting that conditions for lower CO abundance exist beyond the observed limits. The scatter in C i intensity brightness is approximately $-2\text{dex}$⁠, indicating stability, whereas the C i/CO ratio shows greater variability at about $-2.5$ dex. Lower C i/CO ratios at distances $r>200\text{kpc}$ imply less favourable conditions for C i formation due to lower densities and reduced star formation activity. Additionally, time evolution contributes to these trends, highlighting the dynamic interplay between carbon species and the influence of galactic structure on chemical evolution.

The higher C ii emission in the NoAGN model’s core suggests that the absence of AGN feedback helps preserve dense C ii-emitting gas. The overall decrease in C ii emission in outer regions may result from star formation and other galaxy-scale processes. The C i/C ii ratio profile reflects a complex interaction between the AGN and NoAGN cases, with the AGN model initially showing a lower C i/C ii ratio in the central region but a higher ratio in outer regions over time. This suggests that the AGN enhances C i production relative to C ii in these areas due to its influence on ionization and thermal conditions.

In the NoAGN model, column density remains above ${10}^{21}{\text{cm}}^{-3}$ during the early time-scale (t < 6 Myr), with CO, C i, and C ii decreasing as $N_{{\mathrm{H}}_2}$ increases. However, after t = 6 Myr, their intensities rise with $N_{{\mathrm{H}}_2}$ at lower column densities before declining again. This indicates a shift towards more favourable conditions over time. Conversely, in the AGN model, CO consistently decreases across all time-scales, especially after t > 7 Myr, reflecting stronger suppression. C i and C ii also decline with increasing $N_{{\mathrm{H}}_2}$⁠, emphasizing the AGN’s significant impact on the chemical dynamics of carbon species. Additionally, our findings highlight the presence of a ‘CO-dark’ region in the AGN model, where lower J CO transitions [e.g. CO(1–0) and CO(2–1)] are significantly depleted relative to H${}_2$ column density. This depletion is likely a result of the AGN’s mechanical feedback at later time-scale, which disrupts and dissociates CO molecular gas. In contrast, higher J CO transitions may show less pronounced depletion or even an increase in intensity, as the AGN can lead to heating and excitation of the remaining gas.

It is important to note that the specific impact of the AGN on the integrated intensity maps and radial profiles of all molecules and atomic species in this study can vary depending on the properties of the AGN, the molecular gas environment, and the specific stage of AGN activity. Observational studies and modelling efforts are essential for characterizing the AGN’s influence and understanding the underlying physical processes.

ACKNOWLEDGEMENTS

MR and SV acknowledge support from the European Research Council (ERC) Advanced Grant MOPPEX 833460. MR would like to acknowledge SurfSARA Computing and Networking Services for their support (EINF-5315). MR would like to express his heartfelt appreciation to EuroSpaceHub and LUNEX EuroMoonMars Earth Space Innovation for their generous funding and unwavering support.

Data manipulation and analysis were performed using the pandas python package (The pandas development team 2020). pandas is a powerful open-source library that provides efficient data structures and data analysis tools, making it a popular choice for handling structured data in the python programming language.

DATA AVAILABILITY

Data from the simulations are based on the output from publicly available gizmo (Hopkins 2015) code at the following repository: https://bitbucket.org/phopkins/gizmo-public/src/master/, generated on the initial conditions (ICs) described in Section 2. The data underlying this article will be shared on reasonable request to the corresponding author. The simulations in this paper also used the chimes non-equilibrium chemistry and cooling code, which is publicly available from the webpage https://richings.bitbucket.io/.

Footnotes

REFERENCES