The formation of Milky Way-mass disk galaxies in the first 500 million years of a cold dark matter universe

Whether among the myriad tiny proto-galaxies there exists a population with similarities to present day galaxies is an open question. We show, using BlueTides, the first hydrodynamic simulation large enough to resolve the relevant scales, that the first massive galaxies to form are %in fact predicted to have extensive rotationally-supported disks. Although their morphology resembles in some ways Milky-way types seen at much lower redshifts, these high-redshift galaxies are smaller, denser, and richer in gas than their low redshift counterparts. From a kinematic analysis of a statistical sample of 216 galaxies at redshift $z=8-10$ we have found that disk galaxies make up 70\% of the population of galaxies with stellar mass $10^{10} M_\odot$ or greater. Cold Dark Matter cosmology therefore makes specific predictions for the population of large galaxies 500 million years after the Big Bang. We argue that wide-field satellite telescopes (e.g. WFIRST) will in the near future discover these first massive disk galaxies. The simplicity of their structure and formation history should make possible new tests of cosmology.


INTRODUCTION
The search for "primeval" galaxies and the understanding that will come from studying them has long been a guiding principle for much of extra-galactic astrophysics and cosmological research. Theories for the formation of the first galaxies span early models of monolithic collapse of massive bodies (Eggen et al. 1962) to merger scenarios for the hierarchical formation of all structures (Press & Schechter 1974). Disk galaxies are arguably those that have generated the most interest, and were the focus of early theoretical work (Fall & Efstathiou 1980). In the simplest galaxies, the structure of rotating disks should be directly linked to the angular momentum of the surrounding dark matter halo. Finding galaxies for which this is true, in the earliest environments possible would enable powerful tests of the dynamics, gaseous and radiative processes involved in galaxy formation.
In current cosmological models, galaxies form from the gravitational collapse of small perturbations in the matter distribution (Silk & Mamon 2012). This process involves both a hierarchy of merging structures and smooth accretion, so that early galaxies are predicted to be morphologically irregular, clumpy, and compact (Shimizu et al. 2014;Cen & Kimm 2014). This is supported by recent observational data on samples of galaxies at red-shift z = 8 and beyond (Oesch et al. 2010;Ono et al. 2013;Curtis-Lake et al. 2014). The volumes accessible to these studies, both computational and observational are however thousands of times smaller than those that will be probed by upcoming telescopes (Spergel et al. 2015).
The earliest that large rotating disks have been seen observationally is at redshifts z = 2 − 3 (Genzel et al. 2006;Glazebrook 2013), and this has required deep spectroscopy on 10m-class telescopes. At higher redshifts, galaxies in published samples are very compact, with half light radii ∼ 0.5 kpc at z = 8 (Oesch et al. 2010;Holwerda et al. 2014). The fraction of morphologically disturbed galaxies is also known to be high in deep fields (Menanteau et al. 2006). The volumes probed at redshifts z > 3 and particularly the redshifts at the current frontier, z = 8 − 10 are however very small, with a sky areas of 11 sq. arcmin accessible to the Hubble Extreme Deep Field (Illingworth et al. 2013) for example. This has meant that probing the highest luminosity and largest galaxies, has not been possible so far. It is not surprising that observational studies of the structure and morphologies of the highest z galaxies have not revealed large disk galaxies, whether they exist or not.
Theoretically, the situation is similar -it is not known what the prevailing Cold Dark Matter model predicts for the most massive galaxies at redshifts z ≥ 8. Numerical simulations used to make predictions have so far (Jaacks et al. 2012;Cen & Kimm 2014) been limited to volumes 40-100 times smaller than our present work, or else zoomed simulations concentrating on individual galaxies (Pawlik et al. 2011) or overdensities (Yajima et al. 2014). Intriguingly, a dominant fraction of the early forming (dwarf-mass) galaxies simulated so far show evidence of disk-like structure(Romano-Díaz et al. 2011; Figure 1. A sample of disk galaxies selected from the BlueTides simulation at redshift z = 8. We show galaxies identified kinematically to be rotating disks. For each galaxy we show a face-on and side-on view. The colors of top two rows represent the stellar surface density coded by the stellar age (blue to red: 0 ∼ 300 Myear), face-on and side-on. The colors of the bottom two rows represent the star-formation surface density (normalized to 0 ∼ 1). The face-on images have included the effect of dust-extinction. The darker central regions of the disks are mostly caused by the consumption of gas to form bulge stars. There is also an obscuration effect from the higher metallicity (and therefore dust content) of the central regions. Pawlik et al. 2011). It is of prime importance to reach the regime of large galaxies as the early formation of luminous or massive objects (such as first bright quasars, and most massive galaxy clusters) is one of the most potent tests of models of structure formation in hierarchical cosmology. In order to do so, and simulate a statistically useful sample of galaxies, it is necessary to model a volume of order several hundred megaparsecs on a side. This must be done with high enough mass and force resolution to allow studies of structure and morphology (See e.g., Governato et al. (2007)). We have been able to reach this regime, with the BlueTides run, a cosmological hydrodynamic simulation which we have just completed, and we describe the first results below.

THE BLUETIDES SIMULATION
BlueTides was carried out using the Smoothed Particle Hydrodynamics code MP-Gadget (see Springel (2005); Di Matteo et al. (2012) for earlier versions of the code) with 2 × 7040 3 (0.69 trillion) particles on 648,000 Cray XE compute cores of the Blue Waters system at the National Center for Supercomputing Applications (NCSA). The simulation evolved a cube of side length 400 comoving Mpc/h to redshift z = 8, and is the largest cosmological hydrodynamic simulation yet carried out in terms of memory usage, by an order of magnitude (Di Matteo et al. 2012). We refer to Feng et al. (2015) for the physical models employed and other aspects of the simulation in more detail. We also show there that the global star formation rate predicted is consistent with current observations at redshifts z > 8. The volume simulated is approximately 300 times larger than the largest observational survey at these redshifts (Trenti 2012).
Galaxies were selected from the raw particle data (40 TB per output time) using a friends-of-friends algorithm at a range of redshifts from z = 8 − 13. The mass resolution of the simulation is 1.72 × 10 7 M per dark matter particle, so that at redshift z = 8, 36 million galaxies were found containing 64 dark matter particles or more, above a mass threshold of 10 9 M . The BlueTides galaxy luminosity functions at the redshifts of interest here are fully consistent with current observations, and this is shown in Feng et al. (2015). In the present case, we are most interested in the most massive galaxies with at least 10 4 stellar particles, which contain enough particles to study their rotational support. The corresponding stellar mass is 0.8 × 10 10 M or more.
We have made images of all galaxies present in the simulation, by adaptively smoothing the star particle distribution to reveal as much structure as possible given the resolution available. We have also used the star formation rate (SFR) of the gas particles in the simulation to make other images. The SFR can be related to the restframe UV (150 nm) luminosity. The UV luminosity is sensitive to the presence of dust, and there is evidence that the highest luminosity early galaxies are significantly obscured (Wilkins et al. 2013;Cen & Kimm 2014). We account for this dust attenuation using a screening model (Joung et al. 2009). This simple model assumes that the dust attenuation A UV in a galaxy pixel is proportional to the metal-mass density in that pixel. There is only one free parameter, which we normalize by setting A UV = 1 for observed M UV = −21 galaxies (Wilkins et al. 2013). We have taken this approach, rather than performing a more detailed radiative transfer simulation because the current major uncertainties lie in the dust content and attenuation curve (Wilkins et al. 2013).

GALAXY MORPHOLOGY
We turn directly to the morphology of the most massive (M stellar ≥ 0.7 × 10 10 M ) galaxies. We compute the inertia tensor of the star particles and use these to plot views of each galaxy parallel and perpendicular to the major axis. Studying them reveals that a significant fraction are visually disk-like, and surprisingly regular in shape.
In Figure 1 we show images of 10 disk galaxies in the simulation at z = 8, both side-on and face-on for each. These are among those kinematically selected to be disks (see below). We can see that the disk components extend to a diameter of ∼ 4 kpc, with a median half-light radii ∼0.6 kpc for those shown. Although the z = 8 BlueTides disk galaxies have masses comparable to the Milky Way, their sizes are significantly smaller, by a factor of ∼ 5. The disk galaxies also have mean gas fractions of 0.85 at z = 8 and 0.83 at z = 10. The stellar populations in the central regions of each galaxy are significantly older than the outer disk, with a proportion of stars that formed as early as z = 13, i.e. age 300 Myr at z = 8.
Kinematic decomposition of the galaxies lets us know if the disk-like structures seen in Figure 1 are rapidly rotating. We have used the standard technique (Governato et al. 2007) to determine the fraction of stars in each galaxy which are on planar circular orbits, and which are associated with a bulge. We note that this technique is not sensitive to the clumpiness of disks. At redshift z = 8 we find that 70% of galaxies above a mass of 10 10 M are kinematically classified as disks (using the standard threshold of disk stars to total stars (D/T) ratio of 0.2 (Governato et al. 2007). This percentage of disks changes with redshift as shown in Figure 2. In all cases it is much higher than the fraction of observed galaxies of comparable mass which are seen to be disk-like at z = 0, 0.14 ± 0.03 (Buitrago et al. 2013).
We note that the majority of gas-rich disks seen at high redshifts accessible so far to spectroscopy (z=2-3) are clumpy, turbulent and thick (Förster Schreiber et al. 2009) (Genzel et al. (2006) is an exception), and it is natural to expect this to be the case with the BlueTides massive disks, which have high gas fractions. A relevant statistic we have computed is the ratio of circular velocity, V to vertical velocity dispersion σ of our disk galaxies. We have computed this V /σ ratio using information from the disk star particles and plot a histogram of the results in Figure 3. We can see that the values range from Figure 3. Disk galaxy kinematics in BlueTides characterized by the ratio of circular velocity V to vertical velocity dispersion σ. We show histograms of results for this V /σ ratio computed using star particles in the disks of galaxies selected to be disks from their D/T ratio. 3 different redshifts are shown. At redshift z = 10 results for only 3 galaxies are available and they are plotted as individual points V /σ = 2 to 5, with a mean value of 3.08 ± 0.08 measured from all redshifts. The galaxies are therefore definitely disks, having V /σ well above the 0 − 1 expected for nondisk galaxies. This supports our conclusions made using the D/T ratio above. The V /σ values are lower than the V /σ ∼ 10 expected for thin spiral galaxies seen at z = 0, and instead consistent with gas rich turbulent disks from earlier times (Förster Schreiber et al. 2009). Although this is reasonable, confirmation of this aspect awaits future higher resolution resimulations which are able to resolve the fine vertical structure of the disks.

HALF-LIGHT RADIUS
A number of recent studies (Oesch et al. 2010;Ono et al. 2013;Holwerda et al. 2014;Kawamata et al. 2014) have begun to make measurements of galaxy structure and in particular sizes at the highest redshifts. As the observational volumes are small, the magnitude range of observed samples varies from M UV = −21 to M UV = −18. We can therefore compare BlueTides galaxies with these measurements, bearing in mind that the brightest galaxies in the simulation are much rarer than can be observed in current surveys. Using the standard observational software SExtractor (Bertin & Arnouts 1996), we measure the half-light radius (r 1/2 ) of each galaxy. This is shown in Figure 4, where we compare our redshift z = 8 results with observational data from (Kawamata et al. 2014), (Holwerda et al. 2014) and (Curtis-Lake et al. 2014). We can see that the current observations cover the lower end of the magnitude range, and are in the same realm as (but not a good match to) our simulation results, with a mean r 1/2 ∼ 0.6 for M UV ∼ −20.5 → −21.5. Galaxies are therefore extremely compact at these high redshifts. The galaxy sizes in the simulation do appear to be on the low end compared to observations, something which has also been seen at z=3 (Joung et al. 2009) and could be a sign of interesting physics. It has also been shown however(Curtis-Lake et al. 2014) that selection effects, flux cuts and choice of size measurement algorithms can affect measured galaxy sizes by factors of 2 or more. In future work we will explore this in more detail with mock observations. Interestingly, the simulation predicts a positive size-luminosity correlation.

CONCLUSIONS
At redshifts z = 8 − 12 the Universe is about 1000 times denser than the present day. In these very different physical conditions, and in the very short time available, massive disk galaxies are able to form in the cold dark matter model. They share many visual and kinematic characteristics with those that are seen when the Universe was 20 times older. This means that an image one might have of the early Universe containing only small protogalaxies, merging clumps and irregular structures is not appropriate.
An immediate question we can ask is how these disks did in fact form in the simulation. We have tracked particles in the most massive galaxies (including those plotted in Figure 1) back in time to visually assess their formation history. We have found that of the most massive 20 galaxies at redshift z = 8, only one was the product of an obvious major merger between similarly sized progenitor galaxies. A popular picture of gas rich major mergers producing some disks (Springel & Hernquist 2005) therefore does not appear to be relevant here. Given that the distribution of matter around galaxies is known to be highly anisotropic in the CDM model, the infall of gas cannot proceed quasi-spherically as in the classical picture of galaxy formation (Fall & Efstathiou 1980). Instead, smooth infall and accretion build up the disks, with cold gas arriving along the directions of nearby filaments in the density distribution.
These massive disks will be fascinating objects to study with next generation telescopes. The WFIRST satellite (Spergel et al. 2015) will have a field of view 200 times that of the WFC3 instrument on HST, making a planned deep (J=26.7) sky survey covering 2000 square degrees possible. Using the luminosity function in BlueTides and its evolution we predict that this survey should find ∼ 8000 disk galaxies of the type we have kinematically identified as massive disks, about 1 per field of view on average. This compares to a probability of 0.3 for finding a single one of these objects in the current largest area HST survey (BoRG ). The rotation curves will be exciting targets for Thirty Meter class ground based telescopes. Forming from lightly enriched material, these primeval disk galaxies are special objects. Observing them and comparing them with simulations and analytic theory should enable us to make inferences about dark matter, the role of angular momentum and the fundamental principles of galaxy formation (such as tidal torques, and smooth accretion of unenriched material) which are not possible in the later, more complex Universe.