HOW DENSE IS YOUR GAS? ON THE RECOVERABILITY OF LVG MODEL PARAMETERS

For all of our tests we generate artificial lines using RADEX and collision rates from the Leiden Lambda database (van der Tak et al. 2007; Dumouchel et al. 2010; Yang et al. 2010). The specifics depend upon the investigation in question, but the features common to most of the models are as follows.

• 1.  
  The physical free parameters are drawn from a uniform distribution in log space.
• 2.  
  We restrict the thermal pressure, $P/{k}_{{\rm{B}}}$, of the models to be <1 × 10⁸ K cm⁻³. Higher thermal pressures are very unlikely/rare (even ULIRGs only exhibit thermal pressures ∼10⁷ K cm⁻³; Dopita et al. 2005), but also RADEX tends to fail at high thermal pressures, biasing the results.⁴
• 3.  
  The virial parameter, K_vir,⁵ is set to unity (corresponding to virialized clouds), so that for a given ${n}_{{{\rm{H}}}_{2}}$ dv/dr is determined, and not a free parameter. It is, however, kept as a free parameter when attempting to recover input parameters with the MCMC models.

We introduce errors on the lines produced by RADEX. The motivation is to emulate reality, where there is a true, unknown value of the line brightness temperature, but the observed value may be larger or smaller, primarily due to calibration uncertainties. In these tests, we know the true value and generate the randomized values.

For the tests including observation calibration uncertainty, we randomize the flux of each line individually, such that the new, uncertain flux is drawn from a normal distribution with mean equal to the artificial line brightness temperature and with a standard deviation equal to 0.1 times the artificial line brightness temperature. We calculate line ratios using randomized line brightnesses, propagating the uncertainties so that the uncertainty in the line ratios is ∼14%. This is the uncertainty used when calculating χ² in the grids and MCMC traces.

In all but one of the cases we work entirely in line brightness temperature, the exception being when we plot SLED fluxes. This is a further idealization, as it does not take into account the additional errors when observing galaxies due to variations in molecular excitation which can lead to source sizes varying as a function of line transition.

We explore two methods for recovering parameters from line ratios using RADEX: grids and MCMCs. We first adopt the traditional grid-based method before demonstrating the capabilities and limitations of an MCMC approach.

The grid method is, in principle, uniquely simple: a grid of input parameters is generated, and RADEX fluxes and line ratios are generated for each point. Lines or line ratios are then fitted (minimum χ²) to the grid, and the best fit parameters thus recovered. Some authors extend this by generating a full likelihood grid (Ward et al. 2003; Kamenetzky et al. 2011; Rangwala et al. 2011; Papadopoulos et al. 2014). The grid must be trimmed of spurious points where either RADEX or the LVG framework breaks down (due to extremely high optical depths). This is accomplished with appropriate checks of the RADEX output when generating the grid.

We also use an MCMC model to recover the parameters, as in Tunnard et al. (2015b) and similar to Kamenetzky et al. (2014), who adopted a nested sampling approach for two phase modeling of CO emission. Instead of fitting to a discrete grid, the MCMC explores the continuous parameter space, sampling the posterior probability distribution. This also allows for the modeling of much more complex parameter spaces than is easily achievable with the grid method. Furthermore, it has the advantage of naturally producing posterior probabilities and being trivial to marginalize.

The results of the MCMC are a chain of parameter values (called a trace): for a converged MCMC using the Metropolis–Hastings sampler, the results should oscillate about the most probable values (assuming a singly peaked posterior), so that the mean of each parameter over the trace corresponds to the most likely parameter (e.g., Hanson 2001). We present both the means and the parameters corresponding to the minimum ${\chi }^{2}$ set of parameters.

The MCMC method is slower than grids for low dimensional parameter spaces, presenting a significant overhead when running for many realizations of the randomized line ratios. Its true power lies in the ability to simultaneously fit multiple species and gas phases, avoiding the bias inherent in fitting multiple gas phases sequentially or by eye (Kamenetzky et al. 2014). Here, we run all models for 5 × 10³ steps with a burn in of 1 × 10³ steps, based off preliminary tests on convergence rates. This is less than ideal, but necessary to be able to run a large number of trials in a reasonable time.

For the grid recovery models, we generate 10⁵ trials, each with physical parameters drawn randomly from the parameter space and recorded, subject to the pressure upper limit. For the MCMC, only 1000 trials are generated due to the much greater run time of the MCMC. Since we are not directly comparing the grid and MCMC methods, this is not a significant concern.

We present four simple models for testing the precision and accuracy of the LVG approach.

• 1.  
  A general randomized model, which randomly selects parameters from uniform distributions in log space and uses a single molecular species.
• 2.  
  As above, but using a molecular species and its ¹³C isotopologue lines. We explore how well we can recover the isotopologue abundance ratio, and the variation when using a single ¹³C line and when using four ¹³C lines. This model is only run with the MCMC.
• 3.  
  A study of the effects of assuming an erroneous molecular abundance, relative to adopting the molecular abundance as a free parameter.
• 4.  
  A study of the effects of assuming an erroneous isotopologue abundance ratio when including a single ¹³C isotopologue line with four main isotopologue lines, relative to adopting the isotopologue abundance ratio as a free parameter, and two investigations of the effects of background radiation.
• 5.  
  An exploration of the effect of the background radiation blackbody temperature, T_bg, on line emission of HCN and CO.
• 6.  
  A toy model of the effects of varying CMB temperatures on galaxy-integrated CO and HCN emission.

These models represent a wide range of conditions and provide essential estimates of the optimal accuracy of LVG models. All of these model are run with and without the 10% randomization errors to distinguish between inherent and random noise-induced offsets and uncertainties.

For the LVG testing, we fix T_bg = 3 K. For all comparison grids and MCMCs, we adopt the parameter ranges $3\leqslant {T}_{{\rm{k}}}\leqslant 1000$ and ${10}^{2}\leqslant {n}_{{{\rm{H}}}_{2}}\leqslant {10}^{8}$. For the grids K_vir is set to unity, while for the MCMC runs we also adopt $0.1\leqslant {dv}/{dr}\leqslant 1000$, $0.5\leqslant {K}_{{\rm{vir}}}\leqslant 2.0$, and for the models including ¹³C isotopologues $10\leqslant {X}_{{}^{12}{\rm{C}}}/{X}_{{}^{13}{\rm{C}}}\leqslant 1000$.

The MCMC models require initial positions to either be provided or randomly generated. In this work, we generate random initial conditions for each MCMC, with each parameter, θ_i, drawn from from a uniform random distribution centered on the middle of the prior range and with a width 0.3 times the parameter prior range (in log space).

The general model explores randomized T_k and ${n}_{{{\rm{H}}}_{2}}$, with fixed molecular abundances⁶ , T_bg = 3 K and K_vir = 1. The model was run separately for HCN and CO, with and without the 10% randomization errors on the lines fluxes, and the input parameters were recovered using both grid and MCMC methods.

This model serves to explore the ability of LVG models to recover the true physical parameters under the simplest and most ideal conditions.

For HCN, we use only the J = 1–0 to J = 4–3 transitions, as these are the most commonly available lines (although all four lines are available in the literature for only a few sources). For CO, we use the full SLED up to the J = 13–12 transition, as is available for many galaxies but in particular those in the Herschel Comprehensive ULIRG Emission Survey (HerCULES, PI: van der Werf, van der Werf et al. 2010; Kamenetzky et al. 2015; Rosenberg et al. 2015). However, in reality, this full SLED covers a wide range of physical conditions: here they all sample a single-parameter set. Therefore, the derived uncertainties are a lower limit on the true uncertainties when interpreting real observations. Indeed, fitting a single gas phase to the full CO SLED of a galaxy is almost never possible, and even less frequently reliable. We leave an investigation into the effects of reducing a continuum of gas phases across a galaxy to a two- or three-phase model to a future paper.

This model is identical to the Randomized, Single Species case, except that we now only use HCN and include H¹³CN lines. We only run these models with the MCMC as the grid becomes a cube and starts to take a prohibitively long time to generate and compare with the synthetic line ratios. The [HCN]/[H¹³CN] abundance ratio is entered into the model as a free parameter, while we keep ${X}_{{\rm{HCN}}}=2\times {10}^{-8}$ for both line generation and parameter recovery.

We run these tests using both the complete H¹³CN SLED from J = 1–0 up to J = 4–3 as well as just the J = 1–0 line. It is rare enough to have a single H¹³CN observation in a galaxy: to have four lines is unheard of. Nevertheless, we explore the four line case to assess the potential benefit of obtaining additional lines.

It is standard practice with LVG modeling to assume a canonical molecular abundance. As well as simply being frequently necessitated by a paucity of observed lines, this has a physical justification in that LVG models are only sensitive to the combination ${X}_{{\rm{mol}}}{({dv}/{dr})}^{-1}$, not to either individually. This is sometimes then used in the RADEX native form N_mol/Δv. However, since

$$\begin{eqnarray}&&{N}_{{\rm{mol}}}={X}_{{\rm{mol}}}{N}_{{{\rm{H}}}_{2}}=\displaystyle \frac{3.09\times {10}^{18}\;{X}_{{\rm{mol}}}\;{n}_{{{\rm{H}}}_{2}}\;{\rm{\Delta }}v}{{dv}/{dr}},\end{eqnarray} \tag{ 2 }$$

fixing ${N}_{{\rm{mol}}}/{\rm{\Delta }}v$ actually forces certain tacit assumptions about the virial state of the gas.⁷ Due to these complications, we present a simple exploration of whether model parameters are better recovered assuming a canonical, but potentially erroneous, X_mol, or whether better results are obtained by leaving X_mol as a free parameter.

Note that we do not claim that the molecular abundance can be recovered accurately when only a single species is observed. What we are investigating is whether the recovery of T_k and ${n}_{{{\rm{H}}}_{2}}$ is seriously affected by assuming a canonical X_mol.

To this end, we generate HCN lines with physical parameters drawn from the aforementioned parameter ranges, but we now also draw X_HCN from the range 10⁻¹²–10⁻⁴. However, when recovering the parameters with the MCMC, we run two cases: the first where X_HCN is introduced as a free parameter, and the second where we fix X_HCN to the canonical value of 2 × 10⁻⁸.

When only a single ¹³C isotopologue line is available, it is common to assume a canonical molecular abundance and/or isotopologue abundance ratio. However, these values are poorly constrained for CO, let alone less abundant molecular species such as HCN. Furthermore, isotopologue selective chemical effects such as selective photodissociation and isotope charge exchange reactions can lead to isotope fractionation. When coupled with real variations in the elemental isotope abundances, such as in accreted primordial gas or cosmic-ray-dominated star formation paradigms, can lead to large variations in the isotopologue abundance ratios from species to species and galaxy to galaxy (Henkel et al. 2010, 2014; Papadopoulos 2010; Papadopoulos et al. 2011; Ritchey et al. 2011; Roueff et al. 2015). We therefore explore the effects of assuming an incorrect abundance ratio on the recovered T_k and ${n}_{{{\rm{H}}}_{2}}$, using HCN.

The [HCN]/[H¹³CN] ratio is drawn from the range 10 to 1000 (as a uniform random variable in log space), and we recover the parameters assuming [HCN]/[H¹³CN] = 60 and adopting a fixed X_HCN = 2 × 10⁻⁸.⁸ This allows us to test whether assuming a canonical value is reasonable, or whether it is in fact better to introduce [HCN]/[H¹³CN] as a free parameter (as in the Randomized, with Isotopologues model), even when there is only a single H¹³CN line available.

The exploration of the effects of T_bg on the lines is not a test of the abilities of LVG models to recover input parameters. Rather, it is a demonstration of the effects of background radiation on the molecular lines. We generate line fluxes for T_bg = 3, 5, 10, 15, 20, 25, and 30 K. For each T_bg, we generate SLEDs from J = 1–0 to J = 4–3 for HCN and from J = 1–0 to J = 13–12 for CO, with constant density ${n}_{{{\rm{H}}}_{2}}$ = 1 × 10⁴ cm⁻³ (HCN) and 1 × 10³ cm⁻³ (CO), and X_HCN = 2 × 10⁻⁸ and X_CO = 5 × 10⁻⁵, for T_k = 3, 5, 10, 50, 100, 300, 500, and 1000 K. For the CO lines, we reduce the density to 1 × 10³ cm⁻³ since at 1 × 10⁴ cm⁻³ many lines, but in particular the J = 5–4 line, present pathological fluxes (e.g., for T_bg = 3 K and ${T}_{{\rm{k}}}\geqslant 500$ K RADEX fails to converge, ultimately because the column density and associated optical depth becomes too large). This is also more representative of the likely conditions of the CO on galactic scales. We also generate grids to explore the effect of T_bg on the line brightness ratios, with high resolution in T_bg and T_k.

We are also interested in how the CO and HCN line peak brightness temperatures, T_B,peak, are affected by the CMB as a function of the kinetic temperature of the emitting gas. This allows us to investigate the detectability of cold gas reservoirs at high-z. For five redshifts (0, 2, 4, 6, and 8) we run RADEX for the physical conditions above for 1000 values of T_k distributed evenly in log space from 1000 to 3 K. Other physical parameters are kept as for the models above.

We extend the functional models above to a very simple one-dimensional (1D) toy galaxy model in order to explore the effect of CMB temperature on potential galaxy-integrated SLEDs.

The model uses a 1/r temperature profile falling from 100 K at r = 0 down to 10 K at the galaxy edge, and we sample the profile with 10⁴ molecular clouds following an exponential distribution, and then integrate the total emission from all of the clouds in the galaxy. We do not explicitly include dust, but we do set T_CMB as a lower limit on T_k. It is, however, still dust-free, as RADEX does not include dust within the gas: dust can only be added as a source of background radiation incident on the cloud. This toy model allows us to demonstrate the significant effects of the CMB on CO and HCN emission and on their line ratios at high-z.

These toy models are much less complex than the work of da Cunha et al. (2013) and Narayanan & Krumholz (2014), with an emphasis lying somewhere between the two. While da Cunha et al. (2013) focussed heavily on the theory of line excitation as a function of redshift, they did not investigate galaxy-wide trends. On the other hand, Narayanan & Krumholz (2014) focussed almost entirely on detailed modeling of galaxy emission, with little investigation of the effect of the CMB at high redshift (although see their Figure 12). Our approach here is to simply demonstrate the possibility of unifying the two approaches.

The results of the Randomized, Single Species models are presented in Figure 1 in the text and Figures 10 and 11 in the Appendix, where we plot the histograms of the logarithm of the ratio of the recovered over true parameters (hereafter the "deviation"). Numerical results are included in Table 1. For the grid models we also plot the deviation of the product ${T}_{{\rm{k}}}\times {n}_{{{\rm{H}}}_{2}}$, which is more tightly constrained than T_k or ${n}_{{{\rm{H}}}_{2}}$ alone. We record the means and standard deviations of the deviations, as well as the bootstrap-estimated mean and standard deviation of the median of the deviations.

Zoom In Zoom Out Reset image size

Table 1.  Deviations of Recovered Parameters from Grid Models

	Statistical			Bootstrap Median
	${T}_{{\rm{k}}}$	${n}_{{{\rm{H}}}_{2}}$	${T}_{{\rm{k}}}\times {n}_{{{\rm{H}}}_{2}}$	Tk	${n}_{{{\rm{H}}}_{2}}$	${T}_{{\rm{k}}}\times {n}_{{{\rm{H}}}_{2}}$
HCN	−0.2 ± 0.5	0.3 ± 0.8	0.1 ± 0.5	−0.085 ± 0.001	0.194 ± 0.002	0.058 ± 0.001
HCN-noE	−0.02 ± 0.20	0.05 ± 0.37	0.03 ± 0.29	0.0001 ± 0.0003	0.0007 ± 0.0003	0.0038 ± 0.0003
CO	−0.2 ± 0.5	0.2 ± 0.9	0.0 ± 0.9	−0.058 ± 0.001	0.071 ± 0.001	0.024 ± 0.001
CO-noE	−0.1 ± 0.4	0.0 ± 0.5	−0.1 ± 0.7	−0.0120 ± 0.0003	0.0081 ± 0.0004	0.0013 ± 0.0003

Download table as:  ASCIITypeset image

The numerical results for the grid models are shown in Table 1. We include the statistical mean and standard deviations of the deviations, as well as the bootstrap-estimated means and standard deviations of the medians of the deviations. Reassuringly, all parameters are recovered without bias for the no-error cases. There is a trend for both HCN and CO, when line uncertainties are added, to underestimate T_k and overestimate ${n}_{{{\rm{H}}}_{2}}$. This trend is far weaker in the medians, suggesting that it is due to the rare but extreme outliers. There is no sign that the CO model, with nine more lines than the HCN model, obtained improved recovery rates above those of HCN.

The thermal pressure, ${T}_{{\rm{k}}}\times {n}_{{{\rm{H}}}_{2}}$, is better determined than either T_k or ${n}_{{{\rm{H}}}_{2}}$ alone. The effects of increasing T_k and increasing ${n}_{{{\rm{H}}}_{2}}$ are not identical, but they are partially degenerate, leading to the improved recovery of thermal pressure.

The results show that the recovery rates are worryingly low. Without any errors, there is a spread of 0.05–0.10 dex in all of the recovered parameters. That is, even if we knew the line ratios perfectly, in a spatially resolved, virialized gas cloud of uniform density, temperature, and velocity gradient, with a uniform and known molecular abundance, we should not expect to recover parameters with a precision better than ∼10%. Keeping all other idealizations and introducing 10% observational uncertainties, this falls to between 50% and 300%. These are lower limits on the true uncertainties. The accuracies are very good in the no-error cases, but for the randomized cases the systematic offsets can be as large as a factor of 2.

The MCMC results are recorded in Table 2. For the MCMC, we have results for both the trace mean parameters and the minimum χ² parameters. We also record the median of the standard deviations of the traces for each parameter, as these are one measure of the uncertainty in the MCMC mean results. The systematics seen in the grids are also present here, and once again it is readily apparent that they are due to outliers.

Table 2.  Deviations of Recovered Parameters from MCMC Models

	Statistical			Bootstrap Median
	Tk	${n}_{{{\rm{H}}}_{2}}$	dv/dr	Tk	${n}_{{{\rm{H}}}_{2}}$	dv/dr
HCN Best fit	−0.2 ± 0.6	0.2 ± 0.9	0.1 ± 0.6	−0.013 ± 0.003	0.07 ± 0.02	0.07 ± 0.02
HCN Mean	−0.2 ± 0.6	0.2 ± 1.0	0.1 ± 0.5	−0.015 ± 0.003	0.04 ± 0.01	0.03 ± 0.01
HCN-noE Best fit	−0.3 ± 0.7	0.2 ± 1.0	0.1 ± 0.6	−0.0010 ± 0.0003	0.004 ± 0.002	0.009 ± 0.003
HCN-noE Mean	−0.3 ± 0.7	0.2 ± 1.1	0.1 ± 0.6	−0.004 ± 0.001	0.02 ± 0.01	0.031 ± 0.009
CO Best fit	−0.5 ± 0.9	0.5 ± 1.8	0.3 ± 0.9	−0.010 ± 0.003	0.05 ± 0.02	0.11 ± 0.03
CO Mean	−0.5 ± 0.9	0.5 ± 1.8	0.3 ± 0.9	−0.011 ± 0.003	0.08 ± 0.01	0.11 ± 0.02
CO-noE Best fit	−0.6 ± 1.0	1.0 ± 2.0	0.6 ± 1.0	−0.005 ± 0.001	0.06 ± 0.02	0.09 ± 0.03
CO-noE Mean	−0.6 ± 1.0	1.1 ± 1.9	0.6 ± 1.0	−0.018 ± 0.005	0.34 ± 0.07	0.29 ± 0.03

Download table as:  ASCIITypeset image

The MCMC results are particularly interesting. The MCMC appears to have much less Gaussian deviation distributions, being very strongly peaked about the true parameters, but with low levels of extreme outliers. For example, in the no-error case in Figure 11, between 55% and 68% of the T_k values were recovered very nearly perfectly, but there is a low level tail of underestimates extending out to a factor of 300.

The mean parameters of the MCMC traces should be more meaningful than the best fit parameters. This is consistent with the dv/dr deviations for the randomized MCMC model, where the best fit values almost form a top hat distribution, while the mean values are very nearly normally distributed.

The CO results are consistently worse than the HCN results, with a far greater tendency toward extremely low T_k and high ${n}_{{{\rm{H}}}_{2}}$. It is not clear why this is the case, although further examination revealed that the outliers all had very poor MCMC traces due to non-physical synthetic fluxes. These arose even when RADEX converged and optical depths were reasonable: there were specific parameter combinations which appear to produce pathological fluxes (as was found for some combinations in the T_bg exploration). Small variations away from the specific input parameters leads to the recovery of realistic SLEDs. While no method was found to reliably exclude these values from the synthetic fluxes, we do not consider this of great concern: when fitting to real data, the MCMC will find very poor fits for these values and so exclude them, jumping to a nearby point. Since the surrounding phase space will produce good fits, this excludes only an N-dimensional delta-function from the trace, and the results are still reliable. This is supported by the inability of the MCMC to fit these pathological synthetic fluxes: if the responsible phase-volume was non-zero, the MCMC would still be able to "home in" on it. Therefore, this is only a problem for the synthetic fluxes and not for real observations. For the synthetic fluxes it produces the extreme outliers, which are clearly distinguishable post-hoc. The pseudo-Gaussian spread in the deviations is unrelated to this issue, and is due to the many-to-one relationship between physical conditions and line ratios, which is the focus of this work.

These grid and MCMC results demonstrate that even under ideal conditions a 10% calibration uncertainty corresponds to uncertainties in LVG-derived parameters of approximately 0.5 dex/a factor of 3. This is believed to be due to both the slow variation of the line ratios across the parameter space and the many-to-one relationship between physical conditions and line ratios.

The results of the Randomized, with Isotopologues trials with four H¹³CN lines and one H¹³CN line are shown in Figures 12 and 13, respectively, and the numerical results are collated in Table 3.

Table 3.  Deviations of Recovered Parameters from MCMC ¹³C Isotopologue Models

		Four H13CN lines		One H13CN line
		Best	Mean	Best	Mean
Statistical	Tk	−0.1 ± 0.4(0.04)	−0.1 ± 0.4(0.04)	−0.1 ± 0.5(0.1)	−0.1 ± 0.5(0.1)
	${n}_{{{\rm{H}}}_{2}}$	−0.0 ± 0.7(0.14)	−0.1 ± 0.7(0.14)	0.0 ± 0.6(0.2)	0.0 ± 0.6(0.2)
	${dv}/{dr}$	0.0 ± 0.4(0.18)	0.0 ± 0.4(0.18)	0.0 ± 0.4(0.2)	0.0 ± 0.3(0.2)
	$\frac{[{\rm{HCN}}]}{[{{\rm{H}}}^{13}{\rm{CN}}]}$	0.01 ± 0.21(0.03)	0.01 ± 0.22(0.03)	0.05 ± 0.28(0.05)	0.05 ± 0.29(0.05)
Median	Tk	−0.003 ± 0.002	−0.005 ± 0.002	−0.006 ± 0.003	−0.010 ± 0.003
	${n}_{{{\rm{H}}}_{2}}$	−0.019 ± 0.015	−0.025 ± 0.010	0.015 ± 0.016	−0.007 ± 0.012
	${dv}/{dr}$	−0.063 ± 0.016	−0.034 ± 0.009	0.029 ± 0.016	−0.006 ± 0.010
	$\frac{[{\rm{HCN}}]}{[{{\rm{H}}}^{13}{\rm{CN}}]}$	0.002 ± 0.002	−0.001 ± 0.003	0.005 ± 0.003	0.009 ± 0.003

Note. Values in parenthesis are the MCMC standard deviations.

Download table as:  ASCIITypeset image

The figures demonstrate the (well known) power of adding isotopologue lines. Even with the additional free parameter of [HCN]/[H¹³CN], the model places much tighter and more accurate constraints on T_k and ${n}_{{{\rm{H}}}_{2}}$, as well as placing very tight constraints on [HCN]/[H¹³CN]. This is not simply an effect of having more lines to work with since the fits are better than those obtained using the 13 CO lines. The good T_k results are due to the strong degeneracy between T_k and [HCN]/[H¹³CN]. This can be roughly understood as follows: T_k and ${n}_{{{\rm{H}}}_{2}}$ set the primary shape of the SLEDs of HCN and H¹³CN—with this determined, the model can then find [HCN]/[H¹³CN] since, to first order, this corresponds to a scaling of the SLED but not a change in shape.⁹ Put another way: iso-ratio contours for HCN and H¹³CN intersect roughly perpendicularly in ${T}_{{\rm{k}}}-{n}_{{{\rm{H}}}_{2}}$ plots. None of this needs to be entered into the model: it explores the posterior probability space, without needing to understand the underlying relationships. With a single H¹³CN line, the ability to uniquely identify a T_k is lost, again due to the degeneracy between T_k and [HCN]/[H¹³CN].

For the four H¹³CN lines, no-error results are exceptionally good, with the models recovering the correct T_k and [HCN]/[H¹³CN] in almost 100% of the cases, and ${n}_{{{\rm{H}}}_{2}}$ and dv/dr being exceptionally well constrained. There are, however, a few extreme outliers which are responsible for the statistical deviations in Table 3 being so much greater than those of the medians. The single H¹³CN line fits are still very good; although worse than the 4 H¹³CN line fits, they are still better than the HCN-only models, despite the additional free parameter.

The value of multiple ¹³C isotopologue lines is clear. However, despite only being 10–15 times fainter than the ¹²C lines,¹⁰ the acquisition of even 1 H¹³CN line in a local galaxy $(z\leqslant 0.1)$ represents a significant time investment, even with ALMA. The small size of these sources, coupled with the inherent faintness of the ¹³C lines, leads to them being washed out by beam dilution, requiring even more time to observe with single-dish instruments. The one exception to this is ¹³CO, which is sufficiently abundant to emit relatively brightly and be quite easily detectable. In cases where accurate and precise non-LTE gas parameters are desired from molecular observations, it is imperative that ¹³C isotopologue lines be observed.

The results of the third model, with free and fixed X_HCN, are shown in Figure 14 and Table 4. Note in particular the much wider ranges of the x-axes than in the previous plots.

Table 4.  Deviations of Recovered Parameters from MCMC Models For the Free and Assumed X_HCN

		Free XHCN		Fixed XHCN
		Best	Mean	Best	Mean
Statistical	Tk	−0.1 ± 0.7(0.06)	−0.1 ± 0.6(0.06)	0.0 ± 0.7(0.08)	0.0 ± 0.6(0.08)
	${n}_{{{\rm{H}}}_{2}}$	0.5 ± 1.4(0.4)	0.5 ± 1.3(0.4)	0.0 ± 1.2(0.2)	−0.1 ± 1.2(0.2)
	${dv}/{dr}$	0.2 ± 0.7(0.3)	0.3 ± 0.7(0.3)	0.0 ± 0.7(0.2)	0.0 ± 0.6(0.2)
	XHCN	−0.7 ± 2.4(0.7)	−1.2 ± 2.1(0.7)	0.2 ± 2.3	0.2 ± 2.3
Median	Tk	−0.001 ± 0.003	−0.011 ± 0.004	0.044 ± 0.005	0.035 ± 0.003
	${n}_{{{\rm{H}}}_{2}}$	0.32 ± 0.05	0.34 ± 0.05	−0.10 ± 0.04	−0.12 ± 0.04
	${dv}/{dr}$	0.15 ± 0.03	0.18 ± 0.03	−0.05 ± 0.03	−0.05 ± 0.02
	XHCN	−0.63 ± 0.13	−1.07 ± 0.09	0.06 ± 0.14	0.05 ± 0.13

Download table as:  ASCIITypeset image

There is no clear sign of any difference in the recoverability of T_k or ${n}_{{{\rm{H}}}_{2}}$ in either model. Furthermore, the molecular abundance in the free parameter case is essentially unconstrained and although the model is slightly able to recover the true value there is a huge spread in the results as well as a systematic offset of almost −1 dex. It would appear therefore to be perfectly acceptable to assume a canonical molecular abundance. If two species are available, and their abundance ratios are of interest, it should therefore be acceptable to fix one abundance and have as a free parameter the abundance ratio itself, thereby removing a free parameter. Further testing is needed, however, to see whether this result holds when the synthetic line ratios emerge from clouds with ${K}_{{\rm{vir}}}\ne 1$.

The results of the Assumed [HCN]/[H¹³CN] model using one H¹³CN line and with 10% errors are shown in Figure 15 and Table 5.

Table 5.  Deviations of Recovered Parameters from MCMC Models for an Assumed [HCN]/[H¹³CN]

	Assumed [HCN]/[H13CN]
		Best	Mean
Statistical	${T}_{{\rm{k}}}$	−0.1 ± 0.7(0.03)	−0.1 ± 0.7(0.03)
	${n}_{{{\rm{H}}}_{2}}$	−0.3 ± 0.9(0.1)	−0.3 ± 0.9(0.1)
	${dv}/{dr}$	−0.2 ± 0.6(0.1)	−0.2 ± 0.6(0.1)
	$\frac{[{\rm{HCN}}]}{[{{\rm{H}}}^{13}{\rm{CN}}]}$	−0.3 ± 0.6	−0.3 ± 0.6
Median	Tk	−0.042 ± 0.014	−0.032 ± 0.009
	${n}_{{{\rm{H}}}_{2}}$	−0.20 ± 0.04	−0.23 ± 0.04
	${dv}/{dr}$	−0.29 ± 0.03	−0.232 ± 0.026
	$\frac{[{\rm{HCN}}]}{[{{\rm{H}}}^{13}{\rm{CN}}]}$	−0.282 ± 0.028	−0.281 ± 0.029

Download table as:  ASCIITypeset image

The striking result here is that the parameters are recovered much more poorly than in Figure 13, where we ran the same model but without fixing [HCN]/[H¹³CN] to an assumed value. Simultaneously, the MCMC trace reports much smaller errors. In other words, the MCMC fixes very tightly on the wrong point in ${T}_{{\rm{k}}}-{n}_{{{\rm{H}}}_{2}}$ space. That is, we are actually better off adding [HCN]/[H¹³CN] as a free parameter. In doing so, we will not see as dramatic reductions in the parameter space as can be seen by fixing the [HCN]/[H¹³CN] ratio, but we are also less likely to exclude the true parameters.

Fixing [HCN]/[H¹³CN] to a canonical value causes three problems: (1) due to the degeneracy between [HCN]/[H¹³CN] and T_k, fixing [HCN]/[H¹³CN] artificially restricts T_k, and if [HCN]/[H¹³CN] is not close to the canonical value this then offsets T_k. However, since T_k and ${n}_{{{\rm{H}}}_{2}}$ are also partially degenerate, this also offsets ${n}_{{{\rm{H}}}_{2}}$. Then, since the virial state of the gas is dependent upon ${n}_{{{\rm{H}}}_{2}}$ and dv/dr, dv/dr is also affected (for a fixed ${K}_{{\rm{vir}}}$). (2) As well as increased deviations, there is a systematic offset introduced into the physical parameters. (3) An incorrectly assumed [HCN]/[H¹³CN] can force T_k, and hence ${n}_{{{\rm{H}}}_{2}}$, into unphysical regions where RADEX becomes unstable, leading to wildly incorrect parameter recovery.

When using an MCMC there is a further issue, which is that the MCMC will dramatically understate the true uncertainty on the recovered parameters. Fixing [HCN]/[H¹³CN] is a very strong prior, so that the MCMC will show only small oscillations about the mean parameters. This is precisely why ¹³C isotopologue line observations are so powerful: they are extremely good at breaking the degeneracies found in single species models. However, this means that if the wrong isotopologue abundance ratio is assumed, then the recovered parameters will be forced to be incorrect. This effect is so large that you can actually be worse off assuming a canonical isotopologue abundance ratio than if you had no isotopologue lines at all.

The HCN and CO SLEDs for ranges of T_bg and T_k are shown in Figure 2. Each panel corresponds to the indicated T_bg, increasing from left to right, and we plot the line flux in arbitrary units against the J_u of the line transition. All panels share the same vertical scale. The color corresponds to T_k from cold (10 K, sky blue) to hot (1000 K, dark magenta). While it is usually the case that ${T}_{{\rm{k}}}\geqslant {T}_{{\rm{dust}}}\geqslant {T}_{{\rm{CMB}}}$, we include regions where the kinetic temperature is lower than the incident blackbody radiation temperature for completeness.

Zoom In Zoom Out Reset image size

We note that these are the observable fluxes; i.e., the line flux seen in excess of the continuum level, which RADEX provides by default. This "loss of contrast" is the primary cause of the "shrinking" of the SLEDs as T_bg increases.

As expected, the SLEDs fall into absorption when their T_k is less than T_bg. All of the fluxes fall as T_bg increases: this is due to the reduced contrast of the line emission against the hotter background radiation field (e.g., Papadopoulos et al. 2010). The CO SLEDs display fundamentally similar, if offset, behavior to the HCN SLEDs.

This model is, as ever, a gross simplification of reality, where hot background radiation fields on 10 pc scales are usually due to active-galactic-nucleus- (AGN-) or starburst-heated dust, which also emits significantly in the mid-IR, leading to optical pumping of molecular vibrational lines (Aalto et al. 2007, 2015a, 2015b; González-Alfonso et al. 2013). Nevertheless, the low T_bg SLEDs, where the background temperature of 10 K is reached by the CMB at z = 2.3, demonstrate that line ratios indicative of certain behaviors in the local universe are not necessarily applicable in high-redshift galaxies. Similarly, any investigation into correlations between line ratios and galaxy properties must take into account variation in the incident radiation field, as well as variation within individual galaxies.

It is noteworthy that since RADEX presents the line brightness in excess of the background, if there is a hot background field, but it is completely obscured by the foreground gas and dust, then the background temperature should be added back onto the line brightness produced by RADEX. On the other hand, this very specific geometry is probably rare in nature, and may be hard to identify. Furthermore, in these cases it is highly likely that only the outer edges of the very dense clouds are visible, where the gas is not exposed to the full force of the intense background field, further complicating the analysis.

Line ratios are frequently invoked where there are few line observations of a source. In principle, they allow an exploration of the molecular gas conditions and there have been extensive efforts to develop diagnostic molecular line ratios, e.g., as a means of distinguishing between obscured AGNs and obscured starbursts (Kohno 2003, 2005; Meijerink & Spaans 2005; Imanishi et al. 2006; Meijerink et al. 2007; Krips et al. 2008). More recently, Izumi et al. (2013) have attempted to develop diagnostic plots using the HCN(4–3)/HCO + (4–3) and HCN(4–3)/CS(7–6) line ratios. However, these studies are frequently single-dish observations of entire galaxies. There are therefore potential concerns that gas and dust temperature gradients in the galaxies may have significant effects on the integrated line ratios,¹¹ especially if the ratios are dependent on T_bg, which, in active galaxies in the local universe, is dominated by dust emission. If these effects are varying between galaxies, then this may lead to either scatter in true relationships or anomalous relationships where none exist. Furthermore, extending these line ratios to higher redshift, where T_CMB is elevated, could present further problems.

We present two-dimensional (2D) plots of the HCN line ratios as functions of T_k and T_bg in Figure 3. There is a clearly visible degeneracy between T_k and T_bg, especially in ratios with the J = 4–3 line. While the low-J line ratios do not change significantly, the HCN $\frac{(4-3)}{(3-2)}$ line ratio changes dramatically between T_bg = 3 K, T_k = 10 K, and T_bg = 30 K, T_k = 1000 K, rising rapidly.

Zoom In Zoom Out Reset image size

The results of our model exploring the effect of T_bg on ${T}_{{\rm{B,peak}}}$ as a function of T_k are shown in Figures 4 and 5 where we plot the peak brightness temperature for the first nine CO rotational transitions and first six HCN rotational transitions as functions of T_k for a range of redshifts. The models demonstrate how cooler molecular gas, such as the at outer edges of molecular reservoirs or in interclump regions, can become indistinguishable from the CMB. If galaxies at high-z possess kinetic temperature gradients or harbor cold molecular gas environments, the aforementioned effect could lead to a bias in the gas morphologies inferred from low-J CO lines. This is not a new or surprising result; from a purely analytical approach, this is to be expected since the total intensity is:

$$\begin{eqnarray}&&{I}_{\nu }^{{J}_{u}}\propto {S}_{\nu }^{{J}_{u}}(1-{e}^{-{\tau }_{\nu }^{{J}_{u}}})+{B}_{\nu }({T}_{{\rm{CMB}}}){e}^{-{\tau }_{\nu }^{{J}_{u}}},\end{eqnarray} \tag{ 3 }$$

where S is the source function, and therefore that the intensity in excess of the background field (here just the CMB) is:

$$\begin{eqnarray}&&{\rm{\Delta }}{I}_{\nu }^{{J}_{u}}\propto {S}_{\nu }^{{J}_{u}}(1-{e}^{-{\tau }_{\nu }^{{J}_{u}}})+{B}_{\nu }({T}_{{\rm{CMB}}})({e}^{-{\tau }_{\nu }^{{J}_{u}}}-1),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&=(1-{e}^{-{\tau }_{\nu }^{{J}_{u}}})({S}_{\nu }^{{J}_{u}}-{B}_{\nu }({T}_{{\rm{CMB}}})),\end{eqnarray} \tag{ 5 }$$

so that when the gas is emitting approximately thermally, or ${T}_{{\rm{ex}}}\simeq {T}_{{\rm{CMB}}}$, then ${S}_{\nu }^{{J}_{u}}\simeq {B}_{\nu }({T}_{{\rm{ex}}})$, and the line becomes indistinguishable from the CMB (e.g., da Cunha et al. 2013; Wilson et al. 2013, p. 5). This will introduce an additional bias toward detecting hotter, more excited gas as redshift increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The effect is more pronounced for HCN: this is due to the lines being naturally fainter than the CO lines (due to the much lower abundance of HCN cf CO), and so the fall in contrast against the CMB in much more pronounced. There is also a similar effect for dust, with the same associated implications on observability (da Cunha et al. 2013).

The results of our toy model are shown in Figures 6–8. As described in Section 2.10, this model distributes 10⁴ clouds exponentially with radius across an artificial galaxy, with the clouds identical except for their T_k, which is determined by their radius. The galaxy is then placed at a range of redshifts and the galaxy-wide emission integrated. The SLED in Figure 6, left, shows how the flux, relative to the background, decreases with redshift, but also how the peak of the SLED shifts to higher J. These results are consistent with da Cunha et al. (2013).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The right SLED of Figure 6 shows the same data, but with the SLED presented relative to the intensity of the J = 1–0 line. These SLEDS are in turn consistent with Narayanan & Krumholz (2014). Note, however, that da Cunha et al. (2013) and our Figure 6 (left) show that the increase of the SLED with redshift, defined relative to the J = 1–0 line, is primarily due to the much greater reduction in line contrast of the low-J lines. While the increasing T_bg does indeed raise the high-J lines, this effect is small in comparison to the contrast reduction of the low-J lines.

This reduction in the overall SLED is important, as it leads to a higher α_CO,¹² e.g., Equation (1) in Ivison et al. (2011) and Equations (17) and (18) in Bolatto et al. (2013). Since the molecular gas mass of these models is identical at all redshifts, and since $L{^{\prime} }_{\mathrm{CO}(1\mbox{--}0)}$ "de-redshifts" emission (i.e., for two identical SLEDs at z = 0 and z = 2, excluding CMB effects, $L{^{\prime} }_{\mathrm{CO}(1\mbox{--}0)}$ will be identical), any difference in $L{^{\prime} }_{\mathrm{CO}(1\mbox{--}0)}$ in these models is solely due to changes in the CO(1–0) line flux, which is solely due to the CMB. We show these effects in Figure 7, where we plot the line luminosities relative to the z = 0 line luminosities (left) and the line luminosities relative to the z = 0 CO(1–0) line luminosity. For the adopted temperature profile and cloud distribution profile this effect is very large, with 80% of the mass recovered at z = 2 if a standard α_CO is adopted, 60% at z = 4, and by z = 8 only 35% of the mass is recovered.

We note that while Figure 7 (left) shows that the high-J line luminosity increases significantly with redshift, Figure 7 (right) shows that it still only traces a small fraction of the gas (for this toy model). The CO(1–0) line at z = 0 traces all of the molecular clouds, so the fraction of this luminosity recovered is directly comparable to the mass that would be estimated if a standard z = 0 α_CO factor was used. The HCN emission only traces the dense molecular gas, so the fall in the HCN integrated emission is representative of the dense gas that would be missed, not the total molecular mass.

In Figure 8 we present the line brightness temperature CO and HCN line ratios for our toy model galaxy. We present a selection of line ratios only. It is clear that there is very little variation in the galaxy-integrated line ratios with redshift, except for in the $\frac{5-4}{1-0}$ and $\frac{4-3}{1-0}$ line ratios, although in all cases there is greater variation in the HCN ratios than the the CO ratios. This is, as in the T_bg investigation above, due to HCN simply being fainter than CO, so is more strongly affected by the brightening of the CMB. However, these ratios are specific to the galaxy excitation model used here, and are not generally applicable: there is likely to be a set of ratios which vary significantly, but the particular ratios will depend upon the peaks of the galaxy SLEDs, which in turn depends upon excitation conditions. The excitation conditions used here are rather weak, especially compared to the very active star-forming models of Narayanan & Krumholz (2014). In general, however, the galaxy-integrated line ratios appear to be rather robust to changes in redshift, unlike the individual line of sight HCN brightness temperature ratios in Section 3.5.

LVG models are quick and easy to run, and have been used to a varying extent in Galactic, local universe, and high-redshift studies of molecular gas. One of most common and indeed most reliable uses of the models is to generate grids of results, then overlaying either χ² or likelihood contours (e.g., García-Burillo et al. 2000; Rangwala et al. 2011). This graphically illustrates the reasonable parameter ranges, and does not necessarily claim specific values. However, even these methods are susceptible to assumptions which can bias the results in unexpected and unusual ways, especially if canonical abundances or abundance ratios are assumed. For example, assuming an [HCN]/[H¹³CN] ratio of 60 actually imposes strong constraints on the possible T_k range. If the true [HCN]/[H¹³CN] ratio is closer to 300 (such as in NGC 6240: Papadopoulos et al. 2014; Tunnard et al. 2015b), then the assumptions can throw off T_k by over a dex. Furthermore, when multiple lines and species are available, these approaches can be restrictive: they cannot be used to explore higher dimensional parameter spaces.

One great reassurance of the tests presented here is that there is no significant systemic offset in the model deviations, with most trials returning centrally peaked deviation distributions. In almost all models there was a small set of extreme outliers, almost always with an underestimate of T_k and concomitant overestimate of ${n}_{{{\rm{H}}}_{2}}$. These outliers are greatly exacerbated by including regions with a thermal pressure >10⁸ K cm⁻³. This regime is so extreme that its exclusion is unlikely to remove many realistic regions with significant emission, especially when integrated across a galaxy.

The concern arises when model parameters are interpreted as hard and fast values, or as being known better than within a factor of 3. Even in the most optimistic scenarios and with zero observational uncertainties, LVG models cannot be expected to provide answers more accurately than ±0.2 dex. The real uncertainties on the results are most likely much higher, presenting real and significant concerns with any possible precision uses of LVG models.

A topic we have not mentioned until now is that of the collisional rate constants. LVG models require these (temperature dependent) rates which are primarily determined in laboratory experiments. In the models presented here, we should be largely insensitive to any errors in these rates, as we are using the same rates for generating both the original lines and comparison lines (whether they are in grids or MCMC). In real world observations, however, these errors could be significant, especially for the less well studied molecules (such as HNC, a chemically interesting isomer of HCN). Furthermore, RADEX does not extrapolate the rate constants beyond the provided data, keeping them constant instead.¹³ We are not suggesting that this is the wrong approach; indeed, any extrapolation of rate constants must be performed with utmost caution. We are simply noting that this should be appreciated as a limitation of LVG, and indeed any molecular line radiative transfer models, including LIME. Nevertheless, we have no means of quantifying these potential errors with our models, and as such consider the issue beyond the scope of this paper, leaving it as a further reason to consider the model precisions presented here as best cases only.

No discussion of the reliability of LVG models would be complete without mentioning chemistry. The LVG framework does not include chemistry; chemistry can be added to LVG models ante-hoc, post-hoc, or iteratively, but in all cases this is just altering the molecular abundances. In other words, LVG codes, e.g., RADEX, can find excellent fits for completely unrealistic physical conditions and molecular abundances.

To some extent, unrealistic physical conditions may be excluded by the use of appropriate priors. Similarly, a skeptical interpretation of the LVG results is essential, preferably including an examination of the appropriateness of the derived molecular abundance ratios for the physical conditions. A potentially very powerful technique for this would be post-hoc chemical modeling or comparison with lookup tables to verify the results are self-consistent, somewhat in the fashion of Viti et al. (2014). Unlike running chemistry for every step in the MCMC, this method would not add significantly to the run time while providing significantly more physical insight and potentially highlighting interesting results.

We also point to the risks of modeling emission from entire galaxies as single clouds. Not only is this physically unmotivated, it produces unrealistic column densities and masses. Dynamically motivated and self-consistent methods exist which are also physically meaningful, and although they may require prior assumptions. These assumptions are in all cases more reasonable than the tacit assumptions implicit in the single cloud alternative (see appendices in Papadopoulos et al. 2012; Tunnard et al. 2015b).

There are some subtleties associated with assuming a virial state for the modeled gas. Since, for a uniform density cloud, the mass (in kg) is given by

$$\begin{eqnarray}&&{M}_{{\rm{cloud}}}=\displaystyle \frac{4\pi }{3}{\left(\displaystyle \frac{\zeta {\rm{\Delta }}v}{2{dv}/{dr}}\right)}^{3}{n}_{{{\rm{H}}}_{2}}\;{\mu }_{{{\rm{H}}}_{2}},\end{eqnarray} \tag{ 6 }$$

where ${\mu }_{{{\rm{H}}}_{2}}$ is the helium-corrected mass per H₂ molecule, and ζ = 3.08 × 10¹⁸ cm pc⁻¹, or equivalently that:

$$\begin{eqnarray}&&{r}_{{\rm{cloud}}}=\displaystyle \frac{\zeta \;{\rm{\Delta }}v}{2{dv}/{dr}}.\end{eqnarray} \tag{ 7 }$$

The issue is that if a virial state is assumed, then for an observed line width the size of the cloud is fixed. Therefore, if the assumed virial state of the cloud is changed (through dv/dr), and all other variables held constant, then the size of the cloud has also been changed. This subtlety is important both for resolved and extragalactic observations: if the size of the cloud is known from observations then this should be fixed in the modeling. If the size of the cloud is unknown, then there is a degeneracy between the size of the cloud and multiple clouds distributed in velocity space across the galaxy.

Further discussion of this topic can be found in Papadopoulos et al. (2012) and Tunnard et al. (2015b). Here, we simply present a comparison of the peak brightness temperature from RADEX models where K_vir is varied, first without any constraint on the cloud size, and then with the cloud size constrained to be 1 pc, with the results shown in Figure 9. In the unconstrained case, there is less emission in all lines as K_vir increases: this is due to the cloud shrinking, and the column density of CO falling. In the r = 1 pc case, however, there is a pivot around the J_u ≃ 5–6 lines, with emission increasing with K_vir at lower J and decreasing at higher J. Here, the column density is the same but the lines become less optically thick as K_vir increases, so more of the column can be seen, and emission increases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The key point here is that there are strong variations in the line emission as well as the relative line emission with variations of K_vir and dependent upon the underlying assumptions. In situations where the cloud size is known this can be particularly important.

We have presented several models of the potential effects of the background radiation field on LVG model results. These models do not, and cannot, account for dust within the emitting gas, as this is not a feature included in RADEX. They can, however, model gas illuminated by a separate hot dust background, and the effects of the CMB at higher redshift.

As we have noted, it is also important to recognize that RADEX presents line results as line brightnesses in excess of the background field. When dealing with the CMB this is ideal, however, when dealing with a hot dust source that is completely obscured by cooler, intervening dust this is misleading. In these cases, it may be more appropriate for the applied background field to be re-added to the line brightnesses. On the other hand, due to the high column depths and line optical depths of CONs, we are only able to observe the cooler, shielded outer layers in the major isotopologue lines (Aalto et al. 2015a), at which point the background dust has cooled significantly and a high T_bg is not actually incident on the observable regions.

We have demonstrated the well known degeneracy between T_k and T_bg, and demonstrated the effects of a hot dust background on HCN and CO SLEDs (Figure 2). These effects can be huge, including the effects on the line ratios. However, these conditions are likely to be significant in a minority of objects, and more comprehensive modeling of dust included in the gas phase is beyond the scope of this paper. Papadopoulos et al. (2010) argued that if the gas and dust are in thermal equilibrium, then any heating of the dust will increase the line and continuum brightness equally, thereby reducing the line/continuum contrast. While essential when analyzing objects such as Arp 220, this is unlikely to be a concern in most galaxies.

At z = 4, both CO(4–3) and CO(5–4) lie in ALMA Band 3 while the lower J lines, observed by Carilli et al. (2010, 2011) and Hodge et al. (2012), require the JVLA. These studies have revealed a clumpy medium: precisely what is expected if the real conditions are GMC-like clumps embedded in an extended diffuse medium. The low-excitation diffuse medium becomes undetectable against the hotter CMB and only the high-density clumps are seen (see also Casey et al. 2015). The effects on planned blind surveys of CO with the ngVLA and SKA are more complex, but will risk introducing a significant redshift bias on both the low-J CO luminosity function and the H₂ mass fractions, as well as the cosmological ${\rm{\Omega }}({{\rm{H}}}_{2})$ (Casey et al. 2015; Leroy et al. 2015). How significant this effect will be is not clear, as it depends heavily on the kinetic temperature of the molecular gas at high-z. If large fractions of the molecular reservoirs are ≲20 K, then the effects are likely to be very significant as early as z ≃ 3–4. However, if high-z molecular reservoirs are mostly ≥30 K, then we do not expect any significant effect even out to z = 8.

The change in α_CO with redshift highlights the importance of dynamical mass estimates for high-z line observations (Carilli et al. 2010, 2011; Hodge et al. 2012), and presents new challenges for high-z mergers, where dynamical mass estimates may no longer be meaningful. Our model, while very simple, demonstrates that even at z = 2 there are risks of significantly underestimating the molecular masses of galaxies due to the loss of contrast against the CMB, including the potential invisibility of the coldest, yet perhaps in some cases also the densest, regions of the galaxies.

The line ratios most affected by the higher CMB are those including the low-J lines, so that higher J line ratios should remain relatively robust even at redshifts as high as z = 8, depending upon the galaxy excitation.

We include a final model of less general interest, but which is nevertheless important for appreciating the limitations of LVG codes. We explore the case where we have four HCN lines and four H¹³CN lines, and where both X_mol [HCN]/[H¹³CN] are free parameters.

One of the key parameters necessary to predict the dominant heating source in a gas (PDRs, XDRs, MDRs, etc.) is the molecular abundance. While a single molecular abundance is not particularly informative, with reliable estimates of the abundances of several species significant progress can be made (Meijerink & Spaans 2005; Meijerink et al. 2007, 2011; Bayet et al. 2008; Viti et al. 2014). An interesting question therefore is whether LVG models can in fact be used to reliably recover molecular abundances. While we showed in the main paper that with HCN lines alone this is not possible, we investigate whether the addition of H¹³CN improves the situation.

Due to the degeneracy between X_mol and dv/dr, it is not possible to recover the molecular abundance with a single species, unless one is willing to assume a fixed value of K_vir. If K_vir is assumed or known, then for a given ${n}_{{{\rm{H}}}_{2}}$, dv/dr is determined. However, even this is not necessarily enough, due to the degeneracy between T_k and ${n}_{{{\rm{H}}}_{2}}$. Furthermore, the assumption of a fixed K_vir is dubious, especially in luminous infrared galaxies where K_vir is very likely elevated (Papadopoulos et al. 2014). This is before we take into account that even high-resolution interferometric extragalactic observations (in all but the closest galaxies) will confuse multiple gas phases within the beam.

One possible solution appears to lie in isotopologue lines. Here, we present a model with four HCN lines and four H¹³CN lines, with the free parameters T_k, ${n}_{{{\rm{H}}}_{2}}$, dv/dr, X_HCN, and [HCN]/[H¹³CN]. The results are shown in Figure 16. Unfortunately, while T_k and [HCN]/[H¹³CN] are well determined, ${n}_{{{\rm{H}}}_{2}}$, dv/dr, and ${X}_{{\rm{HCN}}}$ are not. The conclusion here is that while LVG models cannot reliably estimate molecular abundances, they can very reliably estimate abundance ratios.

As a continuation of the discussion of the effect of T_bg on the observed line fluxes and line brightness ratios in Section 3.5, we present the CO line ratios in Figures 17 and 18. It is not practical to present all 78 possible line ratios, so we include only the 12 ratios with the J = 1–0 line and the 12 $J+1/J$ line ratios.