Predictions for CO emission and the CO-to-H$_2$ conversion factor in galaxy simulations with non-equilibrium chemistry

Our ability to trace the star-forming molecular gas is important to our understanding of the Universe. We can trace this gas using CO emission, converting the observed CO intensity into the H$_2$ gas mass of the region using the CO-to-H$_2$ conversion factor (Xco). In this paper, we use simulations to study the conversion factor and the molecular gas within galaxies. We analysed a suite of simulations of isolated disc galaxies, ranging from dwarfs to Milky Way-mass galaxies, that were run using the FIRE-2 subgrid models coupled to the CHIMES non-equilibrium chemistry solver. We use the non-equilibrium abundances from the simulations, and we also compare to results using abundances assuming equilibrium, which we calculate from the simulation in post-processing. Our non-equilibrium simulations are able to reproduce the relation between CO and H$_2$ column densities, and the relation between Xco and metallicity, seen within observations of the Milky Way. We also compare to the xCOLD GASS survey, and find agreement with their data to our predicted CO luminosities at fixed star formation rate. We also find the multivariate function used by xCOLD GASS overpredicts the H$_2$ mass for our simulations, motivating us to suggest an alternative multivariate function of our fitting, though we caution that this fitting is uncertain due to the limited range of galaxy conditions covered by our simulations. We also find that the non-equilibrium chemistry has little effect on the conversion factor (<5\%) for our high-mass galaxies, though still affects the H$_2$ mass and Lco by $\approx$25\%.


INTRODUCTION
Molecular gas plays a vital role within the formation of galaxies, as it fuels the formation of new stars.However, one of the most abundant molecules, H 2 , is difficult to detect directly as it does not produce emission at the cold temperatures (10-20K) typical of the dense molecular clouds that contain most of a galaxy's reservoir of molecular gas.For that reason, we have to rely on using other methods to detect and measure the molecular gas content within galaxies.Some methods use IR-based emission, and assume some dust-to-gas ratio to calculate the mass of H 2 , whilst others use emission from the molecular carbon monoxide, or CO.
CO is used to trace the molecular gas as it has its first transition (J=1-0) at ~5.5 K, which can be easily excited at the typical temperatures found within molecular clouds.CO molecules within these molecular clouds can be excited through a combination of collisions with H 2 and radiative trapping, and de-excite spontaneously through collisions and emissions.However, CO is not a perfect tracer for H 2 , as it is more easily dissociated, resulting in a 'dark gas' (Van Dishoeck & Black 1988;Wolfire et al. 2010).This CO-dark molecular gas ★ Contact e-mail: O.Thompson@hull.ac.uk emits CO weakly and it has been estimated that more than 30% of the molecular H 2 gas is CO-dark using gamma ray emission within the Milky Way (Grenier et al. 2005), with evidence of CO-dark gas in various extragalactic studies too (Poglitsch et al. 1995;Madden et al. 1997).Suggestions for tracers of this faint-CO within the Milky Way have been proposed, such as OH and HCO+ (Liszt & Lucas 1996;Allen et al. 2015;Nguyen et al. 2018).More general suggestions that can be used outside the Milky Way have included atomic carbon (Gerin & Phillips 2000;Papadopoulos et al. 2004;Offner et al. 2014;Glover et al. 2015;Li et al. 2018;Clark et al. 2019), [C ii] (Zanella et al. 2018;Madden et al. 2020;Vizgan et al. 2022;D'Eugenio et al. 2023;Ramambason et al. 2024), and in a small range of environments, HCl (Schilke et al. 1995).Due to a significant amount of H 2 being undetectable by CO, we require a way to infer the H 2 mass from the CO we can observe.
We can link observed line emissions from CO to the total H 2 mass using the CO-to-H 2 conversion factor,  co , which is defined as: where  h 2 is the column density of H 2 , and  co is the velocityintegrated CO(1-0) brightness temperature, which is related to the intensity of CO(1-0) emission line.We can also calculate an alternate version of the CO-to-H 2 conversion factor known as  co , which is defined as: where  h 2 is the H 2 gas mass,  co is the CO luminosity from the galaxy, meaning  co is simply a mass-to-light ratio.Bolatto et al. (2013) and references therein showed that in the Milky Way, the  co conversion factor is estimated to be 2 × 10 20 cm −2 (K km s −1 ) −1 ±30%, and the corresponding  co is 4.3M ⊙ (K km s −1 pc 2 ) −1 ±30% .These dust-based approaches are able to estimate the H 2 mass through IR emission from dust, allowing them to calculate the H 2 column density, assuming some dust-to-gas (DTG) ratio.These observations also showed a correlation between metallicity and  co (see Bolatto et al. 2013, Figure 9) where as metallicity increases, the observed  co decreased, with a sharp increase with metallicities below  ∼ 1/3 -1/2  ⊙ .Observational studies are able to use  co and  co to calculate the molecular gas mass from CO observations.Without a good knowledge and understanding of the conversion factor and how it depends on environment, observational studies would be unable to reasonably estimate the molecular gas content of local and distant galaxies.
It can be shown that  co varies within different galaxy populations by looking at the Kennicutt-Schmidt (KS) relation.When applying a constant  co to the KS law, metal-poor galaxies ( < 0.3 ⊙ ) do not line up with the main relation (see also Cormier et al. 2014) and require higher values of  co to line up (see Kennicutt & Evans 2012, Figure 11).For starburst galaxies, we also see a bimodal KS relation when applying a bimodal  co .To achieve a unimodal relation, we require a continuous  co (see Narayanan et al. 2012, Figure 7).We also see the non-universality in  co in ULIRGs, where a Milky Waylike conversion factor would infer a gas mass that would exceed the dynamical mass of the galaxy (Solomon et al. 1997).
Studies have attempted to create a function for the conversion factor dependant on metallicity (Wilson 1995;Schruba et al. 2012;Genzel et al. 2012) due to the relation mentioned prior.However, these prescriptions have a large scatter on the predicted conversion factor, and there are inconsistencies across these studies.For this reason, we require a multivariate function that captures additional dependencies that metallicity alone can not.
The xCOLD GASS survey (Saintonge et al. 2017) observed CO emission from 532 galaxies and then used those observations to calculate various properties of those galaxies, including the molecular gas mass,  h 2 .Rather than use the value stated in Bolatto et al. (2013), this molecular gas mass was calculated by using the multivariate form of  co from Accurso et al. (2017), which depends on the local conditions of the galaxy as follows: where  [C ii] is the observed [C ii] 158m emission line luminosity, Δ(MS) is the offset from the main sequence, and 12+log(O/H) is the metallicity.This model was found by fitting observational data from the Hershel PACS observations of the [C ii] emission line (Poglitsch et al. 2010;Pilbratt et al. 2010), and low-metallicity galaxies from the xCOLD GASS survey.Using observed scaling relations between the  [C ii] / co(1−0) ratio as a function of integrated properties, Bayesian analysis revealed metallicity and MS offset as the only two parameters needed to quantify the variations in the luminosity ratio.Combining these variations with a linear fit to their radiative transfer models gave the formula for a multivariate  co .We discuss the different forms of the Accurso et al. (2017) multivariate conversion factor in more detail in Section 5.2.
Other prescriptions for  co have been derived, relying on other galaxy parameters.Teng et al. (2024) uses dust-based measurements to calculate  co for objects in the PHANGS-ALMA survey Leroy et al. (2021), and then derived a prescription for  co that uses velocity-dispersion.Similarly, Hirashita (2023) derive a prescription for  co that uses the dust-to-gas ratio and the surface area-to-mass ratio of the dust, which reflects the effect of the dust grain distribution.Shetty et al. (2011a,b) investigated  co and its dependence on galaxy parameters, finding a weak dependence of  co on temperature from 20 K to 100 K, and also finding  co increases for surface densities > 100 M ⊙ pc −2 due to saturation of the CO line in the high-surface density regime.
It is important to better understand this relationship between metallicity and  co , especially as we are now observing high-redshift low-metallicity galaxies from the early Universe using ALMA (Fujimoto et al. 2024).Hu et al. (2023) looked at simulating CO from a low metallicity (<0.1Z ⊙ ) galaxy to match the observation of the Wolf-Lundmark-Melotte (WLM) galaxy made by Rubio et al. (2015), and how dust evolution affects the models created.A hybrid method for the chemical evolution was adopted, where the H network was solved on-the-fly and the C chemistry was post processed using astrochemistry 1 .When including dust growth in their models, they were able to reproduce the observed CO luminosity.The conversion factor that was simulated was found to be close to the Milky Way value, although dust-based observations would suggest it should be higher due to the low metallicity.The authors conclude that using a metallicity-dependent  co factor would underestimate the dust-togas ratio in their simulations.
Other computational studies that have looked into modelling CO and the  co conversion factor include Keating et al. (2020), which used the cosmological zoom-in simulations from the FIRE project (Hopkins et al. 2018b) and combined them with a post-processing method assuming chemical equilibrium using the chimes chemistry solver 2 (Richings et al. 2014a,b), allowing them to directly model the line emission from the simulations using the line radiative transfer code radmc-3d (Dullemond et al. 2012) 3 .When comparing the column densities of H 2 and CO in these investigations, it was found that the shielding length in the simulations had to be reduced by one or two orders of magnitude for the simulated column densities to match with observational data.It was also noted that the  co calculated for their Milky Way-like galaxy was an order of magnitude lower than expected from the Milky Way value of 2 × 10 20 cm −2 (K km s −1 ) −1 .This was noted to potentially be due to the resolution of the simulation, and that simulations with higher resolution could help fix this.
Higher resolution simulations, such as the SILCC simulations (Walch et al. 2015), have simulated pieces of galactic discs, about 500 pc × 500 pc × ±5 kpc.These simulations however solve both the H 2 and CO chemistry on-the-fly, rather than post-processing.A subset of these simulations, SILCC-Zoom (Seifried et al. 2017), used a spatial resolution of 0.06 pc to study the formation and evolution of molecular clouds in high resolution.When examining the CO-to-H 2 conversion factor, they found the average  co scattered around 1 -4 x 10 20 cm −2 (K km s −1 ) −1 , agreeing with observational data and Bolatto et al. (2013).It was also found that in these high resolution simulations, the average  co increased over time in agreement with Glover & Clark (2016).Other high resolution simulations include Gong et al. (2020), which used the TIGRESS simulation framework (Kim & Ostriker 2017), simulating sections of a Milky Way-like galaxy at various radii from the galactic centre up to a resolution of 1 pc.These simulations were able to reproduce the relation between metallicity and  co , finding that  co decreases with a power-law slope of -0.8 for the CO (1-0) line.This metallicity scaling of  co agrees with prior work by Feldmann et al. (2012a,b), which combined non-equilibrium chemistry for H 2 and full radiative transfer cosmological simulations with a subgrid CO model.
Richings & Schaye (2016a,b) also found agreement with Glover & Clark (2016), finding that the average  co decreases with time within their low metallicity (0.1 Z ⊙ ) non-equilibrium simulations.However, at ages > 15 Myr, they note that this trend is uncertain, as there are too few clouds of a high age.Within Richings & Schaye (2016a), they compared the effects of non-equilibrium chemistry on their simulations, and found that CO emission and H 2 mass can be both enhanced and suppressed by non-equilibrium chemistry, meaning that the average  co could be affected up to a factor of ~2.3.They also note that usage of a fluctuating UV radiation field would influence these non-equilibrium effects further and may drive additional effects not captured in their current simulations.
In this paper, we use a suite of isolated galaxy simulations that couple the CHIMES chemistry and cooling module to the FIRE-2 subgrid models for galaxy formation (Richings et al. 2022).We also quantify the effects of non-equilibrium chemistry on those results compared to simulations where we assume chemical equilibrium.
The remainder of this paper is organised as follows.In section 2 we describe the methods we used to create and evolve isolated galaxy simulations, how we track the chemistry within, and how we create simulated spectra from the simulations.In section 3 we detail the isolated galaxy simulations we use, their initial conditions and evolution.In section 4 we compare our simulations and simulated spectra to observational data such as xCOLD GASS and dust observations from the Milky Way.In section 5 we look at the conversion factor and its relation to metallicity, as well as prescribing a new fitting for a multivariate function for the CO-to-H 2 conversion factor.We also test whether  co changes when averaged and blurred over different spatial scales to better match our simulated spectra with observations.In section 6 we summarise our conclusions and our main findings.In Appendix A we show results from our resolution tests.

METHODS
Here we detail the methods used to create the simulations, and then from them the simulated spectra that we use throughout this work.We use the simulations from Richings et al. (2022) that were run with their fiducial model.That work also considered model variations with a uniform radiation field, and using a constant dust-to-metals ratio with no dust depletion.However, for this study we use their fiducial model, with local stellar fluxes computed from star particles and a density and temperature dependent model for the dust depletion factors, which we describe in section 2.2 below.

FIRE subgrid models
The simulations within this work were run with a meshless finite mass solver, using the gravity and hydrodynamics code gizmo (Hopkins 2015).These simulations were run using the fire-2 sub-grid physics models (Hopkins et al. 2018b) which were developed as part of the FIRE collaboration4 .
Stars can form from gas particles that become self-gravitating, Jeans-unstable, and are above a density threshold  > 10 3 cm −3 .Details of the star formation algorithm are further described in appendix C of Hopkins et al. (2018b).We assume each star particle is a single stellar population, where we assume a Kroupa (2001) initial mass function.Stellar feedback from each star particle is then calculated with simple fits to the stellar population models from starburst99 (Leitherer et al. 1999).The stellar feedback mechanisms include Type II and Ia SNe, stellar winds from OB and AGB stars, and radiative feedback.
The FIRE-2 models track the enrichment and evolution of the 11 elements that are used in the chimes chemistry solver.Star particles inject metals through SNe and stellar mass loss, the yields we use are as follows: Type Ia SNe yields from Iwamoto et al. (1999), Type II SNe yields from Nomoto et al. (2006), and OB/AGB yields from Van Den Hoek & Groenewegen (1997), Marigo (2001), andIzzard et al. (2004), which are summarised in appendix A of Hopkins et al. (2018b).Type Ia SNe feedback is calculated using Mannucci et al. (2006), for both the prompt and delayed stellar populations.SNe and stellar winds use a feedback system described in appendix D of Hopkins et al. (2018a) (see also Hopkins et al. (2018b)).

The CHIMES chemistry module
We use the chimes chemistry and cooling module to track the evolution of 157 molecules and ions (Richings et al. 2014a,b).This includes the ionisation states of H, He, C, N, O, Ne, Mg, Si, S, Ca and Fe, as well as 20 molecular species; including H 2 and CO.
As chimes is capable of post-processing the abundances assuming chemical equilibrium whilst also solving the non-equilibrium abundances on-the-fly, we are able to make a direct comparison between the two (e.g., Richings & Schaye 2016a;Richings et al. 2022).Calculating these non-equilibrium abundances on-the-fly adds a large computational cost to the simulations (5-10×), so comparing nonequilibrium and equilibrium results allows us to determine what effect the non-equilibrium chemistry has on the simulation A fluctuating time-dependent UV radiation field is used within this work, which follows the propagation of radiation from star particles using a method based on the lebron algorithm used to model stellar radiation pressure in fire (Hopkins et al. 2018b(Hopkins et al. , 2020)).As the simulations are coupled with chimes, this method is modified.The lebron method tracks three stellar fluxes from all particles; the IR, optical, and UV bands.However, for photochemistry, the radiation is tracked in eight separate stellar age bins allowing accurate tracking of the age-dependence of the UV spectra.The radiation from each age bin is then divided into further bands, splitting the non-ionising FUV band from the ionising EUV band, giving 16 stellar fluxes in total for the photochemistry.The optical and IR bands are not required for the photochemical reactions.
Shielding from this radiation field within local gas clouds is based on a Sobolev-like approximation, which uses a shielding length ( sh ) calculated from the density gradient ∇ and the mean inter-particle spacing, ℎ inter , like so: where the first term of this equation accounts for the size of the gas cloud around the particle, and the second term accounts for the size of the particle.We include a factor of 1/2 in this equation so that it matches equation A10 of Gnedin & Kravtsov (2011) in the limit where the density term dominates.The local column density of species  is then given by   =    sh where   is the density of species .Using the methods described in Richings et al. (2014b), the photochemical rates are then suppressed as a function of the local column densities of Hi, H 2 , Hei, Heii, CO, and dust.This equation treats the shielding using an average shielding length, and therefore column density for the local gas cloud.Therefore, this treatment can tend to overestimate shielding, which can lead to higher molecular abundances than expected.
The chimes non-equilibrium chemistry solver also requires the gas-phase abundances of each element in the reaction network to be input.The simulations track the total elemental abundances, but for some elements, a fraction of their total elemental abundance will be locked up in dust grains and therefore will not contribute to the gas-phase chemistry.For that reason, we need to determine the fraction of each element in dust grains, and then reduce the gas-phase abundances to match.
We use the density and temperature-dependent dust depletion model developed by Richings et al. (2022).This model uses observations from Jenkins (2009), which determined the depletion fraction of metals onto dust grains within the solar neighbourhood, by measuring the column densities of 17 metals and Hi along 243 sightlines in the Milky Way.The model also uses observations from De Cia et al. (2016), which expanded on these observations by recording 70 damped Lyman- absorbers observed through quasar spectra, in addition to the original Milky Way sightlines.Jenkins (2009) found that the overall strength of dust depletion,  * , is best fit to the relation:  * = 0.772 + 0.461 log 10 ⟨ H ⟩. (5) This relation can be used to calculate the strength of dust depletion for each gas particle as a function of density, assuming the particles total hydrogen density  H,tot , is approximately equal to the average neutral hydrogen density along the line of sight to the background source in the observations, ⟨ H ⟩. However this assumption can overestimate the strength of the depletion.We also limit  * at higher densities ( H,tot > 3.12 cm −3 ) which corresponds to the strongest dust depletion strength observed in the Milky Way sightlines, as anything higher would be extrapolating to higher depletions than are observed.A limit linked to temperature is also imposed, where above 106 K all metals are in the gas phase due to any dust grains being destroyed through sputtering (Tsai & Mathews 1995).
The depletion factors of individual elements are then linear fits as a function of the parameter  * , with fit coefficients collated from Jenkins (2009) and De Cia et al. (2016), though these fits are uncertain due to limited observations.Jenkins (2009) only contains a handful of carbon depletion measurements.Sofia et al. (2011) finds the gas-phase column densities of carbon are a factor of ≈ 2 than those measured in Jenkins (2009), which would result in a stronger depletion of carbon in dust than expected from the linear fits by that same factor of 2. These depletion factors are used to reduce the gas-phase abundance of each element.We also sum the mass in dust grains of the 17 elements featured in the two observational papers to determine the total dust abundance, and use that to scale the rate of reactions that occur on the surface of the dust grains and any thermal processes that involve dust grains e.g.photoelectric heating.
As this model is based on Milky Way and damped Lyman- absorbers, it does not explicitly follow the formation and destruction mechanics that govern the abundances of dust grains.However, other studies have been able to develop models that are able to capture these processes, allowing for a more detailed evolution of the dust grain population (Asano et al. 2013;Bekki 2015;Hirashita et al. 2015;McKinnon et al. 2018;Choban et al. 2022).

The RADMC-3D radiative transfer code
Once we have the final snapshots from the simulations (described in Section 3.1), we post-process the CO (1-0) emission line using radmc-3d (Dullemond et al. 2012).radmc-3d allows us to follow the emission, propagation and absorption of the spectral line, along with stellar emission.We are also able to follow the absorption, scattering, and thermal emission from dust grains within the simulation.This means we are also able to calculate the total emission including the continuum, as well as the emission from solely the continuum.Obtaining the line emission is then as simple as subtracting the former from the latter.radmc-3d first constructs an adaptive mesh refinement (AMR) grid.We refine the grid until each cell contains no more than 8 star particles and/or gas particles.The abundances of the ion and molecules, as well as their ion-weighted temperatures and velocities, are projected onto the AMR grid.The star particles are also projected onto the grid, and are smoothed, after being split into the eight stellar age bins mentioned in section 2.2.We also include graphite and silicate grains in the radmc-3d simulated spectra, using the abundances of each at solar metallicity from the grain abundances from v13.01 of the cloudy photoionisation code (Mathis et al. 1977;Ferland et al. 2013), and we scale these for each gas particle first by the total metallicity relative to solar for each gas particle, and then by the density and temperature-dependent dust to metals ratio predicted by our dust depletion model.The dust temperature in each cell of the AMR grid is then calculated using stellar radiation.
We then calculate the level populations of ions and molecules within each cell using the local velocity gradient (LVG) method (Castor 1970;Goldreich & Kwan 1974;Shetty et al. 2011a), as this approximates the effects of non-local thermodynamic equilibrium (non-LTE).This method calculates an optical depth for each cell using the gas velocities from neighboring cells, and finds an escape probability for each photon by assuming each photon will eventually be Doppler shifted away from the line centre.Gas can also be excited using this, depending on the temperature.We use atomic data and collisional excitation rates from lamda5 (Schöier et al. 2005) and chianti 6 (Dere et al. 1997;Landi et al. 2013), and we also account for the effects of the cosmic microwave background.The line emission from each cell is then calculated using the level populations.
radmc then produces a 3D data cube in position-position-velocity space, with velocities spanning ±200 km s −1 about the line centre of the emission.This data cube has a spectral resolution of 2 km s −1 and a spatial resolution of 20 pc for each simulated spectra created.
Table 1.Initial conditions of the simulations used within this work.A more complete insight into these is included in Richings et al. (2022). 200, crit is the total halo mass in M ⊙ of the galaxy, which each model is named after,  gas is the gas fraction, and  init is the initial metallicity compared to  ⊙ . 0,tot is the total metallicity at time  = 0 Myr to show the enrichment that occurs during the 300 Myr settling-in period, and  0,gas is the gas-phase metallicity when accounting for dust-depletion using the model described in Section 2.2. 0,tot and  0,gas are mass-weighted averaged over all gas particles in the simulation.

Galaxy formation models
We describe the initial conditions of the simulations in table 1, with a complete insight into these initial conditions available in table 3 of Richings et al. (2022).The simulations range from dwarf galaxy halo masses ( 200, crit = 10 10 M ⊙ ) up to Milky Way-mass ( 200, crit = 10 12 M ⊙ ) where  200, crit is the total halo mass.The total stellar masses were calculated using the abundance matching model from Moster et al. (2013), which has been modified according to Sawala et al. (2015) so the model can account for the inefficiencies of galaxy formation at lower halo masses.We also use two variants of the m3e11 galaxy with a low and high gas fraction, so we can test the effect of gas fraction on our results.Each simulation has a mass resolution of 400 M ⊙ for each gas and star particle.Using the cosmological density parameters for the total matter and baryonic content of the Universe, Ω m and Ω b respectively, we multiply the baryonic particle mass by (Ω m -Ω b )/Ω b to find and set the mass of dark matter particles to 1910 M ⊙ .We apply a constant gravitational softening length of 2.8 pc for star particles, and 1.6 pc for dark matter.For gas particles, we use an adaptive gravitational softening length which we set equal to the mean interparticle spacing, down to a minimum gas softening length of 0.08 pc.For example, at the density threshold for star formation,  H = 10 3 cm −3 , the gas softening length is set to 2.2 pc.
Resolution tests of each simulation were also ran to test how numerical resolution affects our results.Low-resolution versions of each simulation are included with a mass resolution of 3200 M ⊙ , eight times lower resolution than our fiducial model.For high-resolution tests, only the m3e10 galaxy is included due to the extended runtime of the increased resolution combined with the non-equilibrium chemistry running alongside it, and this high-resolution m3e10 was run with 8 times higher mass resolution, at 50 M ⊙ .We include results for these resolution tests in Appendix A.
The initial metallicities7 mentioned in table 1 were set according to the mass-metallicity relation of SDSS galaxies from Andrews & Martini (2013).As these simulations include the injection of metals from winds and SNe (described in Section 2.1), this metallicity will evolve over time.The initial relative abundances between different (2009) throughout this work, where the total solar metallicity is Z ⊙ = 0.0129 metals were assumed to be solar, and the initial He abundance was scaled to be between primordial and solar, according to the total metallicity.
The galaxies are evolved for 800 Myr from the initial conditions.The simulation were first run for a period of 300 Myr with the supernova time-scales drastically reduced, to enable the disk to settle into a quasi-steady state.Within this 300 Myr period, we see an increase in metallicity within our simulations as we can see in Table 1 for both our  0,tot and  0,gas .This is likely due to a sudden burst of star formation within this period, causing the simulation to be injected with metals.
We then used the snapshot at the end of this 300 Myr settling-in period as  = 0 Myr and ran the simulation for a further 500 Myr with the standard models, taking one snapshot every 100 Myr for a total of 5 snapshots per simulation.It is important to mention we do not use the data from the initial 300 Myr period within our analysis or in creating our radmc spectra.
We show the Milky Way-mass galaxy fiducial simulation (m1e12) at 500 Myr in Figure 1, with a Hubble-like image made using the post-processing software FIRE Studio (Gurvich 2022) in the left panel, and the mass-weighted gas temperature distribution in the right panel.

Simulated emission spectra
We run each simulation snapshot through the radmc-3d radiative transfer code mentioned in section 2.3, and we also run the 500 Myr snapshot with the post-processed abundances assuming chemical equilibrium.The final radmc-3d output is given in erg s −1 cm −2 Hz −1 sr −1 , which we can convert into K km s −1 using the equation below: where   is the Boltzmann constant,  is the central wavelength of the emission, and   is the line intensity which is given for each frequency  in the range calculated by radmc from the central line.
We can perform the integral on the right for each position within the 3D data cube we mentioned in section 2.3.After multiplying by the other factors used in the equation, using  = 0.26 cm for the CO (1-0) line emission, allows us to calculate and create a map of the velocity-integrated surface brightness W co like the one we show in Figure 1.
As can be seen in Figure 1, the CO emission follows the spiral structure of the galaxy, as well as within clusters around the arms.For this image, we filter out any surface brightness below 0.1 K km s −1 , but we do not do this for our main analysis and results unless we state otherwise.
As we mentioned in section 1, a portion of the H 2 gas within galaxies is CO-dark and cannot be traced by its emission.Using our simulated emission map, and the simulation data, we can create a trace of this emission against the H 2 column density to show how well our emission traces the gas with CO emission above the 0.1 K km s −1 observational limit.We can also use this to find where gas that would be CO-dark to observers would be in our simulations, and what percentage of each simulation is CO-dark.This fraction of CO-dark gas expected to increase in low-metallicity environments (Amorín et al. 2016;Schruba et al. 2017;Madden et al. 2020), which leads to an underprediction of the H 2 mass when using a constant conversion factor within low-metallicity galaxies.In Figure 1, we show the H 2 column density plot of our Milky Way-mass galaxy m1e12 with the simulated emission we see in Figure 1 traced atop in black outline.This trace is a contour, showing regions within that have velocity-integrated surface brightness > 0.1 K km s −1 .
In Table 2 we show the percentage of CO-dark H 2 with  co < 0.1 K km s −1 , co < 0.25 K km s −1 , and  co < 1 K km s −1 , for each of our galaxy models in their 500 Myr snapshot for the fiducial model with non-equilibrium chemistry, and the model where we assume chemical equilibrium and post-process the abundances.These intensity cuts are based on observations from Leroy et al. (2011) of the Small Magellanic Cloud, which found the SMC to have a 3 intensity threshold of 0.25 K km s −1 , and NGC 6822 to have an intensity threshold of 0.03 K km s −1 .
Recent investigations using the SILCC-zoom simulations found the fractions of CO-dark H 2 gas varying from 15% to 65% when applying an observational detection limit of 0.1 K km s −1 , which increased to 20% to 75% when this detection limit was set to 1 K km s −1 (Seifried et al. 2020).These simulations were in agreement with previous studies by Levrier et al. (2012), which saw that the densest parts of the molecular gas in their simulations were 35% to 40% CO-dark.Another study we can compare to is Richings & Schaye (2016b), which found CO-dark fractions of 76%-86%, in their models with a halo mass matching our m1e11 model.These fractions are not too dissimilar from our calculations for our m1e11 model, although in their work, they used an observational limit of 0.25 K km s −1 , which would explain the slightly higher fraction of CO-dark gas.
Observational limits for the percentage of CO-dark gas tend to vary, with gamma ray emissions showing more than 30% of H 2 is CO-dark (Grenier et al. 2005), and sub-pc resolution observations of OH in the boundary of the Taurus molecular cloud showing the CO-dark gas fractions range from 20% to 80% (Xu et al. 2016).As can be seen in Table 2, our results are in good agreement with both observational data and prior computational work.We also see how the fraction of CO-dark gas decreases as metallicity and gas fraction increase within both our non-equilibrium and equilibrium models, as expected from observational studies (Amorín et al. 2016;Schruba et al. 2017;Madden et al. 2020).

CO to H2 column density relation
In this section we compare our simulation predictions for the relation between the column densities of Previous work in Keating et al. (2020) compared the observed relation between these Milky Way column densities to their simulations, and found their simulations over-produced CO by an order of magnitude at fixed H 2 column densities.To correct this, the shielding length needed to be decreased by a factor of 10 for the column densities to match the observed ratios.The authors suggest that the original discrepancy may have been due to the limited resolution of their simulations.The simulations in Keating et al. (2020) also post-processed their chemical abundances, assuming chemical equilibrium, so we will see if our combination of higher resolution and non-equilibrium chemistry is able to remove this discrepancy.
As can be seen in Figure 2, our non-equilibrium simulations cover the range of CO column densities at fixed H 2 column density from the observational studies.We see how the normalisation of this relation increases as the galaxy halo mass increases.This is due to the increasing metallicity of the simulations as the halo mass increases, which can be seen in prior work (see Hu et al. 2022, Figure 5).As the observations are from the Milky Way, we need to carefully consider which simulations are most appropriate to compare to.The m1e12 simulation has a similar halo mass to the Milky Way, but overproduces CO at a fixed H 2 column density by ≈0.5dex compared to the Milky Way observations.This is likely due to the higher metallicity in m1e12 as can be seen in Table 1.However, the total metallicity of our simulation m1e11 is close to solar metallicity, and reproduces the observed relation between CO and H 2 column densities seen in the Milky Way.We therefore conclude that no further re-scaling of the shielding length (as in Keating et al. 2020) is necessary, and any remaining discrepancies between m1e12 and the Milky Way observations are primarily driven by the metallicity.
In our equilibrium simulations, we see that we predict a higher ratio of CO when compared to H 2 .We therefore conclude that the improved agreement between the observations and our predicted CO and H 2 column densities, without the need for re-scaling of the shielding length, is due to a combination of our high resolution and the inclusion of non-equilibrium chemistry.
Within Appendix A, we show how this column density relation is affected by resolution for our dwarf galaxy m3e10, and see how whilst we match the low-resolution at low column densities, at higher column densities we see that our fiducial model agrees more with our high-resolution version over the low-resolution.
We also compared the  co vs  h 2 relation in the high, fiducial, and low gas fraction variants of m3e11 (not shown) to see if the gas fraction of the simulation has any affect on the resulting ratio of CO and H 2 produced, and found no significant changes.
We further verify our H 2 model by comparing the fraction of molecular hydrogen, f h 2 , against the total neutral gas column density using observations of H 2 absorption lines that have inferred the total neutral gas column density in Figure 3.We collate observations from surveys from the Far Ultraviolet Spectroscopic Explorer (FUSE) of the galactic disk (Shull et al. 2021), the galactic halo (Gillmon et al. 2006), and the LMC and SMC (Tumlinson et al. 2002), as well as measurements of the Milky Way from Wolfire et al. (2008).We group these observations and match them to the simulation that has the closest metallicity to them, similar to that done in Gnedin et al. (2009).These are m1e11 at 0.9 Z ⊙ for the Milky Way observations, and m3e10 at 0.3 Z ⊙ for the LMC observations and SMC observations.We also include m1e12 and compare to the Milky Way observations as whilst m1e11 is closer in metallicity, m1e12 is similar in mass.
We see that whilst there are regions where our simulations match the observations, we also see some disagreement such as for m3e10 at low gas fractions (10 −6 to 10 −5 ).This could be due to how the simulations are measuring the column density compared to the observations, as in the simulations we are simply projecting the column densities onto an image grid, whereas the observations are based on UV absorption lines along a line of sight.This may account for the discrepancies we are seeing between the simulations and the observations.We do see however that our simulations line up with the observed trend of the transition from atomic-to-molecular gas decreasing with increasing metallicity (Gnedin et al. 2009).

xCOLD GASS
Another observational study we can use in comparing our simulations to is the xCOLD GASS survey (Saintonge et al. 2017), which we discussed in section 1.As mentioned before, this survey includes CO measurements from 532 galaxies using the IRAM 30m telescope, measuring the CO (1-0) emission from each galaxy.In this section we will compare our simulation predictions with the directly observed CO luminosity and the inferred H 2 mass from the xCOLD GASS survey.
To calculate the CO luminosity for our simulations, we use an altered version of equation 6 which we use to convert the radmc output into K km s −1 pc 2 .The xCOLD GASS survey calculates the CO luminosity using the observed frequency and luminosity distance (see Saintonge et al. 2017, Equation 2).This luminosity is then corrected to account for the aperture, and includes an error based on the measurement uncertainty, an 8% flux calibration error, and a 15% uncertainty due to the aperture correction.
We plot the 532 xCOLD GASS galaxies against our simulations using the CO luminosity that we calculate using radmc against star formation rate (SFR) in Figure 4. We show results for each 100 Myr snapshot for each galaxy simulation, and we include the results from our non-equilibrium chemistry simulations as well as our equilibrium simulations.Within our simulations, we calculate SFR by averaging the true mass of stars formed over the preceding 50 Myr in the simulations (see Flores Velázquez et al. 2021).xCOLD GASS uses the sum of SFRs from the NUV and a WISE MIR band to probe unobscured and obscured star formation, respectively, which is described in Janowiecki et al. (2017).We also include measurements from Schruba et al. (2012) and Zhou et al. (2021) to fill out the low CO luminosity/low SFR region to allow a comparison to our low-mass galaxies.Readings from Schruba et al. (2012) that are an upper-limit on the CO luminosity are depicted with an arrow indicating they should be lower, similar to those within the xCOLD GASS survey which we mentioned prior.
Whilst our mid-to-high gas-mass galaxies (with high SFR) are consistent with the most CO luminous galaxies seen in the xCOLD GASS observations at same SFR, they are nevertheless ~0.5dex above the average relation seen in the observations.Given the limited sample size spanned by our suite of galaxy populations, it is unclear whether this is a systematic discrepancy or if our simulated galaxies are just extreme examples at fixed SFR.However, despite the high CO luminosities exhibited in our simulations, we have confirmed that our models produce the observed ratio of CO to H2 (based on their column densities, see Section 4.1), and therefore we are also producing an accompanying higher amount of H 2 in these simulations, meaning any calculation of the conversion factor using these galaxies should not be affected by the increased CO within.
Our low gas-mass galaxies, m1e10 and m3e10 lie below the cloud of xCOLD GASS observations.This is due to a lack of observational data in the xCOLD GASS survey at the mass ranges of these low-mass simulations.xCOLD GASS looks at galaxies with stellar masses of M * > 10 9  ⊙ , whereas our low-mass simulations of m1e10 and m3e10 have total stellar masses of 6.6 × 10 6 , and 8.9 × 10 7 respectively.Therefore, these simulations fall far outside of the range of observations from the xCOLD GASS survey, so we would expect them to be below the data in our figure.On top of this, CO emission below 0.2 ⊙ is near undetectable (Leroy et al. 2005;Schruba et al. 2012;Madden et al. 2020), meaning the metallicity range of these galaxies is incredibly hard to probe.xCOLD GASS contains many non-detections of CO at these star formation rates which they set as a 3 upper limit for the CO luminosity by assuming the full-width at half maximum of the Hi emission is equal to the full-width at half maximum of the CO emission for a galaxy.Recent observations with ALMA have been able to detect CO emission at 0.1 ⊙ in the WLM galaxy (Rubio et al. 2015).We can see that measurements from Schruba et al. (2012) are measured at similar SFR as our m3e10 galaxy, with upper-limits on CO luminosity measurements placed above where we calculate the luminosity for our fiducial model to be, but slightly below where our equilibrium chemistry model is.More observations at these low-metallicity ranges will also help us better understand the nature of molecular gas in the low-metallicity regime by lowering these upper bound measurements Observations which are upper-limits on the H 2 column density, and therefore the molecular gas fraction as well are shown with black arrows.The grayscale corresponds to the density of the plotted bins.We see our simulations are able to match the observations as well as the observed transitions from atomic-tomolecular gas.  . co against SFR for each galaxy for each snapshot, with the non-equilibrium chemistry results in black diamonds, and the equilibrium chemistry results in white pentagons.The xCOLD GASS survey points are included in the blue, with the measurements of galaxies with no CO (1-0) emission detected shown in the yellow arrows, which are set as a 3 upper limit on the CO luminosity.We also include measurements from Schruba et al. (2012) in green and Zhou et al. (2021) in red to show readings in the low CO luminosity range to compare to our low-mass galaxies.Our predicted CO luminosities for our mid-to-high mass galaxies are able to match objects seen in the xCOLD GASS survey at fixed SFR, with our high-mass galaxies being ~0.5dex above the average luminosity for the xCOLD GASS survey at fixed SFR.
from the xCOLD GASS survey, as well as those within Schruba et al. (2012).
The chemical modelling in our simulations allows us to directly predict the H 2 mass.In contrast, the H 2 masses in xCOLD GASS are inferred from the CO luminosity, assuming the multivariate  co conversion factor from Accurso et al. (2017), which requires the metallicity and specific star formation rate of the observed galaxies.This Accurso function has a 35% uncertainty within, and when combined with the error on the CO luminosity mentioned earlier gives a great uncertainty on this H 2 mass.On top of that, a portion of the xCOLD GASS galaxies were recorded with a non-detection of the CO (1-0) line, and therefore the values recorded for this are placed as a 3 upper limit on the H 2 gas mass reading.
We plot the H 2 mass from the xCOLD GASS galaxies and our nonequilibrium and equilibrium simulations against SFR in Figure 5, and we highlight those galaxies with readings which are a 3 upper limit on the H 2 mass.We also plot the H 2 mass which we calculate using our simulated CO luminosity and the Accurso function we showed in equation 3.This is to show what an observer would calculate our simulated galaxies H 2 mass at, as the xCOLD GASS survey has done.
As before, our mid-to-high gas-mass galaxies m3e11, m3e11_hiGas, and m1e12 are consistent with observed xCOLD GASS galaxies.As before, we see many non-detections with SFRs similar to our mid-mass galaxies.With CO (1-0) measurements of these galaxies which have their CO luminosities marked as a 3 upper limit by the xCOLD GASS survey, those observations will fall, though there is no lower constraints on the value of these galaxies yet meaning we do not yet know if they will lie within the bounds of our mid-gas mass simulations.However, our simulations line up with the lower edge of the envelope covered by the xCOLD GASS objects, when our calculated CO luminosities for those simulations were on Figure 5.  h2 against SFR for each galaxy for each snapshot with the nonequilibrium fiducial models in black, and the snapshots with equilibrium chemistry shown in white.Each galaxy mass groups by star formation rate, with m1e10 with the lowest SFR and m1e12 at the highest.The blue rings show xCOLD GASS survey points from galaxies where the CO (1-0) line emission was detected.Each calculated  h2 has an error associated with it, which includes the error from the observed CO (1-0) luminosity and the error from the Accurso conversion function used to calculate  h2 .The yellow arrows are a 3 upper limit on the  h2 from the galaxies with no CO (1-0) emission detected.The red hexagons show the  h2 we calculate for each simulation using the Accurso multivariate function we show in equation 2.
We see that the Accurso function overpredicts the H 2 mass of our simulations by an order of magnitude.
the upper edge of the observations.We can also see our calculated H 2 mass using the Accurso et al. ( 2017) multivariate function fitting an order of magnitude higher than our actual H 2 mass for these simulations.Therefore if an observer was to measure our simulated galaxy and its CO luminosity, the calculated H 2 mass would be an order of magnitude too high if they used the Accurso et al. ( 2017)  co factor.As this multivariate conversion factor was used within the xCOLD GASS survey, it is important that it can accurately estimate the H 2 mass of an observed galaxy.For these reasons, we discuss the Accurso function, and a possible replacement fitting further in section 5.At the bottom of these results once again is our low-mass galaxies m3e10 and m1e10.
One possible reason for these discrepancies in  co and  h 2 could be the star formation model used within the simulations, as xCOLD GASS measures CO luminosities and calculates H 2 gas mass equal to what we find in the simulations.However, prior work by Orr et al. (2018) verified the star formation model within the FIRE simulations, as they were able to reproduce a KS-like relation with a slope and scatter that matched observational data.On top of that, a change in the star formation we calculate would not make our results any closer to observations overall, as we require an increase in SFR to match the observations more closely for  co , yet a decrease in SFR to match for  h 2 .Therefore, we can rule out the star formation model as the source of the discrepancies.
Within Appendix A, we analyse whether the numerical resolution has any effect on our simulated CO luminosity and H2 mass with our low and high resolution tests, and we compare these to the xCOLD GASS survey.We find similar results for our low resolution tests to our fiducial model, although the results for the high-resolution test of m3e10 is limited by the bursty nature of the dwarf galaxy models (see Appendix A for further discussion).

CALIBRATING THE CO-TO-H2 CONVERSION FACTOR
In section 1 we discussed the CO-to-H 2 conversion factor,  co , and how we can use it to estimate molecular gas content within galaxies using CO emission.Using our simulated CO emission, we are able to calibrate the conversion factor for each simulation as we have the exact H 2 mass available.

Metallicity relation
We discussed the relation between  co and metallicity earlier i.e. how when metallicity increases, the conversion factor decreases.Bolatto et al. (2013) gave an overview of dust-based observations, showing the decreasing  co with increased metallicity, as well as how  co rapidly increases below metallicities between 1/3 to 1/2 of solar.As we simulate galaxies across a range of metallicities, we are able to test this relation using our simulations.However, our low metallicity galaxies have extremely bursty star formation histories.This bursty star formation can cause the gas fraction of the disc to vary by a large factor (see Richings et al. 2022, Figure 7), which can affect how much dense gas is within the galaxy, which is needed to form CO efficiently.We compare to the dust-based approaches from Israel (1997); Leroy et al. (2011) and Sandstrom et al. (2013).For each galaxy snapshot, we calculate the average  co across the galaxy and plot those in Figure 6 against the dust-based observations for both our nonequilibrium simulations, and those assuming chemical equilibrium.We also plot metallicity-dependent relations for the conversion factor with coloured dashed lines from Schruba et al. (2012); Amorín et al. (2016); Accurso et al. (2017) and Madden et al. (2020).
We can see that both our non-equilibrium and equilibrium simulations follow the trend shown in the dust-based observations.For our low-metallicity galaxies like m1e10, which have highly bursty star formation,  co varies greatly from snapshot to snapshot.However, even with its bursty nature, the calculated  co for our low-metallicity galaxies are still in line with expectations of  co at low metallicities.
We see in Figure 6 how our non-equilibrium and equilibrium results for  co appear similar.In Table 3, we show the non-equilibrium effects on  co and  h 2 when compared to our equilibrium simulations.We then also quantify the non-equilibrium effects on  co .
We see that non-equilibrium chemistry typically reduces  co and  h 2 by ≈ 70-76%, atleast in the mid-to-high mass galaxies.However, as they are both reduced by a similar degree, the resulting  co factor is only weakly affected by a few percent.For our low-mass galaxies, we note that these simulations exhibit highly bursty star formation histories and gas fractions, so it is unclear whether the Figure 6. co against gas-phase metallicity for each snapshot with the nonequilibrium fiducial models in blue, and the snapshots with equilibrium chemistry shown in white.We plot dust-based observational measurements from Israel (1997); Leroy et al. (2011) and Sandstrom et al. (2013), as well as the  co for the Milky-Way with the blue-dotted line.We also plot metallicitydependent relations for the conversion factor with coloured dashed lines from Schruba et al. (2012); Amorín et al. (2016); Accurso et al. (2017) and Madden et al. (2020).We show metallicity when compared to solar metallicity, Z ⊙ .Our simulations reproduce the observed trends between  co and metallicity for both our non-equilibrium and equilibrium models.
non-equilibrium effects suggested in Table 3 for low-mass galaxies will hold generally or may just be stochastic.
This result of seeing only a small fluctuation of  co between nonequilibrium to equilibrium contrasts previous studies by Richings & Schaye (2016a), where the average  co within their simulations could be affected up to a factor of ~2.3.However, these simulations did not include a fluctuating UV radiation field, as well as link the UV field to the star particles, which we would expect to increase the results within non-equilibrium simulations.Our simulations also use different subgrid models, such as those for star formation and stellar feedback.It is possible the inclusion of one of these has driven additional effects, counteracting the non-equilibrium effects in the simulations, leaving only small fluctuations in the  co between nonequilibrium to equilibrium simulations for our mid-to-high mass galaxies.Further study is required to know which is causing the balance between  co in non-equilibrium and equilibrium simulations compared to previous studies.

Multivariate conversion factor
Figure 6 demonstrates that  co is not a constant, and the conversion factor depends on conditions within the galaxy.We see that metallicity is one such factor that causes fluctuations in the conversion factor, but Accurso et al. (2017) demonstrated that other factors are also important by examining observed scaling relations between the ratio of luminosity of [C ii] to CO and local conditions from galaxies from the xCOLD GASS survey and the Dwarf Galaxy Survey (DGS) (Madden et al. 2013;Rémy-Ruyer et al. 2014;Cormier et al. 2014).The factors that accounted for most of the variation in the luminosity ratio were metallicity, and offset from the main sequence, ΔMS.Performing a linear fit to their radiative transfer models, they found the multivariate function we showed in equation 3. The version we show in equation 3 uses the log of the ratio of the luminosity of [C ii] to CO, though there is a version without, instead using a constant in its place which we show below: The substitution of the log of the luminosity ratio was achieved by performing a linear fit using those observed scaling relations to find an equation for the luminosity ratio in terms of metallicity and ΔMS, and then substituting out the luminosity ratio.In Figure 7 we show the observations used to find these scaling relations alongside our simulations.For the [C ii] luminosity, we use the results found in Richings et al. (2022) to calculate the luminosity ratio with our  co calculated within this work.Richings et al. (2022) found that nonequilibrium chemistry only affects their calculated  [C ii] by ≈ 5%, and with our results in table 3 showing only a fluctuation of ≈ 25% for the CO luminosity and 5 % for the conversion factor, we only focus on our fiducial model results when calibrating our multivariate fitting.
The version of the function shown in equation 7 is used by the xCOLD GASS survey to calculate the H 2 mass for each galaxy they observe.As we showed in Figure 5, this multivariate function overpredicts the mass of our galaxy models by up to an order of magnitude.For this reason, we use the method used in Accurso et al. (2017) to find a fitting from our simulations.
Using equation 8 as a baseline with constants , , , and , we show the constants for both versions of the Accurso function, and both versions of our fitting function in Table 4.We see for the luminosity ratio functions, our function sees more of a reliance on metallicity and less on ΔMS.We also see that our function matches the observations with us seeing a negative correlation on metallicity, though we do not see the positive relation on ΔMS due to not probing a full range of ΔMS values within our simulations, which we can see in Figure 7.We see that when including results from Cormier et al. (2015), our low-metallicity observations are within reason.We see that some simulations are also outside the region of the observations, which we will discuss shortly.The dependence of our fitting on the luminosity ratio is also reduced by a factor of ~30, with a greater dependence on the constant .For both functions where  = 0, we see close agreement for  and , though even slight variations in these will cause large variations in the final  co .
We plot the H 2 mass of the xCOLD GASS survey in Figure 8 against the H 2 mass calculated with our alternative fitting function with constants  = 0,  = -2.0362, = -0.0123, = 17.3036.We also plot both our non-equilibrium and equilibrium simulations, as well as the H 2 mass of our simulations calculated using the Accurso function.We can see our fitting estimates the H 2 mass by an order of magnitude below that which was calculated by the Accurso function for the xCOLD GASS survey.We also see our fitting now matches our  co predictions when compared to xCOLD GASS in Figure 5, where our simulations were in the upper edge of the observations.
We caution that the fitting from our simulations is uncertain, as

LCII/LCO
Luminosity ratio against ΔMS Figure 7. Scaling relations of  [C ii] / co against gas-phase metallicity (left panel) and ΔMS.Observations taken from a sample of 23 xCOLD GASS objects in blue, and 7 Dwarf galaxy survey (DGS) objects in red, as well as a sample of local dwarf galaxies from Cormier et al. (2015) in green.We plot our non-equilibrium fiducial models against in black.Whilst the simulations are broadly consistent with the observational data, the limited sample size in the simulations means they do not fully probe the complete distribution covered by the observations.Table 4. Comparison of the constants , , , and  between both variations of the Accurso multivariate function, and both variations of the alternative multivariate function we propose in this work.The two functions where there is a non-zero constant  are the functions where the ratio between luminosity of [C ii] to CO is used, and the two functions where  is set to 0 are those where that ratio has been substituted out.The source column shows whether the equation comes from Accurso et al. (2017)  the limited range of galaxy conditions spanned by our simulation suite does not fully probe the complete distribution of metallicities and ΔMS values observed in the xCOLD GASS survey and DGS, as seen in Figure 7.This may bias the results of these fits.With a greater range of simulations and more observational data, we would be able to see if our simulations truly match the relations between the luminosity ratio, metallicity, and ΔMS.

Variation with spatial scale
We have looked at how  co varies with factors such as metallicity or non-equilibrium chemistry effects, but one factor we also need to look at is spatial scale.Observations are not guaranteed to be as high-resolution as our simulations, so we need to test whether our simulated  co varies largely when adapting our simulated spectra to match observations more closely.We can investigate the variation of  co with different spatial scales through differing ways, the first is by using a Gaussian filter to blur the finer details of our simulations to better match observations.We show how this Gaussian filter affects the H 2 column density of our Milky-Way mass galaxy m1e12 as we increase the standard deviation of the filter, , in Figure 9.We also apply this filter to our velocity-integrated surface brightness  co , which we see in the Figure 8.  h2 against star formation rate.We show the  h2 calculated within the xCOLD GASS survey using the Accurso multivariate function in blue, and compare to the  h2 of the xCOLD GASS objects calculated using our alternative fitting function in purple.We also plot each of our snapshots with the non-equilibrium chemistry in black, the equilibrium chemistry in white, and the mass of our snapshots calculated using the Accurso function in red.
bottom half of Figure 9.We can see how the gas is smoothed out from  = 0 pc, with no blur applied, up to  = 2000 pc, where the galaxy is completely smoothed and all features have been lost.This Gaussian filter is in addition to the gas particles being smoothed over their kernel size with a cubic spline kernel.We show how the average  co changes when the Gaussian filter 'blur' is increased in Figure 10 for our Milky-Way mass galaxy m1e12 for each 100 Myr snapshot.We average the H 2 and  co on a spatial scale that matches the full-width at half-maximum (FWHM) of the Gaussian filter that has been applied.We can see that despite some fluctuations in the  co value, the average value for each snapshot increases by less than 0.5 dex from the initial average value at  = 0  pc where no Gaussian blur is applied.The only exception to this is snapshot 100, which deviates by an order of magnitude when the Gaussian blur is at a FWHM of ~3000 pc.The other deviations may be due to just slight variations in the H 2 gas distribution or  co causing much lower than expected values, which we can see in Figure 9 where  = 40 pc and  = 200 pc, where the minimum values drop to much lower than without a blur by several orders of magnitude.We can see how the scatter on  co decreases on larger spatial scales as well (see also Feldmann et al. 2012b).This is due to averaging on much larger spatial scales, meaning the resolution of our spectra decreases as  and the FWHM of our Gaussian filter increase.

Variation with other physical parameters
In Section 5.2, we showed how the conversion factor is multivariate, but we only show metallicity and offset from main sequence as these values are easily available for an observer, and therefore our function is easily usable to calculate the conversion factor for observers.However, there are other parameters where we see variations in  co , such as molecular gas surface density, Σ h 2 .
Prior work by Shetty et al. (2011a,b) found that line saturation from lines of sight with high surface densities caused  co to rise.Feldmann et al. (2011) also found that on small scales,  co increases with Σ h 2 .However, Narayanan et al. (2011) found the opposite trend, and argued that as Σ h 2 increases,  co decreases due to an increase in velocity dispersion and kinetic temperature caused by more efficient dust-gas thermal exchange at higher densities.This increase in both the velocity dispersion and kinetic temperature caused an increase in CO intensity, and therefore a decrease in  co .This decrease occurs at the same range as the increase noted in Shetty et al. (2011a,b) (see Figure 6 in Narayanan et al. 2011 andFigure 4 in Shetty et al. 2011b).
We show  co against Σ h 2 for our galaxy m1e12 using data from all 5 snapshots taken at 100 Myr intervals in Fig. 11.We see fluctuations in  co against Σ h 2 , with a near constant  co between surface densities of 10 −6 to 100 M ⊙ pc −2 , but a sharp increase for surface densities > 100 M ⊙ pc −2 .This sharp increase aligns with the results  We see how areas with high-surface density and high-FUV flux correspond to a conversion factor similar to the value we find when averaging the  co of galaxy m1e12.
found in Shetty et al. (2011b), as they also found an increase at surface densities > 100 M ⊙ pc −2 .This is because at high surface densities, the CO line is saturated, so the H 2 column density increases, as does the  co value.We can see that we see a similar upturn at 100 M ⊙ pc −2 as seen in Shetty et al. (2011b) which started at a similar gas surface density, which was attributed to the saturation of the CO line.
We can also look at the strength of the far UV flux, which will correlate with star formation to see how  co changes when measuring the strength of the FUV flux.We show  co against FUV flux, and colour each point by Σ h 2 for galaxy m1e12 at time  = 500 Myr in Figure 12.We see how  co changes with both Σ h 2 and the flux of the FUV field.We can also see that when the FUV flux and surface density are at their highest, the conversion factor matches what we measure on average in our simulations, showing how the conversion factor is dominated by areas of both high flux and molecular gas surface density, which corresponds to dense star forming molecular clouds within the galaxy.We also see how in regions with low FUV flux and low star formation,  co has a large scatter.However, as the FUV flux and star formation increase, this scatter is reduced and we see how the conversion factor is compacted around the average value we measure for the system.

CONCLUSION
We have analysed the CO emission from a suite of simulations based on the FIRE-2 subgrid models, and coupled them with the CHIMES chemistry module to evaluate the on-the-fly non-equilibrium chemistry against post-processing the abundances assuming chemical equilibrium.We then take snapshots from these galaxy models, and post-process through the radmc radiative transfer code, creating simulated spectra of the CO emission from our galaxy models.We find that: (i) Our non-equilibrium fiducial models are able to match the relation seen in the Milky Way between the column densities of H 2 and CO.We show this in Figure 2. We note that our Milky Way-mass galaxy m1e12 does not match the observations completely, as the median relation is ≈ 0.5 dex above the observations.However we believe this is due to the large metallicity increase that occurs within the initial 300 Myr of evolution of the model.We also further verify our H 2 model by comparing the fraction of molecular gas to the total neutral gas column density in Figure 3. (ii) Our simulations are able to match the CO luminosities as a function of star formation rate seen within the xCOLD GASS survey, which show in Figure 4. We see our high-mass galaxies are ~0.5 dex above the average for the xCOLD GASS survey, though still match to objects on the upper edge of the data for fixed SFR.For our lowmetallicity galaxies outside the range of the xCOLD GASS survey, we note that they measure galaxies with SFR similar to our low-mass galaxy models and calculate a 3 upper limit on the CO luminosity for those.Future observations may be able to lower these upper limits and show whether our low-mass galaxies truly match.(iii) We plot our simulations against dust-based observations to evaluate the relation between metallicity and  co in Figure 6.Our simulations follow the expected trend even at lower-metallicities, though we are unable to verify due to a lack of observational data due to the difficulties of measuring CO at low metallicity.(iv) We evaluate the multivariate conversion factor from Accurso et al.
(2017), which in Figure 5 we see overpredicts the H 2 mass of our galaxy simulations by an order of magnitude.We follow the method used in Accurso et al. (2017) to find a new fitting to our simulations, which we show in Table 4.However, as we do not probe a full range of [C ii] luminosity and ΔMS ranges, we see our simulations show a strange distribution when compared to a sample of the xCOLD GASS survey and DGS in Figure 7 as some snapshots lay outside the observational data.For this reason, more simulations and/or observations are needed to reduce these caveats with our fitting function.(v) We show how  co is mostly constant when applying a Gaussian filter across differing spatial scales to better match our simulated spectra to observations in Figure 10.These slight variations we see in  co as the Gaussian filter FWHM increases are likely due to the Gaussian filter spreading H 2 gas or  co to regions where we see no gas or CO emission without the filter, which we can see in Figure 9. (vi) We compare  co to the surface density and the FUV flux, and show how  co within our simulations fluctuates with both in Figure 12.We also note that when the surface density and FUV flux is at its greatest,  co aligns with the conversion factor we measure in our m1e12 Milky Way-mass simulations (~10 19 -10 20 cm −2 (K km s −1 ) −1 ).(vii) We quantify the effects non-equilibrium chemistry has on our galaxy simulations in table 3. Whilst  co and  h 2 are reduced in equilibrium simulations by 24% to 30%, in higher-mass simulations these reduction factors are close, meaning  co in high-mass simulations is unaffected (<5%) between equilibrium and non-equilibrium simulations Within this work, we have demonstrated that modelling nonequilibrium chemistry and synthetic emission line observations in hydrodynamical simulations is a promising avenue, allowing us to understand the CO-to-H 2 conversion factor in different environments.Whilst we see little difference in our  co estimates (<5%) for our high-mass galaxies, we see how non-equilibrium chemistry can affect the simulated emission line as well as the molecular gas content of the models.However, this study is limited by our relatively small range of galaxy parameters spanned by our simulation suite.In future work, we will therefore expand the range of our simulations.One such way would be to look at  co in galaxy mergers, as we have only looked at isolated galaxies within this work, as this may reveal more about the nature of the conversion factor in extreme conditions.Collating the merger simulations, as well as our current simulations, could probe a greater parameter space than our current suite of simulations.
Another avenue for future work would be the use of machine learning to predict properties of galaxies and clouds.Recent work by Garcia et al. (2023) designed slick, a package that calculates realistic CO luminosities for clouds and galaxies in hydrodynamical simulations using a model found through random forest machine learning methods.In future work, we will explore how machine learning techniques can be used to model the CO-to-H 2 conversion factor.

APPENDIX A: RESOLUTION TESTS
As mentioned in section 3, we ran resolution tests for our galaxy models.We ran low-resolution versions of each galaxy model, with 8x lower resolution, and a high-resolution version of m3e10 with 8x higher resolution.
We have shown results for how the CO emission traces the H 2 gas with the CO-dark ( co > 0.1 K km s −1 ) H 2 gas percentages in table 2 for our standard resolution simulations with non-equilibrium chemistry and post-processed equilibrium chemistry.We can now look at how these percentages change when we alter the resolution of the simulation in Table A1, where we show the percentage of CO-dark ( co > 0.1 K km s −1 ) H 2 gas within the low-resolution simulations, and the high-resolution version of m3e10.Unlike for our non-equilibrium and equilibrium simulations, we do not run resolution tests for the m3e11 gas fraction variants.
Using these results, we can see how resolution can affect the amount of CO-dark gas within these simulations, as within the m1e10 galaxy we get no CO emission above our 0.1 K km s −1 detection threshold, and so the entire galaxy is classed as CO-dark gas.Within the low-resolution versions of m3e10 and m1e11 we see that we are predicting ~20% less CO-dark gas compared to our non-equilibrium simulations.We also see that there is no general trend of CO-dark gas with resolution for our m3e10 galaxies, which suggests that the differences in these dwarf galaxies are driven partly by stochastic effects.Finally, in the low-resolution versions of m3e11 and m1e12, we see that we nearly match the amount of CO-dark gas within our non-equilibrium simulations in our low-resolution versions, differing by only 2% to 6%.
We can look at how resolution can change our H 2 mass, as well as the luminosity of our simulated CO emission, and the resulting  co calculated from the two.We compare these in Table A2, where we can show the effects of resolution on our simulations for our lowresolution simulations, as well as our high resolution run of m3e10.We average over all 5 snapshots for each galaxy to reduce the effect the bursty nature of the dwarf galaxies has on the results.
In our dwarf galaxy m1e10, we see a reduction on all values when compared to our fiducial model values.However, as we mentioned prior, m1e10 exhibits extremely bursty star formation histories and gas fractions, which we see when we compare our fiducial model to our equilibrium simulations as well.We see in Table A1 that we see no CO emission in our snapshot taken at time  = 500 Myr for m1e10, meaning we are producing little-to-no CO emission at times, causing the exaggerated reductions.In our high-mass galaxies, m3e11 and m1e12, we see that there is only a small fluctuation between our fiducial model and the low-resolution tests for  co , and a reduction of ≈ 12% to 22% for the  co and  h 2 .For our high-resolution simulation of m3e10, we see that  co is enhanced by over a factor of 2, and  h 2 is reduced by ~22%, causing a greater reduction in the final  co down to 35% the value we see in our fiducial simulations.Like m1e10, m3e10 also exhibits bursty star formation, which may contribute to this result.
We can also plot the relation between  co vs  h 2 for our resolution tests in Figure A1, showing how our low-resolution tests compare to the observations as we did for our fiducial models in Section 4.1.We see that our low-resolution models also match the observed relation within the Milky Way between column densities.As our low-resolution version of our Milky Way-mass galaxy m1e12 at time  = 500 Myr is closer to solar metallicity than our nonequilibrium fiducial model, we do see that it is ~0.25dex below our fiducial model.
We also run a version of this plot comparing both our resolution tests for m3e10 to our fiducial non-equilibrium model in Figure A2.We see that for a fixed H 2 column density at  h 2 ~10 20 cm −2 , we produce CO column density an order of magnitude above our fiducial model in our high-resolution version of m3e10.However, at higher H 2 column densities ( h 2 >~10 21 cm −2 ), we see that our fiducial model agrees more with our high-resolution version over the lowresolution.We also see in Table A2 that whilst we are producing a higher ratio of CO to H 2 in this resolution test, that our  co is only affected by a factor of 2.3, and the final  co is affected by a factor of 3, which both are similar to other resolution tests, such as our m1e11 low-resolution test.This smaller difference in the overall CO luminosity is because it is dominated by the high column density regions.
We can also compare our resolution tests to the observed CO luminosities of the xCOLD GASS survey in Figure A3.We see that our mid-to-high mass runs (with high SFR) match observed luminosities of objects with similar SFR in the xCOLD GASS survey.Our low-mass galaxies, and therefore our high-resolution runs of m3e10, lie below the xCOLD GASS survey, as we expect and discuss in Section 4.2.
We also compare our resolution tests against the relation between  co and metallicity in Figure A4, as well as to our fiducial model.We see that our resolution tests match the observed relation between  co and metallicity.
As we see in our results, we see multiple differences in the highresolution version of m3e10 when compared to our low-resolution and fiducial models.However, we believe these differences are driven by the gas content and star formation history of this high-resolution model, which we can see in Figure A3, where our high-resolution m3e10 has a much higher SFR than the low-resolution version.We also note that our non-equilibrium fiducial model and high-resolution test differ by a factor of 2 for the CO luminosity, and the conversion factor.This is to be expected, as we see large fluctuations from snapshot to snapshot for these dwarf galaxies in our simulations, which can be seen in both Figure A3 and Figure A4 where we see large fluctuations within snapshots of the same resolution.Therefore, we also believe part of these differences to be due to the bursty nature of the dwarf galaxies.Further study into the dwarf regime is required to understand this further.
This paper has been typeset from a T E X/L A T E X file prepared by the author. . co against gas-phase metallicity for each snapshot with the low-resolution tests in grey, and the m3e10 high-resolution tests in white.We also plot our non-equilibrium fiducial results in blue hexagons.We plot dustbased observational measurements from Israel (1997); Leroy et al. (2011) and Sandstrom et al. (2013), as well as the  co for the Milky-Way with the bluedotted line.We also plot metallicity-dependent relations for the conversion factor with coloured dashed lines from Schruba et al. (2012); Amorín et al. (2016); Accurso et al. (2017), andMadden et al. (2020).We show metallicity when compared to solar metallicity, Z ⊙ .

-Figure 1 .
Figure 1.Galaxy model m1e12 at time  = 500 Myr shown in various ways Top left: using a mock Hubble image created using post-processing the simulations through FIRE Studio.Top right: The mass-weighted temperature distribution of the gas.Bottom left: the velocity-integrated surface brightness of the CO emission made using radmc.Bottom right: Trace of the CO emission shown compared to the H 2 column density, showing the CO-dark gas.This trace is a contour showing all regions within that have velocity-integrated surface brightness > 0.1 K km s −1 .
CO and H 2 to observations compiled from UV absorption lines within the Milky Way.These observational surveys includeFederman et al. (1990);Rachford et al. (2002);Crenny & Federman (2004);Sheffer et al. (2008), andBurgh et al. (2010).We show the readings from these observational studies, as well as data from the simulations we use in this work within Figure2.The shaded regions indicate the 15th to 85th percentile within each H 2 column density bin, while the solid curve shows the median, taken from time  = 500 Myr within each simulation.The left-and right-hand panels show our simulation results using non-equilibrium and equilibrium abundances, respectively.

Figure 2 .
Figure2.Relation between CO and H 2 column densities for the range of halo masses.The shaded coloured regions indicate the 15th to 85th percentile within each H 2 column density bin, while the solid curves show the median, for each halo mass at time  = 500 Myr.The black points included show observational data from UV absorption lines in the Milky Way fromFederman et al. (1990);Rachford et al. (2002);Crenny & Federman (2004);Sheffer et al. (2008), andBurgh et al. (2010).

Figure 3 .
Figure 3.The molecular gas fraction against the total neutral gas column density for Top: m1e12 at time  = 500 Myr compared to measurements from the Milky Way.Middle: m1e11 at time  = 500 Myr compared to measurements from the Milky Way.Bottom: m3e10 at time  = 500 Myr compared to measurements from the LMC and SMC from FUSE.Milky Way observations collated from Wolfire et al. (2008); Gillmon et al. (2006) and Shull et al. (2021).LMC and SMC observations fromTumlinson et al. (2002).Observations which are upper-limits on the H 2 column density, and therefore the molecular gas fraction as well are shown with black arrows.The grayscale corresponds to the density of the plotted bins.We see our simulations are able to match the observations as well as the observed transitions from atomic-tomolecular gas.

Figure 4
Figure 4.  co against SFR for each galaxy for each snapshot, with the non-equilibrium chemistry results in black diamonds, and the equilibrium chemistry results in white pentagons.The xCOLD GASS survey points are included in the blue, with the measurements of galaxies with no CO (1-0) emission detected shown in the yellow arrows, which are set as a 3 upper limit on the CO luminosity.We also include measurements fromSchruba et al. (2012) in green andZhou et al. (2021) in red to show readings in the low CO luminosity range to compare to our low-mass galaxies.Our predicted CO luminosities for our mid-to-high mass galaxies are able to match objects seen in the xCOLD GASS survey at fixed SFR, with our high-mass galaxies being ~0.5dex above the average luminosity for the xCOLD GASS survey at fixed SFR.
Figure 9: Galaxy m1e12 at time  = 500 Myr, with a Gaussian filter applied, where the standard deviation of the filter, , is increasing from left to right from 0 pc, 40 pc, 200 pc, up to 2000

Figure 10 .
Figure10.Relation between  co and the spatial scale of which the simulation has then been averaged over, which matches the FWHM of the Gaussian filter, for each snapshot taken every 100 Myr of galaxy m1e12 for our non-equilibrium fiducial model.The shaded region shows the 15th to 85th percentile of  co values, and the solid line shows the median, for each snapshot.We see how  co is near constant on all spatial scales, with some slight fluctuations.

Figure 11 .
Figure 11.Relation between  co and Σ h2 for galaxy m1e12 for our nonequilibrium fiducial model.The shaded region shows the 15th to 85th percentile of  co , and the solid line shows the median.We also plot a dotted line at 10 2 M ⊙ pc −2 to show the upturn occurs at a similar range to Shetty et al. (2011b).

Figure 12 .
Figure 12.  co against FUV flux for each grid point with resolution of 20pc for galaxy m1e12 at time  = 500 Myr, where we colour each point by Σ h2 .We see how areas with high-surface density and high-FUV flux correspond to a conversion factor similar to the value we find when averaging the  co of galaxy m1e12.

Figure A1 .Figure A2 .
Figure A1.Relation between CO and H 2 column densities for the range of halo masses of our low-resolution tests.The shaded coloured regions indicate the 15th to 85th percentile within each H 2 column density bin, while the solid curves show the median, for each halo mass at time  = 500 Myr.The black points included show observational data from UV absorption lines in the Milky Way from Federman et al. (1990); Rachford et al. (2002); Crenny & Federman (2004); Sheffer et al. (2008), and Burgh et al. (2010).

Figure A3 .
Figure A3. co against SFR for each galaxy for each snapshot, with the low-resolution tests in grey diamonds, and the m3e10 high-resolution tests in white pentagons.The xCOLD GASS survey points are included in the blue, with the measurements of galaxies with no CO (1-0) emission detected shown in the yellow arrows, which are set as a 3 upper limit on the CO luminosity.

Table 2 .
The percentage of H 2 gas which is considered CO-dark for each intensity cut ( co < 0.1 K km s −1 ,  co < 0.25 K km s −1 ,  co < 1 K km s −1 ) within each galaxy at 500 Myr for the non-equilibrium (non-Eqm) and equilibrium (Eqm) simulations.

Table 3 .
Ratios of CO emission line luminosity, H 2 mass, and  co ,  noneq / eqm , at 500 Myr using the fiducial model.Values highlighted in red and blue correspond to an enhancement and reduction, respectively, of these values when non-equilibrium abundances are used compared to postprocessed equilibrium abundances.

Table A1 .
The percentage of H 2 gas which is considered CO-dark ( co > 0.1 K km s −1 ) within each galaxy at 500 Myr for our fiducial model with non-equilibrium chemistry (Fiducial), the low resolution (low_res) and high resolution simulations (hi_res).

Table A2 .
Ratios of CO emission line luminosity, H 2 mass, and  co ,  res / noneq for the resolution tests compared to the non-equilibrium fiducial model, comparing the average across both simulations.ResType tells us whether we are comparing with the low resolution (low_res) or high resolution (hi_res) version of each model.Values highlighted in red and blue correspond to an enhancement and reduction, respectively.