The X_CO conversion factor from galactic multiphase ISM simulations
Munan Gong Eve C. Ostriker Chang-Goo Kim



1Department of Astrophysical Sciences, Princeton University,
Princeton, New Jersey 08544, USA
2Max-Planck Institute for Extraterrestrial Physics,
Garching by Munich, 85748, Germany; 
munan@mpe.mpg.de
2Center for Computational Astrophysics, Flatiron Institute, New York, NY, 10010


    (J=1-0) line emission is a widely used observational tracer of
    molecular gas, rendering essential the X_ factor, which is applied to 
    convert  luminosity to  mass. We use numerical simulations to 
    study how X_ depends on numerical resolution,
    
    
    
    
    
    non-steady-state 
    chemistry, physical environment, and 
    observational beam size. Our study employs 3D magnetohydrodynamics (MHD)
    simulations of galactic disks with solar neighborhood conditions,
    where star formation
    and the three-phase interstellar medium (ISM) are self-consistently
    regulated by gravity and stellar feedback. Synthetic  maps are
    obtained by post-processing the MHD simulations with
    chemistry and radiation transfer. We find that 
     is only an approximate tracer of . On parsec scales,
    W_ is more fundamentally a measure of mass-weighted volume density, 
    rather than  column density.
    Nevertheless, 
    ⟨X_⟩=0.7-1.0×10^20 cm^-2K^-1km^-1s 
    consistent with
    observations, insensitive to the evolutionary ISM state
    or radiation field strength if
    
    steady-state
    chemistry is assumed. Due to non-steady-state
    chemistry, younger molecular clouds have 
    slightly lower ⟨X_⟩ and flatter profiles of X_ 
    versus extinction than older ones. The -dark 
    fraction is 26-79%, anti-correlated with the average
    extinction. As the observational beam size 
    increases from 1 pc to 100 pc, ⟨X_⟩ increases by 
    a factor of ∼2. 
    Under solar neighborhood conditions, ⟨X_⟩ in molecular
        clouds is converged at a numerical resolution of 2 pc. However, the total 
    abundance and luminosity are not converged even at the numerical
    resolution of 1 pc.
    Our simulations successfully reproduce the
    observed variations of X_ 
    on parsec scales, as well as the dependence of X_ on extinction
    and the  excitation temperature.




§ INTRODUCTION

Molecular clouds are the birth places of stars. In addition, molecular gas is
the dominant ISM component in dense and shielded environments.
Measuring the properties
of molecular clouds is therefore critical to understanding the ISM and 
star formation in the
Milky Way and beyond. However, the most abundant molecule in the ISM,
molecular hydrogen , is not directly observable in emission 
at typical ISM temperatures due to its low mass
and lack of dipole moment. As a result, the second most abundant molecule,
, is often used as an observational tracer for . The standard
technique employs a conversion factor X_ to relate the observed
velocity-integrated intensity of (J=1-0) line emission 
W_ to the  column density N_, 

    N_ = X_ W_.


Although the (J=1-0) line emission is bright and easy to detect 
with ground based radio telescopes, it is often very optically thick.
Many observational studies have measured X_ by deriving the  mass
independently
of CO emission, via 
dust emission or extinction, gamma-ray emission, or the virial theorem
<cit.>. Surprisingly, the
value of X_ only varies within a factor of ∼ 2 for many molecular
clouds in the Milky Way and local disk galaxies. This has motivated the
adoption of a constant standard X_ conversion factor in the literature, 
X_CO, MW=2× 10^20 cm^-2K^-1km^-1s 
<cit.>.

It is important to note that this standard X_
is an average value for nearby molecular clouds on scales of tens of parsecs.
X_ is empirically known to vary both on small scales, and for molecular clouds in
different environments.
One of the earliest studies of X_, <cit.>, suggested that X_
varies by a factor of a few for molecular clouds in the Milky Way, decreasing
with increasing  luminosity. Recent high-resolution observations have found
that X_ can vary by more than an order of magnitude on parsec scales, 
although the averages of X_ over individual
molecular clouds are within a factor of ∼ 2 of
the standard Milky Way value <cit.>.  
Beyond nearby molecular clouds, 
X_ in the Galactic center is a factor of ∼ 4 lower than
the mean value in the disk <cit.>, 
and similar results are found for the central regions in nearby spiral
galaxies <cit.>. 
High surface density starburst regions have X_ significantly below
X_CO, MW <cit.>.
Observations also indicate X_ can be much 
higher than the standard Milky Way value in low metallicity
galaxies <cit.>. 

Theoretical models and numerical simulations have provided insights into the
X_ conversion factor. 
<cit.> constructed spherical cloud models with a photodissociation
region (PDR) code, and suggested that X_ is only weakly dependent on
the incident far-ultraviolet (FUV) radiation field strength, and insensitive
to the small variations in metallicity up to a reduction of metallicity by a
factor of 5 relative to the solar neighborhood. These models rely on simple
assumptions about cloud structure and kinematics. 
To model molecular clouds with more realistic structure,
many numerical
simulations have been carried out to study 3D turbulent molecular clouds 
with self-consistent, time-dependent chemistry and radiation transfer
<cit.>.
<cit.> and <cit.> found
similar cloud-average 
X_ to the standard observed value (with significant variations on
smaller-than-cloud scales).
<cit.> concluded that
X_ has a weak dependence on gas density, temperature, and velocity, 
and the nearly constant X_ is the result of the limited range of
physical properties found in the nearby molecular clouds. 
However, these
simulations consider molecular clouds to be isolated from the large scale
galactic ISM, and their key physical properties such as the average
density and velocity dispersion are set artificially based on 
the initial conditions of the simulations and prescribed turbulent driving.

In recent years, more efforts have
been made to investigate X_ in global galaxy simulations 
<cit.>. 
With resolutions of tens of parsecs, however, global galaxy simulations
cannot resolve 
substructures in molecular clouds, and sub-grid models are generally required to
estimate the  emission. There is no systematic study of the
dependence of X_ on the numerical resolution in the literature.
Moreover, the comparisons between simulations
and observations are often focused on the cloud-average X_.
Despite the rich observational data, little comparison has been made 
regarding to the variation of X_ within molecular clouds on parsec
or smaller scales. Furthermore, as observations of galactic and extragalactic
molecular gas probe a range of scales, it is important to understand how
X_ may vary with the effective area of a radio beam.

In this paper, we present a new study of the X_ conversion factor in MHD
galactic disk simulations with solar neighborhood conditions and 
1-4 resolutions. The
high-density

clouds are
formed and destroyed self-consistently within the turbulent, multiphase, 
magnetized ISM by gravity and stellar feedback.
In our models, the distribution of
 and  is obtained by post-processing the MHD simulations with chemistry
and radiation transfer.

While ideally all dynamics and chemistry would be self-consistent, 
<cit.> pointed out that the gas temperature is not sensitive
to chemistry in the neutral ISM  <cit.>; as a consequence, dynamical simulations may
still represent ISM structure fairly accurately even if they do not include
time-dependent chemistry.

Using our models,
we investigate the dependence of X_ on numerical resolution,
non-equilibrium (i.e.
non-steady-state) 
chemistry, variation in large-scale ISM structure and star
formation rates, and the observational beam size. 
Our analyses also identify the density and shielding conditions that are
required for  and  formation (which differ significantly) in
realistic clouds, and break down the dependence of W_ on microphysical
properties.
Additionally, we perform detailed comparisons with observations of X_
in nearby molecular clouds at parsec scales.

The structure of this paper is as follows. In Section <ref>, we
describe the method of our simulations and the parameters in the
numerical models. In Section <ref>, we show our results and
comparisons with observations. The main findings of this work are summarized in
Section <ref>.



§ METHOD

To investigate the X_ conversion factor in molecular clouds, we
carry out MHD simulations of galactic
disks, and post-process the results from MHD simulations with chemistry 
to obtain the distribution of molecular gas, including  and .
Then we use
line radiation transfer code to model the  emission from molecular clouds.



 §.§ MHD Simulation

The MHD simulation is performed with the TIGRESS
(Three-phase Interstellar medium in Galaxies Resolving Evolution with 
Star formation and Supernova feedback) framework introduced by
<cit.>.
Here we briefly describe the key physics in the simulations, and refer
the readers to  for more extensive descriptions.

The TIGRESS simulations model a kpc-sized region of a galactic disk where
the turbulent, multiphase, magnetized
ISM is self-consistently modeled with resolved star formation and feedback.
The physics are implemented within the Athena code <cit.>.
The ideal MHD equations are solved in a vertically-stratified local shearing box
<cit.>. Self-gravity from gas and young stars are included by
solving Poisson's equation,
while a fixed vertical gravitational potential
represents  the old stellar disk and
the dark matter halo. Sink particles are implemented to represent star clusters,
and feedback from massive stars  
are included based on a population synthesis model 
<cit.>. Both supernovae in star clusters
and from runaway OB stars are included.
The radiative heating and cooling of the gas are assumed to be optically
thin. The heating of cold and warm neutral gas
is from the photo-electric effect on dust grains; in the simulations the
heating rate is time dependent and scales with the instantaneous FUV luminosity
of the star cluster particles.
The cooling rate is obtained
from the local gas density and temperature using a simple cooling function
appropriate for the ionized and atomic ISM (combination of <cit.> and
<cit.>). 

The simulations self-consistently generate a representation of the
turbulent and magnetized 
three-phase ISM. In the fiducial model with solar neighborhood parameters, 
much of the volume is occupied by hot ionized gas, and
most of the mass near the midplane 
is in the warm and cold neutral medium (WNM and CNM), 
similar to the observed ISM in
the Milky Way and nearby galaxies. Although molecular gas is not explicitly
modeled in the TIGRESS simulations, large structures of dense gas
naturally develop, and in reality molecular gas
would form within the regions of the CNM where the gas is 
dense and shielded. We model the
formation of molecular gas by post-processing the TIGRESS simulations with
chemistry and shielding, which is described in detail 
in Section <ref>.

We adopt the fiducial solar neighborhood model of . 
The simulation domain size is
L_x=L_y=1024 and L_z=4096. The initial gas surface density
Σ = 13 M_⊙^-2. The simulation reaches a quasi-steady
state after t≈ 200. The total mass of the gas in the
simulation slowly declines as the gas turns into stars or leaves the
simulation domain as galactic winds. In this paper, we focus on the simulation
during the time frame t=350-420 when the surface density of the gas is
in the range 9 M_⊙^-2 < Σ < 10 M_⊙^-2.

In order to study the effect of numerical resolution on X_, we
consider the simulation with three different resolutions: Δ x=4, 2, and
1. The 4 simulation starts from t=0 with the initial condition
described in , and runs until t=700. 
To save computational time, we use
an “extraction” method to refine the resolution. We use the
output of the 4 simulation at time t=350 as the initial condition
of the 2 simulation, and run that for 70 (until t=420).
Similarly, we use the output of the 2 simulation at t=378 as the
initial condition of the 1 simulation, and run that for 4 (until
t=382). We also reduce the domain size in the z-direction to L_z=2240
for the 2 simulation and to L_z=896 for the 1 simulation. Because
the scale-hight H∼ 100 for the CNM and H∼ 400 for the WNM, the
simulation domain in the z-direction is big enough
to capture most of the mass in the neutral and
molecular ISM.

When refining from a coarser resolution, it takes some time for
the turbulence to cascade down to smaller scales and create finer structures.
The line-width size relation <cit.>,

    v(l) ∼ 0.7 km/s ( l/pc)^1/2,

gives the expected timescale for turbulent cascade
in the dense ISM:

    t_turb(l)=l/v(l) = 1.4  ( l/pc )^1/2.

We only use
the outputs from the 2 and 1 simulations 4-6 after the
extraction from coarser resolution, allowing sufficient time for the turbulence
to develop at the refined resolution.

The density threshold for sink
particle creation, n_thr, also depends on the resolution of the
simulation. A sink particle is created if the cell is at the local
gravitational potential minimum, the flow is converging, and the density of the
cell exceeds the Larson-Penston threshold <cit.>
suggested by <cit.>,

    ρ_thr≡ρ_LP(Δ x/2) = 8.86/πc_s^2/G Δ x^2.

The typical density threshold at the equilibrium CNM temperature is 
n_thr=2956, 927, and 304 cm^-3 for resolutions
Δ x=1, 2 and 4.[This is
    assuming the heating rate of the CNM to be the solar neighborhood value
    Γ = Γ_0 = 2× 10^-26 erg s^-1 
    ( Equation (8)). However, 
    n_thr is insensitive to the change of Γ: 
    n_thr increases by less then a factor of
two when Γ increases by a factor of ten.]



 §.§ Post-processing chemistry

To model the chemical composition of the gas, we have developed a
post-processing module within the code Athena++ <cit.>.
This module reads the output from TIGRESS simulations and performs
chemistry calculations assuming the density and velocity in each
grid cell is fixed. 
We use the simplified chemical network of <cit.>, which focuses
on the hydrogen, carbon, and oxygen chemistry, and gives accurate 
abundances of  and . 
We assume an initial chemical composition of neutral atomic gas,
with all hydrogen in the form of H, all carbon in C, 
all oxygen in O, and all silicon in 
Si. The initial temperature is the same as the output from
MHD simulations. Then we evolve the chemistry, temperature, and radiation field
(see below) for time t_chem=50, 
so that the chemical abundances of the gas reach

steady state. In other words, we do not self-consistently calculate the
time-dependent gas dynamics and chemistry, but instead 
consider the state in which the
chemistry and temperature have reached a equilibrium, consistent
  with radiative heating and ISM structure
as determined by the MHD simulations.
Because gas cooling is not sensitive to the chemical composition, chemistry
has minimal effect on the gas dynamics <cit.>. However, dust
shielding can reduce the gas heating, and lower the gas temperature by a factor
of ∼ 2 in shielded regions of the CNM where molecular gas
forms.



[This typical reduction in temperature
in high density regions (n_H≳ 10 cm^-3)
is found by comparing the initial
temperature output from the MHD simulation 
to the steady-state temperature from
the post-processing chemistry simulation.














]
In return, gas dynamics can also influence the chemical composition.
For example, the timescale for  formation can be longer than the 
turbulent crossing time in the molecular clouds, which may lead to much lower
 abundance than the equilibrium
values <cit.>.  

The temporal dependence of the chemical state and observable CO
  properties are considered in 
Section <ref>.

The heating and cooling of the gas is calculated simultaneously with chemistry,
with the details described in <cit.>. We slightly modify the
parameter Ñ () for  cooling in <cit.> by setting

    Ñ () = n()/max(⟨| v/ r|⟩,
        v_th/l_esc),

where ⟨ | v/ r| ⟩  is the mean (absolute) velocity gradient
across the six faces of each grid cell in the simulation,
v_th=√(2kT/m()) the thermal velocity of  molecules, and
l_esc=100 the maximum length scale for a photon to escape. Using the
maximum of two terms in the denominator of Equation (<ref>) ensures
that there is a minimum probability for the photon to escape when the local
velocity gradient is small, given a maximal molecular cloud size,
≲ 100. This
formalism is consistent with the large velocity gradient (LVG) and escape
probability approximation we adopted in carrying out the synthetic observations 
of  line emission (Section <ref>). 

In order to compute the photoionization and photodissociation rates in the chemistry
network, a radiation transfer scheme is needed to calculate the reduction of
FUV radiation by dust and molecule shielding. We use the six-ray approximation
<cit.>: in each cell, the radiation field is 
calculated by ray-tracing and averaged over six directions along the Cartesian
axes. The incident radiation field is assumed to come from the edge of the
computational domain
along each ray, and has the initial intensity the same as that in the
MHD simulations (the MHD simulations themselves do not include shielding).
The main advantage of this approach is the low computational cost. When
comparing to ray-tracing along many more different angles, the six-ray
approximation gives reasonably accurate results <cit.>.
Because chemistry and radiation depend on each other, we iterate to solve
the chemistry equations and six-ray radiation transfer.





 §.§ Synthetic Observation of  Line Emission 


To model the (J=1-0) line emission, we apply the publicly available 
radiation transfer code RADMC-3D <cit.> with
chemistry and temperature obtained as described in Section
<ref>.[We set a temperature ceiling
of T=200 K for the temperature input, because  only forms within
the CNM where T≲ 100 K, and a high temperature input 
from the WNM and hot gas introduces additional
computational cost for calculating the  population levels in regions 
where the  abundance is essentially zero. We have tested
using a higher temperature ceiling of 1000 K and confirmed that it gives
the same result.]
We select the mid-plane region |z| < 256,
where almost all molecules are found. 
 is assumed to be the only collisional partner with
, and we use a fixed ortho-to-para ratio of 3:1.[The collisional
coefficients for ortho- and para-  are very similar, and we have tested
that a ortho-to-para ratio of 1:1 gives very similar results.] The synthetic
observations are performed along the z-axis, i.e., the
observer is looking at the galactic disk face-on. This avoids cloud blending, as
all molecular clouds form near the mid-plane of the galactic disk. 

The  population levels are calculated by using the LVG and
escape probability approximation, which is implemented in RADMC-3D by
<cit.>. This approximation allows the population levels to be
calculated locally in each cell. The escape probability is

    β = 1-e^-τ/τ,

and the optical depth τ = min( τ_LVG, τ_EscProb). The LVG approximation gives

    τ_LVG = λ_10^3/8πA_10n_/⟨ | v/ r| ⟩
        f_1 ( f_0/g_0/f_1/g_1 - 1 ),

where A_10 is the Einstein A coefficient 
A_10=7.203× 10^-8 s^-1, n_ the number density of 
molecules, g_0=1 and g_1=3 the degeneracy for J=0 and J=1 levels, 
f_0=n_0/n_ and f_1 = n_1/n_ the fraction of  molecules in
J=0 and J=1 levels, where n_0 and n_1 are the level populations, 
⟨ | v/ r| ⟩ the same as that in Equation
(<ref>).
















The optical depth from the escape probability approximation
is set by a typical length-scale L_EscProb, and can be written in the
same from as Equation (<ref>) by substituting ⟨ | v/ r|
⟩ with √(π)v_tot/L_EscProb <cit.>.
Here v_tot is the total velocity dispersion (see below).
We adopt L_EscProb=100, consistent with the
  line cooling in Equation (<ref>). In our
simulation, the velocity gradient is usually relatively large, and in most
cells τ=τ_LVG.

Ray-tracing is performed after the  level populations are obtained. 
In general, the emission line intensity is determined by radiative
transfer <cit.>:

    I_ν = - I_ντ_ν + S_ντ_ν,

where I_ν is the line intensity at frequency ν, S_ν the source
function, and τ_ν the optical depth. τ_ν depends on the line
profile, which is set by the velocity dispersion

v_tot = √(v_th^2 +
v_turb^2). We include a sub-grid 
“micro-turbulent” velocity dispersion according
to the line-width size relation (Equation (<ref>)),

    v_turb = 0.7 km/s ( Δ x/pc)^1/2,

where Δ x is the resolution of the simulation. We also include a
background blackbody radiation field with temperature T_CMB=2.73 K 
from the cosmic microwave background (CMB).

We run RADMC-3D with a
passband from -20 km/s to 20 km/s (wide enough to include all
 emission) and velocity resolution of
0.5 km/s. RADMC-3D produces spectral position-position-velocity (PPV) cubes of the
(J=1-0) line.
We then interpolate I_ν to a finer velocity resolution of
0.07 km/s, and calculate the total (J=1-0) line intensity in each
observed pixel, W_, by integrating I_ν over
all velocity channels that have emission above the detection limit,
T_det=0.4 K.  This approach matches the typical
velocity resolution and sensitivity in observations of nearby molecular
clouds <cit.>.
 We define the
 “-bright” region as pixels with

 W_ > 0.1  K km/s, and calculate
X_ for each pixel in the -bright region. The average X_,
⟨ X_⟩ = ∑ N_/∑ W_ is also calculated only within
the -bright region, similar to the common approach in observations
<cit.>. We define the fraction of -dark ,

    f_dark≡M_(W_<0.1 K km s^-1)/M_H_2,tot.




 §.§ The beam size in synthetic
observations

The default beam size in our synthetic observations is the same as the
numerical resolution in the MHD simulations. In real observations, the beam size (in
physical units) varies depending on the telescope and the distance of the
object. The dust extinction or emission map used to derive  column
densities typically has coarser resolution than the
 map. To analyze the X_ values, the dust map and  map
are smoothed to a common resolution (usually the resolution of the dust map),
which we refer to as the “beam size”.[Note that this is often
called “pixel size” in observations. We use “beam size” to distinguish from the
“pixel size” determined by the numerical resolution of our simulation and
synthetic radiative transfer grid.]
The velocity resolution and sensitivity also vary in observations. We have
compiled the observational parameters from the literature of X_ observations
in the Milky Way and nearby galaxies in Table <ref>. All the
observations listed used  mass estimation from dust extinction or
emission. We also list ⟨ X_⟩ obtained by the observations
when available. 



We investigate the effect of beam size on X_ in Section
<ref>. The adopted parameters and beam sizes
are listed in Table <ref>, which is
designed to match the typical observational parameters listed in Table
<ref>.
The synthetic observations with default beam size are based on the original
model data (the same as the numerical resolution). To create synthetic maps
with larger effective beam, we first smooth out the PPV cubes
produced by RADMC-3D to the desired  map resolution. Then we match the
corresponding velocity resolution from the default 
0.5 km/s in the PPV cubes, by either interpolating to finer or
integrating
to coarser velocity resolution. We integrate over all velocity channels with
emission above the detection limit T_det and obtain a 2D map of W_ at
the corresponding  map resolution. Then both the map for A_V(N_) and
the map for W_ are smoothed to the common resolution of the beam size, for
which X_ is calculated. We note
that the “beam” is square, not circular.

[We have compared
    results for our square beam
    to the result for a circular gaussian beam, and find that it makes very little difference for
⟨ X_⟩.] 
The -bright region, for
which ⟨ X_⟩ is calculated, is defined as pixels with 
W_ > 3 T_detΔ v, where
Δ v is the width of the velocity channel. 





 §.§ Model parameters

We consider three sets of models designed to study different conditions
that may affect X_: the numerical resolution, non-equilibrium
chemistry,
and variation in the galactic environment (ISM structure and ambient radiation
field). The parameters for our models are summarized in Table <ref>. 
Model names denote changes in numerical resolution (RES-1pc, etc.), chemical
evolution time (TCHEM-5Myr, etc.), and simulation snapshot time (T-356Myr,
etc.). Note
that t is the time for the MHD simulation, and t_chem is the time for
the post-processing chemistry, as detailed in Section <ref>.
RES-1pc and TCHEM-50Myr are two names for the same model, used for clarity 
in different sections discussing the numerical resolution or
evolving 

chemistry. To do a controlled study, we set the incident
radiation field strength χ=1 (in <cit.> units, corresponding
to J_FUV = 2.7 × 10^-3erg cm^-2 s^-1)
for all models that intercompare numerical resolution and
non-equilibrium
chemistry (model IDs starting with RES or TCHEM).
In the set of models for studying the variation in galactic environments 
(model IDs starting with T), χ
is obtained from the star cluster particles as described in .





§ RESULTS



 §.§ Convergence study: effects of numerical resolution

In this section, we investigate the effect of numerical resolution on both
chemistry and X_. An overview of the models RES-4pc, RES-2pc and RES-1pc
is shown in Figure <ref>, and the overall properties of the
models are listed in Table <ref>. As the
resolution increases, more small structures and dense gas forms in the
simulations. The locations of molecular clouds are similar in all three models,
but the small scale filamentary structures within the molecular clouds can only be
resolved in RES-2pc and RES-1pc. As we shall show (Section <ref>), 
at least 2 resolution is needed 
to accurately determine the average X_ in molecular clouds
for the Solar neighborhood conditions of the present simulations.








  §.§.§ Molecular Abundances and Dependence of Chemistry on numerical
resolution

As the numerical resolution increases from 4 to 1, 
a larger fraction of mass in the simulations is in
the dense gas. This is quantified by the increase of f_100 (the fraction of
gas with density n>100 cm^-3) with resolution in
Table <ref>, 
and the density distributions in Figure <ref>.
The density distributions are similar at low densities where the gas is well
resolved.  At high densities, the distribution cuts off near the
density threshold for sink particle creation, where the unresolved dense gas
is converted into sink particles in the simulations.
As resolution increases, the density threshold for sink particle creation 
also increases, allowing denser gas to form.



The change of density distribution with resolution affects the chemical
compositions of the gas. As the resolution increases from 4 to 1,
the total  mass stays nearly constant, but the total  mass increases
by a factor of nearly 3 (Table <ref>).[
In Table <ref>, M_ first decreases slightly
when the resolution increases from 4 pc to 2 pc, then increases again at 1 pc
resolution. This non-linear variation of M_ with resolution is actually a
result of temporal variations in the simulations. Because the supernova feedback
from the sink/cluster particles is stochastic,
simulations with the same initial condition
can develop slightly different density structures over time. We compared M_ and
M_ in models RES-4pc and RES-2pc between the time when they have the same
initial condition (350 Myr) and the time of comparision in Table
<ref> (382 Myr). We found that the  mass in both
models are similar (up to ∼ 20% variations), but the  mass increases
significantly (up to a factor of ∼ 3) in the RES-2pc model. The  and
 mass weighted density histograms at different times also show very
similar features to Figure <ref>. Therefore, the
conclusion from Figure <ref> is robust despite the temporal
variations. ]
The reason for this is evident in Figure <ref>:
most  forms in
the density range of n=10-100 cm^-3, which is already well resolved with
4 resolution. However,
most  forms at n ≳ 200 cm^-3, which is not well resolved with
2, maybe even 1 resolution. Using adaptive mesh refinement (AMR)
models, <cit.> found a resolution of ∼ 0.2 is needed
for the  abundance to converge.




The chemical composition depends not only on density, which
affects the rate of collisional reactions, but also on shielding, which
determines the photodissociation rate by FUV photons. Which factor,
density or shielding, is more important in determining the  and 
abundances in realistic molecular clouds with complex structures? Figures 
<ref> and <ref> plot the
probability density distributions (PDFs) of the
 and  abundances versus density and shielding in each grid cell.


We weight the PDFs by n f_ or n f_, so that the color
scale is proportional to the  or  mass in each bin. Simple volume
weighted PDFs will show distibutions centered at very low density and low
molecular abundances, since by volume most gas is atomic.

We quantify the shielding
by calculating the effective extinction A_V, eff
for the photo-electric heating <cit.>,

    χ_PE≡χexp(-1.8 A_V, eff),

where χ_PE is the actual radiation field intensity
obtained from the six-ray radiation transfer.  

As shown in Figure <ref>, the  abundance has a much tighter
correlation with density than with shielding. This is because  self-shielding
is so efficient that the photodissociation rate of  is very
small in most regions that have a significant amount of . 
In the absence of photodissociation by FUV radiation, the 
 abundance is
then determined by the balance between  formation on dust grains,

    H + H + gr → H_2 + gr,

with a rate coefficient k_gr=3.0× 10^-17 cm^3 s^-1
(assuming solar neighborhood dust abundance),
 formation by H_3^+,

    H_3^+ + e → H_2 + H,

with a rate coefficient k_<ref>,
 destruction by cosmic-rays,

    CR + H_2 → H_2^+ + e,

with a rate coefficient k_CR = 2ξ_H (2.3 f_ + 1.5 f_H), 
and  destruction by H_2^+,

    H_2^+ + H_2 → H_3^+ + H,

with a rate coefficient k_<ref>.
Reactions (<ref>) and (<ref>) are also the main pathways for
H_2^+ destruction and creation. Equilibrium of H_2^+ requires

    f_ k_CR = f_H_2^+ f_ n k_<ref>.

H_3^+ is mainly created by reaction (<ref>), and destroyed by
reaction H_3^+ + e, which forms H_2 + H (reaction
(<ref>)) or 3H with a branching ratio of 0.35:0.65.
Equilibrium of H_3^+ requires 

    f_H_2^+ f_ n k_<ref> = 1/0.35 f_H_3^+
        f_e n k_<ref>.

Finally, equilibrium of fully-shielded  
(Equations (<ref>) - (<ref>)) requires

    f_H n k_gr + f_H_3^+ f_e n k_<ref>   = f_ k_CR + f_H_2^+ f_ n k_<ref>
    
        f_H n k_gr + 0.35 f_H_2^+ f_ n k_<ref>   = f_ k_CR + f_H_2^+ f_ n k_<ref>
    
        f_H n k_gr   = 1.65 f_ k_CR.

In the above, each f is the abundance of a given species relative to hydrogen
nuclei. 

Equation (<ref>) can be solved with the conservation of hydrogen
nuclei 
f_H + 2f_ = 1, giving the equilibrium  abundance as a function
of n, plotted as the green dashed line in the left panel of Figure 
<ref>. This agrees very well with the upper limit of f_ in
the simulations. The spread of 
f_ at a given density is due to the incomplete shielding of FUV radiation 
in some regions where destruction of  from photodissociation brings its
abundance lower than that in completely shielded regions.
This can also be seen in the right panel of Figure <ref>: there is a
large spread of A_V, eff at a given f_, and there are many grid cells
with A_V, eff < 1 and significant  abundance.



Contrary to the case of  abundance, which is determined mostly by density, 
the CO abundance is determined by both density and shielding, 
as shown in Figure <ref>.  forms mainly in regions with 
n ≳ 100 cm^-3 and A_V ≳ 1. This agrees very well with
the results from 1D slab models in <cit.>.
The main reason  and  form under different conditions is 
that the self-shielding of  and
cross-shielding of  by  are much less
efficient than the  self-shielding. As a result,  formation is limited by
photodissociation, and  can only form in regions with 
A_V, eff≳ 1 where the FUV radiation field is sufficiently
shielded by dust. Moreover,  formation also requires higher
densities, as C^+ and He^+ formed by cosmic rays destroy 
at lower densities. Figures <ref> and <ref> again 
show that the  mass in our simulations is converged,
but the  mass is not, due to the lack of resolution for very high density gas 
(see also Figure <ref> and discussion). 



Because  and  formation require different conditions,  is only a
very approximate tracer of .  Figure <ref> shows the distribution of 
density n versus the effective extinction
A_V, eff for each grid cell. At a given density, there is a large range of
A_V, eff. We roughly delineate loci where H, , and  form:
 exists in high density regions, and f_ > 0.5 corresponds roughly to
densities n≳ 30 cm^-3.  forms in denser
and well shielded  regions, and 
f_ > 10^-5 roughly corresponds to n≳ 100 cm^-3
and A_V, eff≳ 1. Figure <ref> clearly
shows that a significant fraction of  would not be traced by  emission 
(see f_dark in Table <ref>). As Δ x
decreases from 4 pc to 1 pc, more and more high density gas is
resolved, as also shown in Figure <ref>. Nevertheless, for all
resolutions considered in our models, there is gas in the three different
regimes – atomic, -bright molecular, and -dark molecular.



To validate that we can accurately simulate chemistry in molecular clouds,
we compare the  column densities N_ in our simulations to that in the UV
absorption observations of diffuse molecular clouds. 
Figure <ref> shows the comparison between the simulations and
observations, as well as the result from the one-sided slab model in
<cit.>.
The x-axis of Figure <ref> is the extinction from only :

    A_V(N_) = 2N_/1.87× 10^21 cm^-2.

To avoid
foreground/background contamination, we compare N_ to 
N_ instead of the total column N.[<cit.> discussed
that the dispersion in observations is much smaller when comparing N_ to N_
instead of N.]
Compared to the simulations, the
one-sided slab model gives higher  abundance at A_V(N_)∼ 1.
This is because the six-ray radiation transfer in the 3D simulations considers
extinction of FUV radiation from all directions along the Cartesian axes, 
which is generally lower than the
extinction only along the z-axis, A_V(N_) (that is, A_V, eff≲ A_V(N_)). At A_V(N_)≪ 1 or 
A_V(N_)≫ 1, the  abundance in the one-sided slab model and 3D simulations
are more similar, because either the FUV radiation is only weakly shielded at low
A_V(N_) so that the photodissociation rate is insensitive to A_V(N_),
or else already completely shielded at high A_V(N_) so that the limiting factor
for  formation is no longer photodissociation. The UV absorption
observations can only be conducted in diffuse molecular clouds with A_V(N_)
≲ 1, and there is a lack of observations at higher extinctions. For the
range of A_V(N_) where the observational data are available, the
RES-1pc simulation successfully reproduce the observed range of N_. 
Lower resolution simulations RES-2pc and RES-4pc also show similar average
values (magenta lines) and range (not shown in the Figure) of N_ at 
A_V ≲ 1. At A_V > 1, models with
lower resolutions start to show that the  mass is not resolved at high
densities.









  §.§.§ Dependence of X_ on Numerical Resolution

To understand the relation between physical properties of molecular clouds and
 emission, a helpful reference point is 
the simple uniform slab model for molecular
clouds. In a uniform slab with constant  excitation 
temperature T_exc, Equation (<ref>) can be integrated, giving

    I_ν = I_ν(0) e^-τ_ν +
        B_ν(T_exc)(1-e^-τ_ν),

where S_ν=B_ν(T_exc), the blackbody
radiation field intensity at temperature T_exc, and I_ν(0) is the
initial impinging radiation field intensity at τ_ν=0. 
[In observations, the intensity is often referred to as the value
after background subtraction I_obs = I_ν - I_ν(0). Then Equation
(<ref>) is often written as 
I_obs = (B_ν(T_exc) - I_ν(0))(1-e^-τ_ν).]
The line intensity I_ν is usually measured in terms of the antenna
temperature (also often referred to as the radiation temperature) in radio astronomy:

    T_A(ν) = c^2/2kν^2I_ν.

In the limit of τ_ν→∞,
Equations (<ref>) and (<ref>) becomes 

    T_A = T_0/e^T_0/T_exc - 1,

where T_0=5.5 K = hν_0/k, with ν_0=115.3 Hz, the frequency of the 
(J=1-0) line.

Typically, the (J=1-0) line profile (in terms of T_A and v)
is not too far from a Gaussian profile, and to first order,
the total  line intensity W_ is determined by two parameters: the
peak of the line profile T_peak and the width/velocity
dispersion of the line σ_v. Under the assumption that the line center is
optically thick so that Equation (<ref>) applies, the observed
peak antenna temperature, T_peak, would be directly related to the
excitation temperature T_line,

    T_line≡5.5 K/ln (5.5 K/T_peak+1).

We use the notation T_line for the excitation temperature 
derived from the line profile (T_peak) to distinguish 
from the true excitation temperature in the molecular clouds T_exc. 
Although T_line=T_exc in a uniform slab cloud as long as
the  line center is optically thick, in real molecular clouds and also in
our numerical simulations, the excitation
temperature along the line of sight is not constant, and T_line serves as
an estimate of the excitation temperature where most  emission comes
from. For T_peak≳ 5.5 K, Equation (<ref>) gives 
T_line≈ T_peak. Another important parameter for the 
line, the velocity dispersion, is calculated using
σ_v≡√(⟨ v^2 ⟩_T_A - ⟨ v ⟩_T_A^2),
where ⟨ v ⟩_T_A≡∫ v T_A v/∫ T_A  v 
is the intensity weighted average of velocity, and similarly 
⟨ v^2 ⟩_T_A≡∫ v^2 T_A v/∫ T_A  v. 

The relations between W_ and T_line or σ_v in models RES-1pc,
RES-2pc and RES-4pc are shown in Figure <ref>.
T_line ranges from ∼ 2 K (from the CMB background)
to ∼ 20 K (from the kinetic temperature of dense gas as discussed below),
similar to the range of excitation temperature
observed in Perseus and Taurus molecular clouds <cit.>.
The velocity dispersion spans a relatively narrow range σ_v ≈
1-2 km/s, and the lower limit for σ_v is set by the sub-grid
micro-turbulence velocity in Equation (<ref>). The
observations of nearby molecular clouds have higher resolutions of 
∼ 0.2-0.4, and therefore a slightly lower but still limited range of
velocity dispersions 
σ_v ≈ 0.8-1.5 km/s <cit.>.  
W_ increases with both T_line and σ_v. For a Gaussian
profile with T_peak≳ 5.5 K, W_ = √(2π)T_peakσ_v
≈√(2π)T_lineσ_v. Because the variation in σ_v is small, 
W_ correlates very well with T_line, except for regions where T_line
saturates around 20 K. There is no saturation
of W_, and W_ keeps increasing with increasing σ_v.  



W_ is largely determined by the excitation temperature, and the excitation
temperature in turn depends on the gas density and temperature.
Figure <ref> shows the excitation temperature T_exc and gas
temperature T_gas versus the gas density in each grid cell. 
T_gas decreases with increasing density, as the gas cooling
becomes more efficient, and heating is also reduced by shielding of the FUV radiation
field in dense regions. On the other hand, 
T_exc increases with increasing density, because the collisional
excitation rate of  is proportional to density, 
and because radiative trapping increases in dense regions. 

The lower solid magenta line in Figure <ref> shows the median
T_exc from model RES-1pc as a function of density.  
T_exc only reaches equilibrium with T_gas at 
n ≳ 400 cm^-3, implying that local thermal equilibrium (LTE)
approximation would fail in most regions.


At a given density, T_exc is higher at lower resolutions for two reasons.
First, the velocity gradient | v/ r| is smaller at lower resolutions,
leading to higher τ_LVG and thus lower escape probability β 
and higher T_exc at a given density
(See Equations (<ref>) and (<ref>) below). Second, 
at lower resolutions, 
less high-density gas is resolved, and a larger fraction of the  gas is 
found in lower-density gas (see Fig. 
  <ref>b).  This shifts the distributions of
T_exc and T_gas in Figure <ref> to the left at
lower resolution in models RES-2pc and RES-4pc (dashed and
dotted magenta lines).













In general, thermalization is expected
  for densities above a critical value at which collisional deexcitation
  exceeds spontaneous emission.
For  collisions with , the collisional deexcitation rate is
n_H_2 k_10 for 

    k_10≈ 6× 10^-11( T_gas/100 K)^0.2 cm^3 s^-1,

at 10 K≲ T_gas≲ 250 K <cit.>. The spontaneous emission rate is β A_10, where
Equation (<ref>) gives the escape probability β, so that 

    n_crit = β A_10/k_10.

For T_gas=20 K, Equation (<ref>) gives
A_10/k_10≈ 2.1× 10^3  cm^-3.

With increasing density, the optical depth
τ_LVG increases, leading to
decreasing β (Figure <ref>);
at large τ_LVG, β≈ 1/τ_LVG. 
For model RES-1pc, we fit the average
τ_LVG at a given
density with a broken power-law (magenta line in Figure <ref>):

    τ_LVG   = 2.4× 10^-5 (n/cm^-3)^2.3, 
                 n < 350 cm^-3
    τ_LVG   = 0.21 (n/cm^-3)^0.73,      n ≥ 350 cm^-3.
 

Combining Equations (<ref>) and (<ref>) yields
n_crit∼ 300  cm^-3.  
Thus, in regions where n ≳ 400 cm^-3,
the (J=1) level is expected to be thermalized, and this is indeed
consistent with the median T_exc for model RES-1pc.[
We note that τ_LVG depends on the density and velocity structure,
which is resolution dependent, so the density for thermalization
is not expected to be the same for models RES-2pc and RES-4pc as for model
RES-1pc. In fact, the velocity gradient | v/ r| 
is smaller at lower resolutions, leading to higher average τ_LVG,
and lower density for thermalization in
models RES-2pc and RES-4pc (see Equation (<ref>) and discussions of
Figure <ref>).


]





The dependence of T_exc on n can be understood in a
simplified 2-level system model. The excitation temperature is defined as

    T_exc≡T_0/ln(n_0/g_0/n_1/g_1)  .

With the escape probability approximation, the level populations are given by
<cit.> 

    n_1/n_0 = n_c k_01 + g_1/g_0β A_10
        n_γ^(0)/n_c k_10 + β A_10(1 + n_γ^(0)),

where

    k_01 = g_1/g_0k_10e^-T_0/T_gas,

n_c is the number density for the collisional species, 


and
n_γ^(0) = 1/(e^T_0/T_CMB-1) the background
incident radiation field from the CMB.



If the CMB terms are negligible,
Equation (<ref>) becomes

    T_exc = T_gas/1+ T_gas/T_0ln (1 +
      β A_10/n_c k_10).

For β/n_c small,
T_exc→ T_gas.




The excitation temperature can be estimated as a function of density by
Equations (<ref>), (<ref>), (<ref>) and
(<ref>) (assuming τ=τ_LVG in Equation
(<ref>) and using the average value of T_gas at a given density).



The analytic 2-level system 
approximation for simulation RES-1pc (green line)
agrees well  with the result from radiation transfer
by the RADMC-3D code (lower solid magenta line)
at low and high densities,
while there are differences within a factor of two at intermediate densities n∼
300–1000 cm^-3.









This is because the  rotational levels J=1, 2, and 3 
have energies of 5.5, 16.6, and 33.2 K, all lower or comparable to
the gas temperature. Indeed, there are significant populations in the J≥ 2 levels,
as expected given that T_gas > 5.5 K. The
analytical expression in Equation (<ref>) only takes into account 
the J=0 and J=1 levels, and therefore cannot predict the
excitation temperature very accurately. At low and high densities the
differences are small because the excitation temperature there is determined by 
the background CMB temperature or the gas temperature as the  rotational
levels approach LTE.  Nonetheless, the analytical 2-level
system approximation agrees with the general
trend from the radiation transfer calculations, 
and gives some insight into the relation between T_exc, T_gas and n.
As a further test, we ran RADMC-3D only including 
    the first J=0 and J=1 rotational levels of ,
and found that it can indeed reproduce the analytical result of the 2-level
system model. Figure <ref> shows this comparison.





The relation between W_ and T_line, as well as the
relation between T_exc and density, give rise to the strong correlation between
W_ and the average (mass weighted) density ⟨ n ⟩_M
along the line of sight (Figure <ref> left panel).
Moreover, we found that in the simulations, N_ increases
systematically with ⟨ n ⟩_M (see Figure
<ref>).
This results in a correlation between W_CO and N_ 
(Figure <ref> right panel). Although X_ is measured
in terms of W_ and N_, there is a smaller dispersion in the
correlation between W_ and ⟨ n ⟩_M. This suggests that the 
 emission is more fundamentally a measure of  density than column
density.





The effect of numerical resolution on W_ is already evident in Figures
<ref> and <ref>. As the resolution increases,
more high density gas forms in the simulation, and thus there are
more pixels with high
W_. Numerical resolution also has an effect on X_, as shown in 
Figure <ref>, the histogram of
X_CO,20=X_CO/(10^20 cm^-2K^-1km^-1s) weighted by
W_CO. The average X_ in a certain region can be written as:

    ⟨ X_⟩ = ∑ N_/∑ W_
        = ∑N_/W_ W_/∑ W_
        = ∑ X_ W_/∑ W_.

In other words, ⟨ X_⟩ is simply the W_
weighted average of X_ in each pixel. Therefore, the peak of the histogram in
Figure <ref> roughly indicates the average X_ in the whole
simulation domain. The distributions of X_ in models RES-1pc and RES-2pc
are very similar, with a slightly higher peak in RES-1pc. As a result, 
⟨ X_⟩ is almost the same in RES-1pc and RES-2pc (Table
<ref>). The model RES-4pc,
however, peaks at larger X_ than the higher resolution models, and therefore
has a higher ⟨ X_⟩. This is because the peak of X_
distribution, X_CO,20≈ 0.5, comes from regions with moderatly high
density n≈ 100 cm^-3 and  emission
W_≈ 40 K· km/s, which can only be resolved at a
resolution finer than 2 (see histograms of ⟨ n ⟩_M
and W_ in Figure <ref>). 
Therefore, we conclude that a numerical resolution of at least 2 is needed
in order to resolve the average X_ in molecular clouds for solar
neighborhood conditions.



Finally, we compare the distribution of W_ versus A_V(N_) in model RES-1pc
to observations of the Orion A and B molecular clouds by <cit.>, as
shown in Figure <ref>. Considering the noise level in
the observation, we use a higher 
threshold of W_>1 K· km/s to compare to the -bright region in Orion.
Because most  emission comes from regions with
W_≫ 1 K· km/s, X_ is not sensitive to the W_
threshold.
Our simulations shows a similar distribution of W_ versus A_V(N_) to
that in Orion.  The dispersion of W_ at 
a given A_V(N_) is large, as much as more than an order of magnitude at
low A_V(N_). However, despite the large dispersion
of W_, the average X_ (which is inversely proportional to the slope)
in different A_V(N_) bins is very similar, only varying by a factor of
∼ 2. Similar features are observed in many Milky Way molecular clouds by
<cit.>, and we discuss this in more detail in Section
  <ref>.




There are also some differences between the simulation and observations. 
The average X_ in RES-1pc is a factor of 1.4 lower than that in Orion.







As we shall show based on other analyses and comparisons, the 
typical X_ in our simulations is about a factor of ∼ 2 lower than the
standard Milky Way value;
we discuss possible reasons 

for this discrepency at the end of Section <ref>.
We also note that 
because the observation in <cit.> has a higher
spatial resolution of ∼ 0.2, there are more pixels at A_V(N_)
≳ 4 in the observation in Figure <ref>.
The simulation has more pixels at 
W_≳ 60 K· km/s, a result of the slightly higher velocity
dispersions (see Figure <ref> and discussion). 
In spite of these differences, the
general good agreement between the models and observations indicates
that the simulations can succesfully reproduce the basic physical properties of
observed molecular clouds. 





 §.§ Non-equilibrium Chemistry

The realistic ISM is highly dynamical: turbulence constantly
creates and disperses molecular clouds, and moves gas to environments with
different density, temperature and radiation field strength.
As a result, non-equilibrium
chemistry is likely to be important, especially in low density diffuse gas
where the chemical timescales are long compared to the dynamical timescales. This
is especially an issue for .  Molecular hydrogen
can form in low density gas due
to its effective self-shielding, but its formation timescale, 
t_≈ 10  (n/100 cm^-3)^-1 <cit.>, can be
longer than the dynamical timescale (Equation (<ref>)). Because
 formation chemically relies on the existence of , the  abundance
in molecular clouds can also be far from equilibrium. 
In this section, we carry out comparisons between models with
different t_chem (model IDs start with TCHEM in Table <ref>)
to investigate the effect of non-equilibrium chemistry on X_.

Both  and  abundance increase over t_chem, 
reaching a steady state at t_chem≈ 50, 
as shown in Figure <ref>. Over timescales relevant for
clouds of size ≳ 10 (Equation (<ref>)),
there is a larger increase in the 
 abundance than : From t_chem=5 to 50, the 
abundance increases by a factor of ∼ 3, while  abundance increases
only by ∼ 30%.



The difference in the evolution of  and  abundance 
comes from their different distributions.
As shown in Figure <ref>,



both  and  abundances 
are closer to equilibrium at higher densities, as the rate of
collisional reactions increases with density.
In fact, at a given density in the
  range ∼ 40-400 cm^-3, at 5 the abundance of 
  is closer than the abundance 
  of  is to its final value.
However, in equilibrium most of the  is in gas at
intermediate densities
n≈ 10-100 cm^-3, whereas most  is in gas at high densities 
n≳ 200 cm^-3 (Figure <ref>). This leads to
a shorter timescale for the overall  abundance to reach equilibrium than
. Since the  luminosity also increases much less than the  mass,
this leads to a lower
X_ value at early t_chem (Table <ref>).



Non-equilibrium chemistry also has an effect on the distribution of 
W_ vs. A_V(N_). For model TCHEM-5Myr (Figure <ref>), 
the distribution of the pixels are shifted to the left compared to that in
TCHEM-50Myr (left panel of Figure <ref>). 
This is because W_ is close to equilibrium, but N_ is a factor of ∼ 2
smaller than the equilibrium values, for the same reasons discussed above.
Moreover, the distribution of W_ vs. A_V(N_) in TCHEM-5Myr shows some hints of
a plateau for W_ at high A_V(N_), especially in the binned average
value of W_, which is not present in TCHEM-50Myr. This implies that
younger clouds may not only have lower X_ on average, but also different
distributions of W_ vs. A_V(N_) compared to older clouds. We discuss
this further in Section <ref>.


Note that N_ includes all  along the line of sight, both in high density
clumps where  forms, and in the foreground/background
low density envelopes with only  and no .
Because most 
(in equilibrium) 
lies in these low density envelopes, the fractions of
 in -bright and -dark regions increase by similar proportions
with t_chem, and f_dark stays constant from 
t_chem=5 to t_chem=50 (Table <ref>).





 §.§ Variations in Galactic Environments

Galactic environment fundamentally impacts the molecular content of the ISM.
Supernova feedback creates and destroys molecular clouds, shocks and turbulence
shape molecular clouds in different morphologies, and the
radiation field varies with the star formation activities. Some of these effects
can be seen visually in Figure <ref>. The morphology of molecular clouds
varies from dense concentrated structures (such as in T-356Myr), to more diffuse,
smaller clouds (such as in T-406Myr). The mass and number of young clusters also
changes over time, reflecting the variations in the star formation rate.
To quantify the effect of time-varying galactic environment on X_, we compare 
models with 2 resolution
at different times during the galactic evolution (model IDs start with T
in Table <ref>). As discussed in Section <ref>, 
the average X_ is well resolved with a resolution of 2 in
these simulations. 



A summary of models T-356Myr – T-416Myr is listed in
Table <ref>. In these models, 
M_ and L_ vary by factor of ∼ 3, and the incident radiation field 
strength varies by a factor of ∼ 8. However, 
despite these large variations in the environment, ⟨ X_⟩
stays almost constant, changing only by ∼ 40%. We found no strong correlation
(coefficient of determination R^2<0.4 in linear regression) between 
⟨ X_⟩ and M_, the radiation field strength χ, or 
the average extinction from  in -bright
regions ⟨ A_V ⟩_CO. <cit.> measured
⟨ X_⟩ in individual Milky Way molecular clouds 
using γ-ray observations, and they also found 
no strong correlation between ⟨ X_⟩ and M_ or 
⟨ A_V ⟩_CO.[<cit.> shows a correlation between 
⟨ X_⟩ and ⟨ A_V ⟩_CO with R^2≈ 0.6.
However, this relation is largely driven by one outlier, the Perseus cloud,
which has much lower ⟨ X_⟩ and higher ⟨ A_V
⟩_CO than the rest of the sample. Excluding the Perseus cloud, we
found no strong correlation (R^2≈ 0.3) between 
⟨ X_⟩ and ⟨ A_V ⟩_CO for the rest of their
sample.] 



  <cit.> found a slight anti-correlation of
⟨ X_⟩ and ⟨ W_⟩:
⟨ X_⟩∼ -0.051 ⟨ W_⟩.
We similarly found a slight anti-correlation 
(Figure
<ref> left panel), with
⟨ X_⟩_20 = -0.011 ± 0.005 (⟨ W_⟩/K km s^-1) + 1.0 ± 0.09, where the uncertainties represent the
90% confidence intervals for the fitted slope and intercept. 
The slope of the linear fit is very shallow, 
and ⟨ X_⟩ is not sensitive to the
change of ⟨ W_⟩. 
We note however that <cit.>
focuses on the nearby low mass molecular clouds with
much lower values of ⟨ W_⟩≈ 2-10 K· km/s
than ⟨ W_⟩≈ 10-20 K· km/s 
in the GMCs in our simulations, and therefore may not be directly
comparable to our results.


Large scale galaxy simulations by 
<cit.> found a similar trend that ⟨ X_⟩
decreases with increasing ⟨ W_⟩, although the range of
⟨ W_⟩ is much larger in their simulations as they consider a
wide range of galactic environments. <cit.> found that the
⟨ X_⟩–⟨ W_⟩ relation is caused by the
increase of gas temperature and velocity dispersion at high ⟨ W_⟩, which leads to a faster increase of W_ than N_, resulting
in the decrease of X_. Similarly, we found that the snapshots in our
simulations with higher ⟨ W_⟩ 
also have larger velocity dispersions, 
although the gas temperature is roughly constant in
the  forming regions in our models (see discussion of Figure 
<ref> in Section <ref>). Interestingly, this
is also consistent with the fact that the galactic center molecular clouds 
have larger velocity dispersions and lower X_ compared to the solar
neighborhood clouds. We plan to carry out numerical simulations with



galactic-center-like environments 
in the future to study the variation of X_ in
detail.



Unlike ⟨ X_⟩, the fraction of -dark ,
f_dark, does show
significant variations and a strong correlation (R^2=0.6) with
⟨ A_V ⟩_CO (Figure <ref> right panel). 
Linear regression gives f_dark = -0.31± 0.14 ⟨ A_V
⟩_CO + 1.0± 0.2, where the uncertainties represent the
90% confidence intervals for the fitted slope and intercept.
f_dark increases with decreasing ⟨ A_V ⟩_CO. 
In other words, there is more -dark  in more diffuse molecular
clouds, which is not surprising as  forms in denser gas than .
The same trend was identified in the
simplified spherical molecular cloud model by <cit.>.[The result
from <cit.> shown in Figure <ref> 
is taken from their model with
metallicity Z'=1.9 and incident radiation field χ=10. <cit.>
found that f_dark is not sensitive to Z' or χ in their studies.]
We note that <cit.> uses a slightly different definition of -dark
    , and we use Equation (<ref>) to translate their
definition to ours.
We have also performed an experiment by running the
T-381Myr model only varying the radiation field strength, and found
f_dark stays constant over χ=0.4-3.5, confirming the result from 
<cit.> that f_dark is not sensitive to χ.



Another comparison of X_ with observations is shown in 
Figure <ref>, where the X_ in each pixel is plotted against
T_line. Comparing to the California cloud observed by <cit.>,
our simulations shows a similar slope for the relation between X_ and
T_line at T_line>6 K (the observational data are not
available at lower T_line). However, the 
value of X_ at a given T_line is about a factor of
∼ 4 lower than the observations. One reason for this discrepancy may be
that <cit.> observed (J=2-1) line and
assumed a fixed line ratio of W_(J=2-1)/W_(J=1-0)=0.7, and this ratio is
very uncertain. As discussed below in more detail, 
generally different observations and
also our simulations show similar trends for the variations in X_, but the
absolute value of X_ can differ by a factor of a few. 



Using all of the simulation models, a summary of X_ as a function of
A_V (N_) and comparison with observations is shown in Figure <ref>. 
Because of the large uncertainties in observations of X_ 
at low A_V (N_), we only plot the data at A_V (N_) > 1. 
Both in our simulations and the observations, there is a factor of ∼ 2
variation in X_ over A_V (N_)=1-12.
Simulations with t_chem=50 (RES-1pc,
RES-2pc, T-356Myr – T-416Myr) show a
decrease of X_ at A_V (N_) ≲ 3, regardless of the resolution and
variations in galactic environments. Similar trends can be seen in the
observations of Orion molecular clouds by <cit.> and <cit.>. 
In contrast, the TCHEM-5Myr model shows a flatter profile at
A_V (N_) ≲ 3 and a slight increase of X_ at A_V (N_) > 3.
Interestingly, the California cloud observed by <cit.>
also shows a similar trend. Compared to Orion, the California cloud 
has similar mass and distance, but an order of magnitude lower star formation rate,
and therefore is believed to be much younger <cit.>. 
This has interesting implications that the profile of X_ as a function of
A_V(N_) may be used as an indicator of the age of molecular clouds.

Although the trend for the correlation between X_ and A_V (N_) is
similar in our simulations and observations, there is a discrepancy in the
absolute value of X_. This may be due to systematic
errors in either observations or simulations. One major uncertainty in
observations of X_ comes from the assumptions in deriving N_.
Estimations of  based on γ-ray emission 
systematically give a factor of 2 lower X_ than dust-based methods,
consistent with the value of X_ in this paper 
<cit.>.[The 
 observation by
<cit.> in Figure <ref> is averaged over the molecular
clouds instead of individual pixels in a given A_V(N_) range. Nonetheless, it
indicates the systematically lower X_ in γ-ray observations.] Even
within the dust-based methods, the estimate of X_ in Orion A 
based on dust emission 
is a factor of ∼ 2 higher in observations of <cit.> 
compared to that in <cit.> based on dust extinction.
As another example, the X_ in Perseus measured 
by <cit.> (dust emission) is a factor of ∼ 7 lower than that
in <cit.> (dust extinction).

Several possible factors 
can contribute to the systematics in dust-based observations:
different assumptions of the dust to gas ratio, uncertainties in 
foreground/background subtraction, and different resolutions/beam size
(although the resolution effect is relatively mild, as
noted by <cit.> and discussed in Section
<ref>). <cit.> discussed in detail for the
case of Perseus molecular cloud, that all these
factors can indeed lead to a different estimate of X_. 
  
Sample differences in observations may also play a role. Most
observations of X_ are for nearby low-mass star forming regions, while
most molecular clouds in the Milky Way and our simulations are forming or close
to high-mass stars. The feedback form high-mass stars may lead to slightly
higher velocity dispersions, and lower X_. The only nearby high-mass star
forming molecular cloud is Orion, and it does have a lower value of
X_ compared to the Milky Way average (Figures <ref> 
and <ref>).

For the numerical
simulations, the main uncertainties lie in the assumptions of equilibrium
chemistry and the sub-grid model of micro-turbulence in calculating the 
emission. 


As a further test,
  we produced synthetic observations of
model RES-2pc with half of the fiducial micro-turbulence velocity and no
sub-grid micro-turbulence (only thermal line-broadening on the grid scale), and
found that the values of X_ increase by a factor of 1.4 and 1.8.
Therefore, the uncertainty in sub-grid micro-turbulence may account for part but not
all of the discrepencies in X_ between our simulations and observations.
Future AMR simulations with higher numerical resolution and non-equilibrium chemistry
will be able to provide more insight into these issues. 





 §.§ Dependence of X_ on the Observational Beam Size


Observation of molecular clouds often have different physical beam
sizes/resolutions, which depend on the telescope 
as well as the distance of the object.
In order to investigate the effect of observational resolution on X_, we
smooth the synthetic observations to different beam sizes as described in
Section <ref>.

⟨ X_⟩ increases by a factor of ∼ 2 as the beam size
increases from ∼ 1 to ∼ 100 pc, 
as shown in Figure <ref>.
This is a result of
the -dark . The total  emission remains the same as the beam
size increases, because the detection limits for different beam sizes (Table
<ref>) are generally sensitive enough to detect most of the 
emission. This is not surprising as the sensitivity in observations are
designed to serve the purpose of accurately measuring the  emission.
However, the  emission is smoothed out spatially as the beam size
increases, resulting in a larger area of -bright regions. Although the total
mass of  remains the same, because X_ is calculated only within
-bright regions, a larger area of -bright regions leads to a larger
fraction of  mass accounted for, and therefore an increase of X_. 
This is clearly illustrated in Figure <ref>, showing the
correlation between f_dark and ⟨ X_⟩. 
From beam sizes of ∼ 100 to ∼ 1 kpc, some simulations show
a continued increase of X_ (e.g. T-401Myr), but some simulations with
more diffuse molecular clouds (e.g. T-381Myr) start to have part or all of
their
 emission falling below the detection limits, leading to a non-detection
of W_ or reduction of X_. This suggests that some diffuse molecular
clouds may not be detected with a beam size coarser than ∼ 100 in
extragalactic observations.

In Figure <ref>, 
we plot the observations of ⟨ X_⟩ in 
Milky Way molecular clouds and nearby galaxies (Table <ref>). 
Because of the large uncertainties in the observations (as discussed above, and 
also seen directly in the different ⟨ X_⟩
from two Perseus observations) and dispersions
of ⟨ X_⟩ in different molecular clouds,
we cannot identify any obvious trend for the X_ variation with beam size.
The general range of X_ in the simulations is similar to the
observations.









§ SUMMARY

In this paper, we theoretically model the X_ conversion
factor by post-processing MHD galactic disk ISM simulations with chemistry and
radiation transfer to produce synthetic observations of molecular clouds.
We conduct detailed analyses of the dependence of molecular abundances and
observed line strengths on ISM conditions, and also consider numerical and
observational effects on calculated and measured X_.
Our main findings are as follows:


    
  *  is only a very approximate tracer of . In our
        simulations, most  forms at 
        intermediate densities n≈ 10-100 cm^-3, 
        but most  forms 
        at higher densities n≳ 200 cm^-3 
        (Figure <ref>). The  abundance is
        determined mostly by density, while the  abundance by dust
        shielding (Figures <ref>, <ref>).
        With a 2 numerical resolution,  abundance is
        converged, but  is not. Although there is considerable scatter, 
        the mean relation between the  and 
        column densities in the simulations are in agreement with 
        observations of UV absorption spectra (Figure <ref>).
    
  * For  emission, the high optical depth of the line further
        complicates the observable relation to . 
        On parsec scales, W_ is largely determined by the mean
        excitation temperature of  (Figure <ref>),
        which is in turn determined by the mean gas density. 
        Thus, W_ most directly probes the mean gas density along the line
        of sight. However, for the turbulent clouds in our simulations, the
        mass-weighted mean volume density along a line of sight tends to be
        correlated with column density. This leads to a correlation between 
        W_ and N_ (Figure <ref>).
    
  * A numerical resolution of at least 2 is needed in order to
        resolve the average X_ in molecular clouds for solar neighborhood
        conditions (Figure <ref>). In our simulations with
        environmental conditions similar to the solar neighborhood, we
        found ⟨ X_⟩ = 0.7-1.0 × 10^20 cm^-2K^-1km^-1s, about a factor of 2 lower than the
        estimate from dust-based observations, and consistent with the
        X_ from γ-ray observations. The value of 
        ⟨ X_⟩ is not sensitive to the 
        variations in molecular cloud mass,
        extinction, or the strength of the FUV radiation field 
        (Table <ref>).
    
  * We found the -dark  fraction f_dark=26-79%, which
        has an anti-correlation with the average extinction of molecular
        clouds (Figure <ref> right panel).
    
  * The chemical timescale for  abundance to reach 
        equilibrium is longer than that for  (Figure
        <ref>), because of differences in characteristic
        densities. As a result, younger
        molecular clouds are expected to have lower ⟨ X_⟩
        values and flatter profiles of X_ versus extinction
        compared to older molecular clouds 
        (Figures <ref>, <ref>).
    
  * As the observational beam size increases from ∼ 1 to
        ∼ 100, ⟨ X_⟩ increases by a factor of 
        ∼ 2, due to the decrease of the -dark  fraction (Figures
        <ref>, <ref>).
    
  * Our numerical simulations successfully reproduce the observed 
        variations of W_ on parsec scales, as well as the 
        trends for the dependence of
        X_ on extinction and the  excitation temperature.
        However, the value of X_ in our simulations is
        systematically lower by a factor of ∼ 2 compared to dust-based
        observations (Figures 
    <ref>, <ref>, <ref>). 


The overall
agreement between our numerical simulations and observations of Milky Way
molecular clouds give us
confidence that similar simulations can be used to probe the X_
conversion factors in different environments, such as the Galactic center, 
low metallicity dwarfs, and extreme star-forming systems (ultra luminous
infrared galaxies and high redshift galaxies). 
In a follow-up study, we will investigate the
properties of individual molecular clouds in our simulations.
In the future, we also plan to integrate full non-equilibrium
chemistry with the MHD simulations. 



§ ACKNOWLEDGMENT

This work was supported by grants NNX14AB49G from NASA, and AST-1312006
and AST-1713949 
from the NSF.
We thank the referee for helping us to improve
  the overall quality of this paper,
Jim  Stone, Kengo Tomida and Christopher J. White for
making the code Athena++ available and their help in developing the chemistry
module in the code, Mark H. Heyer for
providing the observational data from <cit.>, Adam K. Leroy for
providing the data from <cit.> and helpful
discussions of comparisons with observations, Shuo Kong for providing the data
from <cit.>, and Simon C. O. Glover for
suggesting an investigation in the effects of the numerical resolution.



§ DEFINITIONS OF -DARK 

In this paper, we define the -dark  as the molecular gas without
 emission along a given line of sight. <cit.> uses a slightly
different definition in their spherical cloud model, 
and they refer the -dark  as the molecular gas
outside of the optical depth τ_CO=1 surface. 
Their definition of -dark  includes the 
along the line of sight in the projected -bright areas on the plane of the
sky, as long as it is outside
of the τ_CO=1 surface (see their Figure 1). In other words, the definition of
<cit.> is in 3D physical space while our definition is in 2D
observational space.

To compare the result from <cit.> to our simulations, we need to
translate their definition of -dark  fraction, denoted by f_DG
(their Equation 1) to our definition denoted by f_dark (Equation
(<ref>) in this paper). Below we derive the relation between
f_DG and f_dark. We refer the readers to Figure 1 in
<cit.> for a useful illustration for this derivation.

From Equation (<ref>), f_dark can be written as:

    f_dark = M_-M_br/M_ = 1 -
        M_br/M_,

where M_ is the total  mass (same as M_H_2, tot 
in Equation (<ref>)), M_br is the mass in -bright
areas on the projected sky. From Figure 1 in <cit.>, 
M_br=M_CO + M_DG, where M_ is the mass with r < R_,
and R_ is the radius of the cloud where τ_ = 1.
M_DG is the mass that lies within R_ in the 2D projected sky,
but outside R_ in the 3D cloud. Compared to the
definition in <cit.>,

    f_DG = 1- M_/M_,

M_DG is the part of the cloud that <cit.> considered to be
-dark, but we do not.

<cit.> assumes
the cloud has a density profile n(r) = n_0 (r_0/r), where n_0 and r_0 are
constants. This gives:

    M_ = ∫_0^R_ 4π m_H n r^2  r = 2π n_0 m_H r_0
        R_^2,

and similarly,

    M_ = 2π n_0 m_H r_0 R_^2,

where m_H is the mass of the hydrogen atom.
M_DG can be estimated by

M_DG≈ 2π R_^2 Σ_DG,
where Σ_DG = m_H∫_R_^R_ n  r = n_0 m_H
r_0 ln(R_/R_). Therefore,

    M_br = M_ + M_DG≈ 2π n_0 m_H r_0 R_^2[1 +
        ln( R_/R_) ].

Equations (<ref>) – (<ref>) then gives the relation
between f_DG and f_dark:

    f_dark≈ f_DG - 1/2(1-f_DG)ln(1/1-f_DG).





§ TEST OF THE RADMC-3D CODE

Figure <ref> shows a test for the RADMC-3D radiation
transfer code.  The level populations of  are solved with only the first
two rotational levels instead of the default 41 levels. The analytical model
uses Equations (<ref>) and (<ref>) to compute T_exc
versus n. Note that because τ_LVG depends on level populations (see
Equation (<ref>)), the average values of τ_LVG in this case 
are slightly larger than that given by Equation (<ref>). 






apj

