The Most Ancient Spiral Galaxy: A 2.6-Gyr-old Disk with a Tranquil Velocity Field

The spiral galaxy candidate (hereafter A1689B11) was first recorded as "source 11" with a photometric redshift of $z=2.9\pm 0.2$ in the strong-lensing analysis of the galaxy cluster Abell 1689 (Broadhurst et al. 2005). Source 11 is gravitationally lensed into two highly magnified images B11.1 and B11.2 and a central image B11.3 (Figure 1). Broadhurst et al. (2005) were the first to point out the spiral feature of A1689B11 and speculated that it would be the highest-redshift spiral galaxy if confirmed. The first spectroscopic redshift of z = 2.5 was provided in Limousin et al. (2007) based on rest-frame UV spectroscopy with Keck/LRIS. The detection of Hα lines at z = 2.54 was reported in the near-infrared (NIR) multislit survey of Yuan et al. (2013b). This work presents the NIR IFS observation of image B11.1 (Figure 1). The imaging and photometric data used in this work are obtained from the Hubble Space Telescope (HST; proposal IDs: 11802, 9289, 11710) and the Spitzer Space Telescope archives (program ID: 20439).

Zoom In Zoom Out Reset image size

AO-aided NIR IFS observations with NIFS (Near-infrared Integral-Field Spectrograph; McGregor et al. 2003) were conducted between 2013 March and 2014 March on the Gemini North telescope under excellent weather conditions (average seeing ∼05, airmass 1–1.6). Our total allocated observational time was 13.5 hr in band 1 (program ID: GN-2013A-Q-23-64). The observation implemented a dithering pattern of "ABAABA" or "AB," i.e., ∼33%–50% of the target exposure time was spent on sampling sky frames in order to facilitate a good sky background subtraction. The observation centered on image B11.1 instead of B11.2 because of the laser-guide star requirement of NIFS. The coordinates of the pointing center are given in Table 1. The field of view of NIFS is 30 × 30, with 29 slitlets, each 01 wide. NIFS delivers a spectral resolving power of R ∼ 5300 in the K band, corresponding to a rest-frame Gaussian velocity resolution of σ ∼ 24 km s⁻¹.

Table 1.  Physical Properties of A1689B11

R.A. (J2000)	13:11:33.336
Decl. (J2000)	–01:21:06.9
Spectroscopic data measurements:
Redshift (zHα)	2.540
Velocity dispersion (${V}_{\sigma }$, mean)	23 ± 4 km s−1
Velocity dispersion (${V}_{\sigma }$, outer disk 1–2 kpc)	15 ± 2 km s−1
Star formation rate (SFRHα)	22 ± 2 M⊙ yr−1
Dust attenuation (E(B–V)Balmer decrement)	0.73
GALFIT/exponential disk best fit:
Inclination (i)	55° ± 10°
Position angle (PA)	−36° ± 6°
Scale length (rs)	1.3 ± 0.4 kpc
Half-light radius (${R}_{1/2}$)	2.6 ± 0.7 kpc
Effective radius (Re)	2.0 ± 0.4 kpc
2D velocity disk model best fit:
Inclination (i, velocity)	51° ± 2°
Position angle (PA, velocity)	−37° ± 2°
Radius rotation (Rt, velocity)	1.7 ± 0.1 kpc
Rotation Velocity (Vrot)	200 ± 12 km s−1
Inferred physical parameters:
Star formation rate (SFRSED)	22 ± 3 M⊙ yr−1
Stellar mass (${M}_{\mathrm{star}}$)	${10}^{9.8\pm 0.3}$ M⊙
Dynamical mass (${M}_{\mathrm{dyn}}$)	${10}^{10.2}$ M⊙
Average SFR surface density (ΣSFR)	0.3 M⊙ yr−1  kpc−2
Average gas surface density (Σgas)	158 M⊙  pc−2
Gas fraction (${f}_{\mathrm{gas}}$)	∼0.18
Toomre parameter ${{\rm{Q}}}_{\mathrm{gas}}$	∼0.7

Note. All values have been corrected for lensing magnifications. PA are defined as follows: PA = 0 when the major axis is up, and positive if rotated counterclockwise. Assumptions for inferred values are described in Section 3.

Download table as:  ASCIITypeset image

The data were reduced using a varied Gemini IRAF package following standard procedures (e.g., Storchi-Bergmann et al. 2009). This IRAF package produces sky-subtracted, telluric-corrected, and flux-calibrated datacubes. Each individual 900 s exposure was spatially aligned based on the brightest Hα spaxels and co-added based on a mean sigma-clipping procedure. The 1σ error datacubes were generated during the sigma-clipping process. We used the star HIP67004 of spectral type A0V for telluric correction and flux calibration. The systematic uncertainties in the flux calibration are estimated to be within 20%. A total of $27\,\mathrm{ks}$ on-source exposures were obtained in the K band with an angular resolution of ∼01, corresponding to a median physical scale of ∼400 pc on the source plane of A1689B11. Our observation reached a 3σ Hα emission line surface brightness depth of 3 × 10⁻¹⁷ erg s⁻¹ cm⁻² arcsec⁻².

The astrometry of the datacube is calibrated by assigning the HST coordinate of the galaxy center to the brightest spaxel in the wavelength-collapsed 2D NIFS image and rotating with the positional angle from the observation. The astrometric uncertainty of our IFS datacube is about one spaxel, i.e., ∼01.

Our K-band observation was originally planned to detect Hα and [N ii] lines. To analyze the IFS spectra, we first collapse the datacube in the wavelength dimension in the vicinity of the Hα line to create an Hα 2D map. We use the raw Hα 2D map to generate an initial 2D mask that flags spaxels with no obvious Hα line detections. We then manually inspect the spectrum of each individual spaxel and refine the mask with three types of visual flagging: (1) significant emission lines, (2) possible emission lines, and (3) no obvious emission lines. The visual inspection of the datacube is necessary in order to reject spaxels that are spurious for low signal-to-noise ratio (S/N) data. We use this visual mask to supervise our subsequent automatic emission-line fitting procedure, i.e., we demand the fitting result to be consistent with our visual mask. For example, an emission-line fitting result of S/N ≥ 5 should have a visual flag of (1), whereas a line detection of S/N $\lt 3$ should be consistent with a flag type of (3). Those with $3\leqslant$ S/N $\lt 5$ should match flag type (2).

Our automatic emission-line fitting procedure involves fitting Gaussian profiles simultaneously to three emission lines in each spaxel: [N ii] λλ6548, 6583 and Hα (Figure 2). The line profile fitting is conducted using a ${\chi }^{2}$ minimization procedure weighted by the inverse of the variance spectrum. The fitting result is very sensitive to the weighting used. The 1σ error datacube generated from the mean sigma-clipping procedure described in Section 2.2 does not provide good weighting to the emission-line fitting, as it results in the S/N of the emission line being overestimated and inconsistent with our visual mask flags. Instead, we use the variance over a sky region devoid of emission lines as the weighting. We select a range of sky spaxels to generate a variance spectrum (green lines in Figure 2) for each spaxel using a bootstrap procedure. The S/N of the emission-line fitting result is consistent with the visual mask. This weighting method is similar to that used in Leethochawalit et al. (2016) and is reasonable because the variance of the NIR data is usually dominated by sky residuals.

Zoom In Zoom Out Reset image size

Hα lines are detected in ∼130 individual spaxels at ≥5σ level (Figure 2, panel (1)). [N ii] lines are not detected above 3σ in either individual spaxels or the integrated spectrum (Figure 2, panel (2)). Yuan et al. (2013b) report an [N ii] line detection on one of the four slits configured with different positional angles on B11.1 and B11.2, implying a spatial variation of metallicities. We extract a mock slit spectrum from our NIFS data based on the slit setup of Yuan et al. (2013b), but no [N ii] line is detected above 3σ. Deeper IFS data from our ongoing OSIRIS/Keck observation on image B11.2 will help to discern [N ii] lines and the spatial metallicity distribution of A1689B11. Figure 2 shows the Hα 2D S/N map and examples of single Gaussian fits to Hα emission lines of individual spaxels. An integrated spectrum from co-added spaxels is also presented in Figure 2 for comparison. Figure 2 highlights the necessity of high spatial resolution in order to distinguish beam smearing from intrinsic line width. It is also apparent from Figure 2 that the Hα lines detected on individual spaxels have very narrow widths (approaching the NIFS instrumental line width) and show systematic offsets in the centroids. The focus of this paper is to report the kinematics of Hα lines (see Section 3).

Abell 1689 is one of the most studied lensing clusters with well-constrained mass models. We use the best-fit model from Limousin et al. (2007) and the software Lenstool⁸ (Kneib et al. 1993; Jullo et al. 2007) to reconstruct the source-plane properties of A1689B11.

We show the lensing configuration and critical lines of A1689B11 in the left panel of Figure 1. A1689B11 is lensed into a three-image system: B11.1, B11.2, and B11.3. Image B11.2 lies in between the tangential and the radial critical line, resulting in a large flux magnification factor ($\mu =12\pm 2$) and a relatively undistorted image (i.e., equally magnified in the major and the minor axis). Our NIFS observation is centered on the less magnified ($\mu =7.2\pm 0.8$) image B11.1. B11.1 is magnified more in the minor-axis direction, making it appear more face-on than the intrinsic image. The geometric magnification for image B11.1 is ∼4 (along the minor axis) × 1.8 (along the major axis). A demagnified central image B11.3 is predicted by the lens model and also observed to be overlapping with one of the brightest cluster members of Abell 1689 (Broadhurst et al. 2005; Limousin et al. 2007).

The statistical error of the lens modeling is estimated from a set of Markov Chain Monte Carlo (MCMC) realizations implemented in Lenstool and is typically ∼10% for Abell 1689. Both multiple images (B11.1 and B11.2) yield consistent source-plane morphologies and positions within 02. Because B11.1 and B11.2 are not close to the critical lines, the systematic errors in the source-plane morphologies are small compared to the case of giant arcs. Our source-plane reconstructed morphologies are robust within the model statistical uncertainties. The NIFS datacube is reprojected to the source plane after being remapped to the HST coordinates using the calibrated astrometry. Emission-line fitting analyses have been carried out on both the image and source plane, and the results are self-consistent. We present both the image-plane and source-plane properties in the following sections.

3.1.1. NIFS Hα Map and Kinematics

In Figure 3 we show the 2D Hα flux intensity map and 2D Hα kinematic maps on the observed frame (image plane), i.e., before lensing reconstruction. To facilitate the comparison with HST images, we mark the extent of HST broadband detections in contours in all panels. We use the HST/F814W band image as an example because it has the deepest exposure among all HST filter observations. We include the morphology of other HST bands in Section 3.2.

Zoom In Zoom Out Reset image size

In order to preserve the spatial resolution offered by lensing and because we have sufficient detections on individual spaxels, we do not attempt to bin the original IFS datacube. The NIFS maps in Figure 3 are unsmoothed and include only data with S/N ≥ 5. Hα emission lines are detected throughout the HST contours; Hα emission lines are relatively stronger near the galaxy core and star-forming (SF) clumps in the spiral arm.

The rest-frame line-of-sight velocity 2D map is derived from the Hα line centroids with respect to the systematic redshift and the velocity dispersion map from the Hα Gaussian line width. We use the median Hα line center to calculate the systematic redshift. The instrumental profile has been subtracted in quadrature from the best-fit Gaussian width to derive the intrinsic line width. Figure 3 shows a systematic rotation and uniformly small velocity dispersion across the disk on the observed plane.

The intrinsic morphology and NIFS measurements are shown in source-plane maps in Figure 4. We use a grid subsampling of ${\text{}}{echant}=20$ on the image plane and ${\text{}}s\_{echant}=10$ on the source plane in Lenstool to optimize the spatial resolution of the source-plane reconstruction. The source-plane HST ACS image in Figure 4 is unbinned and has a pixel scale of 0005 after subsampling. The source-plane NIFS datacube is rebinned adaptively by 5–11 pixels to allow for Hα line fitting of S/N > 5 in each bin. Note that because of the difference in the point-spread function (PSF) on the image and source planes, the improvement in the source-plane S/N is usually less than what the magnification map predicts. On average we achieve a source-plane spatial resolution of ∼150 pc on the HST ACS images and ∼400 pc on the NIFS 2D maps.

Zoom In Zoom Out Reset image size

The source-plane 2D velocity map in Figure 4 clearly shows a velocity gradient consistent with a systematic rotation. The velocity dispersion is enhanced slightly in the kinematic center owing to beam smearing from the rotation. We estimate a maximum beam-smearing effect at the kinematic center to be ∼24 km s⁻¹ based on the velocity map of panel (3) in Figure 4, consistent with the 2D velocity dispersion map of panel (4). We do not see any significant spatial correlation of the Hα intensity map with the velocity dispersion map. Figure 4 confirms that A1689B11 has a uniformly low velocity dispersion across the disk on the source plane. The mean velocity dispersion averaged over all spaxels is ${V}_{\sigma ,\mathrm{mean}}\,=23\pm 4$ km s⁻¹. Excluding the central spaxels that are affected most by beam smearing, the average velocity dispersion on the outer disks is ${V}_{\sigma ,\mathrm{outer} \mbox{-} \mathrm{disk}}=15\pm 2$ km s⁻¹.

3.1.2. Disk Model Fitting

We use an empirically motivated arctangent function from Courteau (1997) to model the 2D velocity field:

$$\begin{eqnarray}&&V(R)={V}_{0}+\displaystyle \frac{2}{\pi }{V}_{c}\arctan \displaystyle \frac{R-{R}_{0}}{{R}_{t}}.\end{eqnarray} \tag{ 1 }$$

The source-plane line-of-sight velocity v_s(R) is related to the intrinsic velocity V(R) by the inclination angle i:

$$\begin{eqnarray}&&{v}_{s}(R)=V(R)\sin (i).\end{eqnarray} \tag{ 2 }$$

The source-plane radius vector R_s is related to the intrinsic radius vector R by

$$\begin{eqnarray}{{\boldsymbol{R}}}_{s}={\boldsymbol{R}}\,\cos (i)\left[\begin{array}{cc}\cos (\mathrm{PA}), & -\sin (\mathrm{PA})\\ \sin (\mathrm{PA}), & \cos (\mathrm{PA})\end{array}\right],\end{eqnarray} \tag{ 3 }$$

where PA is the positional angle. For the definition of PA and i, we use the same convention as GALFIT (Peng et al. 2010), i.e., PA = 0 is to the north (up), and PA = 90 is to the east (left); i = 0 is face-on, and i = 90 is edge-on.

The seven free parameters of the 2D disk model are thus central velocity V₀, inclination i, position angle PA, disk dynamic center ${{\boldsymbol{R}}}_{0}$ (R_0x, R_0y), turnover radius R_t, and asymptotic velocity V_c. Note that V₀ is close to zero if the median Hα line center is a good approximation for the central velocity. The best-fit model is obtained using a ${\chi }^{2}$ minimization procedure similar to that described in Jones et al. (2010b). We use the statistical 1σ errors from the emission-line fitting to compute the ${\chi }^{2}$. The uncertainties of the best-fit parameters are estimated by perturbing the model until the ${\chi }^{2}$ increases by one standard deviation from the best-fit model.

Figure 5 shows the best-fit 2D disk model and the residual maps. The reduced ${\chi }^{2}$ of our best-fit model is 2.1, and the best-fit parameters are as follows: asymptotic velocity V_c = 200 ± 12 km s⁻¹, turnover radius R_t = 1.7 ± 0.1 kpc, inclination angle $i=51^\circ \pm 2^\circ$, and positional angle $\mathrm{PA}=-37^\circ \pm 2^\circ$. The kinematic center is marked in Figure 5 and is consistent with the brightest intensity from both the Hα image and the HST broadband images within 01.

Zoom In Zoom Out Reset image size

We can also use Equation (1) to fit the 1D rotation curve. We extract the 1D rotation velocities using a 08 slit along the major axis with the PA and the kinematic center from the best-fit 2D model (see Figure 5, left panel). The seven-parameter 2D model then reduces to a five-parameter 1D model. In Figure 6, each data point represents one resolution element (01) binned along the slit, and the error bars are the standard deviation of each data bin. We show the best-fit 1D rotation curve and its 1σ variations in Figure 6. We find that the best-fit parameters from the 1D model are consistent with the 2D disk model within statistical errors. Note that we have also tested fitting a 3D disk model such as Di Teodoro & Fraternali (2015) to the datacube. Unfortunately, we do not have sufficient S/N per spaxel for a reliable fit of the 3D model. Our 2D and 1D disk model provides a reasonably good fit to the data. Overall, a regular rotating disk model is sufficient to explain the velocity field of A1689B11. The small residuals on the edge of the disk (Figures 5 and 6) are most likely due to underestimation of the observational noise.

Zoom In Zoom Out Reset image size

The simple 2D disk model does not consider beam smearing and variations of the velocity dispersion. However, as discussed in Section 3.1.1, the disk shows a uniformly low velocity dispersion and the beam smearing does not contribute significantly to the shape of the velocity field. We show the 1D velocity dispersions (red points) along the kinematic major axis in Figure 6. The velocity dispersion in the galactic center is mildly enhanced owing to beam-smearing effects. The dispersion on the outer disk is approaching the thermal broadening threshold (12 km s⁻¹) of the Hα emission line (black dotted line).

We discuss the kinematics of A1689B11 and compare it with other samples in Section 4.

To compare our NIFS data with HST images in more detail, we analyze source-plane morphologies for all HST ACS and WFC3 bands. B11.1 and B11.2 are well detected in all four IRAC bands from Spitzer. We do not gain resolved morphological knowledge from IRAC because of its poor spatial resolution (∼25, larger than the lensed images). Therefore, IRAC data are only used when deriving global properties of A1689B11 in Section 3.3.

We take full advantage of the multiple images from gravitational lensing and combine information from both images B11.1 and B11.2 to derive the best source-plane morphology. B11.1 is covered by six filter bands (F475W, F625W, F775W, F814W, F105W, F140W), and B11.2 is covered by these bands plus the additional two bands of F125W and F160W. Source-plane morphologies for individual HST bands are reconstructed using the same method as implemented for the F814W image described in Section 3.1. The spiral features are more prominent in the HST/ACS optical images (005 resolution), whereas HST/WFC3 images (01 resolution) show mostly the central disk. We include source-plane morphologies for individual bands in the Appendix.

We perform GALFIT (Peng et al. 2010) on the source-plane images using a single exponential disk model. The fitting yields consistent morphological parameters for all bands of HST images and for both B11.1 and B11.2. The mean and standard deviation of the scale lengths (r_s) calculated from all 14 source-plane images are ${r}_{s}=1.3\pm 0.4\,\mathrm{kpc}$, with an inclination angle of $i=55^\circ \pm 10^\circ$ and $\mathrm{PA}=-36^\circ \pm 6^\circ$. The single exponential disk model from GALFIT is in excellent agreement with geometric parameters derived from the NIFS kinematic 2D disk model (summarized in Table 1). For the convenience of comparing with various definitions of radius in literature, we convert r_s to the half-light radius ${R}_{1/2}$ and the effective radius R_e using the empirical relation of ${R}_{1/2}=2.2{r}_{s},{R}_{e}=1.68{r}_{s}$ (Glazebrook 2013). The half-light radius of A1689B11 is therefore 2.6 ± 0.7 kpc. The central area of A1689B11 shows significant Hα emission from our NIFS data and some elongated substructures in the HST ACS images (Figure 4; the Appendix). The central area could be a superposition of a star cluster and a bulge/bar component. We experimented with adding a bulge component to the exponential disk; however, the fitting does not converge for most of the images. We are currently investigating a more sophisticated procedure for the bulge–disk decomposition, and our preliminary result shows a very small bulge component. We will report the full analysis in a future work focusing on the correlation of pitch angle to bulge of this spiral galaxy (T. Yuan 2017, in preparation).

Because gravitationally lensed images cover an extended area, the probability of having contaminated foreground/background sources in the field of view of lensed images is larger than a nonlensed high-z galaxy case. We cross-compare substructures on the source-plane images of B11.1 and B11.2 to reject clumps/knots that may not be associated with the spiral galaxy. We use the F814W HST images that have the deepest exposure and best spatial resolution for the cross-examination. Figure 7 demonstrates our clump identification and rejection procedure. We first manually identify bright clumps in the right panels of Figure 7. There are seven clumps (c0–c4, x1, x2) marked for lensed image B11.1 and nine clumps (c0–c4, x3, x4, x5, x6) for lensed image B11.2. We then predict the source-plane positions and brightnesses of the clumps (c0'–c4', x1'–x6') and cross-examine them in the left panels of Figure 7. Because lensing conserves surface brightness and images B11.1 and B11.2 cover similar intrinsic areas of the lensed galaxy, the resolved source-plane substructures of B11.1 and B11.2 should be identical from the same HST observation.⁹ Source-plane clump brightness and positions that are not consistent within the uncertainties of lens models are considered interlopers. For example, clump x2' is identified initially on B11.1 but is not detected on B11.2 within the lensing position uncertainty (rms = 02); therefore, we reject clump x2' and consider it as a foreground or background source. Similarly, clumps x4' and x5' identified initially on B11.2 are not detected on B11.1 and are rejected. For clump x1' (yellow circles in Figure 7) that is initially identified on B11.1, there is a marginal 3σ detection in the source plane of B11.2; however, the lens model predicts a source-plane brightness that is ∼3 times brighter than what is observed. We therefore reject x1' because of the inconsistency in flux magnification uncertainties (∼10%–20%). Note that clumps/knots that are identified as foreground/background objects in this lensing analysis are at significantly different cosmological distances from the source galaxy and are therefore not satellites. The final combined morphologies of B11.1 and B11.2 are shown in the Appendix.

Zoom In Zoom Out Reset image size

Finally, we use colors of the clumps as an alternative check for substructure identities. We measure aperture photometry for all clumps in B11.1 in Figure 7 using HST broadband images. A reliable (S/N > 5) photometry can only be measured for six clumps across a minimum of three bands: c0'–c4' and x2'. We then carry out spectral energy distribution (SED) fitting for the clumps using the software LE PHARE (Ilbert et al. 2010). We fix the redshift at z = 2.54 and use the stellar population synthesis models of Bruzual & Charlot (2003). We choose the initial mass function (IMF) of Chabrier (2003) and the Calzetti et al. (2000) attenuation law, with $E(B-V)$ ranging from 0 to 2 and an exponentially decreasing star formation history. The best-fit SEDs are shown in Figure 8. Clumps c0'–c4' show similar SEDs, whereas x2' has a very different SED, consistent with our lensing source-plane position analysis that x2' is most likely an interloper. We also show in Figure 8 the best-fit SED for the total photometry of B11.1 (Section 3.3). In total, we confirm five bona fide clumps (c0'–c4') for A1689B11. There are significant Hα emissions detected on the central clump c0' and on clump c2' on the spiral arm; Hα lines are also detected on clump c1' and c3' on the spiral arm (see Figure 4). A detailed analysis on individual clump properties will be reported in a separate work when we combine our ongoing OSIRIS/Keck data analysis on image B11.2. Our preliminary result shows that the sizes (400−600 pc) and surface SF densities ($0.2\mbox{--}0.3$ M_⊙ yr⁻¹ kpc⁻²) of the clumps are comparable to high-redshift SF galaxies and are in the intermediate range of what has been reported for $z\gt 1$ clumpy galaxies (e.g., Jones et al. 2010b; Genzel et al. 2011; Livermore et al. 2015).

Zoom In Zoom Out Reset image size

We estimate the total stellar mass from the SED fitting of broadband photometries from the HST ACS and the Spitzer IRAC data (details in Yuan et al. 2013b; best-fit SED shown in Figure 8). The best-fit stellar mass for A1689B11 is ${M}_{\mathrm{star}}$ = 10^9.8 ± 0.3 M_⊙, and the best-fit extinction value is E(B – V)_stellar = 0.22. The total dust-corrected SFR from the SED fitting is SFR_SED = 22 ± 3 M_⊙ yr⁻¹. All values have been corrected for the lensing flux magnification. The dust-uncorrected SFR from the total Hα fluxes of our NIFS observations is SFR_{Hα nodust} = 3.9 ± 0.4 M_⊙ yr⁻¹. Using the nebular dust extinction $E{(B-V)}_{\mathrm{nebular}}=0.73$ from Yuan et al. (2013b) and the nebular attenuation curve of Cardelli et al. (1989), the dust-corrected SFR from Hα is SFR_Hα = 22 ± 2 M_⊙ yr⁻¹, in agreement with SFR_SED. The total stellar mass and SFR of A1689B11 are consistent with an SF galaxy on the $z\sim 2$ mass–SFR relation (the main sequence) within the ∼0.3 dex 1σ scatter (e.g., Zahid et al. 2012) of the relation. The total SFR of A1689B11 is 10–20 times higher than a typical spiral galaxy of similar mass at $z\sim 0$.

Because we do not have spatially resolved dust attenuation measurements, the following estimation of SFR surface density is indirect and based on a few assumptions. Assuming that the spatially resolved $E{(B-V)}_{\mathrm{nebular}}$ is similar to the global $E{(B-V)}_{\mathrm{nebular}}$ measured from the slit data of Yuan et al. (2013b), the average SFR surface density is Σ_SFR = 0.3 M_⊙ yr⁻¹ kpc⁻². This value of Σ_SFR is in the intermediate range of $z\gt 1$ SF galaxies and is 1–2 orders of magnitude higher than local SF galaxies (e.g., Swinbank et al. 2012; Fisher et al. 2017; Zhou et al. 2017). Assuming a simple Kennicutt–Schmidt (KS) relation (Kennicutt 1998), we convert Σ_SFR into a gas surface density of Σ_gas ∼ 158 M_⊙ pc⁻². Assuming that the surface area of the gas is $2\pi {R}^{2}$, where R is the radius where Hα is detected (∼1.7 kpc), we derive a gas fraction (${f}_{\mathrm{gas}}={M}_{\mathrm{gas}}$/(${M}_{\mathrm{gas}}$+${M}_{\mathrm{star}}$)) of ∼18%. The ${f}_{\mathrm{gas}}$ is in the lower range of SF galaxies at $z\gt 1$ but still significantly higher than local SF galaxies (e.g., Tacconi et al. 2013). We compute the Toomre Q parameter (Toomre 1964) for a gas-dominated disk as defined by $Q\approx \kappa {V}_{\sigma }/\pi G{{\rm{\Sigma }}}_{\mathrm{gas}}$, where κ is the epicyclic frequency of the galaxy's rotation, ${V}_{\sigma }$ is the gas velocity dispersion, and ${{\rm{\Sigma }}}_{\mathrm{gas}}$ is the surface mass density. To compute κ, we assume a Keplerian disk so that $\kappa ={\rm{\Omega }}$, where Ω is the angular frequency. The Ω calculated at the half-light radius is ${\rm{\Omega }}\sim 2\pi /100(\mathrm{Myr})$. Adopting ${V}_{\sigma }=23$ km s⁻¹, we find Q_gas ∼ 0.7 for A1689B11. The small Q_gas is commonly measured in clumpy high-redshift SF galaxies and is consistent with the scenario that large SF clumps form in the violent disk instability (Genzel et al. 2011; Law et al. 2012; Glazebrook 2013; Shibuya et al. 2016). We caution that the intrinsic Σ_SFR can easily change by a factor of two because of the unknown dust attenuation. The systematic errors related to the methodology of deriving ${f}_{\mathrm{gas}}$ and Q_gas are also highly uncertain, and proper calculation requires estimation of the Mach number and molecular gas observations (e.g., Federrath et al. 2017b).

The dynamical mass assuming a rotationally supported disk is ${M}_{\mathrm{dyn}}={{RV}}_{\mathrm{rot}}^{2}$/$G={10}^{10.2\pm 0.1}$ M_⊙, where ${V}_{\mathrm{rot}}$ is the asymptotic velocity V_c and the radius R is the turnover radius R_t from the 2D disk model (Section 3.1.2). The ratio of the dynamical mass to the stellar mass for the inner ∼2 kpc is therefore ${M}_{\mathrm{dyn}}/{M}_{\mathrm{star}}\sim 2.8$. Taking ${M}_{\mathrm{dyn}}$ to be the sum of dark matter and baryonic matter, the inferred dark matter mass fraction within the inner ∼2 kpc (compatible to ${R}_{1/2}$) is ${f}_{\mathrm{DM}}=60 \%$. The value of f_DM is in the typical range of local late-type spiral galaxies of similar V_c and much larger than the massive baryon-dominated clumpy disk galaxies at z ∼ 2 (Genzel et al. 2017; Lang et al. 2017). However, this conclusion is subject to uncertain observational and methodological errors of f_DM. For example, by simply propagating the errors of ${M}_{\mathrm{dyn}}$ and ${M}_{\mathrm{star}}$, we obtain ${f}_{\mathrm{DM}}=60 \% \pm 40 \%$. In addition, in order to decompose the contributions of the baryonic disk and the dark matter halo to the total rotation curve, assumptions on the mass-to-light ratio (M/L) and scale height of the stellar disk have to be made (e.g., Aniyan et al. 2016); both are difficult to constrain for our Galaxy. We therefore caution again that values derived in these two paragraphs are dominated by systematic errors and should be interpreted with some caution.

The morphology of A1689B11 is indicative of a prototype spiral arm but not conclusive. The clumpy morphology is also suggestive of a merger. However, by combining the morphology with the high spatial resolution kinematics presented in this work, we rule out the merger scenario for A1689B11.

Distinguishing mergers from isolated disks is extremely tricky at high redshift. Neither morphological nor kinematic classifications alone can unambiguously exclude mergers. Morphological classifications of mergers rely on footprints of interactions such as bridges, tidal tails, and double nuclei. These features can be easily missed at high redshift because of surface brightness dimming, size evolution, and band shifting (e.g., Hibbard & Vacca 1997; Hung et al. 2015). High-redshift galaxies are clumpier and show more irregular structures than local galaxies, further complicating the visual characterization of mergers (Elmegreen et al. 2007). Kinematic classifications assume that isolated disks exhibit smooth velocity gradients whereas mergers show more asymmetric and chaotic kinematic features (Colina et al. 2005; Dicaire et al. 2008; Shapiro et al. 2008). However, this assumption does not consider post-coalescence mergers, which may also display disk-like kinematics (Bellocchi et al. 2012). A "morpho-kinematic" classification that combines the morphological and kinematic criteria is proposed as a more robust approach (Rodrigues et al. 2017).

The velocity field of A1689B11 is consistent with a rotating isolated disk based on current morphological and kinematic classification schemes. We first use the kinematic criteria of the SINS survey (Shapiro et al. 2008) and derive V_asym and ${\sigma }_{\mathrm{asym}}$ from the 2D velocity and velocity dispersion map. We find ${V}_{\mathrm{asym}}=0.1$ and ${\sigma }_{\mathrm{asym}}=-0.08$, placing A1689B11 in the isolated disk region of the V_asym versus ${\sigma }_{\mathrm{asym}}$ diagram. The unique clump identification method presented in Section 3.2 rejects most of the minor-merger-like features as either foreground or background sources. For clumps that do associate with A1689B11, their velocities are consistent with rotating along with the main disk and their velocity dispersions do not show deviations from the rest of the disk, inconsistent with the scenario that the clumps are mergers. To further rule out the scenario of post-coalescence mergers, we apply the five "morpho-kinematic" criteria of Rodrigues et al. (2017). We find that A1689B11 satisfies all five criteria of an isolated rotating disk: (1) the velocity map has a single velocity gradient (Section 3.1.2 and Figure 5), (2) ${V}_{\mathrm{rot}}/{V}_{\sigma }\gt 1$ (=9–13 for A1689B11), (3) there is a ${V}_{\sigma }$ peak coinciding with the center of rotation (Section 3.1.2 and Figure 4), (4) there is no mismatch between the kinematic and morphological PAs, i.e., ${\rm{\Delta }}{PA}\lt 30^\circ$ (${\rm{\Delta }}{PA}=1^\circ$ for A1689B11; Sections 3.1.1 and 3.2; Table 1), and (5) the rotation center matches the stellar mass center within 04 (≲01 for A1689B11; Sections 3.1.1 and 3.2; Figure 5). Based on these arguments, we therefore exclude mergers as the origin of A1689B11.

Spiral galaxies and irregular galaxies are two distinct morphological classes in the local universe (e.g., Hubble 1926). However, as the general morphology of distant galaxies becomes more chaotic and irregular (e.g., Abraham et al. 1996a, 1996b; Conselice et al. 2005; Elmegreen et al. 2007; Shibuya et al. 2016), spirals and clumpy/irregular galaxies do not have to be mutually exclusive at high redshift. For example, studies focused on spiral morphologies show that spiral structures are highly disturbed and arms are less well developed at $z\gtrsim 0.5$ (Abraham & van den Bergh 2001; Elmegreen et al. 2005). The less well-defined spiral structures at high redshift could be caused by intrinsic evolutions and/or observational biases.

Spiral arms are the main sites of star formation in the local universe. The brightness, sizes, and surface densities of SF regions are much larger at high redshift (e.g., Elmegreen et al. 2005; Jones et al. 2010b; Genzel et al. 2011; Swinbank et al. 2012; Fisher et al. 2014). The size and brightness of the SF clumps in the spiral arm of A1689B11 are comparable to those of $z\gt 1$ clumpy SF galaxies (T. Yuan 2017, in preparation), producing a much clumpier appearance than local spiral arms. In addition, A1689B11 is gravitationally lensed; the reconstructed source-plane morphology is much sharper than a nonlensed case, making it look different from most z ∼ 2 SF galaxies.

To demonstrate the effect of large SF clumps and observational effects such as surface brightness dimming and gravitational lensing on the appearance of spiral arms at high redshift, we manually redshift a local spiral galaxy to z = 2.54 and compare its redshifted morphology with the lensed and nonlensed cases of A1689B11 (Figure 9). The local spiral galaxy template (${\rm{G}}04\mbox{--}1$) is chosen from the DYNAMO (DYnamics of Newly Assembled Massive Objects) sample. The DYNAMO sample is a local analog of turbulent, clumpy disk galaxies at high redshift (Green et al. 2010, 2014; Bassett et al. 2014). ${\rm{G}}04\mbox{--}1$ is one of the few galaxies in the DYNAMO sample that have a spiral morphology. ${\rm{G}}04\mbox{--}1$ has similar half-light radius (${R}_{1/2,{\rm{H}}\alpha }=2.7$ kpc), SF clump size, and brightness to A1689B11 (Fisher et al. 2017). We use the HST Hα narrowband image (pixel size $0\buildrel{\prime\prime}\over{.} 05$) of ${\rm{G}}04\mbox{--}1$, as it best represents the clumpy SF morphology (Figure 9, panel (1)). Panel (2) of Figure 9illustrates the mocked HST/WFC3 (pixel size $0\buildrel{\prime\prime}\over{.} 1$) morphology of ${\rm{G}}04\mbox{--}1$ at z = 2.54. We assume no intrinsic size evolution because the DYNAMO sample consists of compact objects in the local universe, analogous to the sizes of galaxies $z\sim 2$. Because ${\rm{G}}04\mbox{--}1$ is nearly face-on, the morphology of A1689B11 is deprojected in panels (3) and (4) for a better comparison. Panel (3) of Figure 9 shows the best source-plane reconstructed morphology of A1689B11 based on the HST/F814W band image. Owing to gravitational lensing, the effective spatial resolution on the source plane is increased by a factor of 2–3 and the S/N of the image is increased by ∼7. To compare with the nonlensed image of ${\rm{G}}04\mbox{--}1$ at z = 2.54, we show in panel (4) what A1689B11 would have looked like without lensing magnification in the HST/WFC3 IR band.

Zoom In Zoom Out Reset image size

Panels (1) and (2) of Figure 9 underline the importance of separating observational effects from intrinsic evolution of spiral arms. In the local universe, spiral galaxies are divided into three classes based on the number of arms: grand-design (two-arm), many-arm (multiple distinct global arms), and flocculent (multiple less distinct arms) (e.g., Elmegreen et al. 1982; Hart et al. 2017). According to the local image of panel (1) of Figure 9, ${\rm{G}}04\mbox{--}1$ is classified as a many-arm spiral, as is more than half of the local spiral galaxy population (Davis et al. 2014). As a result of surface brightness dimming, ${\rm{G}}04\mbox{--}1$ would be mostly likely classified as a rare one-arm spiral, which are only seen in 14% of local spirals (Davis et al. 2014). It is possible that A1689B11 has more arms like ${\rm{G}}04\mbox{--}1$, but only the longest arm is visible at this redshift and with this magnification.

Panels (3) and (4) of Figure 9 demonstrate the power of gravitational lensing in bringing the otherwise unseen spiral structures into focus. The spiral structure is still detectable in the redshifted image of ${\rm{G}}04\mbox{--}1$, whereas the spiral arm of A1689B11 is barely distinguishable in the unlensed case of panel (4). Note that although ${\rm{G}}04\mbox{--}1$ and A1689B11 share similar galaxy size and SF clump brightness, ${\rm{G}}04\mbox{--}1$ is ∼10 times more massive than A1689B11 and has a large Hα velocity dispersion of ∼50 km s⁻¹. There are noticeable differences in the morphological appearance of panels (3) and (4), implying that the intrinsic spiral structures of A1689B11 are less well developed than those of ${\rm{G}}04\mbox{--}1$. While a detailed comparison of the host properties of ${\rm{G}}04\mbox{--}1$ and A1689B11 is beyond of the focus of the current paper, Figure 9 simply demonstrates that the morphology of A1689B11 is consistent with a spiral galaxy. The effects of gravitational lensing, surface brightness dimming, and larger SF regions at high redshift combine together to make the somewhat unique appearance of A1689B11.

Finally, separating spiral arms from clumpy irregular morphologies requires a quantitative classification scheme. Such a scheme is not available yet at high redshift. For example, there are two geometric parameters that are commonly used in quantifying local spiral arms: the number of arms (or harmonic modes) and the pitch angle. A large fraction of local spirals can be modeled by superpositions of logarithmic spiral functions. Based on this, automatic spiral arm finding and logarithmic function fitting tools have been developed and applied in local galaxy surveys (e.g., Davis et al. 2012; Davis & Hayes 2014; Shields et al. 2015; Hart et al. 2017). Our ongoing effort of testing and adapting these tools at high redshift will help to quantify spiral arms and to separate spiral from irregular structures objectively. We show in Figure 10 a preliminary best-fit logarithmic spiral function created by manually masking pixels that form the spiral arm. We find a large pitch angle ($\theta =37^\circ \pm 2^\circ$) for the arm. We test the automated arm detection and fitting routine of SpArcFiRe (Davis & Hayes 2014) on A1689B11 and notice that the number of arms and pitch angle depend sensitively on the bulge–disk decomposition, lensing PSF reconstruction and masking of noise. We will report the pitch angle analysis in our future work (T. Yuan 2017, in preparation). To conclude, we favor the interpretation of a spiral galaxy with primitive spiral arms developing in A1689B11 instead of a merger or an irregular galaxy.

Zoom In Zoom Out Reset image size

The total SFR (22 M_⊙ yr⁻¹) and stellar mass (10^9.8 M_⊙) place A1689B11 as a typical SF main-sequence galaxy at $z\sim 2$. However, the kinematic properties of A1689B11 are quite "mature" compared to $z\gt 1$ SF galaxies and are more akin to local spiral galaxies.

We compare the gas velocity dispersion and stellar mass of A1689B11 with other samples in Figure 11. We collect Hα gas velocity dispersion measurements on disk galaxies from a few representative IFS surveys at $z\gtrsim 2$: the SINS disks (Förster Schreiber et al. 2006, 2009; Cresci et al. 2009), the KMOS-3D disks (Wisnioski et al. 2015), and gravitationally lensed galaxies (Jones et al. 2010b; Livermore et al. 2015). Only the rotation-dominated (${V}_{\mathrm{rot}}/{V}_{\sigma }\gt 1$) disks from these surveys are shown. The $z\sim 0$ sample is taken from a local reference sample of rotating spiral galaxies with Hα velocity dispersion measurements (Epinat et al. 2010).

Zoom In Zoom Out Reset image size

The local spirals have low Hα velocity dispersions irrespective of their stellar masses. There is a large spread in velocity dispersions of $z\gtrsim 2$ disks, with a median value that is ∼2.5 times larger than local spirals. We also mark the thermal broadening of Hα-emitting gas (10⁴ K) as the minimal velocity dispersion that can be measured from the Hα line. The systematic uncertainty in comparing velocity dispersions across samples could also contribute to the large scatter in high-redshift samples. The major systematics come from the method used to correct for beam smearing, the disk model, and the radius at which the velocity dispersion is measured and the weighting used. The velocity dispersion of A1689B11 is considerably lower than lensed galaxies of similar masses. It is also much lower than the median of all $z\gtrsim 2$ samples. The large ${V}_{\mathrm{rot}}/{V}_{\sigma }$ (∼9–13) of A1689B11 is similar to local spirals and ∼2–5 times larger than the median value of SF galaxies at $z\gtrsim 1$ (e.g., Cresci et al. 2009; Swinbank et al. 2017).

One of the major observational results from IFS surveys of high-redshift ($1\lesssim z\lesssim 3$) galaxies is that the intrinsic gas velocity dispersions as usually measured from optical ionized gas are significantly higher than local SF galaxies (e.g., Förster Schreiber et al. 2009; Law et al. 2009; Lehnert et al. 2009; Genzel et al. 2011; Wisnioski et al. 2015). Popular explanations for the enhanced velocity dispersions include star formation feedback (Lehnert et al. 2009; Green et al. 2014), gravitational instability (Agertz et al. 2009; Ceverino et al. 2010), a combination of star-formation-feedback-driven and gravitationally driven turbulence (Krumholz & Burkhart 2016; Krumholz et al. 2017), a multitude of physical drivers of turbulence (Federrath & Klessen 2012; Federrath et al. 2017a), and cosmological cold gas accretion (Bournaud & Elmegreen 2009; Dekel et al. 2009; Genzel et al. 2011; Genel et al. 2012).

If the ionized gas turbulence of A1689B11 is directly driven by energy ejected from star formation through a supernova, then its SFR surface density suggests a velocity dispersion of ${V}_{\sigma }\sim 40\mbox{--}55$ km s⁻¹ (e.g., Lehnert et al. 2009; Green et al. 2014), depending on how the energy is dissipated into the ISM and assuming a conservative supernova feedback efficiency of 0.25 (Dib et al. 2006; Zhou et al. 2017). In this simple model, the small velocity dispersion of A1689B11 would imply that the supernova feedback efficiency is a factor of ∼2 lower than that of local galaxies, which is difficult to explain.

In the slightly more complicated feedback-driven turbulence model of Krumholz & Burkhart (2016), where ${V}_{\sigma }$ is a function of SFR and Q_gas, the small velocity dispersion of A1689B11 can be reproduced at a ${Q}_{\mathrm{gas}}\sim 0.5$, matching our rough estimation of Q_gas in Section 3.3. The gravity-driven turbulence model of Krumholz & Burkhart (2016) is disfavored because the small ${V}_{\sigma }$ would require an unrealistically large gas fraction ($50 \% \lt {f}_{\mathrm{gas}}\sim 100 \%$). Future direct observations of molecular gas and dust maps of A1689B11 would provide a more robust measurement on Q_gas and gas fraction to distinguish various turbulence driven models.

We also show in Figure 11 the location of the spiral galaxy BX 442 at $z=2.17$ from Law et al. (2012). BX 442 has a high velocity dispersion (${V}_{\sigma }\sim 70$ km s⁻¹) and is thought to have a short-lived spiral triggered by a minor merger (Law et al. 2012). Two conditions must be satisfied according to the minor-merger trigger mechanism. First, the galaxy must be massive enough to stabilize the formation of an extended disk. Second, a nearby merging satellite must be properly orientated and sufficiently massive to excite the spiral feature. A1689B11 is ∼10 times less massive than BX 442 and is not consistent with a merging system as discussed in Section 4.1. Even if the minor-merger scenario works for BX 442, an alternative mechanism is required to account for the existence of A1689B11.

Being an outlier in both the velocity dispersion distribution and morphology (spiral structure) of high-redshift galaxies, A1689B11 provides a few interesting angles to revisit the formation of spiral arms.

Classic theories of spiral arm formation require a dynamically cool and thin disk. In the paradigm of the density wave theory, the amplitude of the induced density wave becomes too weak if the gas velocity dispersion is too high or if the disk is too thick (e.g., Lin & Shu 1964; Elmegreen & Thomasson 1993; Bertin & Lin 1996; Rafikov 2001; Bottema 2003; Sellwood 2014). In the alternative theory of swing amplification or self-gravitational instability, a large velocity dispersion (in the context of a large ToomreQ parameter) dampens the amplification and a thick disk reduces the disturbance gravity that are seeds of the swing amplification (e.g., Julian & Toomre 1966; Elmegreen & Thomasson 1993; Bottema 2003; Sellwood 2014). In the local universe, almost all spiral arms reside in the thin disk (scale height 200–300 pc), with a dynamically cool stellar (${V}_{\sigma }$(star) ∼ 20 km s⁻¹) and gas component (${V}_{\sigma }$(Hα) ∼20–25 km s⁻¹) (Glazebrook 2013). For stellar disks that are dynamically "warm" and have a non-negligible thickness, spiral activity develops mostly in the cold gas component (Bertin & Lin 1996). The exceptions are mergers or bar-induced spirals, where the external driver plays a more important role than internal disk dynamics (e.g., Kormendy & Norman 1979; Bottema 2003; Dobbs et al. 2010).

Both the density wave and swing amplification theory have seen a certain degree of success in explaining local spirals where the model assumptions are easily satisfied. The logarithmic spiral function predicted from the density wave theory has been widely used and confirmed in observations of local spirals, whereas swing amplification is commonly used to explain transient and less regular spirals. The main difference is that density wave theory predicts quasi-steady long-lived arms (wave arm) and swing amplification or its various forms of gravitational instability predicts short-lived reoccurring arms (material arm). It is difficult to constrain the lifetime of the spiral arm observationally and thus hard to distinguish the two theories (e.g., Sellwood 2011). Current simulations seem to prefer either swing amplification or merger in explaining the formation of spiral arms at high redshift (e.g., Law et al. 2012; Fiacconi et al. 2015). For example, self-gravitational instability and fragmentation have been proposed to model the local DYNAMO spiral ${\rm{G}}04\mbox{--}1$ (Section 4.2), which has a high velocity dispersion and clumpy SF regions analogous to high-redshift disks (Inoue & Yoshida 2017). A minor merger is thought to trigger the short-lived spiral structure of BX 442 (Law et al. 2012).

It is reasonable to expect that the quasi-steady spiral arm from the density wave theory is suppressed at high redshift because of the high gas velocity dispersions and geometrically thick disks (Elmegreen & Elmegreen 2006; Cresci et al. 2009). The rarity of spirals at $z\gtrsim 2$ could be partly accounted for by the short-lived arms from either gravitational instability or mergers. On the other hand, the discovery of a dynamically cool disk like A1689B11 could mean that the condition for the classic density wave spirals to develop can exist at $z\gt 2$. The small velocity dispersion of A1689B11 is indicative of a thin disk (Figure 11). Our preliminary fitting of a logarithmic spiral function (Figure 10) is consistent with a density-wave-triggered primitive spiral arm of A1689B11. We speculate that A1689B11 belongs to a population of rare spiral galaxies at $z\gtrsim 2$ that mark the earliest epoch of thin-disk formation (Freeman & Bland-Hawthorn 2002; Freeman 2012; Kraljic et al. 2012; Elmegreen et al. 2017). Note that the spiral structure of A1689B11 would not be visible with current observational capacity without gravitational lensing (Figure 9). Future observations with the James Webb Space Telescope (JWST) will help reveal this population and investigate the earliest onset of spiral arms and thin disks.

Modern theoretical efforts in understanding the origin of spiral arms have focused on semianalytic or N-body simulations in isolated disks (e.g., Wada et al. 2011; Baba et al. 2013; D'Onghia et al. 2013; D'Onghia 2015). Cosmological simulations do not have sufficient resolutions to trace the detail dynamics of spiral arms; however, it is interesting to explore the spiral arm formation from the context of cosmological simulations where the role of environment is included.

Recent cosmological zoom-in simulations (Cen 2014) provide valuable insights into the emergence of the Hubble sequence across cosmic time. Cen (2014) shows the rich physics of cold gas accretion dynamics onto galaxies and evolutionary trends. At $z\gtrsim 2$, the average in situ cold gas accretion streams through the galactic halo virial sphere surface can be characterized by (1) multiple cold streams, (2) high accretion rates, (3) low angular momenta, and (4) high gas densities. Among these four parameters, it is suggested that spiral structures are most sensitive to (1) the number of concurrent streams in the gas accretion. Flocculent spirals only begin to significantly appear at $z\sim 1\mbox{--}2$, when the number of major gas streams is about 2–3, and two-arm grand-design spiral galaxies appear at $z\leqslant 1$ when the number of major cold streams reduces to 1. In the framework of Cen (2014), the average number of gas streams decreases from high to low redshift. However, at each redshift, there is a distribution of in situ environments pertaining to the range of cold streams. It is therefore possible to find "evolved" spiral galaxies at high z and "unevolved" high-redshift analog galaxies at low z at the tail of the environmental distributions.

In the framework of Cen (2014), spiral A1689B11 may be formed in an in situ environment that is characterized by probably one major stream with high accretion rates and gas densities for an extended period of time. Because the number of significant cold streams is correlated with the degree of interactions among galaxies, it is expected that the velocity dispersion of the stellar disk would be positively correlated with the number of cold streams (Cen 2014). This scenario is consistent with accretion energy being one of the main drivers of the turbulence in disks at $z\sim 2$ (e.g., Birnboim & Dekel 2003; Bournaud & Elmegreen 2009; Dekel et al. 2009; Genel et al. 2012).

However, galaxy growth through cold-mode accretion has been challenged by other cosmological simulations (e.g., Nelson et al. 2013; Genel et al. 2014). Whether the Hubble sequence originates from the cold flow process or other mechanisms such as mergers and feedback is highly controversial (e.g., Genel et al. 2015). Moreover, extra care must be given when studying the origin of the spiral arms in cosmological simulations. Spiral features in current cosmological simulations may originate from unphysical perturbations and are sensitive to the resolution and the detailed prescriptions of ISM models. With a larger sample of high-redshift spirals in observations and larger volumes of cosmological zoom-in simulations with subparsec resolutions, we should be able to answer the question of whether the in situ environment plays a critical role in the formation of spiral galaxies.