Dynamical Characterization of Galaxies at z via Tilted Ring Fitting to ALMA [CII] Observations
Until recently, determining the rotational properties of galaxies in the early universe (, Universe age Gyr) was impractical, with the exception of a few strongly lensed systems. Combining the high resolution and sensitivity of ALMA at (sub-) millimeter wavelengths with the typically high strength of the [C ii] 158 m emission line from galaxies and long-developed dynamical modeling tools raises the possibility of characterizing the gas dynamics in both extreme starburst galaxies and normal star forming disk galaxies at . Using a procedure centered around GIPSY’s rotcur task, we have fit tilted ring models to some of the best available ALMA [C ii] data of a small set of galaxies: the MS galaxies HZ9 & HZ10, the Damped Lyman-alpha Absorber (DLA) host galaxy ALMA J0817+1351, the submm galaxies AzTEC/C159 and COSMOS J1000+0234, and the quasar host galaxy ULAS J1319+0950. This procedure directly derives rotation curves and dynamical masses as functions of radius for each object. In one case, we present evidence for a dark matter halo of M. We present an analysis of the possible velocity dispersions of AzTEC/C159 and ULAS J1319+0950 based on matching simulated observations to the integrated [C ii] line profiles. Finally, we test the effects of observation resolution and sensitivity on our results. While the conclusions remain limited at the resolution and signal-to-noise ratios of these observations, the results demonstrate the viability of the modeling tools at high redshift, and the exciting potential for detailed dynamical analysis of the earliest galaxies, as ALMA achieves full observational capabilities.
Rotation curves of galaxies reveal the underlying mass structure, hint at past mergers or gravitational interactions, and allow for simple dynamical mass estimates (Sofue & Rubin, 2001). These profiles of velocity as a function of radius have been regularly measured for local galaxies since Pease (1918) deduced the rotation curve of M31. Subsequent observations of this source (e.g., Babcock 1939; Rubin & Ford 1970) showed an initial rise in velocity near the center of the galaxy followed by a nearly constant velocity out to large radius. This profile type was confirmed in other galaxies (e.g., Burbidge et al. 1960; Bosma 1978), leading to the question of its cause. Solid-body rotation would give a linear relationship between velocity and radius, while Keplerian rotation would result in decreasing velocity at large radius. The sum of the masses of the baryonic components as a function of radius was shown to be insufficient to cause the flat rotation curves, necessitating a “dark matter” component to give rise to the flat rotation curves.
The question of galaxy dynamics in general, and of the dark matter content specifically, of disk galaxies at high redshift is only recently coming into consideration with the advent of new large facilities. The dark matter content of early-type galaxies at high redshift has been inferred from gravitational lensing (e.g., Langston et al. 1990). For late-type galaxies, some progress on determining dynamics at , has been made using integral field units and the H line (e.g. Wuyts et al. 2016). Wuyts et al. find evidence for rotating disks, and that baryons comprise most () of the dynamical mass within the optical disks of galaxies out to . A few cases of galaxy dynamics based on CO observations at have been published, but in all cases the galaxies are the most extreme massive starbursts, and even then, the resolution and signal to noise has been limited (Riechers et al. 2008; Tacconi et al. 2010; Hodge et al. 2012; Tacconi et al. 2013; see Carilli & Walter 2013 for a review).
Tilted ring models have been used as a standard analysis tool to characterize the rotation of galaxies since Rogstad et al. (1974). They consist of a number of rings with variable width, systemic velocity, position angle, inclination, and rotational velocity. By varying each, creating a model velocity field, and comparing it to the observed velocity field, the dynamical nature of the object can be determined. While they allow for the existence of warps (i.e., variations in inclination angle across the disk), counter-rotation, and velocity gradients, their azimuthal averaging wipes out any inherent information regarding spiral arm structure or central bars.
However, the fitting process directly outputs the rotation curve, allowing for dynamical mass calculation and kinematic characterization. These dynamical masses may then be used to test the relationship between the mass of a galaxy and the mass of its central black hole (see Kormendy & Ho 2013 for a review) or to probe the amount of dark matter present (e.g., Richards et al. 2015). In addition, a detection of the flat portion of a rotation curve argues that the underlying structure is a rotating disk, not a velocity gradient due to a merger event or infall/outflow (e.g., García-Burillo et al. 2014; Ueda et al. 2014). We employ the tilted ring modeling herein, to create rotation curves and radial dynamical mass profiles of galaxies.
Excellent rotation curves have been created for local objects using H i data (e.g., Vogt et al. 1996; de Blok et al. 2008; Schmidt et al. 2014; Richards et al. 2015, 2016). These observations have the benefits of high brightness and large spatial extent. However, this tracer does not reveal the dynamics of the central region of the disk, due to the low density of H i in the center of galaxies (Wang et al., 2014). An alternative is H, which traces ionized gas and has been used to observe the kinematics of both galaxies (e.g., Genzel et al. 2006; van Starkenburg et al. 2008; Wright et al. 2007; Epinat et al. 2012; Burkert et al. 2016) and local galaxies (e.g., Persic & Salucci 1995). Additionally, CO (Frank et al. 2016), which traces molecular gas, and [C ii], which traces both H ii regions and colder gas (Sofue & Rubin 2001; Venemans et al. 2016), may be used.
The 158m line of [C ii] () has two pointed advantages in dynamical analysis of early galaxies. First is that the line is typically the brightest emission line at wavelengths from the radio through the FIR from star forming galaxies. And second, [C ii] emission comes from multiple phases of the ISM, including the warm ionized medium, warm and cold atomic medium, and warm dense molecular medium (Pineda et al., 2013). While this generally complicates interpretation of [C ii] luminosity, this wide range of physical conditions makes [C ii] an excellent tracer of kinematics over large areas of the disk. It is also excited by a range of sources, including collisions with electrons, hydrogen atoms, and H (Goldsmith et al., 2012). It can act as a tracer of molecular hydrogen in “dark” clouds deficient in CO (Langer et al., 2010).
It should be noted that [C ii] traces a smaller region than other dynamical tracers (i.e., HI, H, CO; de Blok et al. 2016). Because of this, the outer disk will likely not be observable. While its rest-frame frequency of 1.9 THz (158 m) requires space-based telescopes (or the Stratospheric Observatory for Infrared Astronomy; SOFIA) for study of nearby galaxies, emission from high- objects is redshifted into the ALMA () and NOEMA () bands. With the advent of ALMA and the use of redshifted [C ii], it is possible to extend classical dynamical analysis to the first galaxies, possibly revealing the dark matter content therein.
Here, we present the results of using long-developed tools to fit tilted ring models to a sample of two Lyman break galaxies (LBGs), a damped Lyman- absorber (DLA) host galaxy, two submillimetre galaxies (SMGs), and a QSO host galaxy, all between , as observed with ALMA in [C ii]. In addition, we use a limited form of three-dimensional dynamical fitting to constrain the velocity dispersions and density profiles of one of the above SMGs and the QSO host galaxy. The paper is organized as follows: we detail the observations of the sources in our sample in Section 2; the procedure to fit tilted rings to these sources is explained in Section 3, with the results of this process in Section 4; an alternate fitting procedure is discussed in Section 5; the accuracy of our derived parameters is discussed in Section 6; and we summarize in Section 7.
We will assume (,,h)=(0.682,0.308,0.678) (Planck Collaboration et al., 2016) throughout.
Our sample consists of six galaxies at observed in [C ii] with ALMA. Selection was based on the presence of a velocity gradient that was well modeled by a tilted ring model.
Multiple other sources were considered, but not included in our sample, including HZ6 in the Capak et al. (2015) sample, the SMG in BRI1202–0725 (; Carilli et al. 2013), the SMGs Vd-17871 (; A. Karim et al., in preparation), AzTEC-3 (; Riechers et al. 2014), the QSO host galaxy ULASJ2315+0143 (; Banerji et al. 2017), and the local main sequence galaxy NGC4191 (L. Young et al., in preparation). HZ6 has the third highest [C ii] detection in the same sample as two galaxies that were successfully fit, but our spectral resolution was too low to allow for successful fitting. Similarly, our spatial resolution for the SMG in BRI1202–0725 was too low. The SMGs Vd-17871 and AzTEC-3 showed gradients that were too weak to fit. Both ULASJ2315+0143 and NGC4191 were well fit, but were published separately due to their line emission being CO rather than [C ii].
It may seem strange that the sample of galaxies that we successfully fit (main sequence galaxies, SMGs, and a QSO host galaxy) are qualitatively similar to those that we failed to fit (SMGs and a QSO host galaxy). This shows that while our procedure is applicable to a range of galaxy types, the intrinsic velocity gradient must be strong enough to be observed at a given spatial and spectral resolution. In addition, the rotation must be orderly and disk-like (i.e., not early-stage major mergers or outflow-dominated galaxies).
Below, we detail the observations of each source that was successfully fit.
2.1 HZ9 and HZ10
HZ9 () and HZ10 () were originally selected as LBGs in the COSMOS field (Leauthaud et al., 2007), and were observed in [C ii] and rest-frame-FIR emission in cycle 1 ALMA by Capak et al. (2015). Both HZ9 and HZ10 were detected in rest-frame m emission ( mJy and mJy) and [C ii] ( Jy km s and Jy km s). The star formation rate (SFR) of each source was estimated using SED fits to greybody models based on the m continuum flux density and source intensities in COSMOS and Ultra-Vista images. A range of greybody parameters was assumed from observations of objects, resulting in SFRs of ( M year and M year). [N ii] emission from HZ10 was detected by Pavesi et al. (2016).
ALMA [C ii] observations for HZ9/HZ10 (project 2012.1.00523.S) used 27/25 antennas, consisted of 40/50 minutes of on-source time, and were undertaken on November 14/16, 2013. These band 7 observations resulted in spatial resolutions of and a spectral resolution of km s for the [C ii] line. Our HZ9 data were taken from directly from Capak et al. (2015). For HZ10, we use the updated calibration by Pavesi et al. (2016).
ALMA J081740.85+135138.2 (J0817+1351) is the host galaxy of a high metallicity damped Lyman-alpha absorber (DLA) at . The DLA at was originally detected in Ly- absorption towards the background quasar at a redshift of . These observations yielded an H i column density of (H i)/cm) and a metallicity of [M/H] (Rafelski et al., 2012).
The [C ii] 158 m observations come from Neeleman et al. (2017), who observed J0817+1351 on December 30, 2015 with ALMA in cycle 3 (project 2015.1.01564.S). The total on-source time was 46 minutes with a spatial resolution of and spectral resolution after smoothing of . [C ii] emission was found from the quasar at the redshift of the DLA. In addition to the detection of [C ii], the detected continuum emission at 158 m implies a SFR of M year. The maximum rotational velocity was determined from fitting an exponential disk to the data (see Neeleman et al. 2016), resulting in a dynamical mass limit of M.
Located in the COSMOS field, AzTEC/C159 () is an SMG discovered at 1.1 mm in the AzTEC/ASTE survey (Aretxaga et al., 2011). Smolčić et al. (2015) derived a dust temperature of K, SFR M year, and M. Smolčić et al. determined the dust mass using three techniques: fitting the NUV-mm SED with MAGPHYS ( M), fitting the FIR SED assuming an optically thin blackbody ( M), and fitting the FIR SED assuming the dust model of Draine & Li (2007) ( M). We will use the last dust mass in the following analysis.
Additional observations were made (E. Jimenez-Andrade et al., in preparation) using the VLA in D- and DnC-configurations and NOEMA in D- and C-configurations. These detected emission in CO(5-4) and CO(2-1), resulting in a molecular gas mass M (assuming K km s pc, appropriate for ULIRGs; Solomon & Vanden Bout 2005).
ALMA [C ii] observations for AzTEC/C159 (project 2012.1.00978.S; A. Karim et al., in preparation) used 34 antennas, consisted of 20 minutes of on-source time, and were undertaken on June 15, 2014. These band 7 observations resulted in spatial resolutions of and a spectral resolution of MHz km s for the [C ii] line.
COSMOS J100054+02343 (J1000+0234; ) is a galaxy in the COSMOS field. Originally detected as a V-band dropout (Carilli et al., 2008), followup observations and SED fitting (visible to cm) revealed a possible major merger with a SFR M year and L (Capak et al., 2008). Observations of CO(4-3) and CO(2-1) suggest a gas mass of = M (assuming K km s pc), and a dynamical mass of M (Schinnerer et al., 2008).
ALMA [C ii] observations for J1000+0234 (project 2012.1.00978.S; A. Karim et al., in preparation) used 39 and 36 antennas on June 2 & June 14, 2014, resulting in a combined on source time of 44 minutes. These band 7 observations resulted in spatial resolutions of and a spectral resolution of MHz km s for the [C ii] line.
ULAS J131911.29+095051.4 (J1319+0950;)is a QSO detected in the UKIRT Infrared Deep Sky Survey (UKIDSS) (Mortlock et al., 2009). By fitting an optically thin greybody model with T K and to their FIR SED, Wang et al. (2011) estimate = L. Using this FIR luminosity and an altered form of equation 1 of Bianchi (2013), they find a dust mass of M.
It has been previously observed with cycle 0 ALMA (project 2011.0.00206.S) in [C ii] (Wang et al., 2013), where a velocity gradient was observed. They estimate the dynamical mass to be M at an inclination angle of (from the ratio of the [C ii] major and minor axes), using
where is the full width at half maximum of the [C ii] line, is the inclination angle, and is the radius of the galaxy. Note that this is simply assuming that the [C ii]- emitting material is rotating in a purely circular fashion, and that we will make this assumption as well.
More recently, ALMA [C ii] observations for J1319+0950 (project 2012.1.00240.S; Shao et al. 2017) used 34 antennas in August 2014, resulting in a combined on source time of 36 minutes. These band 7 observations resulted in spatial resolutions of and a spectral resolution of MHz km s for the [C ii] line.
The rotation of J1319+0950 has been fitted by Shao et al. (2017), using a similar method as adopted here. We will compare their results with ours, along with further analysis of the density profile.
3 Tilted Ring Fitting Methods
We began with ALMA [C ii] image cubes of each source. These cubes were converted from frequency to velocity space by assuming zero velocity at the redshifted frequency of the [C ii] line (1900.5369 GHz), where redshifts were taken from previously published [C ii] observations. Velocity fields were made from these cubes and run through the GIPSY task rotcur (van der Hulst et al., 1992) to find the rotation curve of the source. This rotation profile was used in GIPSY’s velfi to create a model velocity field, which we subtracted from a velocity field created from the data to test the goodness of fit.
3.1 Velocity Fields
As discussed by de Blok et al. (2008), multiple methods exist for determining the velocity behavior of a source. The most common is the first moment or velocity weighted image, implemented in CASA as immoments:
This is well suited for high S/N systems with strong lines, as the entire line profile contributes to each integral. However, if the line is weak, noise contributions may render any velocity information unreachable.
The other main methods examine the spectrum of each spatial pixel in an image cube. The velocity may be found from the velocity with the maximum amplitude (i.e., peak velocity field) or from the centroid velocity of a Gaussian fit to the spectral line. While these two approaches are identical in the case of strong lines, the latter is more appropriate for weaker systems. The fit may be be a single or double component fit, or may include an asymmetric term.
Due to the faintness of high redshift sources and the limited (albeit considerable) sensitivity of ALMA, our sources have relatively low S/N, rendering first moment maps useful only as rough indicators of velocity gradients. The best velocity fields were found using the single Gaussian fitter in the AIPS task xgaus. This interactive fitter allowed for rejection of poor fits, resulting in reliable velocity fields.
3.2 Ring Widths
Before we could begin fitting tilted ring models to these velocity fields, the possible dimensions of these rings had to be quantified. Due to the small angular sizes of these sources, we assumed that they could be approximated as flat, warp-less (i.e., constant inclination angle) disks. The simplest model would be a single ring spanning the entire velocity field, but this does not satisfy Nyquist sampling. Instead, we adopt two limits on ring width: an upper limit of the full width at half maximum (FWHM) of the minor axis of the synthesized beam divided by three and a lower limit of the FWHM of the minor axis of the synthesized beam divided by two. The resulting number of possible rings () and the maximum radial extent of each galaxy () are listed in Table 1.
Previous methods (e.g., de Blok et al. 2008; Richards et al. 2015) used a single number of rings to create rotation curves. Since we assume a warp-less disk, we consider the results of combinations of ring numbers. These ring limits are influenced by the range in Table 1 (i.e., to ), but are varied to include all possible permutations (e.g., 2, 3, 4, 2-3, 3-4, 2-4). In this way, we consider six sets of ring numbers for each source.
3.3 Rotation Fitting
The foundations of the approach described in this section are based on the methods of de Blok et al. (2008) and Richards et al. (2015), as outlined in Figure 1. They have been applied to derive rotation curves for The H i Nearby Galaxy Survey (THINGS; de Blok et al. 2008) and as a step in the dark matter mass decomposition of NGC5005 (Richards et al., 2015). Both this work and those upon which it is based assume purely circular rotation (i.e., ), which is considered a valid approximation for well-ordered disks. However, due to our low S/N sources, we will assume a constant inclination and position angle across each disk, rather than allowing for a light variation in each.
The velocity fields from xgaus were fed into the GIPSY task rotcur, which fit a series of tilted rings to the source. Each ring had a set radius () and width, but variable rotation velocity (), inclination angle (), position angle (), center (,), systemic velocity (), and expansion velocity (which we assumed to be zero). These are related to the observed line of sight velocity by
where is the position angle of a point (x,y) in the plane of the observed galaxy (i.e., not in the plane of the sky) with respect to the receding semimajor axis (de Blok et al., 2008):
By direct comparison of the input velocity field and the velocity field created by the model rings, rotcur fits for the best ring parameters.
The results are slightly dependent on the initial estimates of each parameter (central position, , , ). To counteract this, we used two methods of deriving estimates. The first was to obtain estimates of the position angle, inclination, and systemic velocity from the input velocity field by eye, then to run rotcur to verify their plausibility (RC1). The second method found the systemic velocity by fitting a 1-D Gaussian to the [C ii] line profile and the other parameters by fitting a 2-D Gaussian to the [C ii] zeroth moment image (GAUS). When deriving our final curves, we considered all combinations of parameters from each method.
In what follows, a ‘run’ begins with an estimate of central position, , , and , along with a set range of ring widths, and ends with a derived rotation curve and a set of best-fit parameters.
We begin by fixing the ring centers to be identically equal to an estimate (RC1 or GAUS), and setting a minimum and maximum number of rings. Letting all other variables change, we run rotcur for all current permutations of ring number, starting with a run using rings and then repeating with a run using rings, and so on until a run using rings.
The run results are aggregated, and unsuccessful rings are discarded. Success is defined as: , , and , where signifies uncertainty. An average of the systemic velocities of all successful rings is taken, weighted by their fitting uncertainties. Fixing to be this averaged value, we then repeat the fitting procedure, but fitting for and simultaneously. We then fix and to their derived averages.
In the final run, the only variable is . The set of successful rings with the smallest width are chosen as representative of the overall rotation. The uncertainty in rotational velocity for each ring is estimated as the maximum of the uncertainty by rotcur for that ring, the velocity width of a channel, and . In this way, we include the spectral resolution of the input cube, the goodness of fit from rotcur, and the effects of inclination uncertainty. Another common estimate of the uncertainty is the difference in rotational velocities for the receding and approaching halves of the galaxy (de Blok et al., 2008), but we lack the signal to perform this.
In order to test the goodness of fit, we use GIPSY’s velfi task to create a velocity field based on the rotation curve, , , and that were found from rotcur. This field, which is created using the same spatial resolution as our data, is then compared to the field of the original data. The residual velocity values for each pixel are added in quadrature and normalized by the number of pixels, resulting in a normalized error value.
As a summary example, consider HZ9, which had and . By changing the initial estimates for each variable (center position, , , ) between those of our two methods (i.e., by eye and by fits to moment images), we have possible combinations of initial estimates. There are three runs using a single number of rings (2, 3, 4), two runs using a range of two rings (2-3, 3-4), and one run using a range of three rings (2-4), resulting in six possible combinations. By varying both parameter estimates and ring number, we thus consider the results of full runs. Similarly, we considered 96 full runs for HZ10, J0817+1351, AzTEC/C159, and J1000+0234.
4 Results & Analysis
All six objects were well fit by the tilted ring models, with their resulting rotation curves varying in detail. Table 2 shows the resulting redshifts (found from ), position angles, inclinations, and positions of the best model for each source. For J1319+0950, both our results and those of Shao et al. (2017) are given.
For each source, the results are presented from our initial rotcur approach (RC1), fitting 1-D Gaussians to the line profile and 2-D Gaussians to the zeroth moment images (GAUS), and our full exploration of the parameter space (RC2). Note that the position angles are measured counter-clockwise from north, and numbers in parentheses represents uncertainty in the last digit. As a conservative estimate of the uncertainty in redshift, we took the greater of the implied uncertainty from from rotcur and the velocity width of a channel in the initial image cube.
It should be noted that the uncertainties in Table 2 for the position angle and inclination are those outputted by rotcur, not including uncertainties from errors from observational effects, and thus are underestimated. See the Appendix for a discussion of this error. Uncertainties in positions are the cell sizes of the original image cubes, while uncertainties in are half of the ring width.
Note. – Errors may be underestimated (see Appendix). a-See Shao et al. (2017) for details of fitting procedure.
Figure 2 shows the rotation curves derived with rotcur, velocity fields created from the data, and the model velocity fields using the RC2 results. For the best-fit rotation curves, the vertical error bars are the uncertainty in (see Section 3.3). The horizontal error bars represent the width of each ring. The red lines connecting each point are simply interpolations, while the dashed black lines are the better of a linear fit or a fit of the form .
4.1 Dynamical Mass
An estimate of each source’s dynamical mass could be found by assuming circular orbits, as we have done in the above analysis. This requires setting the centrifugal force equal to the gravitational attraction:
where is the velocity at radius and is the gravitational constant. Using the radius and velocity of the largest ring of our models results in the values shown in Table 3. Equation 5 may also be applied to each ring individually, resulting in the mass profiles shown in Figure 3.
In addition to our derived dynamical masses, we include previously found dynamical masses, dust masses, stellar masses, and H masses. The previous dynamical masses were found using equation 1, where the size and inclination angle was estimated from Gaussian fits to the integrated emission, assuming intrinsically circular disks.
Note. – References: a-Capak et al. (2015), b-Neeleman et al. (2017), c-Smolčić et al. (2015), d-E. Jiménez-Andrade et al. (in preparation), e-Schinnerer et al. (2008), f-Capak et al. (2008), g-Wang et al. (2013), h-Wang et al. (2011). All other values are from this paper.
We note that the uncertainties in the above masses are influenced only by the width of the ring and uncertainty in the velocity of the ring. While the velocity uncertainty does include the spectral resolution and rotcur uncertainty in & , there are still systematic uncertainties that we cannot estimate or account for. From an exploration of a model observation at different resolutions and S/N (see Appendix), we estimate a systematic uncertainty of dex, to be added to the uncertainties shown in Figure 3. When this uncertainty is included, our ranges are similar to those already published. While our technique did not result in more precise measurements of the dynamical mass, it should be noted that our procedure is distinct from that used to find masses for these sources previously.
While the velocity gradient observed in the data velocity field is replicated in the model, the two velocity fields show different East-West extents. In particular, while the other sources in our sample show elongation along the axis perpendicular to their isovelocity lines, HZ9 is nearly circular. This difference may imply that HZ9 is a nearly face-on disk that is rotating at a considerable rate, or that the distribution of emission is affected by our low resolution. The restoring beam was , resulting in square cells of . As seen in Figure 2, the source only extends from its center, so few cells were used as input for each ring. To correct for this effect, models must be convolved with the observing beam and compared to the data on a channel-by-channel basis (see Section 5 for implementation of this method).
The difference in spatial extents between the model and data velocity fields is to blame for the disparity between the position angles derived by rotcur and zeroth moment imaging ( deviant). No combination of assumptions on other variables or variations of ring numbers produced a reliable rotation curve using the derived from Gaussian fitting, implying that it does not relate to the rotation. That the inclination angle estimates are similar (within ) may act as marginal evidence that the matter is well organized into a rotating disk. This is supported by the agreement between the estimates of the center of the galaxy and the linear rotation curve.
Similarly to HZ9, the linear rotation curve suggests solid rotation. However, the gradient in the data is not fully reproduced by the model fit. In particular, the southwest region of the data shows a high-velocity area that is not aligned with the eastern velocity gradient.
Observations of Pavesi et al. (2016) show that the centers of [C ii] and [N ii] emission are offset both spatially and in velocity. One possible interpretation of this is an active merger, which we lack the resolution to account for here. This would partly explain why the best fit to the velocity field is dominated by the western velocity gradient of the galaxy, as the east could be galaxy interacting with the main galaxy.
While the effects of observing beam shape are less severe than in HZ9, the model velocity field again does not trace the integrated emission well. However, the inclination angles, positions, and redshifts agree between the rotcur and zeroth moment results.
There is slight evidence for a flat rotation curve, and the model fits the general trend of the velocity field well. While a standard rotating galaxy shows a ‘butterfly’-like velocity field, with isovelocity lines radiating from its center, J0817+1351 simply shows a velocity gradient. The possible ripple in this field could be caused by a warp, but we lack the resolution to confirm this. While the systemic velocity, position angle, and inclination found through RC1 are more applicable than the GAUS estimates, the position from GAUS (only pixels different) returns a better fit.
The slight bend of the rotation curve at large radii may be interpreted as weak evidence for detection of disk-like rotation, but is more likely due to the non-uniform appearance of the northern and southern portions of the input velocity field.
AzTEC/C159 shows the most likely evidence for a flat rotation curve at large radius. The use of a zeroth moment image here is not very helpful, as the integrated emission is nearly circular (). This introduces a large uncertainty in both inclination angle and position angle. However, all derived quantities agree with those found with rotcur. Through the variable combination analysis, the position angle found through the zeroth moment map was found to be inapplicable to the velocity field.
Since we have estimates of the stellar, dust, and molecular gas masses (Table 3 and references therein), we may estimate the total dark matter content () of the inner kpc by assuming
where we are neglecting the mass of the ionized medium () and cold medium (), for which we have no estimates. By inverting this equation, we find that M. The implied baryon fraction ratio (M/M) falls within the low end of the scatter of values found for local spiral galaxies by Richards et al. (2016).
The velocity field of J1000+0234 is unique in this sample in that it features parallel isovelocity lines near its center and two spatially extended regions of high velocity gas at its far edges. That the central lines are parallel suggests that the gas velocity near the center rises linearly with increasing radius, as seen in Figure 2. The best fit model to this velocity field is a steep linear rise near the galaxy center, followed by a flattened rotation profile. The flat portion of the curve produces the “spider diagram”-like radiating isovelocity lines in the lobes, which are not seen in the data.
The spatial distribution of emission and CO excitation level of this source has been interpreted as evidence for it being two galaxies in the midst of merging (Capak et al. 2008; Schinnerer et al. 2008). While this would partially explain the appearance of the velocity field (namely, the lack of radiating isovelocity lines and the small velocity gradient across the northern and southern lobes), the apparent uniformity of the central gradient is evidence for this system being a single rotating galaxy.
Similarly to J1000+0234, J1319+0950 shows a local velocity decrease at the penultimate point in its rotation curve. Due to the non-elliptical shape of each galaxy, the outer rings are fit only to a few data points, instead of the full annuli available to the inner rings. However, as seen in Figure 3, these local velocity minima do not translate to unphysical local dynamical mass minima.
Shao et al. (2017) fit a tilted ring model to this galaxy in a similar fashion, resulting in nearly identical parameters. In addition, they found a dynamical mass of log, which agrees with the mass found through our analysis. The effects of slightly differing inclination angles ( vs ) can be seen in the vertical offset in dynamical mass points in Figure 3.
5 galmod Modeling
Our above approach has several drawbacks. One is that we cannot reliably account for the velocity dispersion () of these objects. rotcur does not allow the user to provide a prior estimate, fix the value, or include its effects in the model. Another limiting property is that we cannot control the density profile of the model. It is possible to estimate a surface density profile from the mass profile. However, this cannot be transformed into a physical mass density without the additional assumption of a vertical mass structure.
Another drawback is that we do not account for beam smearing. This effect, which is strongest for low resolution observations, acts to make strong velocity gradients appear less severe (Begeman, 1987; de Blok & McGaugh, 1997; O’Brien et al., 2010; Kamphuis et al., 2015). As shown in Figure 2 of de Blok & McGaugh (1997), beam smearing alters the nature of the inner rotation curve but leaves the outer flat section relatively unchanged. Thus, while the underlying rotation may have a steeper initial rise, our detection of a signature of dark matter is still valid.
All of these drawbacks, in addition to possible deconvolution effects, can be accounted for by using the GIPSY task galmod. This takes in a radial mass density profile, scale height (z), rotation curve, velocity dispersion, position angle, and inclination angle, outputting a model cube that matches the geometry (i.e., cell size, pointing center, channel width) of an input data cube. If the data cube has arbitrarily high resolution, the model cube may then be convolved with the restoring beam of an observation (e.g., Schulman et al. 1996), and the observed and model cubes may be compared on a channel-by-channel basis, via a PV-diagram, or a global profile. If a sequence of such model cubes were created, spanning a range of density profiles, etc., they could be used to constrain the properties of the observed data. Alternatively, the Tilted Ring Fitting Code (TiRiFiC; Józsa et al. 2007) could be used to fit a model to the cube automatically. In addition to direct comparisons of spectral profiles, the cubes created by galmod may be mock observed using the CASA task simobserve. These observations may be used to test the validity of our spectral profiles and rotation curves.
It should be noted that rotcur has the advantages of speed and simplicity. However, these come at the price of forced assumptions and limited fitting ability. The use of galmod allows for the simultaneous exploration of a large number of variables (e.g., Yim et al. 2014) and can be used to explore the effects of beam smearing (e.g., Frank et al. 2016; Caldú-Primo et al. 2013), but requires manual construction of a set of data cubes with different parameters, rather than automatically fitting for them. While TiRiFiC includes a data cube-based fitting routine, initial applications to our low S/N data resulted in poor fits. Thus, using rotcur to find the basic parameters of a source and then galmod to fit for other parameters while accounting for resolution effects seems to be the ideal path for exploring our datasets.
To explore this avenue, we created a sequence of model cubes based on our data of AzTEC/C159 and J1319+0950, which were best fit by the rotcur process. galmod allows for the specification of a large number of parameters, so we adopted the rotation curves, position angles, central positions, and inclination angles found through our prior analysis to reduce any redundancies in our fitting process. Initial testing showed that results depended weakly on the input scale height and density profile, so these were assumed to be thin (maximum radius; Lang et al. 2017) and flat (). Using these parameters, we created cubes with different velocity dispersions ( to 145 km s).
The CASA task simobserve was then used with these galmod cubes as input to simulate the ALMA cycle 1 observations of A. Karim et al. (in preparation) of AzTEC/C159 and the ALMA cycle 1 observations of Shao et al. (2017) of J1319+0950. It should be noted that the simulated and actual observations had slightly different beam sizes and channel sizes, but are still comparable.
In this way, we created a series of noiseless cubes and generated spectra of each, integrating over the entire galaxy. By comparing the spectra on a channel-by-channel basis with the spectrum of our original data cube, we calculated the goodness of fit.
The resulting goodness of fit for each source as a function of velocity dispersion is shown in Figure 4.
For AzTEC/C159, models with small velocity dispersions were obviously favored. We adopt 15 km s (which is the only minimum) as the velocity dispersion and compare the spectra of the model and data. While the model traces the general profile of the data, multiple features (e.g., the km s and km s components) are not replicated. The km s component is likely noise, as no significant emission is detected in its corresponding channels in the data cube. Instead, multiple low-level, diffuse noise sources are present.
In the case of J1319+0950, the dependence of the fitting error on velocity dispersion shows a significantly different behavior, with a global minimum at 135 km s. The width of the model spectrum is comparable to that of the data, although the model does not replicate the observed asymmetric horn profile, which may only be noise.
We note that a full Bayesian approach (e.g., MultiNest; Feroz & Hobson 2008) is preferable to our discrete exploration of the parameter space, but it is too computationally expensive at this time, given the low S/N of our data (i.e., runs take on the order of a day).
6.1 Accuracy of Fits
Using the method of tilted ring fitting, we are able to provide estimates for the position, position angle, inclination angle, and systemic velocity/redshift of a galaxy as a whole, in addition to its circular rotation velocity as a function of radius. If these parameters are used in the GIPSY task galmod, then we may also constrain the velocity dispersion, scale height, and density profile of the disk. However, since we are using these methods to fit relatively low S/N sources (i.e., with respect to observations of local galaxies that the methods were created to fit), some variables are not precisely constrained. While we will present a detailed analysis of systematic uncertainties and limitations in the process in a future paper (G. Jones, Y. Shao, et al., in preparation), we present our initial analysis in the Appendix.
To summarize, we are confident in our ability to determine the redshift, central position, and position angle of each source. However, the ambiguity between the rotational velocity and inclination of a disk makes the simultaneous determination of both quantities difficult. If the inclination angle could be found using a different technique (e.g., high resolution optical imaging, which is difficult for high-, dusty galaxies), then this ambiguity could be removed, resulting in higher quality dynamical fits. In addition, the quality of the data (i.e., S/N, resolution, channel width) appears to be of less importance than the existence of an intrinsic rotating disk for this work. That is, while very poor quality data will return poor fits, both medium and high quality data will be nearly equally successful.
6.2 Comparison With Other Techniques
The dynamical characterization of galaxies with ALMA and [C ii] is well complemented by recent work with IFU observations of H- emissions of galaxies (Lang et al. 2017; Genzel et al. 2017). The large sample size of Lang et al. allows for stacking of rotation curves, revealing a Keplerian-like decline of rotation speed at larger radii. Similar curves are detected for each of the bright sources in Genzel et al. Comparison with dark matter models suggests a high baryon fraction and pressure support in these outer disks. Since the outer portion of the rotation curves in this study are only tentatively detected, we cannot confirm this.
In addition, rotation curves have already been found for ALMA [C ii] observations of galaxies, using other techniques. De Breuck et al. (2014) fit the velocity field of ALESS 73.1, a SMG, using a Bayesian implementation of the KINematic Molecular Simulation (KinMS) of Davis et al. (2013), a simple arctan rotation curve model, and a dynamical model that accounts for beam smearing and assumes an exponential density profile. All three models returned good fits, resulting in a rotation curve and dynamical mass of M. Since KinMS creates data cubes for mock observations, with given surface brightness, rotation curve, scale height, velocity dispersion, pixel size, channel width, etc., it acts similarly to galmod, but is implemented in IDL.
Even if observations of high redshift sources are not well fit by tilted ring models or general dynamical codes, their dynamical masses may still be estimated using forms of equation 1. In addition to some of the sources in our sample, this widespread method has been applied to the SMGs GN20.2a & GN20.2b (Hodge et al., 2013), the SMG of BRI1202-0725 (Carilli et al., 2013), the LBGs CLM1 & WMH5 (Willott et al., 2015), etc..
Thus, while our technique is not without precedent in the universe, it acts as a suitable complement to those procedures already in place.
We have fit tilted ring dynamical models to [C ii] 158m line observations of two MS galaxies (HZ9 & HZ10), a QSO host galaxy (ULAS J1319+0950), two SMGs (AzTEC/C159 & COSMOS J1000+0234), and a DLA host galaxy (ALMA J0817+1351), all at . Our dynamical analysis is one of the the first attempts to go beyond very simple conclusions based on fitting integrated line profiles for galaxies, using the spatial information inherent in the interferometric observations. The three low-luminosity galaxies (HZ9, HZ10, and ALMA J0817+1351) show linear velocity gradients, consistent with either rotation or other dynamical models, and limited by signal-to-noise and resolution. The three luminous galaxies (ULAS J1319+0950, AzTEC/C159, and J1000+0234) exhibit possible evidence for flattening of the velocity field at large radii, more suggestive of rotating disk galaxies. In the case of AzTEC/C159, evidence for M of dark matter was found.
In addition to fitting models to velocity fields, we created a sequence of model data cubes spanning a range of possible velocity dispersions for AzTEC/C159 and J1319+0950. These high resolution model cubes were then observed using the CASA task simobserve, and then compared to the integrated line profile.By comparing our data to these models, we found evidence for low ( km,s) and moderate ( km s) velocity dispersions.
In an effort to quantify the accuracy of our technique, we tested our code on mock data cubes at two levels of resolution and sensitivity. From this, we learned that high sensitivity and resolution result in better fits, but the effects of beam smearing are still evident. Because of the above uncertainties, our conclusions are broad. For each of the three MS galaxies (HZ9, HZ10, and J0817+1351), the mainly linear rotation curves are likely non-physical and actually caused by beam smearing effects. On the other hand, the two SMGs and the QSO host galaxy in this sample (AzTEC/C159, J1000+0234, J1319+0950), show stronger evidence for a flattened rotation curve. A true interpretation of these results must wait until higher resolution observations are performed, allowing us to truly probe the dynamical nature of these sources. However, these initial results are promising.
While full dynamical characterization or rotation curve decomposition are not yet possible, these first results using ALMA data (taken mainly with short integration times and relatively few antennas) are encouraging. As ALMA attains full capability, high resolution, higher signal-to-noise observations of the [C ii] 158m line will allow for detailed determination of the galaxy dynamical masses at the highest redshifts. In parallel, large new facilities in the near-IR will determine the stellar profile of early galaxies, while future facilities such as the Next Generation Very Large Array, will determine gas density profiles at high resolution (Carilli et al. 2015; McKinnon et al. 2016; Selina & Murphy 2017). Together, these facilities will allow for a full analysis of the baryonic and dark matter components of galaxies out to cosmic reionization.
In an effort to quantify the reliability of our fitting routine, we apply it to a series of simulated observations of a well-studied galaxy. The key to this analysis is that all six of these observations are based the same galaxy, with identical intrinsic dynamics and morphology. Thus, this exploration allows us to test the fitting routine without including scatter from differing source properties.
Appendix A Model Creation
To construct our sky models, we use observations of M51 (NGC5194) in CO() emission taken by Helfer et al. (2003) with the BIMA interferometer between 1997 and 1999. These data were transformed to represent emission from a galaxy with a similar velocity width and flux density as J08171351
The mock datasets are created using two different resolutions: (low, L) and (medium, M). Observations of each were simulated using both 3 hours on source (H3) and 30 hours on source (H30). In total, we here consider four simulated observations, exploring two levels of sensitivity and two levels of resolution. Each mock observation is then folded through our analysis in a manner identical to the observed sources in the above paper.
This source presents a challenge, as it exhibits both variation in inclination angle ( to ; Oikawa & Sofue 2014) and distinct spiral arms (Shetty et al., 2007), both of which are unaccounted for in our code. However, this example is more useful than an idealized test (e.g., using a galmod dataset), as we are able to determine how these non-regular features affect our results at different resolutions and sensitivities. In addition, our models at most extend to kpc in our adjusted frame, or kpc in the frame of M51, which is an area of relatively small inclination variation (Figure 5 of Oikawa & Sofue 2014).
Appendix B Results
Appendix C Discussion
We now consider the effects of resolution and sensitivity on the fit parameters, rotation curves, and mass profiles.
All models yield similar redshifts, which are found by directly converting the best-fit systemic velocity. This parameter, which influences a velocity field as an overall shift, is therefore derivable using even poor data. Note that the first channel in the cube (i.e., crval3) was set as the rest frequency of [C ii] at , so the redshifts here are all less than . Similarly to the redshifts, the central position of each model galaxy agrees.
The position and inclination angles show more variability between datasets. Both high sensitivity position angle estimates agree with from Regan et al. (2001) to within , while the two low sensitivity datasets are discrepant. Thus, our determination of position angle is weakly dependent on sensitivity, but even low sensitivity observations may yield semi-accurate position angles.
On the other hand, the inclination of each dataset varied from to . Regan et al. (2001) state an inclination angle of , while Oikawa & Sofue (2014) use . Since an ambiguity exists between inclination angle and circular velocity (i.e., ), this spread in inclination angles reflects a similar spread in rotation velocity.
We have scaled the rotation curve of Oikawa & Sofue (2014) (their Figure 1) and plotted it over our mock rotation curves (Figure 5). In all cases, the effects of beam smearing are obvious. The inner, steep rise of the scaled M51 rotation curve is not traced, and the initial rise modeled rotation curve is less severe. The only curve that agrees is that of H30M. Discrepancies between the other curves are caused by the scatter in their inclination angles, especially the high sensitivity, low resolution case (). This uncertainty in rotation curve scale extends to an uncertainty in dynamical mass.
If this method was unaffected by the sensitivity and resolution of a given observation, then the mass profiles of Figure 6 would be coincident. Instead, the two low-resolution datasets are outliers, and the large-radius dependence of the other models is disparate. Interestingly, the two medium-resolution datasets are highly similar. However, H3M shows a decrease in dynamical mass with increased radius, which is unphysical and suggests that the model should only be trusted at the inner radii.
We have also applied the technique of Shao et al. (2017) to the six datasets, resulting in mainly the same best-fit parameters. The only parameters that differed between our two analyses were the inclination angles of the two low resolution datasets (i.e., H3L & H30L). Namely, the procedure of Shao et al. finds for H3L (the above procedure yields ) and for H30L (). However, since the low resolution datasets show strange mass profiles, this difference in inclination angles may only signify that the underlying datasets are difficult to model, but the two analysis methods are in agreement.
Appendix D S/N Effects
The quality of each fit may dependent on multiple characteristics, including the spatial resolution of the disk, the S/N of the detection, the number of channels in which line emission is present, and whether or not the source is in reality an ordered, rotating disk. While we cannot control the last property, it is possible to compare the others. In an attempt to identify which observation characteristics contribute to a good fit, we tabulate the fractional resolution of the source (HWHM of the minor axis of the restoring beam), the S/N of the [C ii] detection, and the fractional velocity resolution (i.e., the number of channels with line emission).
In Table 5, we also estimate the ‘Fit Quality’ of each source, where 1 represents the best fit and 4 the worst. Since known rotation curves for our other sources do not exist, their quality values are qualitatively based on the agreement of the model and data velocity fields. Az159 shows excellent agreement, while J0817+1351, J1000+0234, and J1319+0950 show slight deviations. Only portions of HZ9 and HZ10 are well fitted, but HZ9 is also poorly spatially sampled. The modeled sources were simply ordered by how well their rotation curves agreed with that of Oikawa & Sofue (2014).
The spatial and spectral resolution of each source do not show obvious trends with the fit quality. Curiously, both HZ9 and J0817+1351 show similar spatial fractional resolution, but since J0817+1351 has twice the number of cells per HMWM (8.8 vs 4.2 for HZ9), it appears better resolved. From this analysis, it would appear that poor fits will be returned if any one of these characteristics is weak (e.g., poor spatial fractional resolution of HZ9/H30L).
In summary, our analysis is best applied to medium resolution, high sensitivity data of orderly rotating disks. Intermediate resolution allows for characterization of the galaxy’s rotation as a whole, while minimizing under-resolution effects. Higher sensitivity presents more of the outer disk, while improving the results in the easily detected inner disk. The central positions, position angles, and systemic velocities are accurately determined, even in the case of low resolution/sensitivity observations. The ambiguity between rotational velocity and inclination angle introduces an uncertainty in inclination and a scaling factor in the rotation curve and mass profile. This may be corrected for by determining the inclination angle separately (e.g., using axis ratios of resolved observations).
- affiliation: Physics Department, New Mexico Institute of Mining and Technology, 801 Leroy Pl, Socorro, NM 87801, USA; firstname.lastname@example.org
- affiliation: National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA
- affiliation: National Radio Astronomy Observatory, 1003 Lopezville Road, Socorro, NM 87801, USA
- affiliation: Cavendish Astrophysics Group, University of Cambridge, Cambridge, CB3 0HE, UK
- affiliation: Department of Astronomy, School of Physics, Peking University, Beijing 100871, China
- affiliation: Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China
- affiliation: Department of Astronomy, School of Physics, Peking University, Beijing 100871, China
- affiliation: Kavli Institute for Astronomy and Astrophysics, Peking University, Beijing 100871, China
- affiliation: Infrared Processing and Analysis Center, California Institute of Technology, MC 100-22, 770 South Wilson Ave., Pasadena, CA 91125, USA
- affiliation: Spitzer Science Center, California Institute of Technology, Pasadena, CA 91125, USA
- affiliation: Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA
- affiliation: Department of Astronomy, Cornell University, Space Sciences Building, Ithaca, NY 14853, USA
- affiliation: Argelander-Institut für Astronomie, Universität Bonn, Auf dem Hügel 71, D-53121 Bonn, Germany
- affiliation: University of California Observatories-Lick Observatory, University of California, Santa Cruz, CA, 95064, USA
- affiliation: Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany
- See https://casaguides.nrao.edu/index.php/Simulation_Recipes for general procedure.
- Aretxaga, I., Wilson, G. W., Aguilar, E., et al. 2011, MNRAS, 415, 3831
- Babcock, H. W. 1939, Lick Observatory Bulletin, 19, 41
- Banerji, M., Carilli, C. L., Jones, G., et al. 2017, MNRAS, 465, 4390
- Begeman, K. G. 1987, Ph.D. Thesis
- Bianchi, S. 2013, A&A, 552, A89
- Bosma, A. 1978, Ph.D. Thesis
- Burbidge, E. M., Burbidge, G. R., & Prendergast, K. H. 1960, ApJ, 131, 282
- Burkert, A., Förster Schreiber, N. M., Genzel, R., et al. 2016, ApJ, 826, 214
- Caldú-Primo, A., Schruba, A., Walter, F., et al. 2013, AJ, 146, 150
- Capak, P., Carilli, C. L., Lee, N., et al. 2008, ApJ, 681, L53
- Capak, P. L., Carilli, C., Jones, G., et al. 2015, Nature, 522, 455
- Carilli, C. L., Lee, N., Capak, P., et al. 2008, ApJ, 689, 883
- Carilli, C. L., Riechers, D., Walter, F., et al. 2013, ApJ, 763, 120
- Carilli, C. L., & Walter, F. 2013, ARA&A, 51, 105
- Carilli, C. L., McKinnon, M., Ott, J., et al. 2015, arXiv:1510.06438
- Davis, T. A., Alatalo, K., Bureau, M., et al. 2013, MNRAS, 429, 534
- de Blok, W. J. G., & McGaugh, S. S. 1997, MNRAS, 290, 533
- de Blok, W. J. G., Walter, F., Brinks, E., et al. 2008, AJ, 136, 2648-2719
- de Blok, W. J. G., Walter, F., Smith, J.-D. T., et al. 2016, AJ, 152, 51
- De Breuck, C., Williams, R. J., Swinbank, M., et al. 2014, A&A, 565, A59
- Draine, B. T., & Li, A. 2007, ApJ, 657, 810
- Epinat, B., Tasca, L., Amram, P., et al. 2012, A&A, 539, A92
- Feroz, F., & Hobson, M. P. 2008, MNRAS, 384, 449
- Frank, B. S., de Blok, W. J. G., Walter, F., Leroy, A., & Carignan, C. 2016, AJ, 151, 94
- García-Burillo, S., Combes, F., Usero, A., et al. 2014, A&A, 567, A125
- Genzel, R., Tacconi, L. J., Eisenhauer, F., et al. 2006, Nature, 442, 786
- Genzel, R., Schreiber, N. M. F., Übler, H., et al. 2017, Nature, 543, 397
- Goldsmith, P. F., Langer, W. D., Pineda, J. L., & Velusamy, T. 2012, ApJS, 203, 13
- Helfer, T. T., Thornley, M. D., Regan, M. W., et al. 2003, ApJS, 145, 259
- Hodge, J. A., Carilli, C. L., Walter, F., et al. 2012, ApJ, 760, 11
- Hodge, J. A., Carilli, C. L., Walter, F., Daddi, E., & Riechers, D. 2013, ApJ, 776, 22
- Józsa, G. I. G., Kenn, F., Klein, U., & Oosterloo, T. A. 2007, A&A, 468, 731
- Kamphuis, P., Józsa, G. I. G., Oh, S.-. H., et al. 2015, MNRAS, 452, 3139
- Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511
- Lang, P., Förster Schreiber, N. M., Genzel, R., et al. 2017, arXiv:1703.05491
- Langer, W. D., Velusamy, T., Pineda, J. L., et al. 2010, A&A, 521, L17
- Langston, G., Conner, S., Lehar, J, Burke, B, Weiler, K. 1990, Nature, 1990, 344, 43
- Leauthaud, A., Massey, R., Kneib, J.-P., et al. 2007, ApJS, 172, 219
- McKinnon, M., Carilli, C., & Beasley, T. 2016, Proc. SPIE, 9906, 990627
- Mortlock, D. J., Patel, M., Warren, S. J., et al. 2009, A&A, 505, 97
- Neeleman, M., Prochaska, J. X., Zwaan, M. A., et al. 2016, ApJ, 820, L39
- Neeleman, M., Kanekar, N., Prochaska, J. X., et al. 2017, Science, 355, 1285
- Oikawa, S., & Sofue, Y. 2014, PASJ, 66, 77
- O’Brien, J. C., Freeman, K. C. & van der Kruit, P. C. 2010, A&A, 515, A61
- Pavesi, R., Riechers, D. A., Capak, P. L., et al. 2016, arxiv:1607.02520
- Pease, F. G. 1918, Proceedings of the National Academy of Science, 4, 21
- Persic, M., & Salucci, P. 1995, ApJS, 99, 501
- Pineda, J. L., Langer, W. D., Velusamy, T., & Goldsmith, P. F. 2013, A&A, 554, A103
- Planck Collaboration, Adam, R., Ade, P. A. R., et al. 2016, A&A, 594, A1
- Rafelski, M., Wolfe, A. M., Prochaska, J. X., Neeleman, M., & Mendez, A. J. 2012, ApJ, 755, 89
- Regan, M. W., Thornley, M. D., Helfer, T. T., et al. 2001, ApJ, 561, 218
- Richards, E. E., van Zee, L., Barnes, K. L., et al. 2015, MNRAS, 449, 3981
- Richards, E. E., van Zee, L., Barnes, K. L., et al. 2016, MNRAS, 460, 689
- Riechers, D. A., Walter, F., Carilli, C. L., Bertoldi, F., & Momjian, E. 2008, ApJ, 686, L9
- Riechers, D. A., Carilli, C. L., Capak, P. L., et al. 2014, ApJ, 796, 84
- Rogstad, D. H., Lockhart, I. A., & Wright, M. C. H. 1974, ApJ, 193, 309
- Rubin, V. C., & Ford, W. K., Jr. 1970, ApJ, 159, 379
- Schinnerer, E., Carilli, C. L., Capak, P., et al. 2008, ApJ, 689, L5
- Schmidt, P., Józsa, G. I. G., Gentile, G., et al. 2014, A&A, 561, A28
- Selina, R. & Murphy, E. 2017, ngVLA Memo 17
- Schulman, E., Bregman, J. N., Brinks, E., & Roberts, M. S. 1996, AJ, 112, 960
- Shao, Y., Wang, R., Jones, G. C., et al. 2017, arXiv:1707.03078
- Shetty, R., Vogel, S. N., Ostriker, E. C., & Teuben, P. J. 2007, ApJ, 665, 1138
- Smolčić, V., Karim, A., Miettinen, O., et al. 2015, A&A, 576, A127
- Sofue, Y., & Rubin, V. 2001, ARA&A, 39, 137
- Solomon, P. M., & Vanden Bout, P. A. 2005, ARA&A, 43, 677
- Tacconi, L. J., Genzel, R., Neri, R., et al. 2010, Nature, 463, 781
- Tacconi, L. J., Neri, R., Genzel, R., et al. 2013, ApJ, 768, 74
- Ueda, J., Iono, D., Yun, M. S., et al. 2014, ApJS, 214, 1
- van Starkenburg, L., van der Werf, P. P., Franx, M., et al. 2008, A&A, 488, 99
- Venemans, B. P., Walter, F., Zschaechner, L., et al. 2016, ApJ, 816, 37
- Vogt, N. P., Forbes, D. A., Phillips, A. C., et al. 1996, ApJ, 465, L15
- Wang, R., Wagg, J., Carilli, C. L., et al. 2011, AJ, 142, 101
- Wang, R., Wagg, J., Carilli, C. L., et al. 2013, ApJ, 773, 44
- Wang, J., Fu, J., Aumer, M., et al. 2014, MNRAS, 441, 2159
- Willott, C. J., Carilli, C. L., Wagg, J., & Wang, R. 2015, ApJ, 807, 180
- Wright, S. A., Larkin, J. E., Barczys, M., et al. 2007, ApJ, 658, 78
- Wuyts, S., Förster Schreiber, N. M., Wisnioski, E., et al. 2016, ApJ, 831, 139
- van der Hulst, J. M., Terlouw, J. P., Begeman, K. G., Zwitser, W., & Roelfsema, P. R. 1992, Astronomical Data Analysis Software and Systems I, 25, 131
- Yim, K., Wong, T., Xue, R., et al. 2014, AJ, 148, 127