EVIDENCE FOR A CLUMPY, ROTATING GAS DISK IN A SUBMILLIMETER GALAXY AT z = 4

J. A. Hodge; C. L. Carilli; F. Walter; W. J. G. de Blok; D. Riechers; E. Daddi; L. Lentati

doi:10.1088/0004-637X/760/1/11

1. INTRODUCTION

Submillimeter-luminous galaxies (SMGs; Blain et al. 2002) are dusty, gas-rich, high-z galaxies that were revealed in the first extragalactic surveys using SCUBA and MAMBO (e.g., Smail et al. 1997; Hughes et al. 1998; Barger et al. 1998; Eales et al. 1999; Blain et al. 1999; Bertoldi et al. 2000; Greve et al. 2004). Their huge far-infrared (FIR) luminosities (∼10¹³ L_☉) are believed to be primarily driven by intense (∼10³ M_☉ yr⁻¹) star formation (e.g., Alexander et al. 2003, 2005), adding to their already significant stellar masses (e.g., Borys et al. 2005; Hainline et al. 2011). Together with FIR luminous quasars, they are generally thought to evolve into the giant ellipticals we see in the local universe. They are therefore critical for understanding early-type massive galaxy formation.

SMGs are a relatively rare phenomenon, with typical space densities of 10⁻⁵ to 10⁻⁶ Mpc⁻³. Their rarity may be partly because their enhanced star formation rates (SFRs) are necessarily short-lived (<100 Myr; Greve et al. 2005), or else the galaxies would grow too large, too fast. The highest volume density of radio-selected SMGs occurs at z ∼ 2–3, indicating that they peak simultaneously with the peak epoch of star formation (Chapman et al. 2003a). However, a number of recently discovered higher-redshift SMGs (Schinnerer et al. 2008; Daddi et al. 2009a, 2009b; Riechers et al. 2010; Coppin et al. 2010; Capak et al. 2011; Wardlow et al. 2011; Cox et al. 2011) as well as several studies based on statistical arguments (Greve et al. 2008; Penner et al. 2011, and references therein) may indicate that a high-redshift tail does exist and is able to account for the population of old, massive ellipticals already in place at z ∼ 2–3 (Daddi et al. 2009b).

Many SMGs are believed to be starburst-dominated major mergers (e.g., Chapman et al. 2003b; Engel et al. 2010; Narayanan et al. 2010; Hayward et al. 2011, 2012). This would make SMGs the high-redshift analogs of ultraluminous infrared galaxies (ULIRGs) in the local universe. Indeed, there is direct evidence for multiple CO components and/or disturbed kinematics in some SMGs, supporting the merger picture (e.g., Tacconi et al. 2008; Engel et al. 2010; Ivison et al. 2011; Riechers et al. 2011b).

Recently, it has been suggested that other mechanisms that drive extreme SFRs may also be at play. In particular, a scenario has been put forward in which the star formation in massive, high-redshift galaxies is driven by cold mode accretion (CMA; e.g., Kereš et al. 2005; Dekel et al. 2009a, 2009b). CMA-driven galaxies are constantly forming stars at high rates, and the star formation is sustained by smooth infall and accretion of gas-rich material. This process can result in elliptical galaxies because the streams that feed star formation can cause the disk to break up into giant clumps if they have a high enough gas fraction and degree of turbulence. The clumps then potentially migrate inward and merge into a spheroid (Dekel et al. 2009a, 2009b).

The CMA phenomenon has been extended to SMGs by a number of authors (Fardal et al. 2001; Finlator et al. 2006; Davé et al. 2010; Carilli et al. 2010). In this scenario, SMGs are massive galaxies sitting at the centers of deep potential wells and fed by smooth accretion. They can be thought of as super-sized versions of normal star-forming galaxies (SFGs), providing an alternate formation scenario to the merger-induced model. By identifying simulated SMGs as the most rapidly star-forming systems that match the number densities of SMGs, this theory has been successful at explaining some key SMG properties, including their stellar masses and clustering scales (Davé et al. 2010). Observations of the rest-frame optical morphologies of SMGs provide further evidence that they may simply be the most extreme examples of normal star forming galaxies in that era (Targett et al. 2011, 2012).

According to Shapiro et al. (2008), the best way to distinguish between CMA and a gas-rich merger model is through observations of the gas dynamics and distribution. A well-defined disk with a smooth rotation curve would be indicative of CMA, while tidally disturbed gas, with a very high T_B starburst nucleus, would point toward a merger. This is simply due to the redistribution of angular momentum. The large-scale gravitational torques induced by gas-rich major mergers are efficient at removing angular momentum (Barnes & Hernquist 1996), thereby funneling the cold molecular gas into the galaxy's center and producing a nuclear starburst. Indeed, in ULIRGs—the canonical merger scenario—the gas is concentrated in the central kpc (e.g., Downes & Solomon 1998; Bryant & Scoville 1999). The gas-rich, star-forming disks in Arp220 have a size of only ∼100 pc (Sakamoto et al. 2008). In simulations where SMGs result from mergers, their gas reservoirs show slightly larger extents of a few kiloparsecs (Narayanan et al. 2009).

Conversely, the z ∼ 2 massive SFGs which have been proposed to be due to smooth accretion are large disk galaxies, with star formation occurring at large radii or situated in rings (Elmegreen & Elmegreen 2006; Genzel et al. 2008; Elmegreen et al. 2009). Rest-frame UV/optical imaging shows that, unlike low-redshift galaxies, their disks tend to be broken into multiple giant clumps of ∼1 kpc and 10⁹ M_☉ (Elmegreen et al. 2004; Elmegreen & Elmegreen 2005; Förster Schreiber et al. 2006; Genzel et al. 2008). More recent work has also imaged such galaxies in millimeter/CO emission (Tacconi et al. 2010; Swinbank et al. 2010, 2011), finding evidence for clumpy CO emission extending over the disk (e.g., Tacconi et al. 2010).

Observations of the morphology and kinematics of molecular gas in SMGs may therefore shed light on the physics behind the intense star formation. Previously, Carilli et al. (2010) presented Karl G. Jansky Very Large Array (VLA) observations of CO(1–0) and CO(2–1) emission in the z = 4.05 SMG GN20, the brightest SMG in the GOODS-N field (Pope et al. 2006), and one of three SMGs in what appears to be a massive z ∼ 4 protocluster (Daddi et al. 2009b). These low-J transitions are of particular interest because they trace the cold molecular gas thought to make up the bulk of the molecular gas in these systems. Although their high-resolution CO(2–1) imaging showed evidence that the gas in GN20 was well resolved, their observations suffered from some severe spectral limitations. In particular, the observations utilized the old VLA correlator, with a total bandwidth of only 100 MHz (650 km s⁻¹), causing the line profiles to be truncated on both sides. The observations were also done in continuum mode, resulting in no information on kinematics.

We therefore obtained over 120 hr of time on VLA (Perley et al. 2011) to image the CO(2–1) emission in the GN20 field in the B- and D-configurations. The ∼20 hr of lower-resolution D-array data were presented by Carilli et al. (2011). Here, we present the full data set (B+D-configurations) on GN20, providing greatly improved spatial resolution and image fidelity.

We begin in Section 2 by describing our new VLA observations and data reduction of GN20. CO maps and the derived gas mass are presented in Section 3. Section 4 constitutes our analysis, including a dynamical analysis (Section 4.1), the definition and properties of individual gas clumps (Section 4.2), a multiwavelength comparison (Section 4.3), and an analysis of the spatially resolved gas excitation (Section 4.4). We discuss the implications of our results on the nature of GN20 in Section 5, and we end with a summary of our conclusions in Section 6. Where applicable we assume the standard Λ cosmology of H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_Λ = 0.7, and Ω_M = 0.3 (Spergel et al. 2003, 2007). At a redshift of z = 4.05, 1'' corresponds to 7 kpc.

2. OBSERVATIONS AND DATA REDUCTION

We observed the CO(2–1) transition toward the GN20 field as part of VLA key project AC974. The project was awarded 96 hr in B-configuration (baselines up to 10 km) and 28 hr in D-configuration (baselines up to 1 km), for a total of 124 hr. The observations were dynamically scheduled and took place in 2010 March–April (D-configuration) and 2011 February–April (B-configuration).

The CO(2–1) line (rest frequency ν = 230.5424 GHz) is redshifted to ν = 45.655 GHz at z = 4.05, requiring the Q band. The primary beam is ∼1' (FWHM) at this frequency, and the pointing center was chosen to be 10'' west of GN20 so that GN20, the nearby SMGs (and fellow protocluster members) GN20.2a and GN20.2b, and the z = 1.5 galaxy BzK–21,000 would all fall within the 70% sensitivity radius of the primary beam (data presented in J. A. Hodge et al. 2012, submitted). All images have been corrected for the response of the VLA primary beam. We centered the two 128 MHz intermediate frequencies (IFs) at 45.592 GHz and 45.720 GHz, for a total bandwidth of 246 MHz or 1600 km s⁻¹ (taking into account the overlap). Each of the two IFs had 64 channels, resulting in an instrumental velocity resolution of 13 km s⁻¹. Observations were taken in full polarization mode.

We used fast switching phase calibration (Carilli & Holdaway 1999), with a cycle time of five minutes. VLA calibrator J1302+5748 served as the phase calibrator. The quasar 3C286 was used to determine the absolute flux density scale and for bandpass calibration. The B-configuration data were reduced using the Astronomical Image Processing System (AIPS), and the D-configuration data were reduced in the Common Astronomy Software Applications (CASA) package. During the reduction, we discarded data taken during times of poor phase stability. After accounting for calibration overheads and flagging, the total time on source was approximately 50 hr. The data were then combined within AIPS for further analysis.

We imaged the data using the CLEAN algorithm in CASA, and we cleaned down to 1.5σ in a 2'' × 2'' clean box around GN20. A spectrum was extracted using an aperture of the same size as the clean box. The best compromise between resolution and sensitivity for this data set resulted from using Briggs weighting with a robust parameter of R = 1.0. This results in an angular resolution of 0 farcs 19, or 1.3 kpc at the redshift of GN20. We will refer to this resolution as our "native" resolution for the remainder of the paper.

3. RESULTS

3.1. CO(2–1) Detection and Channel Maps

The CO(2–1) spectrum for GN20 is shown in Figure 1. We used a Gaussian fit to the spectrum at 12 MHz (78 km s⁻¹) resolution to derive a peak flux density of 1.35 ± 0.20 mJy and an FWHM of 730 ± 140 km s⁻¹, consistent with Carilli et al. (2011). The integrated flux is 1.0 ± 0.25 Jy km s⁻¹ (see Table 1). The spectrum indicates a redshift for the galaxy of z = 4.0548 (topocentric), consistent with earlier measurements, including high-J transitions (Daddi et al. 2009b). We detect no continuum emission from GN20 at our native resolution (0 farcs 19 with a 3σ limit of 39.5 μJy and at 1 farcs 75 resolution with a 3σ limit of 65.4 μJy.

Table 1. GN20 Observed and Derived Parameters

Parameter	Value
Position (J2000)^a	12^h37^m1189 +62°22'118
z	4.0548 ± 0.0008
S_CO(2–1)	1.6 ± 0.3 mJy (peak 1)
	1.3 ± 0.3 mJy (peak 2)
FWHM_CO(2–1)	300 ± 90 km s⁻¹ (peak 1)
	320 ± 130 km s⁻¹ (peak 2)
I_CO(2–1)	1.0 ± 0.3 Jy km s⁻¹
L'_CO(2–1)	1.6 ± 0.5 × 10¹¹ K km s⁻¹ pc²
M(H₂)	1.3 ± 0.4 × 10¹¹ × (α_CO/0.8) M_☉
M_dyn	5.4 ± 2.4 × 10¹¹ M_☉

Notes. ^aFrom the 1.4 GHz observations of Morrison et al. (2010) at 1 farcs 7 resolution. All other parameters are from the study presented here.

Download table as: ASCII Typeset image

The spectrum shows the clear indication of a double-peaked, or at least flattened, structure, the signature of which is also evident in a position–velocity diagram created at a later point in the analysis (Section 4.1.2). If we instead fit the spectrum with a combination of two Gaussians (Figure 1; dot-dashed line), we find that the two components have peak flux densities of 1.6 ± 0.3 mJy and 1.3 ± 0.3 mJy, and FWHM values of 300 ± 90 km s⁻¹ and 320 ± 130 km s⁻¹, respectively. The integrated flux of the combined components is 1.0 ± 0.3 Jy km s⁻¹, consistent with the result from the single-Gaussian fit above. We will use the two-Gaussian fit to determine the CO luminosity and gas mass in the remaining analysis.

Channel maps of 78 km s⁻¹ width are shown in Figure 2. The left panel shows the emission at our native resolution (0 farcs 19), and the right panel shows the same channels tapered to 0 farcs 38 resolution in order to increase the signal-to-noise ratio (S/N). The two outer channels in each case sample the continuum on either side of the line. The emission first appears in the south of the plotted field, and it shifts to the north with increasing frequency. We will analyze the kinematics of the system in Section 4.1 below.

3.2. Molecular Gas Mass

From the two-Gaussian fit to the spectrum of GN20, we derive a CO luminosity of L'_CO(2–1) = 1.6 ± 0.5 × 10¹¹ K km s⁻¹ pc². This implies a molecular gas mass of M(H₂) = 1.3 ± 0.4 × 10¹¹ × (α_CO/0.8) M_☉ (see Table 1) assuming the standard relationships from Solomon & Vanden Bout (2005). As justified by Carilli et al. (2010), we assumed thermal excitation for the extrapolation from CO(2–1) to CO(1–0), and we used a CO-to-H₂ conversion factor of α_CO = 0.8 M_☉ (K km s⁻¹ pc²)⁻¹, the value typically assumed for ULIRGs and SMGs (Downes & Solomon 1998; Solomon & Vanden Bout 2005; Tacconi et al. 2006, 2008). In Section 4.1.3 below, we use our dynamical mass estimate to put constraints on α_CO for this system.

4. ANALYSIS

4.1. Dynamical Analysis

4.1.1. Moment Maps

Given the high quality of the data, we created moment maps by following the typical approach used for H i and CO observations of nearby galaxies (e.g., Walter et al. 2008; Leroy et al. 2009). This approach is useful for high-resolution data sets such as this (0 farcs 19), where real emission can be resolved out in a moment map with a single, global S/N cut. To address this fact, this technique involves using the data, tapered to a lower resolution (hence recovering more extended emission), as an additional input. For this purpose, we used data tapered to 0 farcs 38 resolution. This lower-resolution data cube is used to create a mask (on a per-channel basis) of the significant emission, including more diffuse emission. These masks are then used to blank the higher-resolution data, thereby retaining real, diffuse emission in the high-resolution data that would have otherwise been blanked. This standard procedure is the best possible method to recover diffuse, low S/N emission while retaining the highest possible spatial resolution.

We used the method described above to create the zeroth moment (integrated intensity) map, shown in Figure 3 (top panel). As the masking process is done on a per-channel basis, the noise in any given pixel in the map is given by $\sqrt{N}\times \sigma _{\rm chan}$ , where N is the number of channels that were integrated for that particular pixel, and σ_chan is the rms noise per channel. We used this information to apply an additional S/N cut to the first moment map (intensity-weighted velocity; Figure 3 lower panel). In particular, we required S/N > 3 for the first moment map, since any remaining unmasked noise will have a large effect on the resulting velocity field. Although some noise is still present in the outskirts, a clear velocity gradient is apparent across the disk.

As a verification of the blanking process described above, we show unblanked velocity-averaged maps of GN20 in Figure 4 at three different angular resolutions. The left panel has been tapered to 0 farcs 38, the middle panel is at the native resolution (0 farcs 19), and the right panel uses Briggs weighting with R = −0.5 to reach 0 farcs 14 resolution (1.0 kpc at z = 4.05). In contrast to the moment analysis described above, these maps were made by simply averaging over 780 km s⁻¹ (the region indicated in Figure 1), without first blanking the noise on a per-channel basis. The gas distribution appears slightly different than that seen in the zeroth moment map because the channels without emission are included in the average. The highest-resolution map confirms that the emission is not concentrated in just one or two strong peaks, but rather spread out over a larger area.

The total flux density (averaged over 780 km s⁻¹) in these velocity-averaged B+D-array maps is 633 ± 67 μJy. For comparison, the total flux densities in the B-array only (extended configuration) and D-array only (compact configuration) maps (not shown) are 365 ± 57 μJy and 813 ± 65 μJy, respectively. The total flux density in the zeroth moment map (800 ± 74 μJy) is consistent with that in the D-array only map, confirming that the recovered emission in the zeroth moment map is real. The blanking process described above therefore allows us to achieve the native resolution while still recovering the diffuse emission present on larger scales.

We estimate the size of GN20's gas reservoir from the zeroth moment map (Figure 3, top). Defining the radius as the maximum radial extent of the resolved CO(2–1) emission in the map, and conservatively assuming an uncertainty in the measurement of ∼30%, we derive a radius of 1'' ± 0 farcs 3, equivalent to ∼7 kpc (± 2 kpc) at z = 4.055. The total diameter of the source in CO(2–1) is therefore ∼14 ± 4 kpc. The large extent of the molecular gas reservoir is not unlike what has been seen in some other (lower-z) SMGs in low-excitation imaging; while SMGs typically have relatively compact distributions in the higher-J CO lines (half-width at half-maximum, HWHM = 2–4 kpc; e.g., Tacconi et al. 2006), recent observations of SMGs in CO(1–0) show more extended gas reservoirs (e.g., Ivison et al. 2011; Riechers et al. 2011a, 2011b; Ivison et al. 2010).

4.1.2. Dynamical Modeling

For the dynamical modeling of GN20, we used the GALMOD task (part of the GIPSY package). GALMOD creates a three-dimensional model using a set of input parameters, and it then convolves the model cube to the spatial/spectral resolution of the data for comparison. We used an input data cube with a spectral resolution of 26 km s⁻¹, and we tapered the data to an angular resolution of 0 farcs 77 as it was found that higher resolutions resolved out too much emission to be usable. The GALMOD task requires a radial profile as input, which (guided by the zeroth moment map) we set as an exponential radial profile with a slope of –0.4. We used a thin disk model and we found that changes in the thickness of the disk (within a reasonable range, < few kpc) did not result in major changes to the model.

As the S/N of our data have necessitated using a lower-resolution (0 farcs 77) cube for the modeling, it is not possible to measure the intrinsic rotation curve of this source directly from the data. Instead, we have assumed a rotation curve as an input and then compared the model to the data in position–velocity space. Position–velocity diagrams are shown for three different input rotation curves in Figure 5. The columns represent the different rotation curves: the left-hand column assumes a flat rotation curve, the middle column assumes a rotation curve that rises linearly from 0''–0 farcs 5 (0–3.5 kpc) to v_max, then is flat at larger radii, and the right-hand column assumes a rotation curve that rises linearly from 0''–1 farcs 0 to v_max (0–7 kpc), then is flat at larger radii. The rows show four different slices through the three-dimensional data: major axis, taken at a position angle of 25° (top row), minor axis (second row), major axis + 30° (third row), and major axis + 60° (bottom row). These comparison plots demonstrate that, while the resolution of the input data cube (due to our limited S/N) make it difficult to constrain the exact shape of the rotation curve, we can constrain the curve to having reached its flat part at least within 0 farcs 5 (3.5 kpc), with a preference for the flat rotation curve.

**Figure 5.** Position–velocity diagrams for CO(2–1) emission in GN20 for four different slices and three different models (shown as contours). The input data cube has a spectral resolution of 077 and a spectral resolution of 26 km s⁻¹. The different slices are shown in the different rows and are (from top): major axis (i.e., a position angle of 25°), minor axis, major axis + 30°, and major axis + 60°. The columns show the different models (contours): a flat rotation curve (left), a rotation curve that rises linearly to 05 (3.5 kpc) then flattens (middle), and a rotation curve that rises linearly out to 10 (7 kpc) (right). The velocities on the vertical axis are relative. Gray scale and thin contours show the observed data. The data contours are 55%–85% in steps of 10%. For all models, the contours are 35%–80% of the peak brightness, in steps of 15%. The left-hand panel shows our best-fit model (see the text for details).
Download figure:
Standard image High-resolution image

We therefore find that the velocity field is fully consistent with a rotating disk with a flat rotation curve (as discussed above). By comparing different models to the data, we find that the best-fit model is a rotating disk with an inclination of i = 30° ± 15°, a maximum rotational velocity of v_max = 575 ± 100 km s⁻¹, and a dispersion of δ = 100 ± 30 km s⁻¹. Note that the error bars are not statistical, but are liberal estimates of the uncertainty determined through careful comparison between model and data. Note also that deriving the dispersion from a spatially and spectrally convolved disk model, unlike other mean-weighted dispersion estimators, is unbiased by beam smearing (Davies et al. 2011). The spectrum predicted by our best-fit model is a good fit to our observed spectrum (Figure 1).

The relatively large value of v_max is due to the fairly small inclination value; while the quantity v_maxsin(i) is well constrained, the two components are difficult to disentangle at our angular resolution. Following standard procedure for the modeling of low to medium resolution H i observations, the final inclination value of 30° was chosen to reproduce the resulting ellipticity in the zeroth moment map. However, we cannot definitively rule out larger values of the inclination (within the quoted error), and therefore lower values of v_max. The uncertainty quoted for v_max folds in the uncertainty in the inclination.

4.1.3. Dynamical Mass and the CO-to-H₂ Conversion Factor

Assuming the parameters from the best-fit model, we derive a dynamical mass for GN20 of 5.4 ± 2.4 × 10¹¹ M_☉. The uncertainty was estimated assuming 1σ uncertainties of 100 km s⁻¹ and 2 kpc on the rotational velocity and radius, respectively. This estimate is based on dynamical modeling, making it more robust than previous estimates for this source (which also relied on higher-J transition lines; Daddi et al. 2009b; Carilli et al. 2010). We will now use this estimate to set limits on the CO-to-H₂ mass conversion factor.

In Section 3.2, we calculated molecular gas masses assuming a CO-to-H₂ conversion factor of α_CO = 0.8 M_☉ (K km s⁻¹ pc²)⁻¹, the value derived for ULIRGs. The mass conversion factor is thought to vary with environment, however, depending on several factors including metallicity, excitation, and interstellar medium pressure (Bolatto et al. 2008, and references therein). Generally, it is thought that it may decrease for objects with large gas surface densities (Downes & Solomon 1998; Scoville et al. 1997; Tacconi et al. 2008), ranging from 0.8 M_☉ (K km s⁻¹ pc²)⁻¹ for ULIRGs up to ∼4.3 M_☉ (K km s⁻¹ pc²)⁻¹ for giant molecular clouds (GMCs) in the Milky Way. (All values stated here include helium.) While the ULIRG value of 0.8 M_☉ (K km s⁻¹ pc²)⁻¹ is often used for SMGs, there is not yet any firm evidence for what the SMG value should be. It is possible that for the more luminous SMGs, it is even lower. Tacconi et al. (2008) find that a Galactic conversion factor is strongly disfavored for their nine SMGs, with the lowest χ² values for their fits resulting from more ULIRG-like values. However, their resolution only marginally resolved their sources, and they rely on CO(4–3) and even higher-order transitions that may be more centrally concentrated (see Section 4.4).

Our derived dynamical mass allows us to put constraints on the conversion factor for GN20. If we assume that GN20 is composed of 100% molecular gas, then we derive a conversion factor of α_CO = 3.3 ± 1.8 M_☉ (K km s⁻¹ pc²)⁻¹ to account for the total mass in the system. This is consistent (within the large uncertainties) with the Galactic value, but it represents the most extreme case. In reality, the total mass will also include contributions from stars, dark matter, and dust. The stellar component has been estimated to be 2.3 × 10¹¹ M_☉ by Daddi et al. (2009b) based on spectral energy distribution (SED) fitting to the Advanced Camera for Surveys (ACS) through IRAC photometry and assuming a Chabrier (2003) initial mass function. (Note that the estimate relies most heavily on the IRAC data, which are coincident with the CO emission.) The contribution from dark matter is largely unknown, but observations of high-z galaxies suggest the value may be roughly 25% (Genzel et al. 2008; Daddi et al. 2010). Using these estimates, and ignoring the contribution from dust, which is thought to make up only a tiny fraction of the total mass (2 × 10⁹ M_☉; Magdis et al. 2011), we derive a conversion factor⁸ of α_CO = 1.1 ± 0.6 M_☉ (K km s⁻¹ pc²)⁻¹. This estimate is in agreement with the limit α_CO < 1.0 derived recently for GN20 by Magdis et al. (2011) using the local M_gas/M_dust versus metallicity relation.

4.2. Molecular Gas Clumps

4.2.1. Clump Search

In order to study the properties of the molecular gas in GN20 in more detail, we have attempted to identify individual gas clumps in the disk of GN20. We used two different methods to identify potential clumps, judging as reliable only those clumps that were independently identified by both algorithms. The first method used the AIPS task SERCH. SERCH is a "matched-filter" analysis employing a Gaussian kernel and used in source finding (e.g., Begum & Chengalur 2005; Kanekar & Chengalur 2005; Aravena et al. 2012), and the algorithm is described in more detail in Uson et al. (1991). We used a cube at our native angular resolution (0 farcs 19) and with a spectral binning of 39 km s⁻¹. We then searched for emission line sources in the three-dimensional data, optimizing the search for line widths between 80 and 200 km s⁻¹ in steps of the channel size (∼40 km s⁻¹). An S/N cut of 4.5σ yielded 11 positive peaks and two negative peaks in a generous 4'' × 4'' region centered on GN20. Based on the number of negative detections, we can expect ∼2 of the positive peaks to be spurious. Ten of the positive peaks lie on the disk of GN20, as defined by the zeroth moment map in Figure 3, and one positive peak and both negative peaks lie off of the disk.

To test the robustness of the clumps identified in GN20 by SERCH, we repeated the search using an independent method. The second method is built upon the principles of Bayesian inference, and we will briefly describe it here. (For a more detailed description, see L. Lentati et al. 2012, in preparation.) This method involves using the multimodal nested sampling algorithm MULTINEST (Feroz & Hobson 2008) to efficiently sample the posterior distribution by fitting a simple parametric model to the data cube described both spatially and spectrally by a Gaussian. The spatial model parameters are fixed to match the clean beam, since the clumps are not expected to be resolved. Spectrally, the model is a Gaussian line profile described by a profile width σ_ν and peak amplitude A_p. The priors on σ_ν and A_p are set as uniform between 20 and 150 km s⁻¹ and 0.0 and 0.4 mJy, respectively, and are considered independent.

We then quantified the probability that the fitted model describes a real clump by using Bayesian model selection (see, e.g., Hobson & McLachlan 2003). The two models considered were:

1.
No clump is present at that position.
2.
A clump with A_p > 0 exists at that position.

We evaluated the evidence associated with the posterior for these competing models, where the evidence is the average of the likelihood over the prior. Specifically, MULTINEST returns the evidence and associated error for each peak in the posterior which can then be compared with the evidence for there being no clump present. The model corresponding to the posterior mode with the greatest probability was then subtracted from the data, and the process was repeated until no sources above a low threshold probability were found. Individual clumps were identified by requiring a minimum distance of 0 farcs 12 and two frequency channels between their centers. If multiple clumps were returned within that distance, the clump with the greatest evidence was chosen. This method returned a list of positive and negative clumps, sorted in order of their probability. The negative clumps were all very low on the list, with 20 positive clumps judged as more probable than the most probable negative clump.

We then cross-matched the lists of clump candidates produced by both SERCH and the Bayesian modeling, requiring the clump centers returned by the two algorithms to agree within a beam radius in order to count as the same clump. Of the ten potential clumps identified on the disk of GN20 by SERCH, six were independently identified by the Bayesian modeling (and these were judged to be six of the seven most probable clumps by the Bayesian modeling). The four SERCH clumps that were not identified by the Bayesian modeling included the two SERCH clumps farthest out on the disk, and two SERCH clumps which were closer than the minimum distance set in the Bayesian modeling. Of the six common clump candidates, one was discarded because of an unnaturally thin profile (i.e., a single spectral channel). The final sample therefore consists of five clumps. These are the most probable clumps that were independently identified by the two different methods, and therefore are the five strongest clump candidates. They are labeled in Figure 6 at their positions on the CO(2–1) zeroth (left) and first (right) moment maps with beam-sized circles (1.3 kpc). As these clumps were identified in the three-dimensional data, they do not necessarily perfectly align with the most significant peaks in the integrated map. We also note that, while these clumps were judged to be the most reliable clumps, they in no way constitute a comprehensive census of the clumps in GN20.

Spectra for these five clumps are shown in Figure 7, where the panel numbers correspond to the clump numbers in Figure 6. Gaussian fits to the spectra are overplotted. The coherent velocity at the position of each clump is shown in each panel, and is based both on the model velocity field at lower spatial resolution (dotted line) and the data velocity field at its native resolution (dot-dashed line). In general, the clumps are within ∼100 km s⁻¹ of the (observed) coherent velocity at their position on the disk. The fact that they are not centered exactly on the observed coherent velocity is due to the presence of other mass along the line of sight. The model velocity field is based on data at a lower spatial resolution (0 farcs 77) and so does worse (generally) at predicting the clump velocities. While these spectra have low S/N, they allow us to estimate the properties of individual gas clumps in a z = 4 galaxy for the first time.

4.2.2. Clump Properties

With this data, we attempt to estimate some basic properties of the molecular gas clumps. We caution that this analysis is based on low S/N data and will need to be verified by even higher sensitivity observations. To begin with, we use the spectral information to determine brightness temperatures of individual clumps without diluting the values by averaging in velocity. Derived brightness temperatures (using the Rayleigh–Jeans approximation) range from 3.2 to 6.2 K and can be read off of the right-hand y-axis in Figure 7. After correcting to the rest frame, values for individual clumps range from ∼16 K up to 31 K (Table 2). For reference, Planck temperatures (i.e., the temperatures calculated using the full blackbody instead of just the Rayleigh–Jeans approximation) are typically ∼20% higher.

Table 2. Molecular Gas Clumps Properties

Clump^a	S_CO(2–1)	FWHM_CO(2–1)	S/N	T_B (Rest Frame)	M(H₂)	% Total M(H₂)	M_virial
	(μJy)	(km s⁻¹)	(σ)	(K)	(× 10⁹ (α_CO/0.8) M_☉)	(× 10⁹ M_☉)	(× 10⁹ M_☉)
1	360 ± 100	90 ± 40	5.3	31 ± 8	4.4 ± 2.2	3% ± 2%	<1.0 ± 0.9
2	270 ± 60	150 ± 50	5.4	23 ± 5	6.0 ± 2.4	5% ± 2%	<3.4 ± 2.2
3	270 ± 60	140 ± 40	5.1	23 ± 6	5.5 ± 2.1	4% ± 2%	<2.9 ± 1.8
4	180 ± 50	160 ± 60	4.5	16 ± 4	4.2 ± 1.9	3% ± 1%	<3.6 ± 2.6
5	240 ± 70	140 ± 50	4.6	21 ± 6	4.8 ± 2.3	4% ± 2%	<2.7 ± 2.0

Note. ^aSee Figure 6.

Download table as: ASCII Typeset image

Since the clumps appear to be unresolved on the scales probed by our observations (1.3 kpc), the brightness temperature values are most likely lower limits. Nevertheless, it is interesting to note that the peak brightness temperatures derived for the clumps are already approaching (or even exceeding) the dust temperature (33K; Magdis et al. 2011). If we instead compare it to the large velocity gradient (LVG; e.g., Scoville & Solomon 1974) gas model fit by Carilli et al. (2010), where the best-fit model consisted of a lower and a higher excitation component, the higher excitation component still has a kinetic temperature of only 45 K. The change in surface area required to make up this temperature difference is small. This may indicate that we are close to resolving the clumps, though we cannot make this statement any more quantitative with the current data set.

We next fit the line widths of the clumps. There are several factors that contribute to broadening the clump line widths beyond their intrinsic values. One factor is the instrumental velocity resolution, which we have removed in quadrature from the measured line widths. Another factor is the rotational contribution from the large-scale velocity structure of the galaxy. To estimate this contribution, we have taken the model velocity field and measured the velocity gradient across a beam-sized area centered on each clump, dividing by two to approximate the weighting of the beam. The resulting line widths are listed in Table 2 and have a median value of 140 km s⁻¹. The errors were derived from the Gaussian fits and assuming a 30% error on the measured velocity gradient.

Using the Gaussian fits to the spectra in Figure 7, we determined the molecular gas masses of the clumps (Table 2). For this, we used the same relations and assumptions as in Section 3.2, including an H₂-to-CO conversion factor of α_CO = 0.8 M_☉ (K km s⁻¹ pc²)⁻¹. We estimate that the clumps have H₂ masses of ∼5 × 10⁹ × (α_CO/0.8) M_☉, or a few percent of the total gas mass each.

Taking the beam size as an upper limit on the clump size, the mass surface densities of the clumps are >3200–4500 × (α_CO/0.8) M_☉ pc⁻². Similar surface densities have been seen in some compact SMGs with double-peaked line profiles (3000–10,000 M_☉ pc⁻²; Engel et al. 2010). GMCs in the nearby universe, on the other hand, are typically several orders of magnitude (in volume) smaller and have surface densities approximately 20 times lower (Solomon et al. 1987). Tacconi et al. (2008) argue that a major merger is necessary to produce such high-mass surface densities as those seen in SMGs. However, Engel et al. (2010) have pointed out that this argument alone is insufficient to distinguish a merger. Indeed, if we assume that the clumps are spherical and roughly the size of the beam (at maximum), then we derive (average) volume densities of only 70–110 cm⁻³, consistent with the mean density of local GMCs (Solomon et al. 1987). The dense cores of GMCs show much higher values, ∼10⁴ cm⁻³ (e.g., Goldsmith et al. 1997), and the same may also be true for the cores of the clumps in GN20. The mean H₂ volume density of GN20 as a whole is only ∼15 cm⁻³.

We estimated the virial masses of the clumps using the isotropic virial estimator (e.g., Erb et al. 2006):

$\begin{equation} M_{\rm virial} = \frac{C\sigma _{\rm v}^{2}R_{\rm g}}{G} (M_{\odot }) , \end{equation} \tag{ 1 }$

where C = 5 is the dimensionless prefactor for a uniform sphere, σ_v is the observed one-dimensional velocity dispersion in km s⁻¹, R_g is the gravitational radius, and G = 1/232 pc (km s⁻¹)² M⁻¹_☉ is the gravitational constant. If we take the beam HWHM as an upper limit on the gravitational radius, we estimate virial masses for the clumps of <1.0–3.6 × 10⁹ M_☉ (Table 2). These masses are all within 1σ–2σ of the observed molecular gas masses, even assuming a low ULIRG-like mass conversion factor of α_CO = 0.8 M_☉ (K km s⁻¹ pc²)⁻¹ and ignoring the possible presence of significant stellar mass. Within the considerable uncertainties, it is possible that the clumps have masses consistent with self-gravitation.

The gas masses derived for the clumps are already on the high side (in general) of the virial mass estimates, even with a low conversion factor and despite the fact that the virial mass estimates are upper limits. This implies that for the majority of the clumps, the gas masses would significantly overshoot the virial masses for α_CO ≫ 0.8 M_☉ (K km s⁻¹ pc²)⁻¹. Assuming the clumps are entirely composed of molecular gas and using M_virial/L'_CO(2–1) as an estimate of α_CO, we find α_CO for the clumps is <0.2–0.7 M_☉ (K km s⁻¹ pc²)⁻¹.

Recently, a couple of other observational studies have looked at molecular gas clumps in z > 1 SFGs. For example, Tacconi et al. (2010) observed molecular gas clumps in a "normal" SFG at z ∼ 1.2. These clumps had gas masses of ∼5 × 10⁹ M_☉, diameters <2–4 kpc, brightness temperatures of ∼10–25 K, and gas volume densities ∼100 cm⁻³, similar to what we get for GN20 if we assume that the clumps are roughly the beam size. Their velocity dispersions, however, were only ∼20 km s⁻¹, causing them to postulate that the "clumps" observed were more likely loose conglomerates of several GMCs. Another recent paper by Swinbank et al. (2011) identified molecular gas clumps in a lensed SFG at z = 2.3, SMM J2135. The clumps in SMM J2135 had gas masses of ∼0.3–1 × 10⁹ M_☉, a median diameter of 200 pc, and gas volume densities of ∼1000–40,000 cm⁻³. To achieve volume densities this high, the clumps in GN20 would need to be up to an order of magnitude smaller than the current beam size.

4.3. Multiwavelength Comparison

Figure 8 shows the zeroth moment CO(2–1) contours at intermediate resolution overlaid on a selection of the multiwavelength data available for GN20, including the Hubble Space Telescope (HST)+ACS z850-band (top left; Giavalisco et al. 2004), the WFC3 140W-band (top right), the WIRCAM K-band (bottom left), and the IRAC 3.6 μm data (bottom right). Not shown is the higher-resolution VLA+MERLIN data (Casey et al. 2009), which was compared to GN20 in CO in Carilli et al. (2010). The object detected in the HST+ACS z850-band data (top left panel) was followed up spectroscopically with the Deep Imaging Multi-Object Spectrograph (DEIMOS) on Keck II and determined to have a redshift of z = 4.055 (Daddi et al. 2009b), consistent with the redshift derived for GN20 from its CO lines.

Assuming that the object is related to GN20, as indicated by the good match in redshift, then the question is, what is causing the clear offset from the CO and radio/(sub)millimeter position also reported in previous work on GN20 (for e.g., Carilli et al. 2010; Pope et al. 2006; Younger et al. 2008; Iono et al. 2006)? We have confirmed that the astrometry of the HST imaging is good to <0 farcs 15 (Daddi et al. 2007), so astrometric error is unlikely. In addition, the counterpart in the WFC3 140W-band imaging is similar in position/morphology. Such large offsets between dust/molecular gas emission and (rest-frame) UV/optical imaging have been seen in a number of SMGs (Chapman et al. 2005; Capak et al. 2008; Riechers et al. 2010). A possible hint in this case comes from the WIRCAM K-band image, which is similarly offset, but shows extended emission in the direction of the radio/(sub)millimeter counterparts. This extended emission is significant at the ∼6σ level. A progression through the observations shows that the peak of emission shifts toward the radio position as we sample longer wavelengths. One interpretation is therefore that the radio/(sub)millimeter emission is offset from the (observed) optical emission because the dust that is presumably associated with the molecular gas heavily obscures the rest-frame UV and optical emission. A significant amount of dust obscuration is not surprising in a starburst of several thousand M_☉ yr⁻¹, the SFR estimated for GN20 (Magdis et al. 2011; Daddi et al. 2009b). Heavy dust obscuration is thought to explain the complete absence of optical counterparts for some SMGs (e.g., Walter et al. 2012). The top left panel in Figure 8 shows that the rest-frame UV emission is well aligned with the edge of the zeroth moment CO contours, giving the clearest indication that the ACS counterpart is physically related but is simply offset due to differential dust obscuration.

4.4. Spatially Resolved Gas Excitation

To investigate the spatially resolved molecular gas excitation, we compare our CO(2–1) data with the CO(6–5) Plateau de Bure Interferometer (PdBI) data from Carilli et al. (2010). The CO(6–5) data are uniformly weighted at an angular resolution of 0 farcs 84 × 0 farcs 67. We tapered the CO(2–1) data to the same resolution and imaged over the same velocity range (800 km s⁻¹).

We then used the AIPS task IRING to calculate the average CO(6–5)/CO(2–1) ratio as a function of radius. The centroids of the CO(6–5) and CO(2–1) emission peaks were consistent within 0 farcs 08, so we used the average of the positions as the center of the concentric annuli. We set the width of the annuli to be 0 farcs 25, i.e., one-third of the synthesized beam, and we specified a major axis position angle of 115° and an inclination of 30° as determined in Section 4.1.2.

The results are shown in Figure 9. The left panels show the CO(2–1) and CO(6–5) emission, using the native beam shape of the CO(6–5) map. The plot on the right is the result of dividing the CO(6–5) emission by the CO(2–1) emission, averaging in elliptical annuli as explained above. Note that the error bars are purely statistical and do not account for the overall flux calibration uncertainties in the maps (which in both cases are of order 10%).

The CO(6–5)/CO(2–1) ratio is an indicator of the excitation state of the gas. The ratio hovers around a value of ∼3, a factor of three lower than expected for thermal excitation. This indicates that the CO(6–5) line is sub-thermally excited, consistent with the CO excitation ladder shown in Carilli et al. (2010). Carilli et al. determined that the data were best fitted with a two-component gas model including a diffuse, low-excitation component and a more concentrated high-excitation component.

With the data at hand, we can go one step further and investigate the radial dependence of the CO ratio. From Figure 9, we see that the ratio is consistent with a constant value out to at least 1'', corresponding to the observed radius of the CO(2–1) emission in the zeroth moment map. Note that we have spaced the radial bins so as to achieve Nyquist sampling, so the first two data points are within the beam. Nevertheless, this result tells us that to first order, the CO(6–5) emission has the same radial profile as the CO(2–1) emission. This would not be the case if, for example, the CO(6–5) emission (which is more closely related to star formation) was more tightly concentrated than the lower-excitation CO(2–1) emission. To illustrate this point, we show two different models (at the same resolution as the data) where the CO(6–5) emission is more centrally concentrated than the CO(2–1) emission. In the first model (dotted line), the exponential scale length of the CO(6–5) emission is two times smaller than that of the CO(2–1) emission. In the second model (dashed line), the CO(6–5) emission has a scale length that is five times smaller. In both cases, the CO(2–1) disk parameters come from our dynamic modeling in Section 4.1.2. The models have been normalized to the peak value in the maps and are meant to demonstrate the clear drop-off in the excitation ratio that would be detectable for a central star-forming event even at this resolution. Instead, the flatness we observe is consistent with the idea that the two transitions are cospatial, and that the star formation is therefore spread out over a large portion of the (14 kpc diameter) disk.

5. WHAT IS DRIVING THE STARBURST?

GN20 has several properties which make it unique as an SMG. It has a higher redshift than the typical SMG, and it is extremely bright—its 850 μm flux density of 20.3 mJy makes it one of the most luminous starburst galaxies known (Pope et al. 2006). In addition, it appears to be part of a massive protocluster of galaxies: there are two other z ∼ 4.05 SMGs just ∼25''/170 kpc away known as GN20.2a and GN20.2b (Daddi et al. 2009b; Carilli et al. 2011); there are 14 z ∼ 4 Lyman break galaxies within 25''; and there is an overdensity of z > 3.5 IRAC-selected galaxies in the field.

The question is: What process is driving the gas accretion and conversion into stars in this example of one of the earliest extreme starburst galaxies? The majority of SMGs are thought to be gas-rich major mergers. There are multiple examples of SMGs with very disturbed kinematics and even showing multiple components, providing direct observational evidence for the merger picture (e.g., Tacconi et al. 2008; Engel et al. 2010; Ivison et al. 2011; Riechers et al. 2011a, 2011b). For example, by identifying likely radio and/or mid–IR identifications to SMGs in the SHADES Source Catalog (Coppin et al. 2006), Ivison et al. (2007) found that significantly more SMGs have multiple robust counterparts than would be expected by chance alone.

In GN20, on the other hand, we observe an extended, rotating gas disk with a velocity field that is consistent with a flat rotation curve. The gas reservoir appears to be fragmented into multiple clumps, with each clump making up only a few percent of the total mass. This morphology stands in contrast to prototypical merging sources like Arp 220 and other ULIRGs, where the molecular gas is typically concentrated in one or two dense, centrally located regions. If the diffuse and clumpy gas components correspond to the two gas excitation states predicted to exist in GN20 with LVG modeling (Carilli et al. 2010), then the clumpy component (presumably given by the higher excitation state gas) only makes up ∼ half of the total gas mass, with the remainder existing as extended, diffuse emission.

Recently, CMA has been suggested as a possible alternative to the major merger formation scenario for some SMGs (Fardal et al. 2001; Finlator et al. 2006; Davé et al. 2010; Carilli et al. 2010). In this scenario, the star formation in these SMGs would be instead fueled by minor mergers and the smooth infall of cold gas along streams. Unlike a major merger, this process can leave the disk relatively undisturbed, especially for a galaxy as massive as GN20 (Davé et al. 2010). However, wet mergers have been shown to rapidly relax into smooth disks, so the presence of ordered rotation alone cannot immediately be taken as an argument against a major merger (Robertson et al. 2006; Robertson & Bullock 2008; Bournaud & Elmegreen 2009).

If GN20 were a major merger where the gas had re-virialized into a rotating disk, one may ask whether it would still be an SMG at this late stage. According to hydrodynamic simulations of major mergers, the submillimeter-bright phase in a (z ∼ 2) major merger-driven starburst should occur as the galaxies approach final coalescence (Narayanan et al. 2009, but cf. Hayward et al. 2011, 2012). During this brief phase (∼0.03 Gyr), the 850 μm flux is significantly boosted and the galaxy may be selected as a luminous submillimeter source. However, according to the simulations, this phase coincides with the stage when the disk-like morphology is most disturbed, with the largest fraction of gas tidally disrupted from disk rotation (Narayanan et al. 2009, their Figure 6). These simulations also show that the gas during the SMG phase is relatively concentrated, with a characteristic radius of ∼1.5 kpc. Such a compact morphology is expected, since the large-scale gravitational torques induced by gas-rich major mergers are efficient at removing angular momentum (Barnes & Hernquist 1996). In GN20, on the other hand, we simultaneously observe (1) an extremely bright (S₈₅₀ = 20.3 mJy) submillimeter flux, (2) ordered rotation, and (3) a very extended (14 kpc) gas disk. Taken together, these observations may be an indication that the star formation in GN20 is fueled by some process other than a major merger. A similar argument was put forward by Bothwell et al. (2010) regarding the SMG HDF132. Of course, this argument rests on the assumption that these simulations provide an accurate representation of merger-driven SMG formation for all SMGs, which may not be the case, especially for notable SMGs like GN20.

On the other hand, there are certain properties of GN20 which may be difficult to explain without a major merger. In particular, the rotational velocity derived from the dynamical modeling is extremely high (576 km s⁻¹). While we cautioned that this is largely dependent on the galaxy's inclination (which is hard to constrain at our S/N), a rotational velocity this large may be hard to reproduce in a spinning disk embedded in its dark matter halo, and the line width could be more easily explained as resulting from the final coalescence stage of a major merger. Also noteworthy is the observation that GN20 is undergoing an extremely dusty massive starburst. This is supported by its high specific SFR (Daddi et al. 2009b), and also by its hot FIR SED (Magdis et al. 2011) and low CO–to–H₂ conversion factor, which both imply a high star formation efficiency. The critical question is: Is it possible to form such an object without relying on a major merger? Locally, at least, the answer is no—such objects are all major mergers (ULIRGs).

Perhaps the real answer lies somewhere in between the two scenarios. Finlator et al. (2006) studied z ∼ 4 B-band dropout galaxies in cosmological hydrodynamic simulations and found two instances of galaxies forming stars at over 1000 M_☉ yr⁻¹, several times that of their 100 Myr averages. These objects both resided in highly overdense regions resembling the cores of massive clusters. The authors found that, while a merger of two massive galaxies was not the dominant mechanism for fueling star formation, there were still signs of interaction with many of the smaller nearby galaxies. GN20 lies in a similarly dense environment—the overdensity of B-band dropouts and IRAC-selected galaxies suggests a protocluster of total mass 10¹⁴ M_☉ centered on GN20 and the nearby SMGs GN20.2a and GN20.2b (Daddi et al. 2009b). It is likely, therefore, that even if there is a dominant rotating disk, enhanced star formation over the disk is being triggered by the gravitational torques from the other galaxies.

6. CONCLUSIONS

We have presented VLA observations of CO(2–1) in the z = 4.05 SMG GN20. These high-resolution data allow us to do detailed imaging of the molecular gas on scales down to 1.3 kpc. We use these data to model the overall gas dynamics and determine the properties of individual star-forming clumps in an extreme starburst galaxy just 1.6 Gyr after the big bang. We summarize our main findings in the following.

The data reveal a clumpy, extended gas reservoir in unprecedented detail. A carefully constructed zeroth moment map recovers emission over an area 14 ± 4 kpc in diameter, and the first moment map shows a clear velocity gradient. A dynamical analysis shows that the data are consistent with a large, rotating disk with a maximum rotational velocity of v_max = 575 ± 100 km s⁻¹, an inclination of 30° ± 15°, and a dispersion of 100 ± 30 km s⁻¹. While it is difficult to determine the exact shape of the rotation curve, a flat rotation curve gives the best fit to the data.

From the dynamical modeling, we derive a dynamical mass for GN20 of 5.4 ± 2.4 × 10¹¹ M_☉. We use this dynamical mass in combination with the stellar mass to put constraints on the CO-to-H₂ mass conversion factor (α_CO), finding α_CO = 1.1 ± 0.6 M_☉ (K km s⁻¹ pc²)⁻¹. This value is consistent with a low, ULIRG-like value, although better data is needed to confirm this result. Our estimate is also in agreement with recent limits set for GN20 using the local M_gas/M_dust versus metallicity relation.

We use two different source-finding techniques to identify five distinct molecular gas clump candidates on the disk of GN20. These clumps have masses of a few percent of the total gas mass and a median line width of 140 km s⁻¹ (FWHM). Their surface densities are >3200–4500 × (α_CO/0.8) M_☉ pc⁻², corresponding to volume densities of >100 cm⁻³ (where the lower limit corresponds to GMCs in the Milky Way). Their peak brightness temperatures (16–31 K) are approaching both the average dust temperature (33K) and the "warm" component of a two-component LVG model, implying that we may be close to resolving the clumps. Within the substantial uncertainties, all clumps are consistent with being self-gravitating masses, and we constrain their CO-to-H₂ mass conversion factor to be <0.2–0.7 M_☉ (K km s⁻¹ pc²)⁻¹.

We compare the observed distribution of CO to the emission seen at other wavelengths. This multiwavelength comparison demonstrates that the molecular gas is concentrated in a region of the galaxy which is heavily obscured in the rest-frame UV and in the optical. The optical counterpart likely constitutes a small percentage of the stellar light.

Using previous observations of the CO(6–5) in GN20, we examine the CO(6–5)/CO(2–1) excitation ratio in elliptically averaged annuli. The excitation ratio is consistent with a constant value over the extent of the 14 kpc disk. This result implies that the star formation is spread out over a large portion of the disk, rather than concentrated in a central-star-forming event.

We discuss our results in the context of different fueling scenarios for GN20. While the presence of an extended, rotating disk during the SMG phase may point toward a process other than a major merger (e.g., CMA), its large line width and extremely dusty starburst would be more easily explained in a merger picture. If the star formation is fueled primarily by CMA, then, due to GN20's location in a dense protocluster environment, it is likely that it is being enhanced by interactions with nearby galaxies.

We conclude by highlighting the importance of low-J CO studies, which trace the bulk of the molecular gas, in studying the formation of massive galaxies throughout cosmic time. Such studies cannot currently be done in a statistical way due to the time required to attain sufficient S/N at high resolution, even with state-of-the-art instruments such as VLA. This fact makes targeted case studies such as this that much more important.

The authors thank the anonymous referee for helpful comments which improved this paper. We thank Glenn Morrison, Mark Dickinson, Bram Venemans, Gabriel Brammer, Desika Narayanan, Benjamin Weiner, and Dario Colombo for useful comments and discussions. C.C. thanks the Kavli institute of Cosmology for their hospitality. D.R. acknowledges funding from NASA through a Spitzer Space Telescope grant. The National Radio Astronomy Observatory is a facility of the National Science Foundation under cooperative agreement by Associated Universities, Inc. This research has made use of the GIPSY package.

EVIDENCE FOR A CLUMPY, ROTATING GAS DISK IN A SUBMILLIMETER GALAXY AT z = 4

Article metrics

Permissions

Author e-mails

Author affiliations

Dates

ABSTRACT

1. INTRODUCTION

2. OBSERVATIONS AND DATA REDUCTION