The Star-Forming Molecular Gas in High Redshift Submillimeter Galaxies

We present a model for the CO molecular line emission from high redshift Submillimeter Galaxies (SMGs). By combining hydrodynamic simulations of gas rich galaxy mergers with the polychromatic radiative transfer code, Sunrise, and the 3D non-LTE molecular line radiative transfer code, Turtlebeach, we show that if SMGs are typically a transient phase of major mergers, their observed compact CO spatial extents, broad line widths, and high excitation conditions (CO SED) are naturally explained. In this sense, SMGs can be understood as scaled-up analogs to local ULIRGs. We utilize these models to investigate the usage of CO as an indicator of physical conditions. We find that care must be taken when applying standard techniques. The usage of CO line widths as a dynamical mass estimator from SMGs can possibly overestimate the true enclosed mass by a factor ~1.5-2. At the same time, assumptions of line ratios of unity from CO J=3-2 (and higher lying lines) to CO (J=1-0) will oftentimes lead to underestimates of the inferred gas mass. We provide tests for these models by outlining predictions for experiments which are imminently feasible with the current generation of bolometer arrays and radio-wave spectrometers.


I N T RO D U C T I O N
Understanding the origin and evolution of active galaxy populations at cosmological redshifts remains an outstanding problem. Of particular interest is a population of submillimetre-luminous galaxies discovered via deep, blind surveys with Submillimetre Common-User Bolometer Array (SCUBA) on the James Clerk Maxwell Telescope (JCMT) at a median redshift of z ∼ 2 ( Barger et al. 1998;Hughes et al. 1998;Chapman et al. 2003a). These Submillimetre Galaxies (SMGs) are empirically defined with the flux limit S 850 μm 5-6 mJy 1 (e.g. Blain et al. 2002), and appear to be a highly clustered population of galaxies forming stars at prodigious rates [star formation rate (SFR) 10 3 M yr −1 ; Blain et al. 2002Blain et al. , 2004Swinbank et al. 2004;Kovács et al. 2006;Menéndez-Delmestre et al. 2007;Valiante et al. 2007;Coppin et al. 2008a]. In addition to undergoing prodigious star formation, some SMGs additionally host heavily obscured active galactic nuclei (AGN; e.g. Ivison et al. 2002;Alexander et al. 2005aAlexander et al. , b, 2008Borys et al. 2005). These (e.g. Frayer et al. 1998Frayer et al. , 1999Genzel et al. 2003;Neri et al. 2003;Greve et al. 2005;Tacconi et al. 2006, and references therein). During this time, a number of pioneering papers have afforded the community a wealth of CO data from z ∼ 2 SMGs. The general molecular emission properties from SMGs can be summarized as follows.
H 2 gas masses. SMGs are incredibly gas rich, with massive H 2 gas reservoirs. Extrapolation from typical CO (J = 3-2) rest-frame measurements from SMGs suggests H 2 masses of the order of 10 10 -10 11 M , making these some of the most molecular gas-rich galaxies in the Universe (Greve et al. 2005;Solomon & Vanden Bout 2005;Carilli & Wang 2006;Tacconi et al. 2006;Coppin et al. 2008b;Tacconi et al. 2008). An analysis of semi-analytic models of SMG formation by Swinbank et al. (2008) found good agreement between the cold gas mass in their simulations and the observed values.
Linewidths. The CO linewidths in SMGs are typically broad with the median linewidth ranging from ∼600 to 800 km s −1 [full width at half-maximum (FWHM); Greve et al. 2005;Carilli & Wang 2006;Tacconi et al. 2006;Coppin et al. 2008b;Tacconi et al. 2008;Iono et al. 2009]. Comparisons to CO-detected quasars at comparable redshifts have had conflicting results. Analysis of literature data by Carilli & Wang (2006) has suggested that SMGs typically have broader linewidths than z ∼ 2 quasars by a factor of ∼2.5 with Kolmogorov-Smirnov (KS) tests indicating that the two are not drawn from the same parent population at the >99 per cent confidence level. On the other hand, new detections by Coppin et al. (2008b) have found rather similar CO linewidths from quasars and SMGs, and KS tests showing that the two arise from the same parent population at the >95 per cent confidence level.
Images. Interferometric CO imaging of SMGs has shown the spatial distribution of most SMGs to be relatively compact (R HWHM 1-2 kpc; Tacconi et al. 2006). Some individual sources have shown rather extended emission (e.g. Genzel et al. 2003). Curiously, images have shown evidence for compact disc-like motion in at least a few sources (Genzel et al. 2003;Tacconi et al. 2006), as well as extended emission in what appears to be interacting/merging galaxies (e.g. Ivison et al. 2001;Tacconi et al. 2006).
Excitation properties. While very few SMGs have been observed in multiple CO lines, the existing data suggest that these galaxies exhibit relatively high molecular excitation conditions. The CO line spectral energy distributions (CO SEDs; alternatively known as CO rotational ladders) are seen to typically turn over at the J = 5 level (∼83 K above ground; Wei et al. 2005Wei et al. , 2007Greve et al. 2005;Hainline et al. 2006). That said, few objects have been detected in their rest-frame CO (J = 1-0) line (Hainline et al. 2006). As H 2 masses are typically inferred from CO (J = 1-0) measurements, the excitation properties of higher lying levels (e.g. line ratios) with respect to the ground-state transition are crucial to constrain for the purposes of deriving molecular gas masses.
While the information provided by CO observations of z ∼ 2 SMGs is indeed invaluable, questions remain regarding the interpretation of many of the various aforementioned observational characteristics of these sources. For example, what is the true relationship between CO linewidths from SMGs and quasars? Is there an evolution in the CO linewidths of SMGs? How reliably can the CO linewidth be used as a dynamical mass tracer in these galaxies? How can observations of higher lying CO lines (in the rest frame) be extrapolated to the rest-frame CO (J = 1-0) luminosity in order to derive an H 2 gas mass (e.g. Hainline et al. 2006)? Can the observed CO emission linewidths, excitation properties and images be explained by a merger-driven scenario for SMG formation and evolution (e.g. 'scaled-up ULIRGs';Tacconi et al. 2006Tacconi et al. , 2008? It is clear that a theoretical interpretation behind the molecular line emission properties may be valuable for elucidating some of the aforementioned issues, as well as providing interpretation for forthcoming observations. Along with providing interpretation for observations, theoretical calculations of CO emission from simulated SMGs can provide direct tests of models of SMG formation and evolution. In Narayanan et al. (2009), we presented a merger-driven model for the formation of SMGs which reproduced the full range of observed 850 μm fluxes from SMGs, the optical-millimetre wave SED, and characteristic stellar, black hole and dark matter masses. Comparing the simulated molecular gas properties of these model SMGs to the extensive data sets in the literature provides a strict test of the models. In this paper, we investigate the CO emission properties from SMGs by combining these SMG formation models with 3D non-local thermodynamic equilibrium (LTE) molecular line radiative transfer calculations (Narayanan et al. 2006a(Narayanan et al. , 2008a. The goals of this paper are to (a) provide direct tests of merger-driven formation mechanisms for SMGs by comparing the simulated CO emission from the models of Narayanan et al. (2009) to observations and, given a sufficient correspondence between models and observations, (b) provide interpretation for existing and future observational data.
This paper is organized as follows. In Section 2, we describe our numerical methods for our hydrodynamic, molecular line radiative transfer and polychromatic SED radiative transfer simulations. In Section 3, we summarize the evolution of the submillimetre, B band and H 2 properties of our simulated galaxies. In Section 4.1, we study the CO morphology, molecular disc formation and CO spatial extents in SMGs. In Section 4.2, we explore the origin of the broad observed CO linewidths. We use these results to investigate the usage of CO as a dynamical mass tracer in Section 5, and analyse the CO excitation in SMGs in Section 6. In Section 7, we provide imminently testable observational predictions and, in Section 8, we summarize.

N U M E R I C A L M E T H O D S
Generically, our methodology involves three steps. First, we simulate the hydrodynamic evolution of galaxies utilizing the smoothedparticle hydrodynamics (SPH) code,   (Springel 2005). We then calculate the molecular line emission properties and SEDs from these simulated galaxies in post-processing utilizing the radiative transfer codes TURTLEBEACH (Narayanan et al. 2008a) and SUNRISE ). This method is summarized in the flowchart presented in Fig. 1. In this section, we present the details of the hydrodynamic and radiative transfer methods.

Hydrodynamics
A major goal of our programme is to understand the evolution of model SMGs while remaining constrained to observations for physical input to the simulations. When data for SMGs are unavailable, we turn to observational constraints from the Galaxy or local starbursts. We consider the formation of SMGs in gas-rich binary galaxy mergers at high redshift. This is motivated by observed radio, CO and optical morphologies of SMGs which appear to show signs of interactions (e.g. Chapman et al. 2003b;Tacconi et al. 2006Tacconi et al. , 2008, as well as theoretical models which demonstrate that mergers serve as an efficient means of triggering nuclear starbursts C 2009 The Authors. Journal compilation C 2009 RAS, MNRAS 400, [1919][1920][1921][1922][1923][1924][1925][1926][1927][1928][1929][1930][1931][1932][1933][1934][1935] Downloaded from https://academic.oup.com/mnras/article-abstract/400/4/1919/1079421 by guest on 28 July 2018 Figure 1. Flowchart summarizing methodology. The galaxies are simulated hydrodynamically using the SPH code,   (Springel 2005). The simulated molecular line emission and SEDs are then calculated in postprocessing utilizing the radiative transfer codes TURTLEBEACH (Narayanan et al. 2008a) and SUNRISE ), respectively. (Barnes & Hernquist 1991Mihos & Hernquist 1994. Moreover, simulations by Narayanan et al. (2009) have suggested that isolated galaxies and minor mergers (mass ratio 1:10) are unlikely to result in a S 850 > 5 mJy SMG.
The hydrodynamic simulations were performed with the SPH code,   (Springel 2005), which utilizes a fully conservative SPH formalism (Springel & Hernquist 2002). The hydrodynamic simulations include prescriptions for radiative cooling of the gas (Katz, Weinberg & Hernquist 1996;Davé et al. 1999) and a multiphase ISM in which cold clouds are considered to be in pressure equilibrium with hot gas (McKee & Ostriker 1977;Springel & Hernquist 2003). This multiphase ISM is implemented such that cold clouds grow through radiative cooling of hot gas, and heating from star formation can evaporate cold clouds . The effect of supernovae on the ISM is treated via an effective equation of state (EOS). Here, we employ the full multiphase EOS where supernovae-driven pressure optimally maintains disc stability (q EOS = 1; for more details, see fig. 4 of Springel et al. 2005).
Star formation occurs following a volumetric generalization of the Kennicutt-Schmidt relation, SFR ∝ n 1.5 (Schmidt 1959;Kennicutt 1998a, b;Springel & Hernquist 2003;Kennicutt et al. 2007), which results in discs consistent with the local surfacedensity scaling relations . The star formation timescale is chosen such that isolated disc models are consistent with the normalization of the local Kennicutt-Schmidt relations. While there are rather few measurements of observed SFR relations at cosmological redshifts, tentative evidence exists that z ∼ 2 galaxies may lie on the present-day Kennicutt-Schmidt relation (Bouché et al. 2007). Indeed, this is theoretically favourable as the relation SFR ∼ n 1.5 may be understood in terms of free-fall time arguments which are redshift invariant.
Black holes are included in the simulations as sink particles which accrete following a Bondi-Hoyle-Lyttleton parametrization according to the local gas density and sound speed (Bondi & Hoyle 1944;2 The main improvement in GADGET-3 over GADGET-2 (described by Springel 2005) is better load balancing on parallel processors. Bondi 1952). The bolometric luminosity of the black hole is set at L bol = Ṁ c 2 , where = 0.1. Feedback from the black hole(s) is modelled such that a fixed fraction of this luminosity (here, 5 per cent) couples thermally and isotropically to the surrounding ISM. This coupling efficiency is set to match the normalization of the local M-σ relation (Di Springel et al. 2005).
The progenitor disc galaxies are initialized with a Hernquist (1990) dark matter halo profile (for details, see Springel et al. 2005) with virial properties scaled to be appropriate for z ∼ 3 (such that they may represent galaxies undergoing a major merger by z ∼ 2; Bullock et al. 2001;Robertson et al. 2006). The progenitors begin with an initial gas fraction of f gas = 0.8, resulting in a gas fraction of f gas ∼ 20-40 per cent by the time the galaxies approach final coalescence. This is comparable to recent measurements by Bouché et al. (2007) and Tacconi et al. (2008) who find tentative evidence for gas fractions in z ∼ 2 SMGs of ∼40 per cent.
Here, we consider the Narayanan et al. (2009) series of 12 models, varying total mass, mass ratio and orbit. We include one additional lower mass merger which will not make an SMG (model SMG13), though will be useful for predictions of lower flux galaxies (Section 7). The initial conditions of the merger simulations are summarized in Table 1. The primary progenitors are initialized with circular velocities ranging from V c = 320 to 500 km s −1 , similar to measurements tabulated by Tacconi et al. (2008). This corresponds to halo masses of the order of M DM ∼ 10 12 -10 13 M , comparable to those inferred by clustering measurements (Blain et al. 2004;Swinbank et al. 2008). We model the formation of SMGs in 1:1 mergers in ∼10 12 , 5 × 10 12 and ∼10 13 M haloes, and 1:3 and 1:12 mergers with a M DM ≈ 10 13 M primary galaxy. We utilize three different initial orbits for the galaxies. For clarity, throughout this work, we primarily focus on the results from a fiducial merger (in Table 1, model SMG10) which produces a relatively average S 850 ≈ 5-7 mJy SMG. 3 The CO results (next section) are generic for all models considered here. There is a dispersion amongst the simulated submillimetre fluxes when varying merger orbit (even at a constant galaxy mass; see fig. 2 of Narayanan et al. 2009 for a direct quantification of this dispersion). The fiducial SMGs studied here lie in the middle of this dispersion for given halo masses and, as such, are typical. See Narayanan et al. (2009) for results regarding the physical properties of these merger simulations.

Molecular line radiative transfer
We utilize the 3D non-LTE molecular line radiative transfer code, TURTLEBEACH, to calculate the CO emission properties from our model SMGs (Narayanan et al. 2008a). TURTLEBEACH is an exact line transfer code that assumes full statistical equilibrium. The details of the code can be found in Narayanan et al. (2006b) and Narayanan et al. (2008a), and we refer the reader to these papers 3 When considering emission properties which are strongly dependent on merger mass, we will, at times, employ the usage of SMG1 as a 'high-mass' fiducial SMG along with SMG10. This allows us to bracket the range of galaxy masses which form SMGs in our simulations. SMG1 forms an SMG ranging from average (S 850 ≈ 5-7 mJy during inspiral) to extremely luminous and rare (S 850 ≈ 15 mJy during final coalescence), whereas SMG10 forms at peak (final coalescence) only an average S 850 ≈ 5-7 mJy SMG. See Narayanan et al. (2009) Table 1. SMGs are ordered with decreasing halo mass. Column 1 is the name of the model used in this work. Columns 2 and 3 are initial orientations for disc 1, Columns 4 and 5 are for disc 2, Column 6 gives the virial velocity of the progenitors, Column 7 gives their halo masses and Column 8 the mass ratio of the merger. Please see Narayanan et al. (2009) for generalized results regarding the physical properties of these galaxies.
Model for details on the specific algorithms used. Here, we summarize the details relevant to this work. We assume that all cold, star-forming gas in the hydrodynamic simulations is neutral. Following the observational constraints of Blitz & Rosolowsky (2006), we model the molecular fraction based on the ambient hydrostatic pressure: where k is Boltzmann's constant and the external pressure P ext is calculated via the fitting formula derived by Robertson et al. (2004): log P = 0.05(log n H ) 3 − 0.246(log n H ) 2 + 1.749(log n H ) − 10.6. (2) In Fig. 2, we show the radially averaged molecular gas fraction as a function of the galactocentric distance for three model SMGs. 4 In order to investigate numerous snapshots at relatively high temporal resolution (5 Myr) for a number of merger models, the hydrodynamic simulations were smoothed to a spatial resolution of ∼160 pc. Resolution tests presented in Narayanan et al. (2008a) which investigated the spectral and spatial invariance of these methods confirm that the resultant CO morphologies, lines profiles and excitation conditions are convergent at these spatial resolutions. In order to properly account for the molecular density gradients 4 We note that this feature of the code is an improvement over the algorithms described in Narayanan et al. (2008a). Specifically, previous works using TURTLEBEACH have assumed that half the neutral gas (by mass) was molecular, consistent with average conditions in local galaxies (Keres, Yun & Young 2003). While these assumptions tying the molecular gas fraction to global observations may reproduce average molecular emission patterns in galaxies and AGN (e.g. Narayanan et al. 2006aNarayanan et al. , 2008aNarayanan et al. 2008c), it does not provide spatially resolved information regarding the molecular content. By tying the molecular gas fraction to the ambient pressure as motivated by observations of GMCs (Blitz & Rosolowsky 2006), we are able to more accurately model the spatial distribution of the molecular gas. We note that even more sophisticated models for treating the neutral gas breakdown exist in the literature (e.g. Pelupessy, Papadopoulos & van der Werf 2006;Robertson & Kravtsov 2008). However, utilizing methodologies such as these becomes computationally infeasible when considering the numbers of simulated galaxies and snapshots modelled in this work. Radially averaged molecular gas fraction as a function of the galactocentric radius for three model SMGs (varying only the orientation angle of the discs) during three different phases. The three model SMGs produce 'average' SMGs with fluxes ranging from S 850 ≈ 2 to 7 mJy (Narayanan et al. 2009). The black lines are the pre-merger phase for each model, the blue lines denote the peak of the starburst (roughly the peak of the SMG phase), and the red lines represent the post-quasar phase. The molecular gas fraction is calculated via a dependency on the ambient pressure as motivated by observations of clouds in local galaxies (Blitz & Rosolowsky 2006). During the peak of the SMG phase (blue lines), the bulk of the neutral gas in the central regions is molecular, owing to rather high gas densities in the galactic nucleus.
which exist within cells of this volume (i.e. dense cloud cores and diffuse cloud atmospheres), we further model the molecular gas following the sub-grid prescriptions of Narayanan et al. (2008a). Specifically, the molecular gas within the ∼160 pc grid cells is assumed to reside in a mass spectrum of giant molecular clouds (GMCs) which are modelled as singular isothermal spheres. The mass spectrum of clouds follows a power law with index α = 1.8 as motivated by observations of local clouds (Blitz et al. 2007), though tests have shown that TURTLEBEACH results do not vary so long as the mass spectrum indices reside within observational constraints (Narayanan et al. 2008b). More details of the implementation of these sub-grid treatments of GMCs may be found in Narayanan et al. (2008a conservatively assume that Galactic abundance patterns hold and model the CO fractional abundance as CO/H 2 = 1.5 × 10 −4 . While molecular abundances in the ISM of high-redshift galaxies are unconstrained, recent measurements have shown that starforming ultraviolet (UV)-selected galaxies, likely progenitors of massive mergers, have solar abundances in their ISM (Shapley et al. 2004). As such, an assumption of Galactic molecular abundances may be reasonable.
We build the emergent spectrum by integrating the equation of radiative transfer through the H 2 gas: where I ν is the frequency-dependent intensity, S ν is the source function, r is the physical distance along the line of sight and τ is the optical depth. The source function is made up of a combination of the emission from dust as well as line emission. Formally, S ν = j ν /α ν where j ν = j ν (dust) + j ν (gas). Similarly, α ν has components from both the gas and dust. The dust radiates as a blackbody, and is assumed to be at the kinetic temperature of the molecular gas. We note that this element of the calculation is inconsistent with the SUNRISE calculations, which formally derives the dust temperature assuming that the dust and radiation field are in radiative equilibrium. Here, Weingartner & Draine (2001) opacities are assumed, though this makes little impact on the CO line flux.
The line source function is dependent on the CO level populations. Therefore, in order to self-consistently calculate the line intensities of CO, the molecular excitation properties must be known. The relevant physical processes in determining the CO excitation are collisions and radiative (de)excitation (e.g. line trapping; Narayanan et al. 2008a). The molecules are assumed to be in statistical equilibrium, and the population levels are calculated considering both the radiation field and collisional processes.
The methodology is an iterative one. To determine the solution to the molecular excitation, the level populations across the galaxy are first guessed at (in practice, we guess a solution near LTE). The molecular gas is then allowed to radiate model photons based on the assumed level populations, and, when a sufficient number of photons have been realized in each grid cell, the mean intensity field is calculated (Bernes 1979). Under the assumption of statistical equilibrium, the radiative excitation rates in combination with the collisional excitation rates give updated level populations, and new model photons are then emitted. This process is repeated until the level populations are converged.
The line transfer takes into account velocity fields. The 3D velocity field across the model galaxy is returned by the GADGET-3 hydrodynamic simulations. The line-of-sight velocity gradient from cell to cell is accounted for in the line transfer via emission and absorption line profiles, both of which are Gaussian in nature. The emission and absorption profiles have their widths determined by the thermal linewidth in the cell, as well as an assumed microturbulent velocity field (here, set at 0.8 km s −1 ). The difference in the frequency centres of the emission and absorption profiles is determined by the line-of-sight velocity difference between the emitting and absorbing clump (e.g. equations 5 and 6 of Narayanan et al. 2006b).
For the models presented here, ∼13 million model photons were emitted per iteration. The mass spectrum of GMCs is considered with a lower cut-off of 1 × 10 4 M and an upper limit of 1 × 10 6 M (consistent with constraints provided by local GMCs; Blitz et al. 2007). The CO excitation was solved for across 10 levels at a time, and the collisional rate coefficients were taken from the Leiden Atomic and Molecular Database (Schöier et al. 2005). The boundary conditions included the cosmic microwave background which was modelled at z = 2.2 to have a temperature of T = 8.74 K.
Finally, we comment that because of spatial resolution limitations in the molecular line radiative transfer, we are forced to consider 8 kpc boxes. Rather than following the centre of mass of the merging galaxy system (which, at times, may have scant little gas), we choose to follow a single progenitor galaxy through its evolution. Of course, as the nuclear discs of the progenitor galaxies overlap (during both first passage and coalescence), our models will include the emission from both galaxies in the pair.

Polychromatic SED radiative transfer
In order to identify when our simulated galaxies would be selected as submillimetre-luminous sources, we simulate the UV through submillimetre continuum photometry using the 3D adaptive grid polychromatic Monte Carlo radiative transfer code, SUNRISE (Jonsson 2006;Jonsson et al. 2006Jonsson et al. , 2009. SUNRISE calculates the transfer of UV through millimetre wave radiation through the dusty ISM. We refer the reader to these work for details on the underlying algorithms as well as numerical tests. Here, we summarize and explain the physical parameters employed in this study.
The radiative transfer is implemented via a Monte Carlo algorithm in which photon packets representing many real polychromatic photons propagate through the dusty ISM, and undergo scattering, absorption and re-emission. Model cameras are placed around the simulated galaxy to sample a range of viewing angles. The emergent flux is determined by the number of photons that escape the galaxy unhindered in a given camera's direction, as well as those scattered into or re-emitted by dust into the camera. For the purposes of the calculations presented here, we used eight cameras placed isotropically around the model galaxy.
SUNRISE is able to handle arbitrary geometries for the sources and dust. The input spectrum includes contributions from stars and black holes (AGN). The AGN input spectrum utilizes the luminositydependent templates of Hopkins, Richards & Hernquist (2007) of unobscured quasars. The normalization of the input spectrum is set by the total bolometric luminosity of the central black hole, L AGN = ηṀ BH c 2 . Again, η is assumed to be 10 per cent (cf. Section 2.1).
The stellar input spectrum is calculated utilizing the stellar populations code, STARBURST99 (Leitherer et al. 1999;Vázquez & Leitherer 2005), where the ages and metallicities of the stars are taken from the hydrodynamic simulations. We assume a Kroupa initial mass function (IMF) (Kroupa 2002), consistent with recent results from observations of z ∼ 2 star-forming galaxies (Davé 2008;Tacconi et al. 2008;van Dokkum 2008). The stellar particles initialized with the simulation are assumed to have formed over a constant star formation history. In order to match the SFR and stellar mass of the first snapshot of the simulation, this corresponds to a star formation history of ∼250 Myr.
The stellar clusters with ages less than 10 Myr are assumed to reside in their nascent birthclouds. SUNRISE models the effects of reddening of the stellar spectrum through these birthclouds utilizing results from the photoionization code MAPPINGSIII (Groves, Dopita & Sutherland 2004;Groves et al. 2008 (Groves et al. 2008). The effect of modelling the obscuration of stellar clusters by H II regions and PDRs is a redistribution of emergent UV light into the far-infrared (FIR)/ submillimetre bands. The H II regions absorb much of the ionizing UV flux and contribute heavily to the hydrogen line emission from the ISM as well as hot-dust emission. The time-averaged areal covering fraction by PDRs (f PDR ) is related to the PDR clearing time-scale as f PDR = exp(−t/t clear ), and is taken to be a free parameter (see Groves et al. 2008). These PDRs absorb much of the non-ionizing UV radiation field and contribute to the emergent polycyclic aromatic hydrocarbon (PAH) and FIR emission. Here, we assume an areal covering fraction of unity. A covering fraction f PDR = 1 translates to a cloud-clearing time-scale longer than the lifetimes of O and B stars (Groves et al. 2008), and the main consequence of reducing this clearing time-scale is to reduce the emergent submillimetre flux at the time of the starburst (Narayanan et al. 2009). The cloud-clearing time-scales in gas-rich galaxy mergers are unconstrained. That said, some constraints may be placed on this value either by educated guesses or by ansätze that are then verified by a comparison of the simulated SEDs to those observed.
The assumption of a clearing time-scale longer than the lifetimes of O and B stars may be a reasonable guess. The centres of local gas-rich major mergers are known to have large molecular volume filling fractions and are well characterized by a uniform, smooth molecular medium (Downes & Solomon 1998;Sakamoto et al. 1999) which may blanket nuclear O and B stars their entire lives.
While an assumption of f PDR = 1 is not the same as a uniform molecular medium, tests have shown that in this limiting case, the submillimetre SED is quite similar to f PDR = 1 (Narayanan et al. 2009).
Similarly, we can take the reverse approach, and assume the ansatz of f PDR = 1, and compare the simulated SEDs to those observed. In Narayanan et al. (2009), we showed that the mean SED in modelled SMGs with a range of masses compares quite well with observed SEDs of SMGs with spectroscopic redshifts. Moreover, these authors found that the peak submillimetre flux was related to the total mass of the galaxy. Utilizing PDR clearing time-scales longer than the lifetimes of O and B stars resulted in average (S 850 ∼ 5 mJy) SMGs with halo masses of the order of ∼5 × 10 12 -10 13 M , consistent with the inferred halo masses of observed SMGs (Swinbank et al. 2008). Similarly, the simulated stellar masses, black hole masses and H 2 masses both compared well with observations and scaled with peak submillimetre flux. Reducing t clear (or, equivalently, f PDR ) would require more massive galaxies to produce relatively average S 850 ∼ 5 mJy galaxies and quickly violate the inferred halo masses from clustering measurements. As such, the assumption of f PDR = 1 may be reasonable. That said, it is important to interpret the results presented in this paper in the context of this assumption for the birthclouds surrounding stellar clusters, especially in the context of mass-dependent results (such as the CO linewidths; Section 4.2).
The dust and radiation field are assumed to be in radiative equilibrium and utilize a method similar to that developed by Juvela (2005). Here, when photons are absorbed in a grid cell, and the dust temperature updated, a new photon with the SED equal to the difference between the SED emitted for the new dust temperature and that from the old one is emitted. This procedure is iterated upon until the radiation field has converged ). The SUNRISE calculations employ 10 million photon packets per iteration. Jonsson et al. (2009) find that the implementation of dusttemperature iteration in SUNRISE recovers the solution to the Pascucci et al. (2004) radiative transfer benchmarks to within a few per cent for UV through millimetre wavelengths. The Pascucci et al. (2004) benchmarks were recovered through the most stringent test cases of τ = 100. We note that the maximum optical depth we see during the SMG phase of these simulations is τ ∼ 75.
The gas and stars initialized with the simulation are assigned a metallicity according to a closed-box model such that Z = (−y ln [f gas ]), where Z is the metallicity, y is the yield=0.02 and f gas is the initial gas fraction (though note that the fluxes during the SMG phase of the model galaxies are not very sensitive to these assumptions). Because the dust properties of z ∼ 2 SMGs are relatively unconstrained, we utilize Galactic observations as input parameters. We assume a constant dust-to-gas ratio comparable to observations of local ULIRGs of 1/50. Tentative evidence suggests comparable dust-to-gas ratios in SMGs (Greve et al. 2005;Solomon & Vanden Bout 2005;Kovács et al. 2006;Tacconi et al. 2006). Simulations parametrizing the dust in terms of a Galactic dust-to-metals ratio of 0.4 (Dwek 1998) give similar results to within 10 per cent. The dust grain model used is the R = 3.1 Weingartner & Draine (2001) dust model including updates by Draine & Li (2007).
We are exploring the dust and photometric properties of z ∼ 2 galaxies in a number of companion papers at various stages of preparation. In this work, we focus on the CO properties of SMGs; as such, we primarily utilize SUNRISE as an informant as to when the model galaxies may have submillimetre luminosities comparable to observed SMGs. We classify our model galaxies as SMGs when they have an observer frame 850 μm flux ≥ 5 mJy which is comparable to a 3σ detection in current wide-field surveys (e.g. Coppin et al. 2006). We model our sources at z = 2.5, and therefore this fiducial criterion corresponds to a flux limit at a rest-frame wavelength of λ = 243 μm.

E VO L U T I O N O F S U B M I L L I M E T R E , B BA N D A N D H 2 P RO P E RT I E S O F S M G S
We begin with a general description of the evolution of the 850 μm, B band and H 2 properties of our model SMGs, and relate these to the evolutionary status of the galaxy merger. In Fig. 3, we present the 850 μm flux, CO linewidth (which will be discussed more in Section 4.2) and observed B-band apparent magnitude (AB magnitude system, modelled at z = 2.5, typical of SMGs; Chapman et al. 2003aChapman et al. , 2005 of our fiducial merger simulation, SMG10 (Table 1). The blue shaded region represents the sightline-dependent range of potential B-band magnitudes. The CO transition plotted is J = 3-2 in the rest frame, corresponding to millimetre-wave observations at z ∼ 2. As a reference for the global morphology, in Fig. 4, we plot the projected gas density at various snapshots for fiducial model SMG10, with the location of the black holes overlaid.
The initial passage of the galaxies induces a starburst of the order of ∼200 M yr −1 . The galaxies form stars at this rate for a few × 10 8 yr as they inspiral towards final coalescence. During this time, the galaxy builds a stellar mass of the order of ∼10 11 M (Narayanan et al. 2009). This elevated SFR allows the galaxy to produce ∼1-2 mJy of flux at 850 μm rendering this galaxy below the nominal 5 mJy detection threshold for SMGs, though detectable with future sensitive instruments.
When the galaxies approach towards final coalescence (T ≈ 0.6 Gyr), tidal torques from the merger drive large-scale inflows (Barnes & Hernquist 1991Mihos & Hernquist 1994. Physically, the interaction of the galaxies spins up the discs as they transfer angular momentum from the orbit of the galaxies to the discs themselves and triggers the growth of bars. Gas that shocks  (Table 1, designed to be an average SMG). All quantities are plotted at redshift z = 2.5. The grey hatched region in the top plot shows when the galaxy is visible as an SMG above S 850 > 5 mJy. The red shaded region in the CO linewidth plot shows the dispersion over 100 random sightlines and the blue shaded region in the B-band plot denotes the dispersion over eight cameras placed isotropically around the galaxy. As the galaxies spiral in towards coalescence, the observed submillimetre flux is relatively low (S 850 ≈ 1 mJy) and the progenitor galaxies are disclike. Consequently, the linewidths are representative of the virial velocity of the galaxies (here, set at V c = 320 km s −1 ). As the galaxies coalesce (T ≈ 0.6-0.65 Gyr), the ∼1000 M yr −1 starburst drives the 850 μm flux to detectable levels (S 850 > 5 mJy). Concomitantly, the B-band flux rises sharply, owing to both the intense starburst and a rapidly growing AGN. The CO line FWHM doubles as two discs enter the simulation box/observational beam (here, set at 8 kpc). In the post-starburst phase, the submillimetre and B-band flux drops, and the linewidths settle towards the rotational velocity of the combined (two-galaxy) system. See main text for more details. on this bar dissipates energy and loses angular momentum, causing an inflow to the central regions (see e.g. Hopkins et al. 2008Hopkins et al. , 2009). The high densities in the nuclear region of the coalesced system give rise to a massive starburst of the order of ∼1000-1300 M yr −1 . The nascent PDRs surrounding young stellar clusters convert the UV flux intercepted from O and B stars into longer wavelength submillimetre radiation. During this starburst event, the galaxy may be selected as a luminous submillimetre source with 850 μm fluxes approaching ∼5-7 mJy (Narayanan et al. 2009). The peak submillimetre flux observed is directly related to the mass of the merging galaxies -galaxy mergers above a total (halo) mass of ∼5 × 10 12 M will produce the nominal ∼5 mJy at 850 μm to be detectable as an SMG. Mergers of a significantly lower mass will have difficulty in producing a strong enough starburst to drive the observed submillimetre flux (Narayanan et al. 2009). This owes to the fact that the submillimetre flux in our model derives largely from the reprocessing of UV flux from the starburst in the birth clouds surrounding stellar clusters (Groves et al. 2008).
Concomitant to the final coalescence starburst, inflows fuel central black hole accretion. A fraction of the accreted mass energy is deposited into the ISM surrounding the central black hole(s), driving a pressure-driven wind. These AGN winds expel much of the obscuring gas and dust in the central regions, allowing the system to be viewed as an optically selected quasar (T ≈ 0.65-0.7 Gyr; e.g. Hopkins et al. 2005Hopkins et al. , 2006Springel et al. 2005, and references therein). The quasar phase and SMG phase are roughly coincident (though the quasar phase may lag the SMG phase by up to ∼20 Myr; Springel et al. 2005). This is in good agreement with observations by Coppin et al. (2008b), who find an overlap in a subset of their z ∼ 2 observed SMGs and quasars. As demonstrated by Fig. 3, however, there is a large sightline-dependent dispersion in potential B-band fluxes during the final merger. Consequently, the same galaxy can show almost 2 mag dispersion based on the observed viewing angle, and not all SMGs will appear as quasars. However, the galaxy is visible as an SMG at all modelled sightlines.
It is important to note that the selectability of our simulated galaxies as quasars is mass-dependent. As discussed by Narayanan et al. (2009), the final black hole mass of the merged system is dependent on the mass of the galaxy. This is similar to results found by Lidz et al. (2006) and Li et al. (2007) who found that quasar luminosity is tied to progenitor galaxy halo mass. Narayanan et al. (2009) found that the final black hole masses of average SMGs (e.g. S 850 ≈ 5-7 mJy) will be a few ×10 8 M , whereas the mergers which produce the most luminous SMGs (S 850 ≈ 15-20 mJy) may make black holes comparable to those seen in quasars (M BH ≈ 10 9 M ). We therefore continue the discussion of the CO properties of our SMGs during the 'quasar phase' as referring to the time period when the simulated dust-attenuated B-band flux peaks (e.g. third panel, Fig. 3), and note that this may not necessarily correspond to a galaxy selectable as a quasar in current surveys. The exact relationship between SMGs and quasars will be discussed in due course (Younger et al. in preparation; Narayanan et al. in preparation).
The starbursts induced by the merger consume significant amounts of gas. However, while gas consumption is high, supernova pressurization of the ISM sustains large molecular gas fractions. Thermal energy input into the hot phase ISM as well as mass increases owing to the stellar mass returned and evaporation of cold clouds increases the ambient pressure on cold clouds (Springel & Hernquist 2003). This increase in pressure increases molecular fractions in the neutral ISM, in accordance with our pressure-based H 2 formation/destruction algorithm (cf. equation 1). Hence, while star formation of course consumes H 2 gas, during starbursts, this effect may be mitigated owing to conversion of H I to H 2 . The H 2 masses during the final coalescence burst are rather high with typical masses of ∼5 × 10 10 M , though they can be as high as ∼3 × 10 11 M (Narayanan et al. 2009).

Molecular disc formation and disruption
The structure of the molecular gas in SMGs, and in particular the (potential) existence of molecular discs, is essential to a thorough  understanding of the evolution of their CO line properties. We briefly outline the key points related to this topic here, though defer a detailed investigation into the survivability of molecular discs in mergers to a future work.
In Fig. 5, we show the CO (J = 3-2) morphology of galaxy SMG10 as a function of time. The temporal evolution depicts the CO morphology as it evolves through inspiral (first row) and final coalescence (second and third rows). The galaxy is most likely to be viewed as an SMG during final coalescence when the SFRs are most elevated (e.g. fig. 1 of Narayanan et al. 2009). Here, this corresponds to T ≈ 0.65 Gyr. Because our CO calculations centre around a single galaxy at all points in its evolution to maximize spatial resolution, both galaxies only appear in the images when the nuclei are within the 8 kpc model box (second row and beyond).
In Fig. 6, we show the relative fraction of gas in rotational motion as a function of time, with the submillimetre flux overlaid as a reference. Gas is considered to be in rotational motion when its circular velocity exceeds 75 per cent of that expected at its radius for Keplerian motion. 5 As such, the relative fraction of gas in rotational motion may be viewed as a measure of the 'disciness' of the molecular gas. Comparisons between Figs 5 and 6 may be made via the time stamps displayed in the panels of Fig. 5.
The galaxies remain relatively disc-like as they inspiral towards final coalescence. Upon final merging, when the galaxy undergoes its luminous SMG phase, the disc is tidally disturbed. Tidal features 5 Varying this fiducial fraction of 75 per cent makes no difference on the relative temporal evolution of the fraction of disc-like gas; lowering or increasing this value simply increases or lowers the normalization of the curve. As such, throughout this paper, this 'fraction' should be taken as relative, and not absolute.  The simulations focus on a single galaxy through its evolution until both galaxies are within a single 8 kpc box (second row and beyond). For reference, the 8 kpc box employed for the CO radiative transfer simulations is shown explicitly with respect to the global morphology in Fig. 4. During inspiral (first row of this figure), the molecular gas is in a disc-like configuration. As the second galaxy enters the box, and they merge (second row), the disc-like morphology is disturbed and a large fraction of the gas is pulled from disc-like motion into relatively radial orbits. This corresponds to the peak of the SMG phase. During this time, extended features and tidal tails may become apparent (middle row). Tidal torquing drives much of the gas towards the central regions, resulting in a relatively concentrated molecular gas spatial extent (third row).
, gas is considered to be in rotational motion when it has a circular velocity of at least 75 per cent of the expected Keplerian velocity at its radius. During the peak of the SMG phase, most of the gas is disturbed from the central rotationally supported disc.
(and, on occasion, double-nuclei) become apparent in the molecular gas morphology during this final-coalescence SMG phase (e.g. T ≈ 0.65 Gyr, Fig. 5). This can be seen more explicitly in Fig. 7, where we plot the CO (J = 3-2) centroid velocity maps of the same snapshots shown in Fig. 5. Soon after the final interaction/SMG phase, the gas yet again re-virializes and a strong (compact) molecular disc re-forms (this is seen in Fig. 6, though the snapshots in Fig. 7 do not extend far enough in time to show this phase). This history of molecular disc formation/disruption throughout the galaxy merger's history is reminiscent of that seen in models of the molecular ISM in z ∼ 6 quasars (Narayanan et al. 2008c), and will play an important role in our understanding of the CO linewidths and usage of CO as a dynamical mass indicator in the sections to come.

Spatial extent of CO emission
With recent advances in (sub)mm instrumentation, high spatial resolution images of the molecular gas in SMGs are beginning to become available (e.g. Tacconi Fig. 5. The inspiral phase (top row) is characterized by ordered disc-like motion. The molecular discs are rapidly destroyed as the galaxies coalesce during the SMG phase (T = 0.63-0.65 Gyr). Note that the disc angle slowly changes throughout the early part of the galaxy's evolution, thus changing the magnitude of the line-of-sight velocities seen. spatial extent of the CO emission from SMGs can be important for determining the dynamical mass from CO linewidths.
During the merger-induced starburst, tidal torquing drives cold gas from the outer disc into the nuclear regions of the galaxy. Consequently, the CO emission becomes rather compact. In Fig. 8, we plot the distribution of characteristic CO (J = 3-2) radii from model galaxies SMG1 and SMG10 over 100 random sightlines during snapshots where the system may be viewed as an SMG (S 850 > 5 mJy). The characteristic radius for CO emission is the flux-weighted standard deviation in the radial distribution of fluxes. Practically, the map is treated as a histogram of fluxes at varying radii. The standard deviation in this histogram is the characteristic radius.
During the SMG phase, the emission is relatively compact owing to the gas funnelled into the central regions from the merger. The characteristic radius averages at ∼1.5 kpc for most sightlines, with some dispersion owing to the evolutionary status and sightline. Some large CO radii are seen from inspiralling discs in massive (M DM ≈ 10 13 M ) mergers which are already SMGs during the inspiral phase (Narayanan et al. 2009).
The CO spatial extents modelled here are comparable to the measurements by Tacconi et al. (2006) of an ∼4 kpc FWHM diameter in SMGs (which corresponds to an ∼0.85 kpc standard deviation in radius if the emission is Gaussian in nature). The distribution of radii additionally is consistent with the σ ≈ 0.6-3.3 kpc range of spatial extents observed from SMGs (Downes & Solomon 2003;Genzel et al. 2003;Tacconi et al. 2006Tacconi et al. , 2008Younger et al. 2008;Iono et al. 2009).

Model results: evolution of CO linewidths
The observed molecular linewidths of SMGs are exceptionally broad, with a median FWHM of ∼600-800 km s −1 , and linewidths exceeding 1000 km s −1 (Greve et al. 2005;Carilli & Wang 2006;Tacconi et al. 2006;Coppin et al. 2008b;Iono et al. 2009). Here, we explore the evolution of CO linewidths in our model SMGs; we show how our merger-driven formalism for SMG formation and evolution may self-consistently explain the observed broad lines from SMGs. Because the CO emission line is essentially a distribution measuring the power at a range of molecular gas line-of-sight velocities, rather than employing any particular fitting methodology (and thus espousing the associated uncertainties), we treat the  8. Characteristic radius of CO (J = 3-2) emission during the SMG phase of fiducial models SMG1 and SMG10. These models were chosen to bracket the range of masses which produce SMGs [ranging from average (5-7 mJy) to luminous (∼15 mJy); Narayanan et al. 2009]. The characteristic radius is calculated using the standard deviation in the flux-weighted radial distribution of fluxes. The distribution is plotted over 100 random sightlines. The characteristic radii match up well with observed size scale FWHM measured by Downes & Solomon (2003), Genzel et al. (2003), Tacconi et al. (2006Tacconi et al. ( , 2008 and Younger et al. (2008). Please note, however, that because the simulations here are performed in isolation, rather than drawn from cosmological conditions, this distribution represents the range of expected CO (J = 3-2) spatial extents from observed SMGs, rather than a true distribution.
line as a distribution of fluxes and utilize the standard deviation of the distribution (σ ) as a measure of the linewidth. For a perfectly Gaussian line, the FWHM of the line would simply be ∼2.35 × σ . While the lines are not perfectly Gaussian, to better compare with observations, we utilize this conversion between σ and FWHM.
In short, the CO linewidths are representative of the dynamics of the system (Narayanan et al. 2008c). Prior to final coalescence, typically only a single galaxy is within the beam, and the linewidths are narrow, representing the rotational velocity of a single galaxy. During the SMG phase, when the galaxies merge, the linewidths roughly double owing to the contribution of emission from both discs. In the post-SMG phase, as the molecular gas relaxes into a new disc (Fig. 6), the linewidths drop by a factor of √ 2. In more detail, the CO linewidths reflect the dynamics of the molecular gas in the galaxy. This owes to the origin of the CO emission lines from the model galaxy. While emission within a given cell (or neighbouring cells) is typically optically thick, globally, the emission is optically thin. When the CO line emitted from a cell escapes the local region, it typically leaves the galaxy unhindered as steep velocity gradients in the throes of the merger shift the absorption profile out of resonance with the emission profile. The CO emission line, then, is essentially the sum of the emission that escapes individual cells containing GMCs at their individual velocities. As such, the linewidths are reflective of the dynamics of the system. The dispersion (σ ) in the CO linewidth corresponds well with the velocity dispersion of the gas along the line of sight.
As the galaxies inspiral towards final coalescence, the molecular gas in the progenitor galaxies of the SMG is relatively virialized ), and thus the linewidths are reflective of the circular velocity of the progenitor discs. As discussed in Table 1 and Narayanan et al. (2009), our model mergers which produce average SMGs (S 850 ≈ 5-7 mJy) were typically initialized with discs with V c = 320 km s −1 , thus resulting in a CO linewidth of σ ≈ 160 km s −1 (which is equivalent to the 320 km s −1 circular velocity in the disc modulated by sin(30 • ) to account for the average inclination angle of the disc). This corresponds to an FWHM of ∼375 km s −1 for a Gaussian velocity dispersion (Fig. 3). Recall that during the inspiral phase, typically only one galaxy is within the 8 kpc beam, which is comparable to most interferometric beam sizes (Greve et al. 2005;Tacconi et al. 2006).
When the galaxies coalesce, both galaxies contribute to the detected line. The velocity dispersion of gas-detected doubles and, consequently, the linewidth roughly doubles. This is coincident with the SMG phase. Here, this causes the linewidths to increase to 2 × V c . Because the FWHM ≈ 2.35 × V c , and accounting for the average inclination angle of the discs, we then arrive at a generalized expression for the CO line FWHM during the final coalescence SMG phase: which is, of course, simply FWHM ≈ 2.35 × V c when an average inclination angle of i = 30 • is assumed. Here, where circular velocities of V c = 320 km s −1 are employed for the initialization of the progenitor galaxies, this results in modelled linewidths of ∼600-800 km s −1 . The dispersion in the modelled CO linewidths owes to both the viewing angle and evolutionary effects. Generically, the SMG phase and quasar phase occur at or around the nuclear coalescence of the galaxies. However, the exact timing of this event is dependent on a number of things, including galaxy orbit, gas content and mass, amongst other factors. As such, the CO linewidths from the SMGs or quasar host galaxies with a merger origin show a large range of linewidths. In Fig. 9, we plot the histogram of CO linewidths from our fiducial galaxy (SMG10) during both its SMG phase and its quasar phase. As before, we nominally define the SMG phase as when the 850 μm flux is 5 mJy and arbitrarily define the quasar phase as when the galaxy is brighter than the 27th magnitude (AB magnitude; apparent magnitude at z = 2.5). Because the SMG phase and quasar phase are nearly coincident, the linewidths are generally broad throughout both phases, with nearly indistinguishable distributions.

Observational comparisons
The distribution of modelled linewidths for SMGs shown in Fig. 9 is in excellent agreement with those presented by Coppin et al. (2008b) and Carilli & Wang (2006). The linewidths are broadly reflective of the mass of the simulated SMGs, and likely signifies a correspondence between the final masses and evolutionary status of our modelled SMGs and those in nature. In this sense, our model faithfully provides a natural explanation for the broad linewidths observed in SMGs.
A comparison between the linewidths of our model quasars and those in nature is more difficult. Coppin et al. (2008b) presented new CO detections of submillimetre-luminous quasars, and found a nearly indistinguishable distribution of CO linewidths from quasars and SMGs, in excellent agreement with the simulations presented here. In both our models, and the observations of Coppin et al. (2008b), the median CO linewidth is ∼600-800 km s −1 , which we view as a general success of our model. However, a literature compilation of linewidths from CO-detected quasars by Carilli & Wang (2006) found a much narrower median in the distribution, with the linewidths showing a median value of FWHM ≈ 300 km s −1 , in contrast to both the observational results of Coppin et al. (2008b) as well as the model results presented here. The source of this discrepancy is not entirely clear. Because the analysis performed by Carilli & Wang (2006)  . Distribution of modelled CO (J = 3-2) linewidths during the SMG and 'quasar' phases of model SMG10. The SMG phase is considered when S 850 > 5 mJy, and the quasar phase is arbitrarily assigned with an apparent B-band magnitude (AB system) cut of 27th magnitude. See Fig. 3 for more details. The phases occur at similar time periods (typically separated by at most ∼20 Myr), and consequently have similar linewidth distributions. The FWHM distribution for these two sources appears to correspond well with the distributions published by Coppin et al. (2008b), though the quasar distribution appears to be discrepant with those published by Carilli & Wang (2006). See text for details regarding the origin of the broad CO linewidths in SMGs and quasars. Again, because the simulations here are not cosmological in nature, the distribution should be viewed as representative of the range of CO linewidths from average SMGs, not the true distribution from a blind survey. While we normally include higher mass models (e.g. SMG1) in our analysis, we refrain here as the extremely broad (FWHM ≈ 1000-1400 km s −1 ) lines from SMG1 take away from the main point of modelled SMGs with average S 850 fluxes producing average CO linewidths. SMGs as massive as SMG1 will be rare (indeed, this is apparent by the relative lack of observed S 850 > 15 mJy SMGs) though, when observed, they will have linewidths FWHM > 1000 km s −1 on average. sample utilized a variety of CO transitions from a non-uniform sample with objects spanning a large range of redshifts from z ∼ 1 to 6. Different CO lines trace different spatial extents, and thus may show a range of linewidths. Similarly, quasars even of a mass similar to those presented in Coppin et al. (2008b) but at a lower redshift would naturally show evolution in their linewidths owing to shallower gravitational potentials. In contrast, the observations by Coppin et al. utilized only CO (J = 2-1) and (J = 3-2) emissions from a sample of quasars in a smaller redshift range (z = 1.7-2.6). In this sense, the comparison between our models and the relatively uniform observations of z ∼ 2 quasars by Coppin et al. seems most appropriate as it best compares to the redshifts and emission lines investigated in this work.

T H E U S AG E O F C O A S A DY NA M I C A L M A S S T R AC E R I N S M G S
In Section 4.2, we saw that during the SMG phase of the model galaxy's evolution, the CO linewidths were larger than expected from ordered disc-like motion gas. In this context, it is interesting to quantify the usage of CO as a dynamical mass tracer in SMGs.
The exact value of CO linewidths in high-redshift mergers as a dynamical mass indicator depends, of course, on the relationship between the linewidth and dynamical mass assumed. Broadly, for disc-like motion, this relationship can be expressed in the form  Figure 10. Sightline-averaged distribution of M dyn /M actual for models SMG1 and SMG10. These models were chosen to bracket the range of masses of galaxies which appear to produce SMGs in our simulations. The dynamical masses are calculated as in equation (5) for the simulated CO (J = 3-2) emission. Generally, using CO linewidths from SMGs for a dynamical mass calculation requires a modulation of a factor of ∼1.5-2 for translating to a true enclosed mass. See text for details.
where σ is the 1D velocity dispersion in the line, i the disc inclination angle, R the spatial extent of the CO emission (here, taken to be the half-width at half-maximum) and k a constant encompassing the relationship between V c and V FWHM (for a given distribution of mass) and R g and R hwhm . Equation (6) holds due to the global optical thinness of molecular line emission across a galaxy.
However, uncertainties exist in all of these conversion factors. Furthermore, it is not clear what fraction of the molecular disc survives during the merger, and whether the inclusion of a sin 2 (i) term is appropriate. Even if it were, prior to the high-resolution imaging capable only by Atacama Large Millimetre Array (ALMA), we can at best assume an average disc angle of i = 30 • . As such, a constructive method for quantifying the usage of CO as a dynamical mass tracer in SMGs is to say and characterize the relationship between V 2 FWHM R HWHM /G and M actual . Here, note that M actual is the total (baryonic and dark matter) mass enclosed in R HWHM and V FWHM is calculated as 2.35× the standard deviation in the linewidth (cf. Section 4.2). In Fig. 10, we plot the ratio of the dynamical mass as calculated in equation (6) to the true enclosed mass (which includes dark matter and baryonic mass). The radius employed is the half-width at half-maximum from the simulated CO emission maps. The ratio of M dyn /M actual is calculated for SMGs which span the range of masses explored here, for all snapshots which satisfy the fiducial criteria S 850 > 5 mJy, and for 100 random sightlines. On average, the ratio of the dynamical mass inferred from equation (6) to the true enclosed mass in the half-mass radius, k , ranges 6 from 1.5 to 1.9. With this, it is then feasible to modulate equation (6) with this given value of k to construct the appropriate relationship between the CO-derived dynamical mass and true mass enclosed in the CO-emitting region.

C O E X C I TAT I O N A N D L I N E S P E C T R A L E N E R G Y D I S T R I B U T I O N
The excitation of CO is an important diagnostic of high-redshift galaxies. First, molecular line measurements from SMGs are typically made in millimetre bands, which correspond to high excitation CO lines in the rest frame at z ∼ 2. Because the inferred molecular gas mass is generally derived via converting CO (J = 1-0) velocity-integrated intensity, understanding the typical CO line ratios (and the average level of thermalization of the gas mass within the telescope beam) is crucial. Assumptions of thermalized line ratios (e.g. LTE, when brightness temperature ratios between levels are unity) between higher lying lines and CO (J = 1-0) may underestimate the molecular gas mass in the event of substantial quantities of subthermal gas. Second, the CO excitation patterns reveal the line(s) of dominant CO power output. As more broad-band molecular line spectrometers become available (e.g. the ZEUS, Zspec, Zpectrometer spectrometers; Bradford et al. 2004;Harris et al. 2007;Hailey-Dunsheath et al. 2008), CO detections of high-redshift galaxies will be critically dependent on observations of the brightest lines. In an effort to aid interpretation of existing data sets, and help guide future observations of SMGs, we investigate the molecular excitation properties and CO line ratios from our model SMGs. We do not attempt to address the longstanding problem of converting the CO (J = 1-0) flux to the H 2 gas mass in galaxies, but rather simply relate the flux density from higher excitation lines to that from the ground rotational transition.

Model results: highly excited CO in SMGs
In Fig. 11, we plot the sightline-averaged model CO SED for models SMG1 and SMG10 (Table 1) which bracket the mass range of galaxies that satisfy the nominal selection criteria S 850 > 5 mJy. The ordinate is the ratio of the velocity-integrated intensity from various transitions compared to CO (J = 1-0), normalized to CO (J = 1-0). The range of emergent CO SEDs from our models is denoted in the grey shaded region and the mean by the solid line. The dashed line shows the expected CO SED for thermalized level populations. The line ratios are modelled for unresolved observations, and include emission from the entire 8 kpc simulation box. To best compare with the few published constraints (next section), we plot the line ratios relative to the ground (J = 1-0) transition. This has the added benefit of providing a direct measure of the ability of higher lying transitions to convert to H 2 gas masses.
As Fig. 11 shows, there is a broad dispersion in potential CO SEDs from our simulated SMG. However, two generic features are evident. First, on average, the CO is quite excited with the SED turnover only occurring at the J ≈ 5-6 level. This is in reasonable agreement with observations of z ∼ 2 SMGs (see the next section). The higher excitation CO SEDs come from SMGs which arise during final coalescence mergers, while the lower excitation SEDs represent the inspiral phase of massive galaxies (e.g. SMG1) which may be seen as SMGs even prior to final coalescence (e.g. fig. 1 of Narayanan et al. 2009).
Second, while the turnover is typically at a relatively high J level, most CO levels are subthermal over the 8 kpc simulated box. The rotational ladder in Fig. 11 shows that the gas is nearly thermalized for the lowest lying lines (CO J = 2-1). Higher lying lines, however, are not thermalized over the 8 kpc model box. While the gas associated with the nuclear starburst in the central ∼1 kpc is warm, dense and thermalized, the outer regions (R ∼ 1-4 kpc) contain a significant quantity of lower density gas which contributes to the SMM 13120 ERO J164502 Figure 11. Model CO SED from fiducial models SMG1 and SMG10. The SED is presented for all snapshots where the galaxy would be detected as an SMG (S 850 > 5 mJy) and the solid line shows the mean SED while the grey shaded region denotes the 1σ dispersion amongst snapshots. The flux densities from J > 1 levels are compared against the CO (J = 1-0) level. The dashed line represents the predicted CO SED for thermalized level populations (LTE). The green diamond shows observational data from SMG SMM 13120+4242 (Hainline et al. 2006) and the red triangles from SMG ERO J164502 (aka HR 10; Andreani et al. 2000;Greve et al. 2003;Papadopoulos & Ivison 2002). To our knowledge, these represent the only SMGs with tabulated multiline data [including a CO (J = 1-0) detection] which are unlensed. The highest excitation CO SEDs come from final coalescence mergers, whereas massive SMGs that are above S 850 > 5 mJy during inspiral may have lower CO excitation. The mean SED displays reasonable agreement with the observations, and suggest that CO levels J 2 may not be thermalized. This shows that caution must be exercised when assuming a brightness temperature ratio of unity to the CO (J = 1-0) line when deriving molecular gas masses. emergent spectrum via line trapping. This lowers the mean excitation conditions observed. As such, Fig. 11 demonstrates that most observations of SMGs which do not focus solely on the nucleus will reveal subthermal CO emission at higher lying (J 3) levels. 7 This is effectively saying that the filling factor of dense gas is relatively low, thus lowering the mean observed CO SED. {Consequently, assumptions of thermalized CO line ratios from unresolved observations of SMGs will typically underpredict the molecular gas mass}.

Observational comparisons
While few tabulated observational constraints exist for relative intensities from SMGs, the excitation seen in Fig. 11 seems to display reasonable agreement with the existing published data. Line ratios for ERO J164502 and ERO J164502 are given by Andreani et al. (2000), Papadopoulos & Ivison (2002), Greve, Ivison & Papadopoulos (2003) and Hainline et al. (2006). These data are represented by the points with error bars in Fig. 11. While the CO SED is normalized to the ground state (thus making this point an unconstraining comparison), the higher excitation detection of two 7 Higher spatial resolution observations probing the nuclei of SMGs should show higher excitation conditions. This lends itself to a direct testable prediction from these models which we outline in Section 7. lines (J = 2-1 and J = 4-3) is entirely consistent with the predictions made here. The intensity from the J = 5-4 line from ERO J164502 is lower than what the models predict here, though we note that this is the lowest excitation SMG known (Papadopoulos & Ivison 2002;Wei et al. 2007). A detailed survey and compilation of literature data has been presented by Wei et al. (2007). While these data are not tabulated, we can make qualitative comparisons. The bulk of these galaxies shows subthermal emission in the higher lying lines with a turnover at the J = 5 or 6 level. Thus, at face value, these observed CO SEDs are comparable with the predictions made in Fig. 11.
A detailed examination of the Weiß et al. observed CO SEDs, however, suggests a potential discrepancy between the observed values and those modelled here. The observed intensity at the J = 5 or 6 level (where the SED turns over) is a factor of ∼15 above that from the ground (J = 1-0) transition; this is in contrast to our models which suggest that the most excited lines will be only a factor of ∼5 greater than the ground state (Fig. 11). This difference may be reconciled with better constraints on the true CO (J = 1-0) emission from SMGs. Recall that the CO SED is typically plotted (both here and in the observational literature) as relative to the CO (J = 1-0) transition. For the bulk of the observed sources in the literature, the J = 1-0 line has not directly observed, but rather inferred from large velocity gradient (LVG) modelling which is a simplified approximate method compared to the 3D non-LTE radiative transfer methods employed here. It may be that the apparent difference in excitation between the observations and theoretical predictions will be reconciled via future detections of CO (J = 1-0) from SMGs.

Testable predictions
In this paper, we have outlined a model for the CO emission from high-redshift SMGs. Our models have provided a natural explanation for the observed H 2 gas masses, CO spatial extents, linewidths and excitation conditions. We have attempted to compare with literature data when available, and generally found a reasonable correspondence between our model results and galaxies in nature.
This all said, while our models appear to provide a plausible match to observed data, the comparisons made thus far are in essence postdictions. An even more powerful test of any theoretical model would be testable predictions of future surveys. In this section, we sketch out potential observable tests of these models which are imminently possible with the latest generation of bolometer and heterodyne receiver technology.

CO linewidths of z ∼ 2 SMGs and quasars
Our model for the increased linewidths from SMGs relies on a temporary increase in velocity dispersions owing to multiple galaxies in the simulation box/telescope beam. Two features of this model are evident. First, the linewidths are relatively low, and representative of the virial velocity of a single progenitor galaxy during the inspiral phase when the observed 850 μm flux may be relatively low (S 850 ≈ 1 mJy; Fig. 3). They increase concomitant to the increase in the submillimetre flux, and are thus broader when the galaxy is in its transient SMG phase, though remain broad in the post-SMG phase when the galaxies of 850 μm flux decrease again. In this picture, one might expect a broad range of linewidths from galaxies with lower (∼1 mJy) 850 μm fluxes: small linewidths corresponding to inspiralling galaxies and large linewidths corresponding to post-SMG phase galaxies. Similarly, lower mass mergers (e.g. model SMG13) at final coalescence may contribute to the dispersion of linewidths on the low flux end.
Second, as discussed by Narayanan et al. (2009), the 850 μm light curve from merging galaxies has a similar shape for a mass sequence of mergers, though scales in normalization. This means that the submillimetre flux curve shown in the top panel of Fig. 3 simply scales upwards with increasing merger mass (indeed this is evident with a direct comparison of the light curve presented in Fig. 3 with that of fig. 1 in Narayanan et al. 2009). This means that in more massive mergers, the inspiral phase corresponds to 5 mJy sources and the peak SMG phase to rarer, more luminous (S 850 10-15 mJy) SMGs. Consequently, there may be a spread in CO linewidths from even 5 mJy SMGs. The most luminous sources, however, only occur at the final coalescence burst of extremely rare ∼10 13 M mergers. As such, S 850 ≈ 10-20 mJy galaxies will have extremely broad (FWHM > 1000 km s −1 ) CO linewidths 8 with a relatively small dispersion.
The aforementioned effects will have the following generic consequence for CO linewidths from high-redshift mergers: the mean CO FWHM will increase with an increasing submillimetre flux and the dispersion will decrease. We show this explicitly in Fig. 12 by plotting the sightline-averaged model CO FWHM versus 850 μm flux for model galaxies SMG1, SMG10 and SMG13 9 in Table 1 for a flux limit of S 850 > 1 mJy, and in the middle panel we show the CO FWHM-1.1 mm relation. This serves as a prediction for both future sensitive 850 μm surveys and 1.1 mm (e.g. AzTEC; Wilson et al. 2008) counts.
In the right-hand panel of Fig. 12, we attempt to compare our model prediction in the left-hand panel with data culled from the recent review by Solomon & Vanden Bout (2005). An important point is that optically selected quasars cannot be included in this sort of comparison as they may preferentially have their molecular discs in a face-on configuration, thus skewing average linewidths (e.g. Carilli & Wang 2006;Narayanan et al. 2008c). This does exclude a number of sources which may be simultaneously quasars and submillimetre bright (e.g. Fig. 3 and Coppin et al. 2008b). The observed data in the right-hand panel of Fig. 12 are generally inconclusive in comparison to the predictions in the left-hand panel. The effects of uncertain gravitational lens magnifications muddy interpretation (sources that are known to be lensed are marked as red squares; Solomon & Vanden Bout 2005). Better statistics with SCUBA2 and AzTEC will provide a direct and imminent test for this aspect of these models.

Predicted CO excitation for high-resolution observations of SMGs
With the advent of broad-band receivers on a multitude of telescopes, CO SEDs of high-redshift objects will soon become available and serve as a test for the excitation predictions outlined in Section 6.
We first note that Fig. 11 itself is a testable prediction. As few multiline surveys are currently published, Fig. 11 serves as a direct comparative for future constraints on the CO SED from SMGs. A 8 Recall that the CO linewidth scales with circular velocity and, consequently, with galaxy mass. 9 These models were chosen to span a broad mass range. Note that model SMG13 is too low in mass to form a detectable S 850 > 5 mJy SMG (Narayanan et al. 2009 Solomon & Vanden Bout (2005). The models are for models SMG1, SMG10 and SMG13. SMG1 and SMG10 were chosen as usual to bracket the range of halo masses employed here and SMG13 added to probe lower fluxes (though note that it is too low in mass to ever form an S 850 > 5 mJy SMG). The brightest SMGs at 850 μm and 1.1 mm are predicted to have a narrow dispersion of CO FWHMs, and typically broad lines whereas lower luminosity objects may have a larger range of linewidths. The mean CO linewidth will increase with an increasing (sub)mm flux. See text for reasoning. From the observational data set, optical quasars have been eliminated. The comparison to the observations is generally inconclusive. This may owe to many uncertain gravitational lens magnifications in the observed SMGs (known lensed sources are marked as red squares; Solomon & Vanden Bout 2005). Because these models probe both the full dynamic range of simulated 850 μm fluxes and CO FWHMs, adding more models does not increase the dispersion in the models plots. crucial component to this will be the direct detection of CO (J = 1-0) emission from SMGs. In the absence of this, relative intensities to this transition will be difficult to interpret, and direct conversion of emission from higher lying transitions to H 2 masses will become unclear. This is currently feasible at a variety of telescope facilities, though interferometers (e.g. the EVLA) may be preferable as baselines may be problematic on single-dish telescopes (Hainline et al. 2006).
Second, we may make predictions for imminent higher spatial resolution observations (e.g. ALMA). The CO excitation conditions presented in Fig. 11 for our model SMGs were averaged over the entire 8 kpc box, thus including emission from subthermally excited gas. While the gas associated with the nuclear starburst is indeed warm and thermalized, including diffuse, subthermal gas had the effect of lowering the mean excitation condition. This is analogous to observations of z ∼ 2 SMGs which include both dense, nuclear gas and diffuse gas in the telescope beam. Consequently, the modelled excitation conditions in Fig. 11 showed excellent agreement with multiline CO measurements from SMGs.
In principle, one can imagine that higher spatial resolution observations of SMGs which focus on the dense, nuclear regions would show higher excitation CO SEDs. Our approach allows us to quantify this effect, as well as make predictions for the next generation of high spatial resolution interferometers (e.g. ALMA). In Fig. 13, we plot the sightline-averaged CO SED for SMG models SMG1 and SMG10 as a function of a decreasing physical beam size. As before, we only consider SMGs with a fiducial selection criteria of S 850 5 mJy. As the observations probe narrower and narrower columns, the observed gas is on average more highly excited and the turnover point of the CO SED moves to an increasingly high Figure 13. Predicted CO line SEDs for models SMG1 and SMG10 during their SMG phases at increasing spatial resolution. The black solid line shows the predicted CO line SED for thermalized populations. As the warmer, denser gas towards the nucleus is probed, the peak and shape of the excitation ladder shift towards higher levels as more of the gas becomes thermalized. High-resolution observations by ALMA of the nucleus of SMGs will typically show higher excitation CO line SEDs than lower resolution observations. rotational number. The lack of complete thermalization in the lower observed transitions (e.g. CO J = 2-1) owes to the fact that diffuse gas is still folded into the observation, even for rather high spatial resolution observations. C 2009 The Authors. Journal compilation C 2009 RAS, MNRAS 400, [1919][1920][1921][1922][1923][1924][1925][1926][1927][1928][1929][1930][1931][1932][1933][1934][1935]

S U M M A RY
We have combined hydrodynamic simulations of SMG formation and evolution with 3D non-LTE molecular line radiative transfer calculations to provide a model for the CO emission properties from SMGs. Our model has shown a number of successes in matching the observed spatial extent of CO emission, CO linewidths and excitation conditions. We utilized these models to understand the origin of these emission properties.
(i) In our model, SMGs originate in major mergers (Narayanan et al. 2009). Strong gaseous inflows drive highly concentrated molecular gas complexes such that the observed characteristic CO radius is of the order of ∼1.5 kpc. The large radius tail of the distribution arises from pre-coalescence galaxies in extremely massive (∼10 13 M ) haloes which are SMGs even during the inspiral phase (e.g. fig. 1 of Narayanan et al. 2009).
(ii) The large CO linewidths from SMGs owe to the fact that they are typically being observed during a transient phase where the gas is highly non-virialized and multiple galaxies are in the simulation box/telescope beam. During interactions, the CO FWHM from mergers is roughly 2.35× the circular velocity of a progenitor (for a Gaussian line; equation 4). Two merging V c ≈ 320 km s −1 discs naturally produce average (S 850 ≈ 5-7 mJy) SMGs with extremely broad linewidths of the order of ∼600-800 km s −1 .
We have additionally been able to provide interpretation regarding the usage of CO as a diagnostic of physical conditions.
(i) The usage of CO linewidths from SMGs as a dynamical mass estimator may overestimate the enclosed mass. Typical overestimates are of the order of M dyn /M actual ≈ 1.5-2.
(ii) The CO excitation in SMGs is high, with the rotational ladder turning over at the ∼J = 5 or 6 level. The level populations at J upper > 2 are typically subthermal, and assumptions regarding brightness temperature ratios of unity with the ground state will lead to underestimates of the inferred H 2 gas mass. The CO (J = 3-2) line from SMGs is typically a factor of ∼3 below the intensity expected from thermalized level populations.
Finally, we have made predictions for this model which are imminently testable with the newest generation of bolometer arrays and wide-band sensitive CO receivers.
(i) The CO linewidths from galaxies which will become SMGs are predicted to be their broadest when the submillimetre flux is the highest. Consequently, the brightest SMGs (S 850 > 15 mJy) are predicted to only have quite broad (FWHM > 1000 km s −1 ) linewidths with a small dispersion in FWHMs observed. Lower flux galaxies (S 850 ≈ 1) can be either normal discs (where the linewidths are predicted to be relatively narrow) or post-SMG phase mergers (where the linewidths will be relatively broad). As such, lower flux galaxies are predicted to have a broad dispersion in observed CO FWHM. Care must be taken to remove both optically selected quasars (which may be biased to have face-on molecular discs) and lens-magnified sources from the sample.
(ii) The intensity from the CO J = 5-4 or 6-5 line, where the CO SED turns over, is predicted to be a factor of ∼5-7 times that of the ground state. While detections at these higher lying lines are becoming routine, future detections of CO (J = 1-0) emission from SMGs will be required to test the predicted rotational ladders.
(iii) The peak and shape of the CO SED will shift towards higher lying transitions as smaller spatial extents are investigated with future high-resolution arrays (e.g. ALMA).