Evidence of a Flat Outer Rotation Curve in a Star-bursting Disk Galaxy at z = 1.6

In the standard ΛCDM model of the universe, the baryonic and dark matter components of galaxies grow in tandem from high to low redshift (e.g., White & Rees 1978). Galactic rotation curves are an important observational probe of these matter distributions in galaxies (see Sofue & Rubin 2001, for a brief review). Through Keplerian inference, rotation curves provide a direct method for measuring the mass enclosed within a galactocentric radius. When combined with a light profile and a mass-to-light ratio, one may estimate the distribution of dark matter in a galaxy (e.g., Carignan & Freeman 1985; van Albada et al. 1985; Kassin et al. 2006; de Blok et al. 2008; Korsaga et al. 2018).

At low redshift (z), disk galaxy rotation curves are seen to rise sharply in the inner few kiloparsecs and then flatten to a near-constant velocity (e.g., Rubin et al. 1978; McGaugh 2016). This is caused by baryons dominating the mass density in the inner regions of galaxies while dark matter dominates in the outer regions (Kassin et al. 2006; Courteau et al. 2007; Dutton et al. 2007, 2013; Barnabè et al. 2012; Cappellari et al. 2013; Martinsson et al. 2013a, 2013b; Courteau & Dutton 2015). Recent observational studies have begun to push galaxy rotation curve measurements to larger radii at z ≳ 1 and while some find flat rotation curves (e.g., Di Teodoro et al. 2016, 2018; Xue et al. 2018), others, intriguingly, find evidence of decreasing outer-galaxy rotation curves (Genzel et al. 2017; Lang et al. 2017). These studies attribute this to a combination of decreased dark matter in galaxies at intermediate redshifts as well as increased pressure support in the outer disks. It is common for low redshift disk galaxies to show an initial decline before flattening in the outer disk (e.g., Rubin et al. 1965, 1978; McGaugh 2016), but the low or negligible inferred dark-matter-to-total-mass fractions in the Lang et al. and Genzel et al. samples imply that these initial declines will not begin to flatten if detected to higher radii.

Some studies find lower dark matter fractions in galaxies at z ≳ 1 than at z = 0 (Förster Schreiber et al. 2009; Burkert et al. 2016; Price et al. 2016; Wuyts et al. 2016) while others find no evolution (Di Teodoro et al. 2016; Stott et al. 2016; Yuan et al. 2017; Tiley et al. 2018). Many of these works rely on inverting scaling relations with significant scatter, such as the Kennicutt–Schmidt relation (Kennicutt 1998), in order to estimate gas masses. This may be the source of disagreement. There is much greater agreement among many of these studies when comparing the stellar to the total masses of the galaxies (See also Pelliccia et al. 2017).

Different simulations also make different predictions about the shapes of rotation curves at z > 1. Teklu et al. (2018) find 38% of disk galaxies with M_⋆ > 5 × 10¹⁰ M_⊙ at z = 2 in the cosmological hydrodynamic simulation, Magneticum Pathfinder (K. Dolag et al. 2019, in preparation) have declining rotation beyond ∼4 kpc (as estimated from Figure 4 in Teklu et al. 2018). In contrast, the IllustrisTNG cosmological hydrodynamic simulations (Lovell et al. 2018) show initially declining rotation curves at z = 2 that gradually flatten in the outer disk, similar to low-z disk galaxies. It is clear that more observations of the outer rotation curves of high redshift galaxies and their dark matter fractions are needed to puzzle out the evolution of dark matter fractions with time.

In this paper we present the serendipitously acquired rotation curve for a known dusty star-forming galaxy (DSFG) at z = 1.555, detected in Hα to a high galactocentric radius of 14 kpc. The galaxy exhibits a flat rotation curve, well-fit by an arctangent profile and a dark matter fraction at the half light radius typical of disk galaxies at z = 0. The paper is organized as follows. Section 2 details the observations and data reduction, Section 3 discusses the kinematic modeling, Section 4 details the dynamical and stellar mass estimations, and Section 5 presents our conclusions. Throughout this work we adopt a Planck cosmology with H₀ = 67.7 km s⁻¹ Mpc⁻¹ and Ω_Λ = 0.6911 (Planck Collaboration et al. 2016). We assume a Chabrier initial mass function (IMF; Chabrier 2003).

The data were obtained with the Multi-Object Spectrometer For Infra-Red Exploration (MOSFIRE; McLean et al. 2010, 2012) on Keck I as part of an observational campaign to follow up a sample of 450 μm and 850 μm identified DSFGs (for sample details, see Casey et al. 2013, 2017). DSFGs comprise the most intense star-bursting galaxies in the universe (see reviews by Blain et al. 2002 and Casey et al. 2014).

Observations were made on 2013 December 31 with the full width at half maximum (FWHM) of seeing between 08 and 09, as measured by the telescope focus routine, MIRA. We adopt an average seeing of 085 throughout this work. The galaxy, DSFG850.95, was observed in the H band for a total integration time of 1920s with the slit width set to 07. A 15 ABBA nod pattern was used and the slit orientation was randomly placed with respect to our target due to position angle optimization driven by multiplexing. Additional observational details are described for the full parent sample in Casey et al.(2017).

The data were reduced with the MOSPY data reduction pipeline.⁵ Figure 1 plots the reduced 2D signal-to-noise spectrum of DSFG850.95 showing Hα, [N ii], and continuum emission. An unrelated foreground galaxy at photometric redshift of 0.7, COSMOS 2004241 (Capak et al. 2007), sits a projected ∼35 to the south of DSFG850.95 and is seen in faint continuum emission in Figure 1 at spatial extent −30 kpc. It does not affect the analysis of DSFG850.95's rotation curve.

Zoom In Zoom Out Reset image size

The left panel of Figure 2 shows H-band imaging (rest frame ∼6400 Å), MOSFIRE slit, and ALMA 870 μm dust continuum contours of DSFG850.95. The H-band FWHM of seeing was comparable to that of the MOSFIRE observations and is shown as the dashed circle. The green ellipse shows the H-band half light radius. There is a misalignment between the position angles of the H-band exponential disk fit and the slit of 23° ± 4°, which should not strongly affect our measured radial velocities. We will discuss this further in Section 3.1. The right panel of Figure 2 shows Hubble imaging in the F814W filter, at a rest-frame effective wavelength of ∼3200 Å. The centroid of UltraVISTA K-band emission was used for the slit targeting position (rest-frame ∼8600 Å). It shares the same positional centroid as the H-band image in the left panel of Figure 2; however, the slit was slightly offset from the intended target position. The seeing on the night of the MOSFIRE observation is much larger than this offset, and is large enough to cause emission originating from positions outside the slit to be probed by the slit. This offset, therefore, does not strongly affect our analysis. Additionally, the rest-frame UV emission, which is expected to originate in the same regions as the Hα emission is well-probed by the slit, though it is more highly extincted than the Hα emission. The offset between the HST and dust emission centroids does demonstrate the importance of dust obscuration in this DSFG, however. Figure 3 shows another look at the HST I-band imaging in blue with ALMA 870 μm dust continuum in red.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We extracted 1D spectra from the MOSFIRE data using the iraf⁶ routine apall. Apertures of extraction were placed on each pixel across the galaxy with a radius equal to half the average seeing. Spectra were extracted with variance weighting. Radial velocities and velocity dispersions were measured from the extracted spectra using the splot fitting routine in iraf. We fit Hα, [N ii]λ 6548, and [N ii]λ 6583 simultaneously with fixed spacing between spectral features and tied but variable emission line full widths at half maximum (FWHM). Pixels contaminated with telluric line emission were omitted from the fits and errors in the best-fit parameters were estimated using 1000 Monte Carlo perturbations of the data by typical sky noise. Figure 4 shows the position–velocity and position–dispersion diagrams measured from these fits. The velocity curve (upper panel), corrected for galaxy inclination angle (see Section 3.1 for details), rises steeply in the inner few kiloparsecs and then flattens to a near-constant velocity in the outer regions of the galaxy. The dispersion curve (lower panel) is centrally peaked due to the effect of beam smearing.

Zoom In Zoom Out Reset image size

We model the rotational velocity, V(r), and velocity dispersion, σ(r) of DSFG850.95 following the method employed by, e.g., Courteau (1997) and Weiner et al. (2006), modified to allow for the parameterization of a possible decline in rotational velocity at high galactocentric radii. The model assumes the intrinsic rotation curve follows a modified arctan function of the form:

$$\begin{eqnarray}&&V(r)=\displaystyle \frac{2}{\pi }{V}_{{\rm{a}}}\arctan \left(\displaystyle \frac{r}{{r}_{{\rm{t}}}}\right)+{cr},\end{eqnarray} \tag{ 1 }$$

where V_a is the asymptotic velocity, r_t is the knee radius, r is the distance from the center of the galaxy, and c is the outer-galaxy slope which allows us to quantify the degree to which the rotation curve increases, flattens, or decreases.

We also measure an isotropic ionized gas velocity dispersion for the entire galaxy, σ₀, which is taken to be the vertical offset of a Gaussian fitted to the position–dispersion data in the bottom panel of Figure 4:

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)=A\exp \left(-\displaystyle \frac{{(r-m)}^{2}}{2{w}^{2}}\right)+{\sigma }_{0},\end{eqnarray} \tag{ 2 }$$

where Σ(r) is the value of the fitted Gaussian at each r, A is the amplitude, r is the distance from the center of the galaxy, m is the central position, w is the width of the Gaussian, and σ₀ is the the isotropic velocity dispersion. We choose σ₀ as the representative dispersion for DSFG850.95 because it is the value the velocity dispersion approaches in the outer galaxy where beam smearing has less of an effect. Beam smearing is well known to artificially increase the observed velocity dispersion toward the centers of galaxies due to the rise in the rotation curve (e.g., Bosma 1978; Weiner et al. 2006; Di Teodoro & Fraternali 2015), but has less of an effect at higher galactocentric radii.

We correct the observed velocity dispersions for the intrinsic instrumental dispersion prior to fitting using the equation: $\sigma =\sqrt{{\sigma }_{\mathrm{obs}}^{2}-{\sigma }_{\mathrm{inst}}^{2}}$ where σ is the velocity dispersion reported in Figure 4, σ_obs is the FWHM of the emission line divided by 2.355 and σ_inst is the combined instrumental and spectral seeing dispersion, taken to be the average width of several telluric lines measured across the unreduced spectrum (σ_inst = 39 ± 2 km s⁻¹). The reported error on σ_inst is the standard deviation in the mean of the measured telluric line widths and are folded into the errors in velocity dispersion in Figure 4. Measured velocity dispersions that are less than the instrumental dispersion are unreliable because small errors on the observed dispersion correspond to large changes in the inferred physical dispersion (Weiner et al. 2006). Dispersions measured below the instrumental seeing are plotted as upper limits at 2σ significance in the bottom panel of Figure 4 and are excluded from the dispersion curve fits, though their central wavelengths are still fitted for the velocity curve.

We fit the position–velocity and position–dispersion data with the pymc3 python package for Markov chain Monte Carlo (MCMC; Salvatier et al. 2016) so that we may characterize the posterior probability distributions of our fit parameters. The light red curves in the top panel of Figure 4 show 50 randomly chosen MC fits and the best-fit curve is plotted in dark red. The best-fit value of asymptotic velocity is $| {V}_{a}| =310\,\pm \,6$ km s⁻¹, the knee radius is r_t = 017 ± 002 or 1.5 ± 0.2 kpc, and the outer slope is c = −0.5 ± 1.0 km s⁻¹ arcsec⁻¹, which is consistent with a flat slope. The best-fit velocity curve is still seen to be rising at the radius of the outermost data point. There is no evidence of a velocity turnover with declining rotation. The average velocity of the flat portion, taken to be the average of the data points between radii of ∼6–14 kpc, corresponding to 1.2–2.8 disk scale lengths, is V_flat = 285 ± 12 km s⁻¹.

The bottom panel of Figure 4 shows the measured velocity dispersion data with 50 randomly selected MC fits overplotted in gray with the best-fit curve in black. The best-fit dispersion, σ₀, is 48 ± 4 km s⁻¹, typical of mildly turbulent disk galaxies at z ∼ 1.5 (e.g., Wright et al. 2007; Law et al. 2009; Contini et al. 2012; Epinat et al. 2012; Wisnioski et al. 2015; Simons et al. 2017).

3.1. H-Band Light Profile Fitting

Using the 2D image fitting code, GALFIT (Peng et al. 2002), we fit a Sérsic profile (Sérsic 1963) to H-band imaging of DSFG850.95 from the UltraVISTA survey (McCracken et al. 2012, see Figure 2, above).

The fit reveals a Sérsic index of n = 1.29 ± 0.03, implying an exponential disk light profile. Adding a second Sérsic component to fit a bulge causes significantly worse fits. The H-band half light radius is r_1/2 = 8.4 ± 0.1 kpc, or using the relation r_d = r_1/2/1.68, the disk scale radius is r_d = 5.0 ± 0.1 kpc and r_2.2 = 11.0 ± 0.2 kpc. We divide the velocities in the top panel of Figure 4 by sin(i) to correct for galaxy inclination, assuming an intrinsic edge-on axis ratio of 0.25 (Wuyts et al. 2016). From our best-fit axis ratio, b/a = 0.668 ± 0.008, we determine the galaxy has an approximately edge-on inclination of i = 87° ± 2° using the equation (e.g., Wuyts et al. 2016):

$$\begin{eqnarray}&&\cos i=\sqrt{\displaystyle \frac{{(b/a)}^{2}-{0.25}^{2}}{1-{0.25}^{2}}}.\end{eqnarray} \tag{ 3 }$$

Table 1 summarizes the physical characteristics of DSFG850.95 from GALFIT and our other analyses.

Table 1.  Physical Characteristics of DSFG850.95

Property	Value
z	1.555
Vflata	285 ± 12 km s−1
σ0	48 ± 4 km s−1
Vflat/σ0	6
ΔP.A.b	23° ± 4°
nc	1.29 ± 0.03
id	87° ± 2°
r1/2	8.4 ± 0.1 kpc
r2.2	11.0 ± 0.2 kpc
rmaxe	14.1 kpc
SFR	${280}_{-90}^{+110}$ M⊙ yr−1
M⋆	(3.8 ± 3.0) × 1010 M⊙
${M}_{{{\rm{H}}}_{2}}$	(8.88 ± 0.03) × 1010 M⊙
Mdyn(14 kpc)	(2.7 ± 0.3) × 1011 M⊙
Mdyn(R1/2)	(1.9 ± 0.1) × 1011 M⊙
fDM(R1/2)	0.44 ± 0.08

Notes.

^aAverage of flat rotation on both sides of DSFG850.95. ^bPosition angle offset between slit and H-band semimajor axis. ^cSérsic index. ^dInclination angle. ^eMaximum detected radius of Hα.

Download table as:  ASCIITypeset image

There is a misalignment of 23° ± 4° between the position angles of the H-band semimajor axis and the slit. Large misalignments would result in an underprediction of the maximum rotational velocity of a galaxy; however, Weiner et al. (2006) shows that the average decrease in observed rotational velocities for a slit misalignment of ∼23° is small—approximately the same as the error on our measurement of V_flat.

The H-band imaging, together with the rotational velocity and velocity dispersion analyses in Section 3 suggest DSFG850.95 is a rotationally dominated disk galaxy at z = 1.555 that is similar to massive disk galaxies at low redshift.

Our radial velocity measurements as a function of radius, corrected for galaxy inclination and ionized gas velocity dispersion support allow us to measure the total mass enclosed as a function of radius in the galaxy using the equation (e.g., Erb et al. 2003; Swinbank et al. 2004; van Starkenburg et al. 2008; Epinat et al. 2009; Di Teodoro et al. 2018):

$$\begin{eqnarray}&&{M}_{\mathrm{dyn}}(r)=\displaystyle \frac{{V}_{c}^{2}(r)r}{G},\end{eqnarray} \tag{ 4 }$$

where G is the gravitational constant, r is the galactocentric radius, and V_c is the circular velocity, given by:

$$\begin{eqnarray}&&{V}_{c}=\sqrt{{V}^{2}+2{\sigma }_{0}^{2}\displaystyle \frac{r}{{r}_{d}}},\end{eqnarray} \tag{ 5 }$$

where V is the inclination-corrected rotational velocity, σ₀ (measured to be 48 ± 4 km s⁻¹) is the isotropic ionized gas velocity dispersion, and r_d is the exponential disk scale length measured from H-band imaging. The median difference between the rotational and circular velocities is 2%, since σ₀ is small compared to V.

The top panel of Figure 5 shows the dynamical mass as a function of radius. The total dynamical mass enclosed within the H-band half light radius is (1.9 ± 0.1) × 10¹¹ M_⊙, and within 14 kpc is (2.7 ± 0.3) × 10¹¹ M_⊙. The horizontal dashed line and shaded region in the top panel of Figure 5 marks the maximum dynamical mass and its associated error. The bottom panel of Figure 5 shows the fraction of dynamical mass enclosed by radius.

Zoom In Zoom Out Reset image size

Next we compare the dynamical mass at the half light radius with stellar and gas masses. Casey et al. (2017) uses MAGPHYS with the high-z extension (da Cunha et al. 2008, 2015) to fit UV through submillimeter ancillary photometry from the COSMOS collaboration (Capak et al. 2007; Laigle et al. 2016) and finds a stellar mass of ${M}_{\star }=({3.8}_{-1.0}^{+0.4})\times {10}^{10}$ M_⊙. We follow the procedure of Hainline et al. (2011) for similarly selected DSFGs at similar redshifts to derive a more conservative error than Casey et al. We take the error to be half the difference between the stellar mass estimated with an instantaneous burst history and that estimated with a continuous star formation history. Thus, we arrive at a total stellar mass of M_⋆ = (3.8 ± 3.0) × 10¹⁰ M_⊙. Table 1 summarizes all the estimated and measured parameters of this galaxy.

Following the prescription of Scoville et al. (2016) we use our ALMA 870 μm dust continuum observations to estimate a molecular gas mass of (8.9 ± 0.4) × 10¹⁰ M_⊙, all of which is contained within the half light radius. We calculate the dark matter fraction at the half light radius using the equation:

$$\begin{eqnarray}&&{f}_{\mathrm{DM}}({R}_{1/2})=1-\displaystyle \frac{{M}_{\star }({R}_{1/2})+{M}_{\mathrm{gas}}}{{M}_{\mathrm{dyn}}({R}_{1/2})},\end{eqnarray} \tag{ 6 }$$

where M_⋆(R_1/2) is half of the total stellar mass, M_gas is the gas mass, and M_dyn(R_1/2) is the dynamical mass at the half light radius. We find f_DM(R_1/2) = 0.44 ± 0.08. Figure 6 shows the dark matter fraction as a function of circular velocities for a few studies in the literature (Martinsson et al. 2013a, 2013b; Bland-Hawthorn & Gerhard 2016; Genzel et al. 2017). DSFG850.95 has a dark matter fraction closer to that of the Milky Way than the galaxies in Genzel et al. (2017), though it has a higher rotational velocity than the Milky Way.

Zoom In Zoom Out Reset image size

We conclude by noting that our finding of a flat rotation curve may not be inconsistent with the conclusions of Lang et al. (2017). The majority of disk galaxies in cosmological hydrodynamic simulations by Teklu et al. (2018) exhibit flat rotation curves with only 38% hosting declining rotation curves. Lang et al. (2017) split their data out into two subsamples according to the degree of pressure support, one with V_rot/σ₀ > 6.3 and the other < 6.3. They define V_rot as the intrinsic, inclination-corrected peak rotation velocity and σ₀ as the intrinsic velocity dispersion, both equivalent to our V_c and σ₀. The galaxies with less pressure support (V_rot/σ₀ > 6.3) have flatter outer rotation curves than those with more pressure support, supporting the conclusion that much of the observed velocity decrease in their sample is due to increased pressure support. Our rotation curve, with a V_c/σ₀ ∼ 6 is consistent with their stacked curve of galaxies with V_rot/σ₀ > 6.3. A statistical sample of galaxy rotation curves detected well beyond 10 kpc is needed to determine the fraction of flat or declining rotation curves at high redshift.

Recent observations of disk galaxy rotation curves have found declining rotational velocities beyond ≈2.2 exponential disk scale lengths, implying intermediate redshift disk galaxies have lower dark matter fractions than z = 0 galaxies as well as increased velocity dispersion support.

In this paper we present a clear counterexample from MOSFIRE observations of Hα and [N ii] emission in the z = 1.555 disk galaxy, DSFG850.95, showing a flat rotation curve between ∼6 and 14 kpc (1.2–2.8 disk scale lengths), well-fit by an arctangent function. At the H-band half light radius we estimate the dark matter fraction to be 0.44 ± 0.08 based on UV-through-submillimeter photometry and 870 μm dust continuum observations. Ground-based H-band imaging reveals a Sersic index of n = 1.29 ± 0.03 and an edge-on inclination angle. Taken together, these data suggest DSFG850.95 is a massive, rotationally supported disk galaxy in place by z = 1.555 (4.1 Gyr after the big bang) with a flat rotation curve and a dark matter fraction at the half light radius, similar to low redshift disk galaxies.

The authors thank the anonymous referee for a constructive report that significantly strengthened the arguments presented in this paper. The authors also wish to thank Michael Boylan-Kolchin, Michael Petersen, and Justin Spilker for useful discussions. P.M.D. and C.M.C. acknowledge financial support from NASA Keck PI Data Awards: 2012B-U039M, 2011B-H251M, 2012B-N114M, 2017A-N136M, the University of Texas at Austin College of Natural Sciences, and NSF grants AST-1714528 and AST-1814034. This paper makes use of ALMA data ADS/JAO.ALMA#2015.1.00568S. ALMA is a partnership of ESO (representing its member states), NSF (USA), and NINS (Japan), together with NRC (Canada) and NSC and ASIAA (Taiwan) and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO, and NAOJ. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.