Modelling molecular clouds and CO excitation in AGN-host galaxies

We present a new physically-motivated model for estimating the molecular line emission in active galaxies. The model takes into account (i) the internal density structure of giant molecular clouds (GMCs), (ii) the heating associated both to stars and to the active galactic nuclei (AGN), respectively producing photodissociation regions (PDRs) and X-ray dominated regions (XDRs) within the GMCs, and (iii) the mass distribution of GMCs within the galaxy volume. The model needs, as input parameters, the radial profiles of molecular mass, far-UV flux and X-ray flux for a given galaxy, and it has two free parameters: the CO-to-H2 conversion factor $\alpha_{CO}$, and the X-ray attenuation column density $N_H$. We test this model on a sample of 24 local ($z \leq 0.06$) AGN-host galaxies, simulating their carbon monoxide spectral line energy distribution (CO SLED). We compare the results with the available observations and calculate, for each galaxy, the best ($\alpha_{CO}$, $N_H$) with a Markov chain Monte Carlo algorithm, finding values consistent with those present in the literature. We find a median $\alpha_{CO} = 4.8$ M$_{\odot}$ (K km s$^{-1}$ pc$^{2}$)$^{-1}$ for our sample. In all the modelled galaxies, we find the XDR component of the CO SLED to dominate the CO luminosity from $J_{\text{upp}} \geq 4$. We conclude that, once a detailed distribution of molecular gas density is taken into account, PDR emission at mid-/high-$J$ becomes negligible with respect to XDR.


INTRODUCTION
The molecular gas in galaxies is the fuel for star formation (SF) and is affected by active galactic nuclei (AGN) feedback processes (Fabian 2012;Kennicutt & Evans 2012;Cattaneo 2019), such as accretion (Combes et al. 2013;Izumi et al. 2016;Farrah et al. 2022b) and/or outflows (the review by Veilleux et al. 2020, and more recently García-Bernete et al. 2021;Ramos Almeida et al. 2022;Alonso Herrero et al. 2023).Most of the molecular gas resides in Giant Molecular Clouds (GMCs) that are the birthplaces of stars and, for this reason, influence the overall chemical and dynamical evolution of a galaxy (Hoyle 1953;Tan 2000;Ballesteros-Paredes et al. 2020).
From observational and theoretical works (see Elmegreen & Scalo 2004, for an extensive review), it is well established that GMCs have an internal hierarchical structure constituted by clumps resulting from the supersonic compression produced by the turbulence (McKee & Ostriker 2007;Hennebelle & Falgarone 2012).Furthermore, turbulence in GMCs provides global support against the gravitational collapse (Federrath & Klessen 2013).Several studies have collected catalogues of GMCs through the galaxy structure, finding that their ★ E-mail: federico.esposito7@unibo.itdistribution in mass and size is well described by a power-law function (Fukui & Kawamura 2010;Heyer & Dame 2015;Brunetti et al. 2021;Rosolowsky et al. 2021).
The molecular gas in GMCs is mainly observed through the carbon monoxide (CO) rotational lines.The CO(1-0) luminosity is the most used proxy for the total molecular gas mass in galaxies, which can be obtained assuming the so-called CO-to-H 2 conversion factor (Solomon & Vanden Bout 2005;Bolatto et al. 2013, and references therein).Moreover, the CO Spectral Line Energy Distribution (CO SLED), i.e. the luminosity of CO rotational lines as a function of the rotational quantum number , is a very powerful diagnostic for the physical conditions of the molecular interstellar medium (ISM, Narayanan & Krumholz 2014;Kamenetzky et al. 2018;Vallini et al. 2019;Valentino et al. 2021;Wolfire et al. 2022;Farrah et al. 2022a).

MODEL OUTLINE
We build a physically motivated model1 that decribes the interaction of the radiation from AGN and star formation with the molecular gas in galaxies.This model takes the Vallini et al. (2017Vallini et al. ( , 2018Vallini et al. ( , 2019) ) works -focused on the far-infrared (FIR) and molecular line emission from GMCs in high- galaxies -as starting point.In those works, single GMCs were modelled as collections of clumps in a log-normal density distribution, their CO SLEDs were computed with radiative transfer calculations, and a galaxy was filled with a uniform distribution of identical GMCs.The aim of the present Figure 1.A sketch of the workflow followed in this paper.First, we define a synthetic GMC with four physical variables (, ,  0 , M), which define the density distribution of its molecular clumps (as described in Section 2.1).We combine Cloudy results to generate the CO SLED of each clump and of the synthetic GMC (see Section 2.2).Having collected the four radial profiles for a real galaxy ( mol ( ),  mol ( ),  0 ( ),  X ( )), we filled it with the GMC distribution described in Equation ( 4), we associate every GMC with the corresponding incident fluxes, and eventually we have the CO SLED of the studied galaxy.
work is to extend such analysis, including different GMCs in single galaxies (following a mass distribution), and illuminating them with a differential radiative flux (following observed radial profiles).
Figure 1 shows a sketch that summarizes its modular structure, which from the sub-pc scales of clumps within GMCs (see Sec. 2.1) progressively zooms-out to the kpc scales of the gas distribution within galaxies (see Sec. 2.3).More precisely, Section 2.1 deals with the analytical description of the internal structure and mass distribution of GMCs, Section 2.2 outlines the RT modelling implemented to compute the CO emission, whereas in Sec.2.3 we present our assumptions concerning the molecular gas distribution on kpc-scales and the resulting total CO emission from galaxies.

GMC internal structure and mass distribution
We characterize GMCs with four physical parameters: the total mass (), the GMC radius (), the mean density ( 0 ) and the Mach number (M).The density structure of GMCs (see e.g.Hennebelle & Falgarone 2012) is well described by a log-normal distribution (Vazquez-Semadeni 1994;Padoan & Nordlund 2011;Federrath et al. 2011;Vallini et al. 2017).The volume-weighted probability distribution function (PDF) of gas density , in a supersonically turbulent, isothermal cloud of mean density  0 , is     where  = ln(/ 0 ) is the logarithmic density, with mean value  0 = − 2  /2.The standard deviation of the distribution,   , depends on M and on the  factor, which parametrizes the kinetic energy injection mechanism (often referred to as forcing) driving the turbulence (Molina et al. 2012):

Name
We assume  = 0.3, which is the value for purely solenoidal forcing (Federrath et al. 2008;Molina et al. 2012).
We define a GMC as a spherical collection of clumps, whose densities are distributed following the PDF in Equation (1).Accounting for the presence of clumps within the GMCs is one of the fundamental features of this work, as dense clumps emit the bulk of the CO luminosity ( crit > 10 4 cm −3 for  ≥ 2, Carilli & Walter 2013).We then randomly extract clumps with  ≥  0 from the PDF and for each of them we calculate its Jeans mass,   = (4/3) 3    , and the Jeans radius: Here   = √︁   /(  ) is the sound speed, that depends, after we set the adiabatic index to  = 5/3 and the mean molecular weight to  = 1.22,only on the gas temperature  and the clump number density .We assume a fixed temperature of  = 10 K for all our clumps, which is the typical value for dense clumps in molecular clouds (Spitzer 1978;Hughes et al. 2016;Elia et al. 2021).We proceed with the clumps extraction until we fill the entire GMC mass (with a 10% tolerance).Following Vallini et al. (2017), we distribute the leftover mass through the whole GMC, accounting for what we define intraclump medium (ICM).
We consider 15 different GMC models, with masses in the range  = 10 3 − 10 6 M ⊙ (0.2 dex steps), and a fixed surface density Σ = 170 M ⊙ pc −2 (McKee & Ostriker 2007).This translates into GMCs radii in the range  = 1 − 40 pc, and mean number densities,  0 = 3/(  4 3 ), in the range  0 = 10 2 − 3 × 10 3 cm −3 .The Mach number is set to M = 10 for all the GMCs in this work (Krumholz & McKee 2005).In Figure 2 we plot the distribution in density of the clumps within the three GMCs listed in Table 1: different GMCs (shown with different colors) have a different peak in the clumps density distribution, since they have a different mean density  0 .The ICM appears, in Figure 2, in single clumps at low-end densities.
The observed GMC mass distribution in galaxies is well approximated by a power-law (see Chevance et al. 2023, for a recent review).In particular, we follow Roman-Duval et al. (2010) and Dutkowska & Kristensen (2022) and we set the mass distribution as follows:  4) down to  = 10 3 M ⊙ , because we expect the majority of GMCs in galaxies to be small and with low masses (e.g.Fukui & Kawamura 2010;Colombo et al. 2014;Rosolowsky et al. 2021).
We verified that the choice of different power-law exponents for Equation (4), namely 1.39, 1.89 and 2.5, corresponding to the lowest and highest values in Roman-Duval et al. (2010) and the highest value from Chevance et al. (2023) does not greatly affect our results.
In Table 1, we list the two extremes in mass, labelled as "Tiny" and "Huge", respectively, and the "Median" GMCs, extracted from the distribution.

Radiative transfer modelling
We use Cloudy (version 17.02, Ferland et al. 2017) to model the CO line emission from a 1-D gas slab with constant density  and irradiated by a stellar or AGN incident flux.In what follows we label a Cloudy run as "PDR-model" if the incident far-ultraviolet (FUV) photons flux is produced by a stellar population.As it is a standard in PDR studies (Wolfire et al. 2022) the FUV flux ( 0 ) is parametrized in terms of the Habing field 1.6 × 10 −3 erg s −1 cm −2 (Habing 1968), which equates to  0 = 1.Conversely, we label a run as "XDR-model" if the incident flux  X (in the 1 − 100 keV range) comes from an AGN.
For PDR models we adopt the Spectral Energy Distribution (SED) of the incident radiation from the stellar population synthesis code Starburst99 (Leitherer et al. 1999), assuming a continuous star formation model with age  = 10 Myr and solar metallicity, and log( 0 ) = 0 − 6.In the XDR models, we use log[ X /(erg s −1 cm −2 )] = −1 − 4, and the SED is set with the Cloudy command table xdr, that generates a truncated AGN Xray SED,   ∝ ℎ −0.7 , in the 1-100 keV range (see Maloney et al. 1996, for details).3).The clump fluxes have been converted to solar units.Left and right panels are PDR and XDR models, respectively.In the upper panels the orange, purple, and brown lines represent the mean CO SLED for log  = 2.5, 3.5, 4.5 cm −3 , respectively, with the incident flux ( 0 for PDRs,  X for XDRs) left free to vary over the ranges 10 0 − 10 6 and 10 −1 − 10 4 erg s −1 cm −2 , respectively.In the lower panels the yellow, green, and blue lines represent the median CO SLED for log  0 = 0.0, 2.75, 6.0 (bottom-left panel) and for log  X = −1.0,1.5, 4.0 ergs −1 cm −2 (bottom-right panel), respectively, while the volume density is left free to vary over the range 10 0 − 10 6.5 cm −3 .In all panels the solid lines represent the mean CO SLED, while the dashed lines with the same colors show the minimum and maximum CO SLEDs at the given density or incident flux.If the minimum luminosity of a given SLED is out of the plot, then the median curve is plotted with a downward arrow to highlight the presence of very low luminosities.
The second free parameter in the PDR and XDR models is the gas density .To cover the density range of clumps and ICM in GMCs (see Section 2.1), we consider a grid of Cloudy runs spanning log(/cm −3 ) = 0 − 6.5.The ranges of ,  0 and  X have logarithmic steps of 0.25 dex.All our Cloudy runs assume solar metallicity for the gas, and elemental abundances from Grevesse et al. (2010).Moreover, we account for the Cosmic Microwave Background ( CMB = 2.7 K), the Milky Way cosmic ray ionization rate,  = 2 × 10 −16 s −1 (Indriolo et al. 2007), and we set the turbulence velocity  = 1.5 km s −1 (see Pensabene et al. 2021).
Cloudy solves the RT through the gas slab by dividing it into a large number of thin layers.Starting at the illuminated face of the slab, it computes the cumulative emergent flux at every layer.While in this work we are interested in the CO lines from CO(1 − 0) to CO(13 − 12), Cloudy also computes other molecular and atomic lines (e.g.HCN, HCO + , [CII]) that we store in our database.We compute the emergent line emission up to a total gas column density of log( H /cm −2 ) = 24.5:this choice allows us to fully sample the molecular part of the irradiated gas clouds, typically located at  H > 2 × 10 22 cm −2 (McKee & Ostriker 2007).The CO emission from a GMC is obtained by interpolating the Cloudy outputs at the density   and radius   for each -th clump and then summing up all the clumps luminosities.

The CO SLED of single clumps and GMCs
In what follows, we first discuss the resulting CO SLED from clumps of different densities within a given GMC, noting that the luminosity of each -th clump is  = 4 2    (  ), where   (  ) is the flux computed by Cloudy.This is shown in Figure 3.
As shown in the upper panel of Figure 3 the CO SLEDs for XDR clump models are overall brighter at a given gas density, with respect to the CO SLEDs of PDR clump models.Leaving the incident flux free to vary (upper panels of Figure 3) has a limited effect in the CO cooling energy output of PDR models (top-left panel).Varying the incident flux seems more important in XDR models (top-right panel), since the derived CO SLED luminosities can decrease by more than four orders of magnitudes (as pointed by the downward arrows in the top-right panel of Figure 3).
Leaving the density free to vary while fixing the incident flux, however, could make the CO SLED luminosity range over several orders of magnitude for any possible value of incident flux, both for PDR and XDR models (bottom panels of Figure 3).The very  1.The dashed black lines represent the CO SLED of the tiny (upper panels), median (middle panels), huge (lower panels) GMC models, while the grey shaded area highlights the variation in the GMC CO SLED produced by changing the incident flux.Left panels are for PDR models (i.e. the incident flux range is 10 0 ≤  0 ≤ 10 6 ), right panels are for XDR models (i.e. the incident flux range is 10 −1 ≤  X /(erg s −1 cm −2 ) ≤ 10 4 ).The colored lines highlight the contribution of the clumps within each GMC, as a function of their density .
low-luminosity SLEDs correspond to the lowest densities computed (log /cm −3 = 0).It is interesting to note that the highest  X values, which correspond to clumps closer to the AGN, return very faint CO SLEDs: this is because high X-ray photon fluxes lead to the CO dissociation (see Wolfire et al. 2022, for a recent review of XDR processes).
As detailed in Section 2.1, for a given modelled cloud we have a list of extracted clumps.Once the contribution from the various clumps in the GMC is accounted for, by summing up their CO luminosity, the global CO SLED of a GMC depends on the incident flux only.
In Figure 4, we show the CO SLED for the three GMCs listed in Table 1.The variation in the GMC CO SLED introduced by spanning the whole  0 and  X ranges is highlighted by a grey shaded area.We note that varying  0 has a limited effect on the predicted CO SLEDs (left panels in Figure 4), as opposed to varying  X (right panels).
In Figure 4 we plot the contribution of clumps of different densities to the global GMC CO SLED.PDR models (left panels of Figure 4) are dominated by clumps with density  = 10 4 , 10 5 , 10 6 cm −3 in the low ( ≤ 4), medium (5 ≤  ≤ 7) and high- ( ≥ 8) transitions, respectively.In the XDR case (right panels of Figure 4), instead, the very high-density clumps ( ∼ 10 6 cm −3 ) never dominate the CO SLED: their contribution is always at least one order of magnitude lower than clumps with  = 10 4 − 10 5 cm −3 .However, as can be noticed in Figure 3, the very high-density clumps can sustain a high XDR emission only when a large incident flux ( X ∼ 10 4 erg s −1 cm −2 ) is present.This is because at high volume and, thus, column densities ( H =   ∝  1/2 ), the molecular gas is efficiently shielded from external radiation, so that only a large incident flux can penetrate the gas.The same argument applies to PDRs.

Galaxy radial profiles
In this Section we describe how we distribute the GMCs throughout the galaxy volume, and how we associate to each GMC an incident flux ( 0 or  X ) depending on the position within the galaxy.To do so, we need the radial profiles of the molecular mass  mol (), the molecular volume  mol (), the FUV flux  0 (), and the X-ray flux  X ().Spatially resolved observations of nearby galaxies have shown that low- CO emission is well described by an exponential profile along the galactocentric radius  (Boselli et al. 2014;Casasola et al. 2017Casasola et al. , 2020)), with scale factor  CO = 0.17 25 , where  25 is the radius of the galaxy at the isophotal level 25 mag arcsec −2 in the B-band.The CO(1 − 0) emission can be converted into the molecular mass by adopting a CO-to-H 2 conversion factor  CO (e.g.Bolatto et al. 2013, for a review).We adopt the Milky Way value of  CO = 4.3 M ⊙ (K km s −1 pc 2 ) −1 (which includes the helium contribution) as our initial fiducial value, but we will discuss the effect of relaxing this assumption in Section 5.1.We thus obtain the following cumulative molecular mass radial profile: where  CO,tot is the CO(1 − 0) luminosity of the whole galaxy in L ⊙ ,  mol is in M ⊙ .This molecular gas mass is also confined within a volume radial profile  ().We approximate the galaxy as a disc, with half-height  CO = 0.01 25 (Boselli et al. 2014;Casasola et al. 2020), in its external part, while in its inner part ( < 1.5  CO ) we use the equation for a sphere to model a bulge: The profile  () changes at  = 1.5CO to ensure a smooth transition between the spherical and the disc-like regions.Equations ( 5) − (6) imply that the molecular gas volume density  mol () increases towards small radii.The gas surface density Σ() is instead roughly constant in the central part (i.e. at  ≲  CO ), and decreases exponentially as Σ mol () = Σ mol (0) −/ CO .
In our model the molecular gas -distributed according to Equations ( 5) − (6) -is irradiated by a FUV or X-ray flux that depend on the galactocentric radius.Following E22, we model  0 () with a Sersic function (Sérsic 1963): which is characterized by three parameters: the effective radius   , the shape parameter  and the normalization  0 (  ).We assume a minimum  0 = 1 at every radius, which is equal to the Milky Way interstellar radiation field (ISRF).The X-ray profile,  X (), is derived from the intrinsic X-ray luminosity  X , in the 1 − 100 keV range.Given the importance of the attenuation of the X-ray flux due to obscuring gas, usually measured in terms of gas column density  H (), we use the following formula (Maloney et al. 1996;Galliano et al. 2003): where  22 () =  H ()/10 22 cm −2 .
From the derived radial profiles of the molecular mass and volume, we can easily compute the column density radial profile as Whenever the gas column density  H () < 10 22 cm −2 , due to the marginal impact on the X-rays attenuation (e.g.Hickox & Alexander 2018), we use the classical definition  X () =  X /(4 2 ).We term this "Baseline model", as it is also possible to fit the average  H value from the observed CO SLEDs (see Section 4).We will discuss in detail the effect of  H () in Section 5.2.

THE DATASET
In E22, we gathered a wealth of observational data for a sample of 35 local ( ≤ 0.15) AGN-host galaxies.For each source we collected: CO line luminosities from CO(1 − 0) to CO(13 − 12), mostly coming from Herschel observations; the optical radius  25 (from the HyperLeda 2 database); the total IR luminosity  IR (8 − 1000 m), mostly from IRAS observations; the intrinsic X-ray luminosity  X (2 − 10 keV) and the obscuring column density of the X-ray photons,  H, X-ray , from a variety of X-ray observatories; and the FIR flux (observed by Herschel) from which we derived the FUV field radial profile,  0 () (i.e. the three parameters of Equation 7:  0 (  ),   , and ).Any change in the data with respect to E22 is reported in Appendix A.
First of all, we select a sub-sample of the 35 galaxies analyzed in E22.We adopt the following selection criteria: (i) the availability of the CO(1−0) detection, necessary for the derivation of the molecular gas mass of a galaxy; (ii) the availability of both incident flux radial profiles, i.e.  0 () and  X (); (iii) an apparent optical radius  25 ≤ 250 ′′ , since larger objects would have a  CO larger than the CO beams; (iv) at least three CO detections (i.e.without counting the upper limits) between CO(1 − 0) and CO(13 − 12), to ensure we have enough data points to make a meaningful comparison with the model output.These selection criteria reduce our sample to 24 active galaxies.We list, in Table 2, the  CO , the median  0 (), the  X and the  H, X-ray of each galaxy in the sample.
The parent sample was selected considering the availability of Herschel data, and being their intrinsic X-ray luminosity  X ≥ 10 42 erg s −1 in the 2 − 10 keV range.The sample of 24 objects considered in this work is composed by nearby ( ≤ 0.062) moderately powerful ( X ≤ 10 44 erg s −1 ) AGN-host galaxies, but at the same time they display a wide variety of physical properties.Four galaxies of our sample (IRAS 19254−7245, IRAS 23128−5919, Mrk 848, and NGC 6240) are clear mergers, while other six (IRAS F05189−2524, IRAS 20551−4250, Mrk 231, NGC 34, NGC 1275, and UGC 5101) show some signs of interaction as tidal tails.The rest of the sample has a spiral-like morphology.Half of the galaxies (12/24) are luminous infrared galaxies (LIRGs: 10 11 ≤  IR < 10 12 L ⊙ ), while 5 are ultra-LIRGs (ULIRGs:  IR ≥ 10 12 L ⊙ ).Almost all of the galaxies (22/24) have their nuclear activity classified as Seyfert, with the exceptions of Mrk 848 and IRAS 20551−4250 which are classified as lowionization nuclear emission-line regions (LINERs).NGC 1275 (also known as 3C 84, Perseus A) is the only radio-galaxy of our sample, and it is also the central dominant galaxy in the Perseus Cluster.Two of the furthest galaxies (Mrk 231 and I Zw 1) are also commonly classified as quasi-stellar objects (QSOs).We refer to E22 for further details on the galaxy sample and data collection.The CO SLEDs of all the objects in our sample are shown in Figure 5.

RESULTS
The model takes into account the radial profiles defined in Equations (5) − (9) to spatially distribute the 15 GMCs described in Section 2.1 in a real galaxy.The profiles for all the galaxies in the E22 sample are presented in Figure B1.We divide these profiles 2 http://leda.univ-lyon1.frinto logarithmically-spaced radial bins.Based on some tests, we adopt a 0.05 dex bin size, starting from  min = 10 −3 kpc, and up to  max = 2 CO .First, since for every bin we know the molecular mass and the molecular volume, we extract random GMCs from the distribution (Equation 4), until we fill up the entire mass and/or volume of the bin.Then, we associate each GMC, in each radial bin, with the matching incident fluxes  0 (),  X (), and finally produce the expected CO SLED of each galaxy.

The different models
The molecular mass of a galaxy is highly dependent on the choice of the CO-to-H 2 conversion factor  CO , while the X-ray flux can be strongly attenuated by gas clouds between the AGN and the GMCs which we parametrize in terms of a hydrogen column density profile  H ().
The model discussed in Section 2 assumes a default Milky-Way value of  CO = 4.3 M ⊙ (K km s −1 pc 2 ) −1 , and the intrinsic  H () profile from Equation (9).We label this "Baseline model".
For each galaxy we search also for the "Best-fit model" following the Bayesian Markov chain Monte Carlo (MCMC) method: we use  Notes.The molecular masses  mol are calculated within a 2  CO radius and by using the best-fit  CO (except for NGC 7172, for which we used  CO = 4.3). 0 is the median value of the  0 ( ) profile. H,X is the column density derived from the X-ray SED.The best-fit  (1-0) PDR ,  (1-0) XDR ,  (4-3) PDR , and  (4-3)  XDR are the fractions to the total CO(1 − 0) and CO(4 − 3) luminosities, respectively, due to PDR and XDR emission, respectively.Best-fit  CO is in units of M ⊙ (K km s −1 pc 2 ) −1 .For both  CO and log  H , we report the median values of the marginalized posterior distributions, with 1 width as errors.
the  2 likelihood function to fit the observed CO SLED and determine the posterior probability distribution of the model parameters, i.e.  CO and  H , with uniform prior distributions  CO = [0.43,43] and  H = [10 22 , 10 25 ] cm −2 .In this model,  H has a constant value at every radius, and it acts as an average  H seen by the GMCs.To run the MCMC algorithm, we use the open-source Python package emcee (Foreman-Mackey et al. 2013), which implements the Goodman and Weare's Affine Invariant MCMC Ensemble sampler (Goodman & Weare 2010).In this way we are able to fully characterize any degeneracy between our model parameters, while also providing the 1 spread of the posterior distribution for each of them.To be able to include also upper limits in the likelihood function, we follow the approach described in Sawicki (2012) and Boquien et al. (2019), splitting the  2 formula in two sums: where the first sum contains the detections   and their errors   (  covers the first 13 lines of the CO SLED), the second one containing the 3 upper limits  ul,  (whereas we used the 1 upper limit as the measured error   ), and in both sums   are the model values.
Finally, we produce a third set of synthetic CO SLEDs by keeping the default  CO = 4.3 M ⊙ (K km s −1 pc 2 ) −1 but using the column density  H, X-ray derived from the absorption of X-rays along the line of sight (Brightman & Nandra 2011;Ricci et al. 2017;Marchesi et al. 2019;La Caria et al. 2019, and Salvestrini in prep.).We label these CO SLED predictions " H, X-ray model".Also in this case we set  H, X-ray constant at every radius.

Best-fit model results
We run emcee with 10 walkers exploring the parameter space for 10 4 chain steps.The chains have been initialized by distributing the walkers around  CO = 4.3 M ⊙ (K km s −1 pc 2 ) −1 and the median  H () for each galaxy (if larger than 10 22 cm −2 , 10 22 cm −2 otherwise).For each walker, we first run the algorithm with a burn-in chunk of 5000 steps which we later discard.The procedure gives a mean autocorrelation length  = 65 for both  CO and  H , and a mean acceptance fraction of 0.59.In the following we report the median values of the marginalized posterior probability distributions for the two parameters, with a 68% confidence interval as errors.The MCMC code did not converge for NGC 7172, due to the dominance of upper limits (7, with only 3 detections): in the subsequent para- graphs and sections, we only consider the remaining 23 galaxies of our sample.
The best-fit procedure returns median values 0.5 ≤  CO < 9.5 M ⊙ (K km s −1 pc 2 ) −1 , as shown in Figure 6.The median CO-to-H 2 conversion factor for our galaxy sample is  CO = 4.7 +2.4  −2.8 , where the lower and upper errors are the 16 th and 84 th percentiles of the sample distribution, respectively.This value is comparable to that of the Milky Way, and the range agrees well with the available literature (e.g.Bolatto et al. 2013;Leroy et al. 2015;Mashian et al. 2015;Accurso et al. 2017;Seifried et al. 2017;Casasola et al. 2017;Dunne et al. 2022).
The best-fit values of  H are shown in Figure 7 against  H, X-ray .We find values 10 22.0 ≤  H ≤ 10 24.1 cm −2 with a median of log( H /cm −2 ) = 22.1 +1.6 −0.1 .By comparison, the 21/24 galaxies for which we have  H, X-ray have a median log( H, X-ray /cm −2 ) = 23.7 +0.5  −1.8 .The maximum fitted obscuration of our sample is in TOL 1238-364 ( H = 10 24.1 cm −2 ), which also has the maximum  H, X-ray = 10 25 cm −2 .As shown in Figure 7, however, the  H values, derived from the fit of the CO SLED and from the fit of the X-ray spectrum, do not correlate.This, anyway, does not constitute a critical issue of our modelling, as will be explained in Section 5.2.
We find a degeneracy between the  CO and  H sampled values in almost all the galaxies (see e.g. Figure 8 for NGC 3227).This is not unexpected, since it is known that  CO depends on the optical depth (Bolatto et al. 2013;Teng et al. 2023), which in turn depends on  H . Increasing  CO in a galaxy means more molecular mass within the same volume, hence  H has to increase accordingly.The galaxies for which the sampled parameters do not show degeneracy are also among the ones for which our model works worse (e.g.Mrk 231, NGC 4151, NGC 4388, see Appendix C).In Appendix C, together with observed and modelled CO SLEDs for all the galaxies, we also plot the relative (i.e.normalized to the best-fit model) residuals between observed data and best-fit model.We find relative residual values ≥ 2 for at least two detected CO lines in 4 galaxies: IRAS 20551−4250, Mrk 231, NGC 4151, and NGC 4388.Apart from NGC 4388, which has a very peculiar CO SLED shape that our model is not able to reproduce, for the other three galaxies our model fails to reach such high luminosities for the high- lines.This may be due to an additional source of excitation (other than FUV and X-ray flux), as shocks or cosmic rays.In fact, IRAS 20551−4250 and Mrk 231 are known to host powerful CO outflows, with velocities up to 500 and 700 km s −1 , respectively (Lutz et al. 2020).NGC 4388 has also a detected CO outflow, reaching 150 km s −1 (Domínguez-Fernández et al. 2020), while for NGC 4151 we only have a detected H 2 outflow travelling at 300 km s −1 (May et al. 2020).It is tempting to intepret these results as evidences of AGN mechanical feedback on CO excitation, but this is beyond the purpose of this work.
The best-fit  CO and  H , together with the 1 spread of their posterior distributions, are listed, for each galaxy, in Table 2, together with the molecular mass within a 2 CO radius (recalculated with the best-fit  CO ).

Three examples of modelled galaxies
In Figure 9 we show the comparison between the observed CO SLED and that resulting from the "Baseline", "Best-fit", and " H, X-ray " modelling, where the "Best-fit" model is calculated with the ±1 values of the ( CO , log  H ) posterior distributions.We show three sample galaxies: MCG+04−48−002, NGC 3227 and NGC 7130, chosen to be representative (due to the spread of their best-fit parameters) of the results obtained by means of our fitting procedure, with best-fit median values  CO = [7.32, 1.98, 3.59] M ⊙ (K km s −1 pc 2 ) −1 , respectively, and  H = [10 23.96 , 10 22.53 , 10 22.08 ] cm −2 , respectively.In Appendix C we gather the observed and theoretical CO SLEDs for the whole galaxy sample.
MCG+04−48−002 (top panel of Figure 9) has very similar " H, X-ray ", "Baseline" and "Best-fit" CO SLEDs, due to their similar  H values.The effect of a high  H on the CO SLED is visible especially at high .The effect of changing  CO (by a factor of ∼ 80%) is evident by comparing the " H, X-ray ", and "Best-fit" CO SLEDs.
In NGC 3227 (middle panel of Figure 9) we can notice the difference between the "Baseline" and "Best-fit" CO SLEDs due to different  CO and  H : a higher  CO would boost the luminosity of all the CO lines, but a higher  H decreases the high- lines, exposing a typical PDR bump in the low-s (cfr.Figures 3 and 4).The " H, X-ray " SLED is very bright due to its low  H, X-ray , which boosts the XDR emission.
NGC 7130 (bottom panel of Figure 9) has instead three very different modelled CO SLEDs.In the "Baseline" and " H, X-ray " models we can clearly see the PDR bump at low-, due to the high  H which absorbs the X-rays.The higher  H, X-ray makes its SLED to dramatically decrease towards the high- lines, where the XDR emission is dominant.The "Best-fit" CO SLED has instead a lower PDR component due to its low  CO = 3.59, which better reproduces the observed CO SLED.
Figure 11.Galactocentric radius at which the luminosity of a given CO transition reaches 90% of the total value for our sample of galaxies, divided by their  CO .The filled circles, linked with solid lines, mark the  90% for every galaxy, color-coded according to their  X ( CO ), while the brown shaded area represents the values between 16% and 84% of the radii distributions for our sample.

PDR vs. XDR emission
In Figure 10 we plot the relative importance of PDR and XDR emission as a function of  for our galaxy sample as resulting from the "Best-fit" model.It is remarkable the fact that from CO(4 − 3) going upwards the CO luminosity is almost all due to XDR emission, even after taking into account the effect of X-ray flux attenuation with a best-fit  H .This confirms what is shown in Figure 3 for single molecular clumps, and in Figure 4 for single GMCs, i.e. the XDR models dominate the overall CO luminosity, with the exception of the low- lines.These results are extensively discussed in Section 5.3.

The CO SLED radial build-up
Since we modelled the GMC distribution and their CO emission as a function of galactocentric radius, we can study the spatial distribution of different CO lines.In Figure 11 we plot the radius at which the luminosity of a given CO transition reaches 90% of the total value for our sample of galaxies,  90% .It is immediate to see that there is a separation between low- lines, dominated by PDR emission (Figure 10), which emit up to several kpc (i.e. through the whole star-forming galaxy disc), and mid-and high- lines, dominated by XDR emission, which emit most of their luminosity within  = 1 kpc.The spread in radii for a given CO line is due to the different sizes of the studied galaxies, but also, at least for the mid-/high-, to the different incident  X : a higher  X indeed produces a larger XDR emission (as in Meĳerink & Spaans 2005).

DISCUSSION
In this section we discuss the implications of the derived values of the CO-to-H 2 conversion factor  CO and the gas column densities  H for our galaxy sample.Finally, we discuss about the relative importance of PDRs and XDRs to the CO luminosity.

The CO-to-H 2 conversion factor 𝛼 CO
The best-fit procedure returns reasonable values for  CO .We find a median value of  CO = 4.8 +2.4  −2.8 M ⊙ (K km s −1 pc 2 ) −1 , which is similar to the Milky-Way adopted value.Less than a third (8/24) of our sample has a ULIRG-like For some targets we have independent measures of  CO from the literature to compare with.By modelling the gas dynamics, Fei et al. (2023) found, for I Zw 1,  CO = 1.6 ± 0.5 M ⊙ (K km s −1 pc 2 ) −1 , which is lower than our result ( CO = 6.6 ± 1.7 M ⊙ (K km s −1 pc 2 ) −1 ); however, their analysis is limited to the central 1 kpc, whereas our model extends up to 2 CO (i.e.5.6 kpc for I Zw 1).By comparing gas and dust emission, García-Burillo et al. (2014) found values in the range 0.43 − 1.43 M ⊙ (K km s −1 pc 2 ) −1 for the central 700 pc of NGC 1068 (whereas we find  CO = 5.4 +0.9 −0.5 M ⊙ (K km s −1 pc 2 ) −1 up to a diameter of 4.9 kpc).In NGC 6240, the quiescent gas and the outflowing one have different  CO = 3.2 ± 1.8 M ⊙ (K km s −1 pc 2 ) −1 and  CO = 2.1 ± 1.2 M ⊙ (K km s −1 pc 2 ) −1 , respectively (Cicone et al. 2018): our value for this galaxy ( CO = 2.0 ± 0.1 M ⊙ (K km s −1 pc 2 ) −1 ) is closer to the outflowing gas, which is where the bulk of the molecular mass is, according to Cicone et al. (2018).The circumnuclear region (up to a 175 pc radius) of NGC 7469 has a  CO between 1.7 and 3.4 M ⊙ (K km s −1 pc 2 ) −1 according to Davies et al. (2004).Our model gives a higher result ( CO = 5.3 +0.5 −0.7 M ⊙ (K km s −1 pc 2 ) −1 ) but it is derived for the molecular gas extending up to a 4.7 kpc radius.Surprisingly, we find a positive, moderate (Pearson coefficient  = 0.45,  = 0.03), correlation between  CO and the total IR luminosity  IR , as shown in the first panel of Figure 12.The value of  CO is predicted to decrease for higher gas temperatures, hence for starburst galaxies as (U)LIRGs (Bolatto et al. 2013).However, high gas and column densities are expected to increase  CO : for example, Teng et al. (2023) found  CO to increase in the central region of local barred galaxies due to the high CO optical depth, which they find to be a more dominant driver than the gas temperature.A limit of the present analysis is the fact that we used a single  CO for every radius, while it is likely to vary with different local conditions, as metallicity, temperature, density and optical depth (Bolatto et al. 2013;Teng et al. 2023).Another important point to make is that our (U)LIRGs (17/24 galaxies have  IR ≥ 10 11 L ⊙ ) are AGN, so that a fraction of the IR luminosity is due to AGN activity rather than SF (e.g.Nardini et al. 2010;Alonso-Herrero et al. 2012).For example, the largest value we derive is  CO = 9.41 M ⊙ (K km s −1 pc 2 ) −1 for Mrk 231, the only quasar of our sample, with  IR = 10 12.5 L ⊙ , the highest of our sample.
We also find  CO to correlate with the effective radius  ,70 of the 70 m emission ( = 0.47,  = 0.02, second panel of Figure 12).These radii result from the Sersic fit on the Herschel maps analyzed in E22.No significant correlation ( = −0.12, = 0.58) shows up instead between  CO and  CO : this means that the optical and molecular radii ( CO = 0.17  25 ) do not affect the CO-to-H 2 conversion factor, while it seems influenced by the size of the star-forming region (traced by the 70 m emission in our analysis).We do also find  IR and  ,70 to strongly correlate between each other ( = 0.72,  = 8 × 10 −5 ), so it is possible that only one of the two correlations is actually meaningful.There is no significant correlation ( = 0.13,  = 0.56) between  CO and the median  0 listed in Table 2, or the  0 ( ,70 ), so it seems  0 is not the main driver of the  CO variation.
Some clues can come from the way we implemented  CO in our model, i.e. as a normalization factor of the "Best-fit model".A bright IR galaxy has also bright CO emission at every  (e.g.Kamenetzky et al. 2016).However, there is complex interplay between the combined effects of  CO and  H in our model.For example, MCG+04−48−002 (top panel of Figure 9) has a high  CO = 7.32 M ⊙ (K km s −1 pc 2 ) −1 , even though the "Baseline model", made with a lower  CO , is brighter for some CO lines, and this is due to the high  H that we derived for this galaxy, which suppresses the XDR emission.We do not find ( = 0.28,  = 0.20) a convincing correlation between the two best-fit parameters ( CO ,  H ), but we do find (third panel of Figure 12) an anti-correlation ( = −0.43, = 0.04) between  CO and  X measured at a distance  CO from the AGN, and obscured with the best-fit  H : a high  X is indeed expected to boost the XDR CO emission, hence lowering  CO .
In the last panel of Figure 12, we show a positive, moderate ( = 0.44,  = 0.04), correlation between  CO and the PDR fractions to the total CO(1 − 0) modelled emission.For some galaxies (e.g.MCG+04−48−002, which has a high  (1−0) PDR = 0.81 -see Table 2) this is due to the combined effect of  CO and  H : a high  H suppresses the XDR emission, so that a high  CO is needed to boost the otherwise low PDR emission and increase the likelihood.
We stress that the characterization of a class of objects (local AGN as in our study) or even more of individual galaxies by constraining galaxy properties such as  CO , is in line with some recent studies (e.g., Ellison et al. 2021;Casasola et al. 2022;Salvestrini et al. 2022;Thorp et al. 2022).These works highlight the heterogeneity of properties of the local galaxies, in addition to provide observational constraints to theoretical models and for high-redshift studies.

The X-ray attenuation column density
Taking into account the 1 uncertainty, we find non-negligible (i.e. H > 10 22 cm −2 , see Equation 8) X-ray attenuation for most (14/23 galaxies) of our sample.The median is log( H /cm −2 ) = 22.1 +1.6 −0.1 , and the range [10 22 , 10 24 ] cm −2 .By comparison the column density derived from the X-ray SED, for 20/23 galaxies, has median log( H, X-ray /cm −2 ) = 23.7 +0.5  −2.0 and range [10 20.5 , 10 25 ] cm −2 .Although our best-fit  H approximately has the same range of  H, X-ray , the two values do not correlate (Figure 7).This uncorrelation can be due to the different locations of the X-ray absorbers: in the case of  H, X-ray , the X-ray flux is attenuated by gas clouds along the line of sight only, while the  H we derive is due to clouds between the AGN and the GMCs in the whole galaxy volume (i.e.considering multiple lines of sight).
Lately, more complex models of X-ray absorption are leading to different estimates of  H , measuring the obscuration along the line of sight plus a reflective component from the nuclear torus (Yaqoob 2012;Esparza-Arredondo et al. 2021).Several works targeting local AGNs show that these different  H , both estimated from the X-ray SED, fail to correlate between each other (Torres-Albà et al. 2021;Sengupta et al. 2023).
García-Burillo et al. ( 2021) compared instead different estimates of  H , from X-ray absorption, and from high-resolution (7-10 pc) ALMA CO observations, for a sample of 19 local Seyferts, finding similar median values but with a large dispersion around the 1:1 line.As our estimates of Figure 7, they also find less absorbed galaxies (with  H,X-ray < 10 22 cm −2 ) to systematically lie above the 1:1 line, and vice versa.This could be expected if, in less absorbed galaxies,  H,X-ray is dominated by gas clouds smaller than our spatial resolution, which is dictated by the diameter of our smallest modelled GMC (3 pc, see Table 1).
It is believed that, at least in local AGNs, most of the obscuration is due to the nuclear torus rather than gas on galactic scales (see Hickox & Alexander 2018 for a recent review, and Gilli et al. 2022 for high-redshift systems).However, the orientation, filling factor and opening angle of the torus are all important quantities that affect the measurement of  H, X-ray .With our analysis we give a way to measure the attenuation of AGN obscuration on the molecular gas in all the possible directions.

PDR vs. XDR emission
As shown in Figures 3 and 4, the CO luminosity is dominated, for our modelled clumps and GMCs, by XDR emission, at least at  > 3. We find the same by fitting the observed CO SLEDs (Figure 10).The fact that high- lines need X-rays to be radiatively excited is well known and studied (Bradford et al. 2009;van der Werf et al. 2010;Rosenberg et al. 2015;Kamenetzky et al. 2017), however this is the first time a detailed study of CO emission, coming from different gas sizes (clumps, GMCs), produces this same result.The majority of the studies have found, so far, CO excitation to be consistent with PDRs up to CO(6 − 5) (van der Werf et al. 2010;Hailey-Dunsheath et al. 2012;Pozzi et al. 2017).One notable exception is the work by Mingozzi et al. (2018), where their fiducial model for the CO SLED of NGC 34 has a XDR component which starts dominating the luminosity at  upp = 3.
All these studies used a combination of two or three PDR and XDR components, each with a single density and incident flux.To be able to reproduce PDR emission consistent with observations at mid-/high-, these models need very high gas densities ( H ≳ 10 5 cm −3 ), more typical of high-redshift, turbulent star-forming regions (Calura et al. 2022).The mass fraction of such dense clumps within the modelled GMCs is ∼ 2%, being most of the molecular mass in the range [10 2 , 10 4 ] cm −3 , where only XDRs are able to emit significant high- CO emission (see also Figure 4).This also confirms the findings of E22, i.e. that only unlikely extreme gas densities PDR models can reproduce the CO emission.
As a caveat, we note that our model does not include (i) variations in the cosmic ray ionisation rate (CRIR), and (ii) the effect of shocks.It is known that both high cosmic-ray fluxes (Vallini et al. 2019;Bisbas et al. 2023) and shocks (Meĳerink et al. 2013;Mingozzi et al. 2018;Bellocchi et al. 2020) influence the CO SLED, especially the high- lines.Further developments of our model may be able to include a varying CRIR and the treatment of large scale shocks, as in the case of galaxies with winds and outflows (García-Burillo et al. 2014;Morganti et al. 2015;Fiore et al. 2017;Speranza et al. 2022) or mergers (Ueda et al. 2014;Larson et al. 2016;Ellison et al. 2019).

CONCLUSIONS
We have presented a new physically-motivated model for computing the CO SLED in AGN-host galaxies, which takes into account the internal structure of GMCs, the radiative transfer of external FUV and X-ray flux on the molecular gas, and the mass distribution of GMCs within a galaxy volume.The assumptions are that the molecular clumps in a GMC follow a log-normal density distribution, the clumps size is equal to their Jeans radius, and that the GMC masses, in a galaxy, follow a power-law distribution.The model can predict the CO SLED of a galaxy from the radial profiles of molecular mass  mol (), FUV flux  0 () and X-ray flux  X ().We also set two free parameters: the CO-to-H 2 conversion factor  CO , and the X-ray attenuation column density  H .We have used Cloudy to predict the CO lines emission of clumps and GMCs, through a combination of PDR and XDR simulations, and by taking into account the observed radial profiles.
To test the validity of our model, we extracted, from the galaxy sample of E22, a sub-sample of 24 X-ray selected AGN-host galaxies with a well-sampled CO SLED.The sample represents a good variety of local ( ≤ 0.06) moderately luminous (10 42 ≤  X [erg s −1 ] ≤ 10 44 , 10 10 ≤  IR [L ⊙ ] ≤ 10 12.5 ) AGNs.We compared the CO SLEDs produced by our model with the observed ones, and we selected the best-fit  CO and  H with a MCMC analysis for each galaxy.Our main results are here summarized: • the best-fit procedure returned, for our galaxy sample, median values of  CO = 4.8 +2.4  −2.8 M ⊙ (K km s −1 pc 2 ) −1 and log( H /cm −2 ) = 22.1 +1.6 −0.1 , where the lower and upper errors are the 16 th and 84 th percentiles, respectively; we find  CO consistent with the Galactic value, and  H that do not correlate with  H, X-ray ; • XDRs are fundamental to reproduce observed CO SLEDs of AGN-host galaxies, particularly for  upp ≳ 3.
We argue that the larger importance of X-rays to the CO luminosity, with respect to other works in the literature, is mainly due to the fact that the majority of molecular gas mass in a galaxy is at  H ∼ 10 2 − 10 4 cm −3 .This conclusion comes from a physicallymotivated structure for the molecular gas, rather than from singledensity radiative transfer models.
The model predictions can be used to estimate  CO in AGN-host galaxies, by constraining both the relative and absolute intensity of several CO lines.They can also be useful to estimate the effective attenuation of X-ray photons on the molecular gas in a spatiallyresolved way, rather than only on the line-of-sight (i.e.,  H, X-ray ).
Moving forward, the model is able to predict the emission of other molecular species, as HCN, HCO + , H 2 , H 2 O.These molecules have been increasingly used in detailed studies of local AGN (Butterworth et al. 2022;Eibensteiner et al. 2022;Huang et al. 2022).The model is also built to exploit a wealth of observed galaxy data, and compare different molecular lines spatially resolved by ALMA (as done with NGC 1068 in García-Burillo et al. 2019;Nakajima et al. 2023) with our model predictions.

Figure B1.
Radial profiles of our sample of 24 galaxies, sorted by their total molecular mass when measured with a Milky Way  CO .The four panels show, from top to bottom, the molecular mass  mol ( ) (recalculated with the best-fit  CO ), the molecular volume  mol ( ), the incident fluxes  0 ( ) and  X ( ), and the "Baseline model" column density  H ( ). 0 ( ) and  X ( ) are in  0 and erg s −1 cm −2 units and they are plotted with dashed and solid lines, respectively.Colors are the same as Figure 5.

APPENDIX A: UPDATE OF OBSERVED FLUXES FOR OUR SAMPLE
Reviewing the data contained in E22, we have found a better source for one CO line luminosity in one galaxy.The CO(2 − 1) of NGC 1275 was taken, in E22, from Salomé et al. (2011), which was converted to 303 +8 −8 × 10 2 L ⊙ .Since this flux was observed, in Salomé et al. (2011), from a very small nuclear region of NGC 1275, we wanted to improve the measurement by finding another observation with a larger beam.Lazareff et al. (1989) observed the CO(2 − 1) emission with the IRAM-30m telescope, which has a full-width at half maximum (FWHM) of 10. ′′ 5, finding a larger CO(2 − 1) luminosity of 15 +3 −3 × 10 4 L ⊙ .Throughout this paper we use this value, since it is the one with the larger beam we can find for this galaxy; in this way the CO(2 − 1) beam is almost equal to the projected  CO of the galaxy:  CO = 10.′′ 9.

APPENDIX B: RADIAL PROFILES OF GALAXIES
In Figure B1 we show the radial profiles of the molecular mass  mol (), the volume  mol (), the incident fluxes  0 () and  X (), and the intrinsic column density  H ().  mol () has been recalculated with the best-fit  CO .The radial profiles all follow the Equations (5) − (9).

APPENDIX C: OBSERVED AND MODELLED CO SLEDS FOR THE WHOLE GALAXY SAMPLE
Figures C1 -C23 show the "Baseline model", and " H, X-ray model" for each galaxy of the sample with orange and brown lines, respectively.The observed CO SLED is plotted with a black line, with vertical bars representing measurement errors, and downward arrows indicating censored data points (i.e.upper limits).We randomly select, from the posterior distributions of the model parameters  CO and log  H , 100 different modelled CO SLEDs, plotted in red.The blue lines and pink shadings, are the modelled CO SLEDs with the median and 1 spread of the model parameters posterior distributions, respectively; we call "Best-fit model" the one with the median values of  CO and log  H .For each galaxy, the bottom panels of Figures C1 -C23 show in black the relative residuals, calculated as the difference between observations and best fit model, divided by the latter; the best-fit model is the blue line (fixed at 0), with the 1 spread in pink shading.

Figure 2 .
Figure 2. Histograms of clumps density distributions within the three GMCs listed in Table1.The shapes are set from Equation1, and mainly depend on each GMC mean number density  0 .
) Equation (4) has been derived by fitting the observed GMC mass distribution in the Milky Way with a power-law, albeit only in the mass range  = 10 5 − 10 6 M ⊙ due to observational limits that make it hard to sample the distribution towards the low-mass end (Roman-Duval et al. 2010).Dutkowska & Kristensen (2022) extrapolate the relation down to  = 10 4 M ⊙ , but in this work we extrapolate Equation (

Figure 3 .
Figure 3. Synthetic CO SLEDs from CO(1 − 0) to CO(13 − 12) for spherical molecular clumps of radius  =   (Equation3).The clump fluxes have been converted to solar units.Left and right panels are PDR and XDR models, respectively.In the upper panels the orange, purple, and brown lines represent the mean CO SLED for log  = 2.5, 3.5, 4.5 cm −3 , respectively, with the incident flux ( 0 for PDRs,  X for XDRs) left free to vary over the ranges 10 0 − 10 6 and 10 −1 − 10 4 erg s −1 cm −2 , respectively.In the lower panels the yellow, green, and blue lines represent the median CO SLED for log  0 = 0.0, 2.75, 6.0 (bottom-left panel) and for log  X = −1.0,1.5, 4.0 ergs −1 cm −2 (bottom-right panel), respectively, while the volume density is left free to vary over the range 10 0 − 10 6.5 cm −3 .In all panels the solid lines represent the mean CO SLED, while the dashed lines with the same colors show the minimum and maximum CO SLEDs at the given density or incident flux.If the minimum luminosity of a given SLED is out of the plot, then the median curve is plotted with a downward arrow to highlight the presence of very low luminosities.

Figure 4 .
Figure 4. CO SLED of the three GMCs listed in Table1.The dashed black lines represent the CO SLED of the tiny (upper panels), median (middle panels), huge (lower panels) GMC models, while the grey shaded area highlights the variation in the GMC CO SLED produced by changing the incident flux.Left panels are for PDR models (i.e. the incident flux range is 10 0 ≤  0 ≤ 10 6 ), right panels are for XDR models (i.e. the incident flux range is 10 −1 ≤  X /(erg s −1 cm −2 ) ≤ 10 4 ).The colored lines highlight the contribution of the clumps within each GMC, as a function of their density .

Figure 5 .
Figure5.Observed CO SLED of our AGN sample, composed of 24 galaxies.The order of the objects depends on their molecular mass within a radius of 2 CO (blue is low, red is high, with the range being [5 × 10 7 − 9 × 10 10 ] M ⊙ when calculated with a Milky-Way  CO ).Downward arrows are upper limits.

Figure 6 .
Figure 6.Histogram representing the best-fit CO-to-H 2 conversion factor  CO , selected through the minimization of the  2  , for our sample of 24 galaxies.The dashed line is for the Milky Way value (Bolatto et al. 2013), while the shaded hatched area highlights the values associated to ULIRGs (Pérez-Torres et al. 2021).

Figure 7 .
Figure7.A comparison of the column densities  H from our fitting procedure (on the -axis) and the  H, X-ray derived by modelling the absorption of Xrays (on the -axis).The colour of each circle (i.e. of each galaxy) depends on the intrinsic X-ray luminosity  X (2 − 10 keV).The dashed line is the 1:1 line.

Figure 8 .
Figure 8. Corner plot showing the marginalized posterior distributions of  CO (in M ⊙ (K km s −1 pc 2 ) −1 ) and  H (in cm −2 ) for the galaxy NGC 3227.The contours represent (1, 1.5, 2)  levels for the 2D distribution.The best-fitting parameters and the 16 th and 84 th percentiles are plotted with solid black lines and square and dashed magenta lines, respectively.

Figure 9 .
Figure 9.The modelled CO SLEDs of 3 galaxies (from the top to the bottom panel: MCG+04−48−002, NGC 3227 and NGC 7130) of our sample, in physical units of  ⊙ .For each panel, the black line is the observed CO SLED (with downward arrows indicating censored data points), and the orange, brown, red, and pink lines are the modelled CO SLEDs: our baseline model without fitting any free parameter and with a negligible  H ( ), the baseline model with a constant  H,X-ray derived by modelling the absorption of Xrays, 100 MCMC modelled SLEDs randomly picked from the parameters posterior distributions, and the best-fit model covering the 1 spread of such distributions, respectively.The -axis represents the upper rotational quantum number  of each CO line, from CO(1 − 0) to CO(13 − 12).

Figure 10 .
Figure 10.Bar plot showing the relative percentage to the total modelled luminosity of PDR (pink side) and XDR (brown side) emission for each analyzed CO line, from CO(1 − 0) to CO(13 − 12).

Figure 12 .
Figure 12.Scatter plot of best-fit  CO (on the -axis) against, from left to right, the total IR luminosity, the effective radius of the 70 m maps, the PDR fraction of the total CO(1 − 0) modelled emission, and the X-ray flux at  CO , for the whole galaxy sample.The dashed black line is the regression fit between the quantities on the  and  axes, and the brown lines are 100 bootstrapped fits, which highlight the confidence interval of the regression fit.The Pearson correlation coefficient  and the -value appear on every panel on the bottom right.

Table 1 .
Main properties of the smallest, median and largest GMCs in the mass distribution described in Section 2.1.The columns are the mass, radius, mean number density and number of extracted clumps.

Table 2 .
Properties and results on the sample of 24 AGN.