Introduction

Studying the properties of molecular clouds is crucial to understanding star formation (Kennicutt & Evans 2012). The dominant constituent of molecular clouds is molecular hydrogen (H 2 ), which is a perfectly symmetric molecule, rendering it largely invisible at the typical temperatures of molecular clouds. While observable in ultraviolet absorption against background sources, it can only be detected via emission in environments where the gas is excited to temperatures above a few hundred Kelvin. The second dominant species is neutral helium, which remains inert in molecular clouds. Therefore, observational studies of molecular clouds largely rely on tracer species, namely, emission from dust and molecules. The most important of these tracers is carbon monoxide (CO; Bolatto et al. 2013). CO has a relatively high abundance, canonically [CO/H 2 ]≈10 -4 (Hollenbach & Tielens 1999), making it the most abundant molecule after H 2 . The small dipole moment allows its rotational transitions to be easily excited at the cold temperatures of molecular clouds. A crucial CO observable is the J=(1-0) rotational transition at a rest frequency of 115.27 GHz.

It is common for the emission of the lowest rotational transition of CO to be used to measure the total molecular gas (Bolatto et al. 2013). This is encoded in the CO-to-H 2 conversion factor, X CO , and the related quantity α CO . X CO is defined as

where N(H 2 ) is the H 2 column density in cm -2 and W CO (J = 1 -0) is the CO flux in K km s -1 . The fiducial Milky Way (MW) value is X CO,MW =2×10 20 cm -2 (K km s -1 ) -1 (Bolatto et al. 2013). This conversion factor has been used to estimate gas mass in local, resolved studies of MCs and the molecular gas mass in high-redshift galaxies (e.g., the COLDz survey; Riechers et al. 2019). A significant number of studies, both observational and theoretical, have been devoted to measuring, modeling, or applying X CO . Prior work shows it varies with density, metallicity (Bell et al. 2006;Shetty et al. 2011;Lagos et al. 2012;Narayanan & Hopkins 2013;Glover & Clark 2016), cosmic-ray (CR) ionization rate (CRIR; Bell et al. 2006;Wolfire et al. 2010;Bisbas et al. 2015;Clark & Glover 2015;Glover & Clark 2016;Remy et al. 2017;Papadopoulos et al. 2018), and the radiation field (Bell et al. 2006;Wolfire et al. 2010;Shetty et al. 2011;Lagos et al. 2012;Narayanan & Hopkins 2013;Clark & Glover 2015;Glover & Clark 2016;Gaches & Offner 2018a;Gong et al. 2018). Previously, Gaches & Offner (2018a) found that far-ultraviolet radiation feedback from forming stars can reproduce the higher X CO values measured toward diffuse starforming clouds in the outer galaxy.

Traditional one-dimensional photodissociation region (PDR) models have long predicted that neutral carbon will exist only in a thin transitional layer between ionized carbon and CO (Hollenbach & Tielens 1999). However, observations show that forbidden line emission from neutral carbon covers similar spatial extents as CO (e.g., Ikeda et al. 1999;Kulesa et al. 2005;Lo et al. 2014). It is posited that forbidden line emission from neutral carbon is a good tracer of the gas mass (Papadopoulos et al. 2004;Offner et al. 2014;Glover et al. 2015;Glover & Clark 2016). Synthetic observations of hydrodynamic simulations show that X CI has a smaller dispersion than X CO within a molecular cloud and is a better tracer in low-metallicity gas that tends to become CO-dark (Offner et al. 2014;Glover et al. 2015;Glover & Clark 2016). Observational studies using X CI as a tracer of gas mass performs as well as X CO (Lo et al. 2014). X CI is defined analogously to X CO :

where W(CI) 609 μm is the integrated flux of the  P P 3 1 3 0 transition at 609μm. Gaches et al. (2019, hereafter Paper I) presented a modified astrochemical code that includes CR attenuation in situ. Paper I included CRs accelerated by accreting, embedded protostars and CR attenuation in one-dimensional astrochemical models of molecular clouds. We used the code to study the impact of changing the CR spectrum due to differing galactic environments and the effects of embedded CR sources for a subset of species including CO, HCO + , and N 2 H + and tested various prescriptions for constraining the CRIR. We found that ions are enhanced and neutrals are depleted in dense gas due to embedded CRs. Carbon chemistry is substantially altered depending on the assumed CR model: CRs produced by embedded sources create a significant reservoir of atomic carbon, mostly neutral, in dense gas. Embedded CRs reduce the amount of CO in clouds and warm the gas to over 30 K. In this paper we investigate the impact of the above effects on X CI and X CO . In Section 2 we describe the methods used in this paper. In Section 3 we present the results and discuss the implications for observations.

Methods

We use the same astrochemical models from Paper I and summarize the methodology here. See Paper I for further details.

We generate synthetic protoclusters assuming the tapered turbulent core (Offner & McKee 2011) accretion model following the method described in Gaches & Offner (2018a). We directly sample from the bivariate protostellar mass distribution using the method of conditional probabilities. Each molecular cloud is described by a gas surface density and number of constituent protostars, Σ cl and N * , respectively. We only consider models where the star formation efficiency ε g ≡M * /M gas 50%.

We calculate the accelerated proton spectrum due to accretion shocks for each star in the protocluster. CR protons are assumed to be accelerated via diffuse shock acceleration (DSA; reviewed by Drury 1983;Melrose 2009) near the surface of the protostar. DSA predicts a power-law spectrum in momentum space, j(p), with an injection momentum, p inj , set by the shock gas temperature and a maximum energy constrained by collisional energy losses and upstream diffusion (Gaches & Offner 2018b). The CR flux spectrum is

where j 0 is the normalization constant calculated from the total shock energy and efficiency, and a is set by the shock compression factor. We find that the maximum proton energy is typically between 1 and 10 GeV (Gaches & Offner 2018b).

We attenuate the CRs by the gas surface density out of each protostellar core, Σ core =1.22Σ cl , following Padovani et al. (2009). We assume the CRs within the core free-stream outward since shallower attenuation produces too much CR heating in the core (see Gaches & Offner 2018b). CRs may also be attenuated by the accretion flow (Offner et al. 2019), although we do not include this in the model. The total number of CRs escaping into a natal molecular cloud embedding a protocluster is the sum of the CRs accelerated by the individual protostars and then attenuated into the surrounding gas:

. We embed the protoclusters in the center of one-dimensional molecular clouds with a density profile, n(r)=n 0 (R/r) 2 . We set the outer density to n 0 =100 cm -3 , and the radius, R, is determined by the total column density, μm H N(R)=Σ cl . We utilize a modified version of the PDR astrochemistry code 3D- PDR (Bisbas et al. 2012) described in Paper I, 6 which includes CR attenuation in situ. The astrochemistry code uses CR spectra at the surfaces of the gas model as inputs, rather than a global CRIR. It is not known exactly how CRs transport through molecular clouds. Therefore, we consider two different transport regimes: diffusive (1/r) and rectilinear (1/r 2 ). We use the two external CR spectra from Ivlev et al. (2015): a model that extrapolates the Voyager 1 data, , and one that attempts to account for modulation from interstellar gas, .

We also consider the impact of FUV radiation and we irradiate the external surface of the molecular cloud with the normalized interstellar radiation field described in Draine (1978). We model the chemistry with the gas-phase UMIST12 network (McElroy et al. 2013), which includes 215 species and ∼3000 reactions. The network does not include gas-grain reactions, freezeout, or any desorption processes. We do not include the grain-assisted recombination proposed in the reduced network presented by Gong et al. (2017). We explored the impact of grain-assisted recombination for C + and He + on our results and found no significant changes in the CO-to CI-to-H 2 conversion factors. The inclusion of grain chemistry will be investigated in future studies. We include a model following the canonical setup: a low ionization rate with no attenuation, denoted as LNA. Models denoted with H or L utilize the  and  external spectra described above. Models without embedded sources are denoted with NI, while those including sources in the diffusive or rectilinear regimes are denoted DI or RI, respectively. We consider the six different

Note. L/H denotes a low/high external CR spectrum, NI denotes no internal sources of CRs, DI denotes internal sources with a=1 (diffusive transport), RI denotes internal sources with a=2 (rectilinear transport), and NA denotes no internal sources or CR attenuation.

CR models listed in Table 1: LNA, LNI, LRI, LDI, HNI, and HDI. 3D-PDR calculates the CO line-integrated emissivity, ò, for the J-ladder from J=0 to J=41 and the CI 307 and 609 μm emissivities assuming non-local thermodynamic equilibrium and using an escape probability method to account for the line opacity. We calculate the line flux from the emissivity:

where c is the speed of light, k b is Boltzmann's constant and ν is the line frequency. This definition of integrated flux assumes that the interstellar medium is entirely optically thin. We calculate the H 2 column density from the astrochemical models

where x(H 2 ) is the abundance of H 2 and n H is the gas density. Finally, we compute X CO using Equation (1).

Results and Discussion

We present the results from the astrochemical models on the CO-to-H 2 and CI-to-H 2 conversion factors here. A more general discussion on the astrochemical impact of CRs accelerated within protoclusters is presented in Paper I.

Effect of CRs on X CO

Figure 1 shows the CO-to-H 2 conversion factor as a function of cloud surface density, Σ cl , and star formation efficiency, ε g , for four of the CR models in Table 1. We plot X CO normalized to the fiducial MW value X MW =2×10 20 cm -2 (K km s -1 ) -1 (Bolatto et al. 2013).

The behavior of X CO changes significantly with the assumed CR model. X CO varies only as function of surface density for the models without internal sources, LNI and HNI. There is a 0.2 dex offset in X CO between models using the high and low external CR spectrum for Σ cl <3 g cm -2 owing to increased temperatures at low extinction in the model HDI. The decline in X CO at higher surface densities is the result of a larger turbulent line width because of two cooperating effects. First, there is a higher temperature due to the increasing importance of turbulent heating. Second, the turbulent line width produces brighter, but still optically thick, CO emission.

In the models with CRs that attenuate diffusively, LDI and HDI, X CO becomes a sensitive function of the star formation efficiency, losing much of the dependence on surface density. X CO is reduced by up to 0.5 dex due to embedded sources with the lowest values occurring for the highest star formation efficiencies. It is important to emphasize that CRs from embedded sources do little to reduce the overall amount of H 2 (Paper I). However, they cause two effects that act to decrease X CO . First, while they reduce the amount of CO in deeply embedded regions of the cloud, they cause an enhancement of CO in low extinction gas due to an increase of HCO + and, following the formation of OH through + H 3 , the OH formation pathway becoming important (Bisbas et al. 2017). Second, the increased CRIR leads to higher kinetic temperatures making the CO emission brighter overall.

Some prior work has investigated the effect of star formation on X CO . CR and FUV feedback from star formation external to the molecular cloud can be modeled by scaling their intensity linearly with the star formation rate (SFR; Papadopoulos 2010). This is motivated by the relationship between the supernova rate and the SFR and implicitly assumes that CRs are mainly accelerated in supernova shocks. Clark & Glover (2015) used these relations to model how the SFR affects X CO in simulated molecular clouds. They found that X CO increases with the SFR if the cloud properties remain fixed. The increase of X CO with SFR is very weak if the density of the cloud scales with the SFR. Bisbas et al. (2015) modeled the effect of enhanced CRs on the [CO/H 2 ] ratio, comparing different environments. They show that [CO/H 2 ] decreases substantially with an increase in the CRIR. By construction, these models only account for variations in the external CR flux and neglect CRs accelerated within protoclusters due to accretion, jets, or stellar winds.

Effect of CRs on X CI

Forbidden line emission from neutral carbon is a possible tracer for molecular gas, as discussed above. Figure 2 shows X CI as a function of surface density, Σ cl , and star formation efficiency, ε g . X CI shows the same qualitative trends as X CO , although it is more sensitive to the CRIR: a spread of 1.2 dex in X CI and 0.5 dex in X CO for the LDI model. The canonical model, LNA, which has no attenuation, exhibits a maximum value of X CI 4×10 21 cm -2 (K km s -1 ) -1 . Models using the high, external CR spectrum, HNI and HDI, exhibit a 0.2-0.8 dex reduction in X CI compared to the low spectrum models, LNI and LDI, respectively.

The increased CRIR throughout the cloud in the high models and those with internal sources causes atomic carbon to exist outside a thin transition layer. Atomic carbon is formed in the dense gas through the destruction of CO by He + :

He CR He e ,

with neutral carbon forming from recombination of C + . Neutral carbon is also the result of direct dissociation of neutral molecules, such as CO, by CR protons and CR-generated photons. This enhancement leads to a reduced X CI . Embedded sources cause X CI to decrease by over an order of magnitude across two orders of magnitude increase in the star formation efficiency.

Neutral carbon emission is easily observable at high redshifts due to the line shifting to millimeter wavelengths. Starburst galaxies have higher SFRs producing extreme environments and more CO-dark gas (Wolfire et al. 2010;Glover & Clark 2016). Thus, at high redshifts and in galaxies undergoing starbursts, CI may become an optimal tracer of molecular gas.

Statistical Trends

Figure 3 statistically summarizes the impact of the various CR models on X CI and X CO . The violin plots show the distribution of the logarithmic difference between X i as calculated with the canonical model, LNA, and each of the CR models in Table 1 using the clouds across the Σ cl -ε g space as samples. These distributions represent the impact on X CO when CR attenuation or embedded sources are neglected.

We find very little deviation in X CO when attenuation is included in quiescent models without internal sources. Comparison to the star-forming and extreme CR model without internal sources, HNI, shows that X CO will be overestimated by 0.15 dex in calculations using the often-assumed CRIR of ζ≈10 -17 s -1 . CRs from embedded sources, which propagate via diffusion, decrease X CO for all clouds. Furthermore, there is a substantial spread due to variation with the number of protostars, N * . The high model with internal sources, HDI, logarithmic difference with the canonical model exhibits a dispersion of 0.3 dex, similar to the spread derived from MW observations (Bolatto et al. 2013). If CRs from embedded sources transport as r -2 there is no impact on X CO because the CRIR is lower and dominated by the CRs originating from external sources rather than internal ones.

The X CI distributions in the right panel of Figure 3 show much greater sensitivity to the CR model assumptions. All models differ significantly from the often-assumed canonical model in X CI . X CI decreases by 0.5 dex for the high model with no internal sources, HNI, and massive and inefficient starforming regions. In the case of a "Quiescent" CR environment, CRs from embedded sources have a larger impact on X CI . The inclusion of CRs from embedded sources in star-forming and extreme environments, represented by HDI, reduces X CI by nearly a dex compared to the canonical model.

Comparisons to Galactic-scale Observations

The hatching in Figure 1 denotes different CR environments: "Quiescent" regions with z á ñ < -- 10 s

x 16

1 , "Star-forming" regions with z < á ñ < ---10 10 s

x 16 15 1 , and "Extreme" regions Figure 1. Color shows logX CO /X MW where X MW =2×10 20 cm -2 (K km s -1 ) -1 as a function of gas surface density, Σ cl , and star formation efficiency, ε g . White shaded cells show regions where X CO is consistent with Milky Way observations, -0.3logX CO /X MW 0.3. The hatched regions indicate different cosmic-ray environments, where we define z á ñ x , the spatially averaged CRIR, z á ñ < - 10

x 16 as "quiescent," z < á ñ < --10 1 0

x 16 15 as "star-forming," and z á ñ > - 10

x with z á ñ > -- 10 s

x 15

1 , where z á ñ x is the spatially averaged CRIR. These labels are motivated by observational surveys that show the majority of pointings through diffuse gas have 10 -16 <ζ<10 -15 s -1 . Low A V observations where ζ>10 -15 s -1 are primarily sightlines toward the galactic center (Indriolo & McCall 2012;Indriolo et al. 2015).

There have been numerous observational studies measuring X CO in different environments within the MW and other galaxies (see Bolatto et al. 2013 and references within). Remarkably, in the MW and many of the Local Group galaxies, X CO is relatively constant on kiloparsec scales. The consistency of X CO in the MW and Local Group can be explained by similar molecular cloud properties due to star formation feedback (Narayanan & Hopkins 2013). There is a general trend in starforming galaxies of low values of X CO toward the center and larger values in the outer disk (Sandstrom et al. 2013).

The white shading in Figure 1 shows where X CO is consistent with the MW average value and spread. Models without embedded sources, LNI and HNI, are only consistent with the MW value for Σ cl <0.2 g cm -2 and Σ cl <0.6 g cm -2 , respectively. Models with high surface density and low star formation efficiency, similar to clouds in the Galactic Center, exhibit a decreased X CO compared to clouds with Σ≈1 g cm -2 . The introduction of embedded sources increases the agreement with the MW X CO . Clouds with star formation efficiencies greater than 2% in the low model with internal sources, LDI, have X CO =(2.5 ± 1)×10 20 cm -2 (K km s -1 ) -1 , consistent with the MW value. High models with the internal sources show X CO =3±1.5×10 20 cm -2 (K km s -1 ) -1 . Thus, CRs accelerated during the star formation process act to regulate X CO and reduce variation.

Starburst galaxies tend to have lower values of X CO (Downes & Solomon 1998;Papadopoulos & Seaquist 1999;Papadopoulos et al. 2012;Salak et al. 2014). Our models show that X CO always decreases toward regions with more extreme CR environments. Environmental changes, which occur in higherredshift galaxies due to enhanced supernova rates, will also decrease X CO and X CI . In starburst galaxies, which have high star formation rates, this decrease could be compounded by CRs produced during the star formation process.

Summary

We found in Paper I that the inclusion of CR sources, specifically accreting protostars, embedded within molecular clouds and CR attenuation make the CRIR vary throughout the cloud. In this paper, we investigate the impact of different external CR fluxes and the inclusion of embedded CR sources on the CO-to-H 2 and CI-to-H 2 conversion factors. We find that differences in the CR flux caused by changes in the external environment and embedded star formation alter X CI significantly and X CO by factors of a few. However, external environmental changes alone reduce X CO only by 0.2 dex, within the measured spread of X CO in the MW (Bolatto et al. 2013). The difference in X CI is more pronounced: it declines by an order of magnitude for the lowest surface density environments. The inclusion of embedded CR sources removes the strong dependence of X CO and X CI on surface density and reduces the conversion factors by 0.6 and 1.2 dex, respectively. Embedded sources act to regulate X CO and reduce variation as a function of gas surface density and star formation efficiency. Clouds in the low model including embedded sources, LDI, with star formation efficiencies greater than 2% are consistent with the observed MW value and spread of X CO,MW =2×10 20 ±0.3 dex cm -2 (K km s -1 ) -1 . Observations of the CRIR in diffuse gas in the MW show that the average CRIR z á ñ » - 10 16 , which is represented by our models with a high surface CR spectrum. Models with this CRIR and ongoing star formation are consistent with the observed MW value for all regions with star formation efficiencies greater than 1%. Our models reproduce the trends of a decreasing X CO toward more extreme CR environments, such as those observed in the Galactic center, the high-redshift universe, and starburst galaxies. Our results motivate the inclusion of CR physics and the possibility of CR feedback from internal sources when modeling X CO and X CI .