ASTROCHEMICAL CORRELATIONS IN MOLECULAR CLOUDS

We use a turbulent hydrodynamic simulation to model a typical Milky Way cloud. The simulation was performed using the Adaptive Mesh Refinement (AMR) code orion (Truelove et al. 1998; Klein 1999). It was previously discussed in Offner et al. (2013) and Offner et al. (2014) as simulation Rm6, so we only briefly summarize its properties here.

The simulation represents a piece of a typical local molecular cloud. The domain is 2 pc on a side and contains 600 M_☉. The gas temperature is 10 K and the 1D gas velocity dispersion is 0.72 km s⁻¹ such that the cloud satisfies the linewidth–size relation (e.g., McKee & Ostriker 2007). The gas turbulence is initialized by adding random velocity perturbations with wave numbers k = 1.2 for two box crossing times without self-gravity.

Once self-gravity is turned on collapse proceeds. The base grid is 256³ but four additional AMR levels are inserted to ensure that the gas obeys the Jeans criterion with a Jeans number of 0.125 (Truelove et al. 1997). On the highest level (Δx_min = 0.001 pc), sink particles ("stars") are added when the gas exceeds the maximum density (Krumholz et al. 2004). The output we analyze in this study is at 1 global free fall time at which point ∼18% of the gas resides in stars.

The chemical abundances were computed using 3d-pdr, a three-dimensional photodissociation code coupled with a full chemical network (Bisbas et al. 2012). Offner et al. (2014) presented an analysis of the line emission based upon the abundances of molecular hydrogen, atomic carbon, and carbon monoxide (run Rm6_1.0_12_1f_a). In this study, we perform an additional analysis of the most common astrochemical species in the chemical network. Here, we describe the astrochemistry calculation procedure and refer the reader to Bisbas et al. (2012) and Offner et al. (2014) for additional details.

We irradiated the simulated cloud by an isotropic 1 Draine FUV field, where "1 Draine" is the standard interstellar radiation field (Draine 1978). 3d-pdr uses the hydrodynamic densities and assumed external field to compute the temperature and abundance distribution for points in the cloud with densities 200 ⩽ n ⩽ 10⁵ cm⁻³. Below n = 200 cm⁻³ we consider the gas as ionized, using limiting conditions on abundances and gas temperature, whereas above n = 10⁵ cm⁻³ we assume it is fully molecular. We do not calculate the chemistry in those two regimes. 3d-pdr computes the radiation field using a resolution of 12 healpix rays (Górski et al. 2005), emanating from each grid point. The input grid is the density field of the hydrodynamic calculation resampled to a resolution of 256³/12, i.e., a new grid comprised of every twelfth data point. (See Offner et al. 2013 for a discussion of spatial resolution convergence.)

3d-pdr employs the UMIST2012 chemical database (McElroy et al. 2013), which includes 215 species and approximately 3000 reactions. The calculation includes the formation of H₂ on dust grains following Cazaux & Tielens (2002, 2004), photodissociation of H₂ and CO and self-shielding. The initial elemental abundances are [He] = 1.0 × 10⁻¹, [C] =1.41 × 10⁻⁴, [O] =3.16 × 10⁻⁴, [Mg] =5.1 × 10⁻⁶, [S] =1.4 × 10⁻⁶, and [Fe] =3.6 × 10⁻⁷ (Asplund et al. 2009), which are similar to values estimated for local molecular clouds. We adopt a cosmic ionization rate of 5 × 10⁻¹⁷ s⁻¹, which is the average value found in the Milky Way. 3d-pdr solves the chemical network in equilibrium, with the final time parameter representing the time allowed for the chemistry to come to equilibrium. We analyze the calculation after advancing the chemistry to equilibrium at 10 Myr.

The 3d-pdr calculation does not consider shock chemistry or dust-grain chemistry, which includes the freeze-out of species such as CO onto dust grains, surface reactions (with the exception of H₂ formation), or release of grain mantle species into the gas phase by means of evaporation, photodesorption or desorption (e.g., Viti et al. 2004). While these processes may impact abundances under certain conditions, we expect them to have a minimal impact on our results. For example, turbulent intermittency in the form of strong, thin shocks within the diffuse gas may enhance the emission and abundance of tracers such as CH⁺, H₂, and CO (Falgarone et al. 2005, 2009). However, excitation by FUV photons likely dominates the populations of the lowest energy states, which are what we study here. Dust grain freeze-out primarily affects gas with densities ⩾10⁴ cm⁻³; here, only 1% of the volume of our model cloud has $n_{\rm H_2} \ge$ 10⁴ cm⁻³ (see further discussion of the impact of molecular freeze-out in the Appendix). Other grain processes, such as photodesorption, may impact the species abundances over a larger range of densities. Namely, Guzmán et al. (2013) find that photodesorption is needed to reproduce the measured enhanced abundance of H₂CO in the Horsehead PDR region. In contrast, they find that grain-chemistry is not required to model the H₂CO abundance within a UV-shielded dense core, which is instead well-fit by a pure gas-phase model. This suggests that the higher UV radiated outer regions of our cloud may require consideration of additional processes, at least with respect to H₂CO. However, we note that grain chemistry involves a high-degree of uncertainty and is sensitive to the local UV field, grain properties, and cosmic-ray flux, which make it difficult to apply conclusions from case-studies in the literature to our particular conditions. We quantitatively examine the impact of potential effects of dust–grain chemistry and shock chemistry in the Appendix.

In this study, we essentially assume a "one-way" coupling between the hydrodynamics and chemistry. Performing the chemistry by post-processing the simulations allows us to consider much larger chemical networks that would otherwise be computationally impossible when evolving the chemistry in parallel. However, the chemistry is not coupled to the hydrodynamics, so, although 3d-pdr calculates the gas temperature due to UV heating and cooling, it does not influence the gas evolution, and consequently the hydrodynamics and chemistry are not fully consistent. However, even for the warmer (≈100 K) gas, motions are dominated by turbulence rather than thermal broadening, so we do not expect large differences. For further discussion of the chemistry modelling, see Offner et al. (2013). Even though the chemistry considered here is not time-dependent, the approach exhibits good agreement with Glover & Clark (2012). While their chemistry network is considerably simpler than ours, the equilibrium timescales they find are the same order of magnitude of the times at which we evolve the chemistry to in this study, and are the same order of magnitude of the free fall timescale of the modeled cloud.

We carry out synthetic observations for the 16 difference species in Table 1 and compare our theoretical results to observations. We study these species because they are commonly used tracers. We use the Leiden Atomic and Molecular Database (LAMBDA) ⁶ for the reaction rates and cross sections for the different molecules. When performing these calculations, we use the collisional partners defined in the LAMBDA database files, mainly H₂, H, and He, assuming most of the H₂ is para-H₂.

To compute the emission, we use the radiative transfer code radmc-3d ⁷ with the Large Velocity Gradient (LVG) approach (Shetty et al. 2011), a radiative transfer method that does not assume local thermodynamic equilibrium (LTE). This method computes the molecular level populations given the density, velocity and temperature fields. We use the velocity information from the orion calculation, while 3d-pdr calculates the temperature and abundance information. Before performing the radiative transfer, we interpolate all data to a 256³ resolution. We calculate the synthetic spectra for velocities within ±3 km s⁻¹ of the line center. While 3d-pdr computes the level populations for all levels in the LAMBDA data, we only analyze the emission from ground level transitions.

We use a constant microturbulence value of 0.1 km s⁻¹ to account for unresolved turbulence. We adopt a "Doppler catching" parameter d_c = 0.025, which forces an interpolation of the velocity field between cells if there is a jump greater than 0.025 times the local linewidth. We use collisional excitation and line data from the Leiden atomic database for all 16 species (Schöier et al. 2005).

Since the main goal is to study the structure of the tracer species, we do not include the dust continuum in the emission calculation. We neglect heating and UV feedback due to embedded protostellar sources. However, since the cloud is forming low-mass stars, radiative feedback likely has a small impact on the emission. We convert the line emission into a brightness temperature T_b using the Rayleigh–Jeans approximation:

$$\begin{equation} T_b = \frac{c^2 I_{\nu }}{2\nu _i^2 k_b}, \end{equation} \tag{ 1 }$$

where I_ν is the specific intensity and ν_i is the frequency of the line transition.

In this study, we use the Spectral Correlation Function (hereafter denoted as SCF), first introduced by Rosolowsky et al. (1999), to study the spectral structure of the emission cubes. We define a position–position–velocity (PPV) cube as spectral line data consisting of two spatial dimensions and one velocity dimension. Likewise, a position–position–position (PPP) cube is data consisting of the density information in all three spatial dimensions. The SCF is sensitive to the temperature and sonic Mach numbers and weakly sensitive to the magnetic field strength (Padoan et al. 2003; Yeremi et al. 2014). Previous studies have only focused on the SCF of ¹³CO emission, in both simulated and observed molecular clouds. We calculate the SCF for all 16 species emission and density cubes to study how the SCF changes for the same cloud but different observational tracers.

There are several different functional forms for the SCF. We use the form given in Padoan et al. (2003), where the SCF is a function of length scale, l and denoted by S(l):

$$\begin{equation} S(l) = \left< 1 - \left< \sqrt{\frac{\sum _v | O({\bf r}, v) - O({\bf r + l}, v)|^2}{\sum _v |O({\bf r}, v)|^2 + \sum _v |O({\bf r+l}, v)|^2}} \right>_{{\bf r}} \right>_{|{\bf l}| = l}. \end{equation} \tag{ 2 }$$

Here r is a two dimensional position on the image plane, l is the offset vector with length l, and O(r, v) is any PPV spectral data set. By definition S(0) ≡ 1. The SCF is defined to be bounded between 0 and 1, with 0 indicating no correlation. Padoan et al. (2003) found that in driven turbulence simulations, the SCF can be analytically fit by a power law on small scales. They fit the SCF for ¹³CO (1–0) maps of several observed and simulated molecular clouds and demonstrated that the parameter correlations can be used as theoretical model tests. Likewise, we fit the SCF on small length scales using a power law:

$$\begin{equation} S(l) = S_0 l^{\alpha}, \end{equation} \tag{ 3 }$$

where S₀ is the value of the SCF at l = 1 pc. We fit the power laws in log space for length scales between 3× l_min, where l_min corresponds to either the simulation resolution, or the beam size, and l =1 pc corresponding to half the maximum length scale. This was to remove beam size effects from the SCF power-law fit. Figure 1 shows the SCF for four different tracers. The tracers follow a tight power-law behavior for small values of l, as illustrated by the black lines fits. However, Figure 1 also shows that at some length scale, the SCF function appears to flatten and become noisy.

Zoom In Zoom Out Reset image size

When comparing SCFs for different spectral maps, it is useful to define a quantitative scalar value that describes how similar two SCFs are to one another. We use the distance metric between two SCFs as defined by Yeremi et al. (2014):

$$\begin{equation} d_{{\rm SCF}} = \sqrt{\sum _l [S_1(l) - S_2(l)]^2}. \end{equation} \tag{ 4 }$$

Using this distance metric provides a quantitative measure of how similar the spectral structures are between different species, resolutions, and sight lines. Yeremi et al. (2014) showed that this metric is sensitive to global hydrodynamic parameters, although not to the magnetic field strength. Since the distance metric is a 1D statistic we cannot compute a chi-squared value, i.e., a probability measure of uncertainty. In Section 3.2 we will describe an example of using the distance metric to quantitatively comparing SCFs.

Figure 1 shows that the SCFs of species which trace dense gas transition away from power-law behavior at larger scales and shows that the power-law description holds only up to some scale l ≪ L where L is the cloud size. At larger length scales, the SCFs start to flatten for all species. For example, the CO SCF is a continuous power law up to nearly 2 pc, while the SCF for NH₃ starts to flatten around 0.4 pc. This flattening indicates where the gas is not dense enough to thermally emit. CO remains correlated at large scales throughout the cloud, with a high surface filling fraction of f_s = 0.82, and samples both diffuse and dense structures. CO has an low volume filling fraction, f_v = 0.07 since it is still contained mostly in dense environments with emission coming from diffuse regions because of its substantially smaller critical density. The emission from high-density regions is constrained by its high optical depth. In our study, the surface filling fraction, f_s, is the fraction of the area containing the brightest 95% of the emission. The volume filling fraction, f_v, is defined from the 3D abundances, being the fraction of the volume containing the top 95% of the mass. However, NH₃ becomes quickly uncorrelated, because it traces only compact emission, i.e., dense cores, illustrated its much smaller volume filling fraction, f_v ∼ 0.005, and small surface filling fraction, f_s ∼ 0.25.

Figure 2 displays the SCF slope as a function of offset, S₀ for 16 synthetic emission cubes, where S₀ is the value of the SCF at 1 pc. While we define the SCF in terms of a slope and an offset, we show in Figure 2 that these two parameters are tightly correlated. Therefore, henceforth, we only discuss the slope of the SCF. Padoan et al. (2003) also found a strong correlation between the slope and the offset for simulated clouds.

Zoom In Zoom Out Reset image size

A statistic describing cloud structure is only meaningful if it does not vary strongly as a function of cloud viewing angle. Here we verify that this is the case.

We calculate the emission for nine different lines of sights through the cloud in 45° increments about the x-axis for all 16 tracers. We find that the SCF slopes change by only ≈10% over all viewing angles. This variation, though small, is useful as a benchmark to define when two tracers are truly similar. To define an effective uncertainty, we calculate the maximum distance between any two SCFs for all lines of sight. Figure 3 shows this maximum distance for six different species. We find that there is a power-law relation between the average slope, 〈|α|〉, and the maximum distance, d_max. This relation now provides a quantitative way of identifying homologous tracers. We define two tracers as "complementary" if their SCF distance for a given line of sight is less than the maximum d_max of either tracer. Therefore, species A and B are complementary if

$$\begin{equation} d_{\rm SCF}(B \Leftrightarrow A) \le {\rm max}(d_{\rm max}(A), d_{\rm max}(B)). \end{equation} \tag{ 5 }$$

Figure 3 also suggests that diffuse species are more likely to have complements because they have shallower slopes. We can see this by comparing the SCFs of CO and NH₃. If we fix the CO offset and alter the slope by the maximum ± 10% error, the total rms difference is geometrically d_max. Since the SCF slope is shallower for CO, the magnitude of the SCF remains higher. In contrast, if the NH₃ SCF is altered in a similar manner, the magnitude of the SCF becomes smaller much faster as a result of the steeper power law. Therefore, the rms difference between the SCFs is much smaller than what it would be for CO. As such, there is an expectation that d_max should decrease with increasing slope.

Zoom In Zoom Out Reset image size

Figure 4 shows the SCF slopes for all the species in this study with their line-of-sight scatter. The species are seemingly separated into three distinct groups, though only two of the groups show a large jump in slope. Based on these groupings, we will refer to diffuse tracers as those with α ∼ −0.3, intermediate tracers as those with α ∼ −0.5, and dense tracers as those with α ∼ −0.75. CO and C both have the same SCF slope within 1 σ of each other and within 2 σ of all the other diffuse tracers except C⁺, where σ is the line-of-sight slope error. The larger errors for C⁺ seem to correspond to different diffuse structures being superimposed at different sight lines giving somewhat different structures. Note that while the species whose SCFs are within the line-of-sight error limits are homologous (meaning they trace similar density regimes), they are only complementary if they also satisfy the distance metric criteria (Equation (5)).

Zoom In Zoom Out Reset image size

The simulated molecular cloud that we study has no magnetic fields. While Yeremi et al. (2014) find that the SCF only weakly depends on the magnetic field, the presence of a strong magnetic field could create asymmetries in the gas distribution. Large asymmetries could in turn increase the SCF dependence on the viewing angle.

We use 16 common astrophysical tracers to investigate how the SCF depends on species. We generate synthetic emission maps for the lowest energy state transitions shown in Table 1. Figures 5 and 6 show integrated emission maps in position–position (PP) and position–velocity (PV) space, respectively. The various tracers span different ranges of position and velocity space depending on their abundance and excitation. For example, astronomers commonly use N₂H⁺ to trace dense gas, so it is unsurprising that N²H⁺ exhibits very compact emission in both figures. The velocity plots indicate structure that may be hidden by projection. Dense cores stand out in both maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We find that the filling fractions in PP and PV are similar: tracers with compact spatial emission are also compact in velocity space. This is essentially Larson's relation (Larson 1981), which states that small structures should have small velocity extents. CO has a high surface filling factor (f_s ∼ 0.82) and traces both high- and low-density regions, with a relatively high abundance across the entire spatial region, exhibited by its spatial emission distribution. CO is chemically connected to the high-density tracer H₂CO through a number of reactions. H₂CO is photodissociated to form CO when the density becomes low enough that it is no longer shielded from the external UV field. CO can get turned into H₂CO through intermediaries, such as HCO⁺ in both gas phase reactions and dust grain chemistry. As expected then, the H₂CO emission has a much smaller surface filling fraction of f_s ∼ 0.41 and is resides mostly in dense environments, requiring a relatively high density to be excited (n_c ∼ 10⁵ cm⁻³).

Optical depth indicates the degree of transparency. Tracers that are optically thin, such as NH₃, have emission that reflects the underlying density structure more accurately. Figure 7 shows the optical depth of each line transition at the line center calculated by radmc-3d. The figure shows that C, CO, CN, and C+ all have high optical depths. This is a result of having a lower critical density with a relatively high abundance. The intermediate density tracers, CS, SiO, HCN, HNC, and HCO+, are similar, with the gas only being optically thick towards the highest-density regions. High-density tracers, NH₃, H₂CO, H₂CS, and N₂H⁺, remain optically thin throughout almost all of the entire volume. NH₃ is an exception because of its fairly low critical density, allowing it to be excited down to lower densities, though still relatively optically thin except in the densest regions of the filaments. OH⁺ suffers from a very low abundance, and it is not easy to excite, so its emission appears very diffuse.

Zoom In Zoom Out Reset image size

Figure 8 displays two different integrated maps with spectra at two different points. It shows that the CO emission is mostly dominated by gas motions. The red star spectrum shows a single feature with a width of ∼0.7 km s⁻¹, which is consistent with the characteristic turbulence velocity. The CO white star spectrum is broad in part because it is a superposition of multiple density features along that particular line of sight. Similarly, the NH₃ white star spectrum shows a superposition of several density features at different velocities along the line of sight. However, the NH₃ red star spectrum, which is along a rather diffuse line of sight, shows a much narrower feature with a width of 0.3 km s⁻¹. The average sound speed of a cold molecular cloud is approximately 0.2 km s⁻¹, indicating that the NH₃ in the red star region is undergoing purely thermal motion rather than dominated by turbulence.

Zoom In Zoom Out Reset image size

Beam resolution is one of the most important factors in observations and it impacts the apparent gas structure and mean optical depth. A larger beam averages out the emission within the beam size, lowering the overall optical depth. We convolve the emission maps with a Gaussian beam to simulate a realistic resolution observation. We place the simulated cloud at a distance of 1 kpc to establish the angular size. Larger beams significantly blend the dense cores and diffuse gas structure. Interestingly, we find that some species' spectral structures are artificially similar at one resolution and then deviate significantly at some lower resolution. To quantify this, at each resolution we calculate the SCF for each species and then compare their distance metrics. Figure 9 shows the evolution of the SCF slopes for all 16 species as a function of spatial resolution. The overlapping regions indicate where tracers at a particular spatial resolution have similar emission distributions. Note that this only impacts the emission. Therefore, overlapping regions do not indicate similar density distributions, unless the SCF slopes are similar at good (near simulation) resolutions. Thus, Figure 9 quantifies which tracers are statistically similar and useful for studying particular densities and size scales.

Zoom In Zoom Out Reset image size

At the highest resolution (i.e., no smoothing), C, CO, CN, and C⁺ are all very similar tracers since they all trace diffuse gas. The positive correlation between C and ¹²CO has been observationally studied for years (e.g., Plume et al. 1999; Ikeda et al. 1999, 2002; Kulesa et al. 2005; Shimajiri et al. 2013), although historically C was theoretically predicted to only exist in a PDR surface layer (Hollenbach & Tielens 1999). Papadopoulos et al. (2004) predicted that C should be more prevalent in clouds than previously predicted by 1D models due to a combination of nonequilibrium processes and clumpy cloud morphology. Glover et al. (2014) and Offner et al. (2014) both qualitatively demonstrated the similarity between the C and CO distributions in 3D PDR calculations. Our SCF comparison demonstrates that CO and C spectral cubes are quantitatively similar.

At higher densities, species such as SO and NH₃ also appear very similar. As expected, HCN and HNC structures appear nearly identical in spectral space. At lower spatial resolution, several species intersect: NH₃, N₂H⁺, SO, and H₂CO. This occurs as larger beams increasingly blend dense, compact emission.

We note that some species have a non-monotonic dependence on resolution where the SCF slope decreases until ∼15–20'' and then increases again. The decrease in the slope is due to overlapping Gaussian structures creating artificial cores. Figure 10 shows the emission spectral structure evolution with beam size visually for NH₃ emission in both PP and PV space. Figures 11 and 12 illustrate these trends as a function of beam size using the distance metric. The change in the slope occurs because species with significant compact emission (those which are generally optically thin) have their emission smoothed on larger scales, making their emission more extended and thus appear similar to other more optically thick species. Smoothing also somewhat affects the more diffuse tracers C, CO, CN, and C+, which experience blending on larger scales. This can give the appearance of a more core-like structure where these smoothed regions overlap. At lower density, the emission smooths into am even more diffuse looking component.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Section 3.3 we demonstrate that many tracers produce similar SCFs. However, similarity between the SCFs of two emission maps does not guarantee that the underlying densities are also similar. In this section, we compare the SCFs of the emission and density data directly. We obtain a PPV cube based on the gas density (hereafter denoted by PPV_ρ) by taking the simulated density and velocity cubes and constructing a PPV cube where each spectral bin contains the number density (calculated as N_i = N_H · n_i) within a given velocity bin. Here, "V" is the velocity vector projected along a given line of sight. For the PPV_ρ cube, the (i, j, k)^th voxel contains the total mass along some line of sight through the point (x_i, y_i) within a particular line-of-sight velocity range Δv_k. We compare the emission and density SCFs in Figure 12. This figure shows the distance between the emission SCFs and PPV_ρ SCFs for each species. For C, CO, CN, and C+, PPV_ρ does not match the emission structure well. In fact, as the spatial resolution of the PPV cube becomes coarser, agreement between density and emission worsens. High opacity tracers exhibit poor correspondence between emission and density.

In contrast, species with optically thin emission have very similar density and emission SCFs, as indicated by the darker cells in Figure 12. For example, CO has a low critical density, n_c ∼ 2000 cm⁻³, and a relatively high abundance. It becomes optically thick as the gas density approaches 10⁴ cm⁻³. Figure 7 demonstrates that CO is very optically thick throughout most of the simulation box. High optical depth effectively flattens the perceived distribution, i.e., as the gas becomes optically thick the emission no longer traces higher-density regions. On the other hand, N₂H⁺ is optically thin throughout the domain and has a small distance (d < 0.1) between its emission and density SCF. N₂H⁺ has a high critical density of 10⁵ cm⁻³ and has a much lower abundance than CO.

Differences between the true density and the emission arises from a combination of chemistry, which changes the abundance, and excitation, which impacts the line shape. For example, OH⁺ has a low critical density of ∼4000 cm⁻³ but its emission remains optically thin due to its very low abundance. However, HCN is optically thin because it is only excited at gas densities above 10⁶ cm⁻³ despite its modest abundance: [HCN]/[H₂]  ∼10⁻⁷.

We now compare our results to observations of two local molecular clouds (MCs), Perseus and Ophiuchus. We use the ¹²CO data from the COMPLETE survey (Ridge et al. 2006), which were taken using the Five College Radio Astronomy Observatory (FCRAO). The ¹²CO observations have an angular resolution of 46''. Since the Perseus and Ophiuchus clouds are approximately 250 pc and 150 pc away (Ridge et al. 2006) corresponding to physical scales of 0.06 pc and 0.03 pc, respectively, we can compare at the two length scales represented in our simulation. We calculate the ¹²CO SCF for the full extent of both clouds. We also divide Perseus into two parts and compute the SCF of each half. Table 2 displays the resulting SCF slopes and power-law fit errors. We plot the full SCFs in Figure 13.

Zoom In Zoom Out Reset image size

In order to properly compare the model to the observations, we perform several procedures on the observational data. First, we smooth the data with a Gaussian beam corresponding to the FCRAO 46'' beam resolution. Then, we regrid the data so that each pixel has the same spatial pixel scale of the two different clouds. Third, we apply a detection limit where we remove pixels with emission less than 0.01 × 〈T_b〉. Finally, we add noise corresponding to Gaussian thermal noise with a standard deviation of σ_T ∼ 0.3 K to model the noise in the Perseus COMPLETE data. We did not perform biased sampling, which may affect the longest length scales of the SCF due to the irregular observational stencil.

Table 2 gives the values of the slopes and the associated error of the power-law fit. The model cloud SCF initially has a slope of −0.29. However, after matching to the observational resolution the SCF slope changes by quite a bit to −0.13 and −0.19 for the Ophiuchus and Perseus resolutions, respectively. Furthermore, the fit error indicated in the table shows that the power-law fits are all well constrained. Figure 13 shows that the model cloud and the Perseus cloud have similar SCFs when taking into account observational biases, such as beam size and pixel resolution. This is indicated by the similar power-law fits shown in the bottom of the figure. The Ophiuchus SCF is steeper than the mode SCF at the Ophiuchus pixel and beam resolution, though it is similar to both the Perseus SCF and the model SCF at Perseus resolution. The flattening of the SCF at small scales is due to the beam resolution. While the emission is still correlated, the turbulence has been resolved out within the beam area. There is still some difference between the simulation and observations which may be due to other biases. However the model and observations match within the 10% error mentioned in Section 3.2.

In Figure 1, the SCFs for all of the tracers appear to flatten out at some large scale, with the exact scale seemingly changing depending on the tracer. Padoan et al. (2003) also found that their SCFs from simulated clouds traced by ¹³CO flatten, but the SCFs for local clouds, also traced by ¹³CO, steepen. We were not able to replicate this discrepancy with our simulations, though we conclude that the steepening effect is most likely due to some observational bias at large length scales.

A fundamental question of this work is: what does the SCF reveal about the chemistry of various tracers? Using Figures 9 and 4 we group the SCF slopes into three emission categories: diffuse, intermediate, and compact. C is an example diffuse tracer which is excited down to lower densities of ∼800 cm⁻³ and is abundant throughout the entire volume. HCN is excited at fairly high densities and it becomes quite abundant once it is shielded from the UV field, so it traces up to intermediate scales. A common dense core tracer is N₂H⁺ which also gets excited at higher densities, where it is also shielded from photodissociation. Commonly studied tracers such as CO, C⁺, HCN, NH₃, and N₂H⁺ confirm that the SCF slope reflects the expected optical depth, filling fraction, and critical density of the emission. While it is impossible to compute the true density distribution from the SCF alone, the slope does indicate the scale the species traces. For instance, shallower slopes, such as that of CO, show that the gas remains correlated on larger scales.

CO is a species of particular astrophysical importance. Over the past several decades, CO has become the most prominently utilized cloud mass tracer and has received significant theoretical attention. Our results show that the SCF of CO traces diffuse regions, as indicated by its shallow slope. Our results also show that CO traces the gas up to high densities, although proper treatment of dust grain depletion could change this. Since CO has a low critical density (see Table 1), the lowest level transition is easily excited throughout most of the cloud. This allows the cloud to produce emission in both the lower-density environments and the dense cores (where the emission saturates due to optical depth). For this reason, CO also has a high volume and surface filling fraction. Figure 7 confirms that CO has the highest average optical depth, which is due to its low critical density and higher abundance.

However, the SCF analysis identifies several tracers that exhibit similar emission. We find that C and CN are homologous to CO. The correlations between these species makes chemical sense since all three depend upon the abundance of neutral Carbon. The UMIST chemical network shows that both CO and CN form rapidly in the original diffuse environment. Their behavior at higher densities though is quite different. CN is photodissociated much more rapidly towards lower densities than CO because it cannot self shield. CN is depleted through reactions forming more complex molecules faster than CO. This results in a factor of several orders of magnitude difference between the abundances of the two molecules. Therefore, it is reasonable to expect CO and CN to be very similar in low-density environments. In regions of high mass star formation, CN could be used as a proxy for CO, which is optically thick, or of the surrounding diffuse gas.

Another common molecular gas tracer is HCN, which is frequently used to trace gas with densities above 10⁶ cm⁻³. The SCF slope implies it also traces intermediate density regions in filaments, i.e., n ∼ 10³–10⁴ cm⁻³, although the emission level may be too low to detect. Recent work by Forbrich et al. (2014) confirms that HCN and HNC are good tracers of dense environments over a wide range of extinctions. The optical depths we calculate for HCN are off by a factor of three since we do not model the HCN fine structure for the ground state transition.

Finally, there are a variety of species that trace dense cores, including NH₃ and N₂H⁺. Figure 6 confirms the trend that species tracing more diffuse regions have shallower slopes than species that typically trace high-density regions. Here, the regions the species trace are clear in the velocity information. Diffuse gas tracers have significant emission in a broad range of velocities, as illustrated by the horizontal bands in species CO, C, and CN in Figure 6. Species which trace higher densities show no diffuse component. Instead, the emission is located in clumps or filaments, which exhibit a smaller range of velocities.

While the density is important, the gas must also have a high enough temperature to excite the transition. We define a characteristic line temperature T, where hν_i = kT_i, where ν_i is the line frequency. All the species have line temperatures below 10 K except C, C⁺, and OH⁺, which have line temperatures of 23 K, 91 K, and 44 K, respectively. The average gas temperature in diffuse regions is around 100 K, so even these species are easily excited from their ground states. However, OH⁺ is mostly observed in absorption. This is due to its low abundance and higher excitation making measurements from absorption in the dust continuum easier than trying to detect the very small emission signal (as noted by the multiplicative factor of ∼10⁶ in Figures 5 and 6).

Complex molecules, such as NH₃ require higher densities to form, and photodissociate rapidly at lower densities where the UV field is higher. Simple light diatomic molecules such as CN form in diffuse regions. An exception is CS, which appears only in intermediate density regions. It photodissociates faster than CO but slower than CN. However, due to the much lower initial abundance of sulfur, it forms only where the sulfur is concentrated. Most of the tracers in this study that have a shallow slope also tend to have a very high optical depth. The only exception is OH⁺ which has a very low abundance.

The results from Figures 9 and 11, as well as our definition of complementary species, suggest sets of homologous species. These are groups of species whose spectral structure is very similar in PPV space, indicating that, especially for the optically thin tracers, they should trace similar density regimes.

Recent theoretical studies, such as Glover et al. (2014) and Offner et al. (2014), found that C is a good alternative tracer to CO. C has several advantages, including a ground state transition at 609 μm and a lower optical depth than CO. These recent studies challenge the idea that C traces only the surface of the PDR. Our study confirms this picture by showing that C and CO have very similar SCFs, with CO being complementary to C. Our study also predicts that CN is an alternative tracer to CO and C. CN (1–0) has a similar transition frequency with an optical depth around an order of magnitude less than CO, although its slope is slightly flatter indicating that it traces even lower-density regions, and it has a similarly high filling fraction. The slightly flatter slope is expected; CN is destroyed faster than CO in higher-density environments due to its role as a reactant in reactions forming more complex molecules. Since CN and CO form very quickly and depend on the relatively high-abundance of C, they both have very high filling fractions. Their surface filling fractions are both over 90%, and their volume filling fractions are greater than 0.3. In fact, Table 1 shows that CN is both more surface filling and volume filling than CO. The cosmic-ray induced photodissociation rates (i.e., CX + CRPHOT → C + X) given by the UMIST2012 network are R_CN  =  1.4  ×  10⁻¹³ s⁻¹ and R_CO = 7.5 × 10⁻¹⁶ s⁻¹. Cosmic rays penetrate further into the cloud than the external UV radiation, indicating the CN will be destroyed by cosmic rays inside the cloud faster than CO. Indeed, CO has an abundance several orders of magnitude greater. Singly ionized carbon traces more diffuse regions but is not as abundant in higher-density regions where it combines with O to form CO. The C⁺ (1–0) transition is in the infrared, which can be observed using space-based instruments, such as the Herschel Space Observatory (e.g., Langer et al. 2014) and the GREAT Spectrometer on Stratospheric Observatory for Infrared Astronomy (SOFIA).

A commonly used high-density tracer is N₂H⁺, which has a ground state transition at 93 GHz. Our study predicts several homologous tracers to N₂H⁺ including H₂CO, H₂CS, NH₃, and SO. All of these species exist in similar environments and have surface filling fractions between 0.01 and 0.25. However, some of these tracers only form the high-density cores which such as N₂H⁺, with a volume filling fraction of 0.21. H₂CO has a higher volume filling fraction than N₂H⁺ (f_v = 0.37) which indicates that it traces a larger fraction of the core gas. We find that H₂CO and SO have higher brightness temperatures indicating that they should be easier to detect than species with fainter emission such as H₂CS and N₂H⁺. NH₃ and SO show significant emission in filaments, but N₂H⁺ is brightest in dense "cores." NH₃ has a relatively low critical density (n_c = 1991 cm⁻³), so it is excited at lower densities than high density tracers like N₂H⁺. However, it is not as bright as other low critical density tracers like CO (n_c ∼ 2000 cm⁻³) because of its low abundance outside of dense regions. Some of these correlations could change with the inclusion of dust grain chemistry. This would lower the abundances of the higher-density tracers, such as NH₃, H₂CO, H₂CS, and SO, and reduce their emission. However, these molecules only begin to deplete for H₂ number densities ⩾10⁷, greater than the maximum density in this simulation. Detailed treatment of gas–grain chemistry is beyond the scope of this work.

Both CO and HCN, as well as other high-density tracers, may be affected by depletion at high densities. CO depletion is extreme at high-densities with the mean abundances declining by several orders of magnitude (see Bacmann et al. 2002; Bergin & Tafalla 2007; Hocuk et al. 2014). HCN also freezes out, but at higher densities. In order to assess the impact of depletion on our CO and HCN results, we adopt crude treatments for dust freeze-out and recalculate the optical depth and SCF. We set the abundance of CO to 0 where the H₂ number density exceeds 10⁴ cm⁻³. For HCN we set a limit of [HCN]/[H₂] = 10⁻⁸ in cells where $n_{\rm H_2} \ge 10^4$ cm⁻³. We compare the SCF slopes of the tracers with and without depletion and find that the change in slope is less than 10%. Setting a maximum HCN abundance of 10⁻⁸ relative to H₂ reduced the HCN optical depths to τ < 10 across the entire cube, with the typical optical depth being ∼1–4. Furthermore, since in these high densities the HCN gas will remain optically thick, it should not drastically impact the SCF correlations. Overall, we conclude that the exclusion of dust depletion does not significantly affect our results due to the small percentage of the simulation volume at high densities (1%).

When CO starts to freeze-out onto dust grains, several molecules, such as N₂H⁺, will see an increase in their abundances. This occurs because H$_3^+$ reacts with CO to form HCO⁺ which is the main destruction mechanism for H$_3^+$. As the amount of CO decreases there is a surplus of H$_3^+$. In high-density environments, H$_3^+$ can be destroyed to form N₂H⁺ (see Aikawa et al. 2001; Bergin et al. 2002; Rodgers & Charnley 2004) by the following mechanism:

$$\begin{equation} {\rm H}_3^+ + {\rm N}_2 \rightarrow {\rm N}_2{\rm H}^+ + {\rm H}_2. \end{equation} \tag{ A1 }$$

In order to test whether this effect has an impact, we adopt a simple approximation in which the amount of N₂H⁺ increases at high densities. Since CO starts to deplete at H₂ densities above n(H₂) ⩾ 10⁴ cm⁻³, we multiply the N₂H⁺ abundance by a factor of 4 where the gas fits this criteria. Since the emission is mostly optically thin, this merely scales the emission by a constant factor. Similar to the case for CO and HCN depletion, we find that this mechanism has no statistical impact on our results. It is noteworthy that the enhancement of N₂H⁺ is not universal, with studies such as Tafalla et al. (2002, 2004) showing either no N₂H⁺ increase or showing N₂H⁺ depletion.

Since our hydrodynamic simulation does not resolve shock fronts to the necessary resolution to calculate post-shock densities and temperatures, well-known shock tracers such as SiO will not be properly modelled. In shocks, changes in the density and temperature can lead to a large enhancement of the abundance of SiO. Furthermore, at higher temperatures, the abundance of SiO can be increased by sputtering from dust grains. This mechanism is the ejection of Si and SiO from a grain surface following a high enough energy impact of gaseous species. At lower temperatures, the sputtering rate is expected to be small (Draine & Salpeter 1979). Even though sputtering from shock chemistry, and other shock effects, will affect the abundance, the inclusion of shock chemistry is beyond the scope of this work. Finally, the increased brightness of the SiO emission is an artifact of using a nondepleted Si abundance, which is over an order of magnitude larger that observed depleted abundances (Morton 1975).