On the distribution of the Cold Neutral Medium in galaxy discs

The Cold Neutral Medium (CNM) is an important part of the galactic gas cycle and a precondition for the formation of molecular and star forming gas, yet its distribution is still not fully understood. In this work we present extremely high resolution simulations of spiral galaxies with time-dependent chemistry such that we can track the formation of the CNM, its distribution within the galaxy, and its correlation with star formation. We find no strong radial dependence between the CNM fraction and total HI due to the decreasing interstellar radiation field counterbalancing the decreasing gas column density at larger galactic radii.However, the CNM fraction does increase in spiral arms where the CNM distribution is clumpy, rather than continuous, overlapping more closely with H2. The CNM doesn't extend out radially as far as HI, and the vertical scale height is smaller in the outer galaxy compared to HI with no flaring. The CNM column density scales with total midplane pressure and disappears from the gas phase below values of PT/kB =1000 K/cm3. We find that the star formation rate density follows a similar scaling law with CNM column density to the total gas Kennicutt-Schmidt law. In the outer galaxy we produce realistic vertical velocity dispersions in the HI purely from galactic dynamics but our models do not predict CNM at the extremely large radii observed in HI absorption studies of the Milky Way. We suggest that extended spiral arms might produce isolated clumps of CNM at these radii.


INTRODUCTION
The interstellar medium (ISM) exists in a complex, multi-phase, form from hot ionised gas, to cold molecular star-forming clouds (e.g. Mc-Kee & Ostriker 1977;Draine 2011;Tielens 2005;Klessen & Glover 2016). Matter cycles through these phases as stars formed in the cold dense molecular clouds release energetic feedback and momentum into their surroundings, which then influences the subsequent evolution of the galaxy by setting the equilibrium disc structure and depletion time.
A crucial component in this phase matter cycle is the cold neutral medium or CNM. This consists of the neutral atomic hydrogen (H ) with temperatures around 100 K (Kulkarni & Heiles 1987;Dickey & Lockman 1990) that makes up the bulk of the neutral gas in galaxies, alongside the warm neutral medium (WNM) which has temperatures of order 10 4 K. These atomic phases exist together in pressure equilibrium such that they can be considered as a two-phase medium (Field et al. 1969;Wolfire et al. 2003;Bialy & Sternberg 2019). It is from the CNM that gas is compressed and cooled to form individual molecular clouds where stars are born (e.g. McKee & Ostriker 2007; ★ E-mail: rowan.smith@manchester.ac.uk Girichidis et al. 2020). The CNM is consequently a gateway and a pre-condition for star formation in galaxies, and determining its distribution is crucial for any theory of galaxy evolution. In this paper, we seek to use numerical models to investigate the broad trends of where the CNM is located in galaxies, what sets this distribution, and how it corresponds to star formation.
The phases of the ISM are usually explained in terms of pressure equilibrium (Field et al. 1969;McKee & Ostriker 1977). As summarised by Ostriker & Kim (2022), the midplane pressure is in vertical dynamical equilibrium with the weight of the ISM. This midplane pressure then determines the balance between hot and twophase (warm+cold) gas such that they have median pressures within 50% of each other. Ultimately this leads to a near linear relationship with star formation (Σ SFR ∝ 1.2 ) due to the importance of feedback in setting this disc pressure. Similarly, the existence of the cold neutral gas phase should also be related to the local pressure environment. Wolfire et al. (2003) have shown that there is a minimum thermal pressure for the CNM phase to coexist with the WNM, which has a value of th / ≈ 3000 K cm −3 at the solar circle. As the pressure falls with galactocentric radius, this leads to a predicted limit of < 18 kpc for the CNM in the Milky Way.
One of the key ways of observing the CNM is through deep surveys of the absorption of the 21 cm line of neutral hydrogen against radio continuum sources since the emission is dominated by the warm phase. Using the Millenium Arecibo 21 cm absorption line survey, Heiles & Troland (2003) found that in the solar neighbourhood the CNM component makes up about 40% of the neutral gas, whereas more recently Murray et al. (2018) estimated the CNM was 28% of the total H but that 20% of the H was in a thermally unstable phase. Dickey et al. (2009) combined multiple 21 cm emission surveys of the galactic plane to deduce that the CNM in the Milky Way was similarly located to the WNM. Moreover, recent work with the Australian Square Kilometer Array Pathfinder (ASKAP) by Dickey et al. (2022) found that in the inner Milky Way, the CNM has a similar scale height to the molecular gas, but in the outer galaxy the CNM and WNM are well mixed and maintain a constant CNM/WNM ratio out to radii of at least 40 kpc. Similarly, Strasser et al. (2007) found CNM at large galactocentric radii in spiral arms in the outer galaxy. Soler et al. (2022) investigated the filamentary structure in the H emission toward the Galactic disk using a Hessian matrix method. The identified filamentary structures correspond to roughly 80% of the H emission and most likely consist of CNM material due to their higher density. The mean scale height of the filamentary H , was lower than that of the total H in the outer galaxy, suggesting that the CNM and WNM have different scale heights in this regime. These results are at odds with the aforementioned pressure equilibrium models, which suggest that in the low-pressure environment of the outer galaxy, the CNM phase should become increasingly less prevalent (Wolfire et al. 2003).
Another way of probing the CNM phase is through the [C ] 158 m emission line, as this is the main coolant of the diffuse ISM. As the emission is sensitive to the density of the gas, that from the WNM is a factor of ∼ 20 less bright than that from the CNM. This means that [C ] 158 m emission from diffuse regions can be assigned to the CNM phase. This exercise was done by the GOTC++ Herschel/HIFI survey published by Pineda et al. (2013). In contrast to H absorption studies these authors found that the CNM column density decreased more rapidly with galactocentric radius than that of the WNM, and consequently, the fraction of the atomic gas in the cold phase was much lower in the outer galaxy (∼ 20%).
The distribution of the neutral atomic medium is something that has not just been a focus of observational studies, but is also a topic of interest for theoretical numerical studies. In cosmological simulations, resolving the scale height of the H gas is a key challenge for reproducing Milky Way type galaxies (e.g Hopkins et al. 2018). Previously, the H scale heights of simulated galaxies were of order of kiloparsecs (Bahé et al. 2016;Marinacci et al. 2017), however recent FIRE (Gensior et al. 2022) simulations have more realistic heights of ∼ 100 pc in the galactic centres, rising to ∼ 800 pc in the outer galaxy. The authors postulate that the solution to faithfully reproducing galaxy scale heights comes from the inclusion of a realistic multiphase medium with cold gas and small scale stellar feedback. Such models typically omit a detailed chemical treatment of the gas, and even FIRE has a mass resolution of over 6 × 10 4 . To obtain a finer description of the CNM, therefore, simulations of gas in an isolated galaxy are needed.
One approach is to simulate gas in stratified boxes (e.g. Hennebelle & Iffrig 2014;Walch et al. 2015;Girichidis et al. 2016;Rathjen et al. 2021). However, these studies have mainly focused on molecular gas and star formation rather than H . One prominent example is the work of Kim et al. (2013) which showed that the star formation rate surface density varied almost linearly with the midplane pressure set by the weight of the ISM. As earlier stated, the local pressure balance of the ISM is also theorised to be extremely important in setting the CNM fraction (Wolfire et al. 2003).
Unfortunately, stratified box simulations typically only cover an area of a few square kiloparsecs of an idealised disc and are therefore unable to investigate the full galactic distribution. One approach is to embed a co-rotating high-resolution box within a galaxy simulation as is done in the Cloud Factory simulations of Smith et al. (2020). These models, and stratified boxes from the FRIGG project (Iffrig & Hennebelle 2017), have been used to interpret the orientations of the H filaments identified in The H /OH/Recombination-line survey of the inner Milky Way (THOR, Beuther et al. 2016), as reported in Soler et al. (2020). However, the simulations of Smith et al. (2020) only reach a high resolution in a 3 kpc box in the star-forming disc. To investigate the true galactic distribution of the CNM we need the entire galaxy to be included. We, therefore, turn to an updated version of the isolated hydrodynamic galaxy models presented by Tress et al. (2020) and Tress et al. (2021), which reach parsec resolution or better in the cold gas across the galaxy while including full hydrogen chemistry.
The paper is structured as follows. We first outline the details of the galaxy models and how they are analysed in Section 2. Then we describe the radial and vertical distribution of the CNM that we find in Section 3. In Section 4, we investigate how this relates to the local pressure conditions in the galactic disc, and the effect of spiral arms. Finally, in Section 5 we give our conclusions.

Isolated Galaxy Simulations
Our simulations are carried out using the A code (Springel 2010) with custom physics modules to describe star formation and cold dense gas. For a full description of our numerical setup, see Tress et al. (2020). However, we briefly summarise the major features and differences from previous work here.
The models of Tress et al. (2020Tress et al. ( , 2021 consist of two simulations of a galaxy disc consisting of dark matter halo (6 × 10 11 M ), bulge (5.3 × 10 9 M ), stellar disc (4.77 × 10 10 M ), and gas disc (5.3 × 10 9 M ). The first model is isolated and develops large scale spiral structure but with breaks and bifurcations, whereas the other is perturbed by a fly-by from a companion galaxy, inducing the formation of strong spiral arms. For our analysis, we focus on an updated version of the Isolated case as this best lends itself to radial averaging. However, in Section 4.4 we will revisit the latter Interacting model.
The original simulation assumed a constant interstellar radiation field consistent with the solar neighbourhood value 0 = 1.7 in Habing units (Draine 1978;Habing 1968). However, this is not a good description of the outer galaxy where the field will be lower due to the low star formation rate. To account for this we computed the time-averaged radial star formation rate surface density profile from 280 to 320 Myr in the original Tress et al. (2020) model during the steady state period of the galaxy. We then fit an exponential function to the star formation density and scale it such that, at the radius where the star formation rate surface density is equal to the solar neighbourhood value, it takes the value of 0 . The final expression takes the form where is the interstellar radiation field, the galactocentric radius in units of kpc, and 0 the radius where the star formation surface density matches the value of the solar neighbourhood. Figure 1 shows the average radial star formation rate density profiles from the original Figure 1. The star formation rate density as a function of radius from the isolated galaxy simulation of Tress et al. (2020) between 280 − 320 Myr. The blue line shows the exponential function used to represent the scaling of the interstellar radiation field. models, and the radial dependence of the new varying interstellar radiation field. (We checked subsequently that the star formation rate was not significantly different in the new runs with varying field.) Shielding from the field is computed assuming a constant shielding length of 30 pc using the T C algorithm . We also assume that the mean free path of the FUV photons in the Milky Way is much less than the scale on which the SFR density is varying, this is a reasonable assumption out to galactocentric radii of 10 kpc in the Milky Way (Wolfire et al. 2003) but may be less valid in the outer galaxy. Nonetheless, while this is an approximation, it is far better than assuming a constant value for the interstellar radiation field. Cosmic ray ionisation is also included at a constant rate of = 3 × 10 −17 s −1 for atomic hydrogen, with the rates for other chemical species (H 2 , He, etc.) scaled appropriately. However, in the atomic gas (i.e. when < 1), photoelectric heating dominates over cosmic ray heating by an order of magnitude or more (see e.g. Figure  10 in Wolfire et al. 2003).
Heating and cooling of the gas is computed simultaneously with the solution of the chemical rate equations. We use the latest version of the cooling functions used by Clark et al. (2019).A temperature floor of 20 K is imposed on the ISM to prevent anomalously low temperatures occurring close to our resolution limit. For our chemistry, we adopt the approach first used in A by Smith et al. (2014), namely the NL97 network of  which utilises the network for hydrogen chemistry first presented in Glover & Mac Low (2007a,b). A simplified description of CO formation and destruction is also included in this network, based on Nelson & Langer (1997). However, we do not analyse CO in this work.
For simplicity, we assume that the gas has solar metallicity throughout our simulated galaxy and that the dust-to-gas ratio is the same as in the solar neighbourhood. In reality, the Milky Way has a metallicity gradient of approximately −0.037 dex kpc −1 (Arellano-Córdova et al. 2020), and similar results are found in other nearby spiral galaxies (see e.g. Kreckel et al. 2019). Therefore, in our simulation, the gas in the outer galaxy ( > 12 kpc) is roughly 2-3 times more metal-rich than we would expect to be the case in a real Milky Way-type spiral galaxy. However, this will have only a minor effect on the temperature of the atomic ISM, since the dominant cooling and heating processes have the same dependence on metallicity, provided that the dust-to-gas ratio scales linearly with the metallicity (Wolfire et al. 1995).
The resolution of the simulations depends on two criteria. Firstly we set a target gas mass of 300 , which means that by default will refine or de-refine the grid so as to keep the masses of all of the grid cells within a factor of two of this value. However, on top of this we require that the Jeans length is resolved by at least four resolution elements to satisfy the Truelove criteria and avoid artificial fragmentation (Truelove et al. 1997;Greif et al. 2011;Federrath et al. 2011). This leads to a mass resolution of ∼ 10 M between densities of 10 −22 < < 10 −21 g cm −3 , which equates to spatial scales of a parsec or smaller. We are therefore confident that the CNM is well resolved.
Star formation is modeled via sink particles (Bate et al. 1995;Federrath et al. 2010), which are non-gaseous particles that represent collapsing regions of gas that will form small (sub)clusters of stars. These are formed by checking if regions of gas with a density greater than 10 −21 g cm −3 are unambiguously bound, collapsing, and accelerating inwards. Only if these criteria are met will the gas be replaced with a sink particle, which can then accrete additional mass that falls within a radius of 2.5 pc of the cell if it is gravitationally bound to it. As star formation is still inefficient at these scales we assume a 5% star formation efficiency (Evans et al. 2009) and associate a stellar and gas fraction to each sink.
Using the model of Sormani et al. (2017), we sample the IMF and associate supernovae with the massive stars as described by Tress et al. (2020). For each supernova, we calculate an injection radius, which is the radius of the smallest sphere centred on the supernova that contains at least 40 grid cells. If the injection radius is smaller than the expected radius of a supernova remnant at the end of its Sedov-Taylor phase, we inject thermal energy from the supernova; otherwise, we inject momentum (e.g. Gatto et al. 2015). Mass is returned with each supernova explosion such that when the last supernova occurs the gaseous component of the sink is exhausted. The sink is then turned into a star particle. To account for type Ia supernovae, we also randomly select a star particle every 250 years and create a supernova event at its position (based on the star formation history of M51, as quantified by Eufrasio et al. 2017, which is similar to our model). Figure 2 shows our updated Isolated galaxy simulation at a time of 300 Myr, where we carry out our analysis. Note, that we performed a smooth averaging procedure on the original simulations to investigate if choosing a single snapshot in time biased the result. We found that the radial profile followed an identical trend radially but with a smoother shape. However, the averaging process masked real fluctuations due to spiral structure and so chose instead the snapshot analysis.

Data Analysis
In order to explore how the CNM varies as a function of galactocentric position and height above the midplane, we bin the simulation output in two different ways. Firstly, we bin the data in cylindrical coordinates in order to study whether the CNM fraction varies with galactocentric radius. We use 60 bins linearly from an inner radius of 0.1 kpc to an outer one of 15 kpc. In the Isolated galaxy simulation, the star formation is mostly contained within a radius ∼ 10 kpc and so both the inner and outer disc is included. To examine the impact of features such as spiral arms we also bin the data in angle using 60 bins from 0 to 360 degrees. For this binning scheme we include all the material within a vertical extent of 1.5 kpc of the disc midplane.
Secondly, we use a vertical binning scheme where we bin the data in radius and vertical extent rather than angle. The same radial interval is used, but the -direction is now binned from −1.0 kpc to 1.0 kpc relative to the galactic midplane using 60 bins.
In each bin we record both the total H mass, and the H mass with a temperature of 200 K or lower, which is chosen to be consistent with the CNM features identified in absorption by Heiles & Troland (2003). (Although it should be noted that the majority of the gas in the CNM phase is much colder than this.) The CNM fraction is taken to be the mass of cold H with < 200 K divided by the total H mass in the bin. The H surface density Σ HI is then simply the total H mass divided by the area of the bin. The exact choice of temperature threshold chosen for the CNM definition makes only a small (10 − 20%) difference to the results. For example, when testing we found the total CNM fraction for the galaxy to be 0.164 with a threshold of 150 K, 0.192 with 200 K, and 0.227 for 300 K. The 200 K threshold is consequently a reasonable compromise, and has the advantage of being consistent with the observational literature.
In order to investigate the origin of the CNM fraction we calculate the mid-plane thermal pressure th in units of the Boltzmann constant, , using the below sum for all cells in the bin with | | < 50 pc in the radially binned data, where is the cell number density, the cell temperature, and the weighting variable, which can be either mass or volume. Note, that we investigated changing the threshold in at which material was included and found it made little difference to the average below a threshold of 100 pc. Following the approach of Kim et al. (2013) we exclude dense regions with number density > 50 cm −3 from the average. Above such densities regions may become self-gravitating, at which point their pressure would no longer be reflective of the global disc conditions. We also calculate the vertical turbulent pressure by taking the sum over all cells in the midplane defined above, where is the cell density, is the vertical velocity, and is the weighting variable of the cell. We divide by for comparison with the thermal pressure.
To investigate the scale height of the H and CNM we use our vertically binned data and fit a Gaussian as was recently done by Bacchini et al. (2019) for THINGS galaxies, using the equation where is the number density in the bin, 0 the number density at the midplane, the distance from the midplane, and the scale height. The fit is performed using the curve_fit routine. As a second check of the vertical extent we also calculate the massweighted position and its dispersion, in each of our radial bins to investigate local variations.
As discussed in Section 2.1 the simulation uses the sink particle method to track the amount of mass going into bound collapsing regions. We use these to find the growth rate of the stellar mass in each sink. To calculate this we compare the analysed time snapshot with the previous snapshot created by the simulation. If a sink has no precursor we assign it as newly created sink mass, but if it does, then we take the difference between the stellar masses to calculate how it has grown. We then divide by the time between snapshots to get the rate. The sinks are then binned spatially in the same manner as the gas mass and added to get the total star formation rate Σ SFR from all the sinks in the bin.  Figure 2 shows a top-down view of the CNM column density alongside similar maps of the full atomic hydrogen distribution and the H 2 distribution. The CNM closely follows the distribution of the molecular gas. For the inner 12 kpc of the galaxy disc we calculate the area covering fraction of material with surface density greater than = 10 19 cm −2 (chosen from a visual inspection of Figure 2) for the three species. H covered the majority of the disc (92%), but the CNM phase covers only a small fraction (5.5%) of the disc surface. The molecular gas covers an even smaller fraction (3.4%) of the disc surface and a visual inspection of Figure 2 shows that the CNM coincides more closely with the H 2 than the total H . This suggests that the CNM is both intermixed with, and extends beyond the molecular clouds. Figure 3 shows the probability density distribution of the H surface density within this region. For comparison we select the pixels where the H 2 and CNM surface density was above our chosen threshold of = 10 19 cm −2 and plotted the probability density of H surface density in these pixels. Both the CNM and molecular hydrogen are found at similar H surface densities. This is in good agreement with recent results from 3D dust mapping that show that nearby clouds are surrounded by extended CNM (Zucker et al. 2021).

Overlap with Molecular Gas
The full H distribution is much more extended than the CNM and reaches out to 15 kpc beyond the star-forming disc. The morphology of the H also changes between the star-forming disc and the outer regions. In the outer disc the H is less filamentary and fills the surface area smoothly without voids, unlike the inner regions where it follows the arms more closely. In Figure 4 we show the CNM fraction throughout the galaxy model at the same time as shown in Figure 2. This shows that the differences between the CNM and H distribution are not simply due to the total mass of gas being larger in the H component and thus making it appear more extensive and space-filling. Instead the relative abundance of CNM to the total H , changes inside and outside of spiral features.   Figure 5 shows the CNM fraction of the gas in the galaxy disc after the radial binning procedure described in Section 2.2. The CNM fraction peaks at small values of and then gradually decreases until values of ∼ 0.8. The distribution is similar to the Arecibo survey of Heiles & Troland (2003) but, as in that survey, a single mean value does not describe the variation in the CNM well. Figure 6 shows how the H surface density and CNM fraction vary with galactocentric radius. As our simulated galaxy is not designed to be an exact analogue of the Milky Way we will adopt the following terminology. The 'inner galaxy' is where < 2 kpc, the 'star-forming disc' is at 2 < < 9 kpc and the 'outer disc' is at > 9 kpc. The total H surface density peaks in the inner galaxy then declines radially until it is below 0.1 M pc −2 in the outer disc. The HI surface density has multiple peaks in the star-forming disc which correspond to spiral arms. The warm and cold neutral mediums contribute to this surface density in different ways. In the inner galaxy the CNM and WNM (where here we include all H with > 200 K in the WNM) are well mixed and follow a similar distribution, albeit at lower surface density Figure 6. The radially averaged gas surface densities (top) and CNM fraction (bottom). The CNM fraction decreases more steeply with galactocentric radius than the total HI surface density.

Dependence on Galactocentric Radius
for the CNM. In the outer disc there is CNM only out to a radius of approximately 12 kpc in contrast to the full H distribution, which continues to beyond 14 kpc. In contrast the molecular hydrogen falls below surface densities of 0.1 M pc −2 at a radius of 8 kpc and so is confined to the star-forming disc. Note that the H 2 surface density is a lower limit as some mass will be inside sink particles. Dickey et al. (2009) report significant CNM in the outer disc of the Milky Way out to very large galactic radii, however we only find CNM in our model at disc radii < 12 kpc, as opposed to < 14 kpc for H . We will discuss this further in Section 4.4.
The CNM fraction is approximately constant in the disc at a value of roughly = 0.2, apart from a peak at the galactic centre. This trend of constant with galactocentric radius is in agreement with the absorption studies of Dickey et al. (2009Dickey et al. ( , 2022, but is in tension with the GOTC+ survey of Pineda et al. (2013) who see a clear radial decreasing trend in the CNM fraction derived using the C line. Figure 7 shows the angular dependence of the CNM fraction at increasing radii. The contours denote the H surface density and take the form of diagonal stripes in this phase space due to the spiral arms. Unsurprisingly, there is a low CNM fraction outside the spiral arms where the H column density is low and the gas temperature is hotter. However, the CNM fraction remains non-negligible between 10-12 kpc where the spiral arms are less defined.
When comparing to observations it is useful to understand the temperature distribution of the gas. Most obviously because the CNM, WNM split is a division in temperature, but also because the CNM temperature, cool , will determine the observed spin temperature that is directly measured from absorption features. For example, Strasser et al. (2007) found that the values of cool derived from absorption features in the outer galaxy showed no strong dependence on galactic Figure 7. The CNM fraction calculated in radial and angular bins. The contours show the HI column density with contour levels of 5, and 10 M pc −2 respectively. The CNM fraction is higher in the spiral arms but is still substantial at 10 kpc where the arms are less defined. radius. In Figure 8 we show the mass-weighted mean CNM temperatures in our model. In agreement with observations we see a mostly flat trend. There is some hint that in the star forming disc the temperature is ∼ 10% lower than the outer disc but this is only a small variation. In the star-forming disc, the right hand panel of Figure 8 shows that the temperatures are cooler in the spiral arms.

Dependence on Vertical Extent
In the top panel Figure 9, instead of binning with angle we bin in height above the disc. Material with a significant CNM fraction pervades all the H contours, even at the lowest surface density of 0.1 M pc −2 . To investigate this further, in the bottom panel of Figure  9 we plot how the scale-height varies with radius in the total H compared to the H in the CNM phase. The scale height is determined by fitting to Equation 4 as described in Section 2.2. We find typical H scale heights of 100 pc in the star forming disc, which then flares in the outer galaxy as expected from observations (e.g. Yim et al. 2014;Bacchini et al. 2019;Randriamampandry et al. 2021;Soler et al. 2022) The distribution is noisy due to the spiral structure but the two populations follow each other closely. However, at 12 kpc they diverge due to the CNM disappearing from the gas phase. Unlike the H , the CNM scale height does not flare.
To investigate whether this behaviour is due to averaging we investigate another metric which does not depend on a fit. Figure 10 shows radial profiles and angular maps of the mass-weighted mean position and its dispersion. While the mean position is always close to the 0 position in the simulation there are variations of up to 50% of the scale height in Figure 9. The CNM's mean vertical position closely follows that of the full H distribution, but the dispersion shows very different behaviour. In the galaxy centre both the total H and the CNM have very low scale heights (∼ 50 pc). However, while the dispersion of the total H gas steadily rises, the CNM remains tightly confined below values of 100 pc. Note that this is in good agreement with the findings of Dickey et al. (2022) for the inner galaxy and star-forming disc. The CNM distribution, therefore, seems to consist of discrete 'clumps' of CNM which are at a variety of vertical heights.

Local Pressure Environment and Star Formation
To investigate how the local pressure environment corresponds to star formation and how this connects to the CNM fraction, we plot in Figure 11 the mass-averaged total (thermal and turbulent) pressure, , vs the local star formation rate surface density Σ SFR . The grey dashed line shows the empirically determined scaling found from the simulations of Kim et al. (2013) Σ SFR = 2.1 × 10 −9 / 10 4 1.18 (5) where Σ SFR is in units of M pc −2 yr −1 and / in units of cm −3 K. This follows the same trend as the scaling of our models confirming that the total midplane pressure plays a major role in setting the global star formation rate. The normalisation of our model is different, however it should be noticed that our treatment of star formation, using sink particles, is different from the star particle approach of Kim et al. (2013) who assume a constant efficiency per free fall time and this may account for the discrepancy. However, the observational results of Sun et al. (2020) are in good agreement with Kim et al. (2013) so further investigation may be needed. We also considered the volume-weighted pressure relation and found it had similar scaling behaviour but there was a discontinuous jump in Σ SFR at high stellar densities. This is most likely due to the small volume fraction of dense star-forming regions. For this reason we consider only mass-weighted pressures going forward. Regardless of how calculated, the dispersion in the relationship decreases at large star formation rate surface densities as high pressures are needed to create sufficient dense gas.
In Figure 12 we investigate how the total mass-weighted pressure corresponds to the CNM and star formation. The left panel shows the mass-weighted total pressure vs the CNM surface density. The dashed line shows a line with gradient of 0.6, which matches the upper envelope of the data. There is a large amount of scatter, but the CNM column density increases with pressure as would be expected.
The middle panel of Figure 12 shows the star formation surface density plotted against the CNM surface density in a manner analogous to a Kennicutt-Schmidt relation (Kennicutt 1998). Observationally, the star formation density scales with a power of ∼ 1.4 against the total gas surface density (H + H 2 ), as shown, for example, in de los Reyes & Kennicutt (2019). However, when compared to only molecular gas, the scaling is linear down to at least 0.1 solar metallicity (Bigiel et al. 2008;Whitworth et al. 2022). Using the S P linear regression feature we fit a power law to the data in logspace and found a slope of 1.56 with a standard error of 0.06 best described the data. That this value is closer to the Kennicutt total gas exponent scaling with star formation than the molecular value, suggests that the CNM is not gravitationally bound and actively forming stars.
The right-hand panel of Figure 12 shows how the CNM fraction corresponds to the local star formation rate surface density. There is a large scatter, particularly at low CNM fractions, however, the trend is clear that higher values of correspond to higher star formation rate surface densities. In the regions of our galaxy model with the most star formation, the CNM fraction is as high as 80%.
We can now investigate how the pressure contributes to the largescale trends discussed in Section 3.2. Figure 13 shows the massweighted pressure as a function of the galactic radius where the  CNM fraction > 0. As expected the pressure steadily falls from the inner galaxy to the outer, albeit with local variations. However, strikingly the CNM vanishes from the gas phase at radii of about 12 kpc, which corresponds to where the local pressure falls below 1000 K cm −3 .

Dependence on Radiation Field
The original Isolated galaxy modelled in Tress et al. (2020) had a constant UV field and so it is possible to use this as a comparison to investigate the importance of the interstellar radiation field in setting the CNM distribution of the galaxy. Figure 14 shows the radial dependence of the CNM surface density, fraction, and vertical dispersion for the constant field of magnitude 0 . Most strikingly the CNM only extends to a radius of 10 kpc due to the higher field at large radii. This is the same extent as the star-forming disc in these models, and so the CNM no longer extends into the outer galaxy. The surface density of the HI and CNM are flatter as a function of radius than the surface distribution with the varying field shown in Figure  6. Intriguingly, while the radial extent of the CNM changes (from 12 to 10 kpc) the H and H 2 extent remains largely unchanged (at 14, and 8 kpc respectively). Clearly, the CNM is far more sensitive to the local radiation field in this regard compared to the other phases.
With a constant interstellar radiation field the CNM fraction , shown in the middle panel of Figure 14, decreases radially rather than remaining flat in the disc as was the case in Figure 6. The vertical dispersion of the CNM, however, remains similar with and without the varying field with the CNM having a dispersion of around 100 pc throughout the disc in both cases with no flaring.
To explore the origin of this behaviour we plot the mean total midplane pressure as a function of radius for the model with a constant field in Figure 15. The pressure now falls below a value of 10 3 K cm −3 at a galactocentric radius of 10 kpc as opposed to 12 kpc in the varying case. Both of these radii exactly correspond to where the CNM disappears from the gas phase, once again confirming the important role of midplane pressure in setting the CNM distribution.

Turbulence in the Outer Galaxy
However, pressure alone is not a full explanation for the distribution of the CNM in the gas phase. Turbulent velocity fluctuations are important for pushing gas out of equilibrium and creating local over-densities where the CNM can form. The vertical velocity dispersion observed in H in spiral galaxies is observed to be both The correspondence between the local mass-weighted pressure and the CNM surface density Σ . Middle: The surface density of the CNM vs the local star formation rate surface density, analogous to a Kennicutt-Schmidt relation for the CNM. The grey line shows a linear fit to the data, which has a similar exponent to that derived observationally for total column density. Right: The CNM fraction vs star formation rate surface density, which shows that higher CNM fractions correspond to greater star formation rate surface density. Figure 13. The galactocentric radius and mass weighted pressure of the points shown in Figure 12 for all points where > 0. The CNM distribution ends at ∼ 12 kpc, which is where the pressure falls below values of 10 3 K cm −3 . highly turbulent, and varying radially from values of roughly 12 to 15 kms −1 in the central parts, to ∼ 4 − 6 kms −1 in the outer parts e.g. (Dickey & Lockman 1990;Kamphuis & Sancisi 1993;Rownd et al. 1994;Meurer et al. 1996;de Blok & Walter 2006). Such velocity dispersions can be theoretically explained via a number of mechanisms that include supernova feedback (Dib et al. 2006), MRI instability (Piontek & Ostriker 2005, 2007, or accretion from the galactic halo (Klessen & Hennebelle 2010). However, none of these effects are operating in our model. In the outer galaxy there is no star formation and hence no supernova feedback. Our model is purely hydrodynamic, and it is isolated with no accretion of material from the intergalactic medium.
In Figure 16 we plot the mass averaged vertical velocity dispersion in radial bins for H and the CNM. For both phases the velocity dispersion is higher in the star-forming disc where we have supernova feedback, however, it remains substantial and of similar magnitude to observations in the outer galaxy. Therefore we must conclude that galactic dynamics alone are enough to match the observed velocity dispersion. Consequently, we are confident that we are not underestimating the CNM fraction in our models due to a lack of turbulence.

CNM and Spiral Arms
In our model, the CNM extends to midway in the outer disc but does not cover its full extent. However, the absorption studies of Dickey et al. (2022) find CNM in the entirety of the outer disc of the Milky Way, out to radii of 40 kpc. Kinematic distances are hard to estimate and absorption measurements give information along single sight lines. Still, it seems clear that in the Milky Way the CNM can exist beyond the radii predicted in our models.
One possibility is that spiral arm structure is able to concentrate gas in isolated pockets in the outer disc. For example, Strasser et al. (2007) observe CNM in the outer Milky Way, but specifically target spiral arms. Our fiducial model is of an isolated galaxy without strong spiral structure in the outer regions. Therefore to test the importance of spiral arms in forming the CNM we turn to the other model described by Tress et al. (2020), where a close encounter has triggered the formation of well-defined spiral arms to create an 'M51 analogue'. Figure 17 shows the gas distribution of the Interacting model at the same time as the isolated model analysis. Note that the interstellar radiation field is constant in this model so we also compare it to the model discussed in Section 4.2 rather than just our fiducial case. Figure A1 in Appendix A shows the column density maps for the Isolated model with a constant interstellar radiation field.
The atomic hydrogen distribution in the outer galaxy is now very different from the Isolated constant field case. Where before the distribution smoothly filled the surface of the outer disc, now the atomic hydrogen is concentrated in the extended spiral arms. The lower panels of Figure 17 show that both H 2 and the CNM now can be found in the outer galaxy. For instance, a particularly large clump of CNM can be seen in the lower left corner. Once again the CNM closely follows the molecular hydrogen, but with a higher column density over a slightly larger area than the Isolated case. Figure 18 shows the binned CNM fraction of the interacting galaxy. In the spiral arms, the CNM now persists to large radii albeit at a low level. Figure 19 shows the radially averaged profiles for these maps. When averaged over an annulus of the galaxy disc the surface density of the CNM is negligible in the outer galaxy and follows the same radially decreasing trend as Figure 6.
Finally, the lower panel of Figure 19 shows the mean dispersion in the z position, ( ), of both H and the CNM. In the inner galaxy, the CNM and H have increasing dispersion, suggesting that they are both well mixed vertically in the disc. However, beyond 6 kpc, where the dense gas is mainly in spiral arms, ( ) drops for the CNM but continues to rise for H . This further suggests that the CNM in the extreme outer galaxy is confined to distinct clumps of gas where the local pressure has been enhanced by spiral arms.
Certainly, the scenario presented here is only one possibility for how CNM might exist in the extreme reaches of the outer galaxy. This paper focuses on a pair of models, and is not a parameter study. It explores two extremes, one model with weak spiral structure, and another with strong, and consequently brackets a range of possible CNM distributions. Magnetic fields could add additional support to molecular cloud envelopes which could then contain a substantial amount of CNM as well as CO-dark molecular gas (e.g Smith et al. 2014). Another possibility is that cold gas could be accreted onto the outer galaxy from the circumgalactic medium such as in the Magellanic Stream of our own Galaxy (Fox et al. 2014). Alternatively, the discrepancies that we see between Milky Way observations and our model may be due to projection effects. Our analysis has been from an external perspective, in order to better quantify the variation with galactocentric radius and vertical height above the disc midplane within the Milky Way the CNM distribution has to be recovered from sight lines containing emission and absorption from gas at multiple different radial positions within the disc.

CONCLUSIONS
We have investigated the distribution of the CNM in new highresolution simulations of an isolated disc galaxy based on those presented by Tress et al. (2020) but with a radially varying interstellar radiation field. Our A simulations use a variety of custom modules allowing us to follow the chemical and star-forming evolution of the dense gas. These included time-dependent gas chemistry, gas self-shielding from the ambient UV field, sink particles to represent star formation, and supernova feedback. The resulting models allow the distribution of the CNM to be followed down to sub-pc scales in the dense gas. Most of our analysis focuses on an isolated galaxy model, however, we later compare it to a model with a constant interstellar radiation field, and another with identical mass which has been perturbed by a close encounter to generate strong spiral arms that extend to the outer galaxy.
Our conclusions are as follows:   (i) No single value describes the CNM fraction (the fraction of the atomic gas with T<200 K relative to the total H ) everywhere in the galaxy. Values range from = 0 to 0.8 with high values being less likely.
(ii) The CNM is not uniformly distributed in the star-forming disc, but follows a clumpy distribution. A comparison of the column density maps shows that it overlaps more closely with the H 2 distribution than the total H . This is particularly true in spiral arms where the CNM fraction is clearly enhanced.
(iii) The CNM extends into the outer galaxy, beyond the starforming disc but does not extend as far out as the full H distribution. In our model, the radial surface density profiles of H , the CNM, and H 2 remain above a value of 0.1 M pc −2 out to radii of 8, 12, and 14 kpc respectively.
(iv) The CNM fraction remains approximately constant with galactocentric radii in agreement with measurements with H absorption. This effect is a consequence of the falling interstellar radiation field compensating for the decreasing column density at large radii. In our comparison model with a constant interstellar field the CNM fraction decreased with galactocentric radius.
(v) The vertical distribution of the CNM is clumpy in the starforming disc with a scale height of around ∼ 100 pc and no flaring. An analysis of the dispersion, ( ), suggests that the CNM is more localised in the -direction than the overall H distribution.
(vi) The star formation rate in the galaxy is well correlated with the total midplane pressure in logspace as suggested by Kim et al. (2013). We find that the CNM column density also correlates well with pressure as expected from models of the ISM such as Wolfire et al. (2003). We find no CNM in our isolated galaxy discs beyond the point where the total pressure / drops below 1000 K cm −3 . (vii) The 'CNM Kennicutt-Schmidt' relation (i.e. the star formation rate surface density plotted against the CNM column density) has a scaling of ∼ 1.5, which is more similar to the relation seen for total column density than the linear relation seen for molecular gas. This suggests that while the formation of the CNM is a precursor to star formation, it is not the predictor of gravitationally collapsing star-forming regions that the H 2 is. However, the CNM fraction does increase with star formation density, with the highest CNM fractions always associated with active star formation.
(viii) The vertical velocity dispersion of the H in our models is in good agreement with observations of H in nearby galaxies. In the star-forming disc a major contribution to this is supernova feedback, but even in the outer galaxy where there is no star formation galactic dynamics is sufficient to drive velocity dispersions between 5 − 10 km s −1 .
(ix) Spiral arm features in the outer galaxy may give rise to isolated clumps of CNM at extremely large radii, beyond where it is expected in our more symmetric isolated models. This may be an explanation for the CNM seen at extremely large galactic radii in the Milky Way by H absorption studies (Dickey et al. 2022).