Rapid filamentary accretion as the origin of extended thin discs


 Galactic outflows driven by stellar feedback are crucial for explaining the inefficiency of star formation in galaxies. Although strong feedback can promote the formation of galactic discs by limiting star formation at early times and removing low angular momentum (AM) gas, it is not understood how the same feedback can result in diverse objects such as elliptical galaxies or razor thin spiral galaxies. We investigate this problem using cosmological zoom-in simulations of two galaxies forming within 1012 M⊙ haloes with almost identical mass accretion histories and halo spin parameters. However, the two resulting galaxies end up with very different bulge-to-disc ratios at z = 0. At z > 1.5, the two galaxies feature a surface density of star formation ΣSFR ≃ 10 M⊙ yr−1 kpc−2, leading to strong outflows. After the last starburst episode, both galaxies feature a dramatic gaseous disc growth from 1 to 5 kpc during 1 Gyr, a decisive event we dub ‘the Grand Twirl’. After this event, the evolutionary tracks diverge strongly, with one galaxy ending up as a bulge-dominated galaxy, whereas the other ends up as a disc-dominated galaxy. The origins of this dichotomy are the AM of the accreted gas, and whether it adds constructively to the initial disc angular momentum. The build-up of this extended disc leads to a rapid lowering of ΣSFR by over two orders of magnitude with ΣSFR ≲ 0.1 M⊙ yr−1 kpc−2, in remarkable agreement with what is derived from Milky Way stellar populations. As a consequence, supernovae explosions are spread out and cannot launch galactic outflows anymore, allowing for the persistence of a thin, gently star-forming, extended disc.


INTRODUCTION
Galaxies have been observed with a variety of shapes and morphologies ranging from ellipticals to spirals. Although it is yet unclear what are the most important formation mechanisms that determine the final morphology of a galaxy, it is thought that both large-scale processes, e.g. galaxy interactions and gas accretion, and small-scale processes such as star formation and stellar feedback play important roles. The effect of these different mechanisms interact in a non-linear way, making it difficult to disentangle what fundamental physical process really governs their evolution. Recent models of galaxy formation (Stinson et al. 2006;A g e r t ze ta l .2013 ;Teyssier et al. 2013;Hopkins et al. 2014;C a p e l oe ta l .2018) have demonstrated that supernovae-driven feedback has a particularly strong impact on the formation of galaxies (for a recent review, see Naab & Ostriker 2017). It is believed that supernovae (SNe) explosions release enough energy and momentum in their surrounding medium to eject large quantities of gas from the galaxy into its host halo. Other feedback mechanisms, such as stellar winds and UV radiation from massive stars contribute as well, although the momentum provided is smaller than the terminal momentum acquired after the Sedov-Taylor phase of SN explosions (Agertz et al. 2013).
Furthermore, star formation at small scales has been shown to be a crucial and delicate ingredient: The right amount of stars formed in ⋆ E-mail: michael.kretschmer@physik.uzh.ch the right conditions, at the right location and over the right time-scale will affect the fate of the galaxy. Moreover, since stellar feedback and star formation are closely related to one another, self-regulation mechanisms are likely to emerge and set the characteristic timescales of the global star formation and gas depletion rates (Hopkins et al. 2014(Hopkins et al. , 2018b;A g e r t z&K r a v t s o v2015, 2016; Muratov et al. 2015;Anglés-Alcázar et al. 2017;Semenov, Kravtsov & Gnedin 2017).
The study of galaxies similar to our own Milky Way has proven particularly challenging (Navarro & Steinmetz 2000;Abadi et al. 2003a;Governatoetal.2004). They are massive galaxies, with deep potential wells, but their overall star formation efficiency, although probably the largest in the entire galaxy population, remains very low, of the order of 20 per cent of the available baryons in the halo (Persic & Salucci 1992;Fukugita, Hogan & Peebles 1998;B e n s o n et al. Balogh et al. 2001;C o l ee ta l .2001;Baugh 2006;Guo et al. 2010;Mosteretal.2010;Behroozi, Wechsler & Conroy 2013).
On the contrary, present-day Milky Way mass galaxies are, broadly, either (1) quiescently star-forming disc galaxies, featuring thin gaseous discs with low levels of turbulence and no sign of strong galactic outflows or (2) gas-poor, almost non-star-forming (quenched) elliptical galaxies (Kormendy et al. 2009). How galaxies transition between the early, turbulent, outflow driven, evolution at z>1, to settling into large and thin galactic discs, seemingly undisturbed by vigorous feedback, is not yet understood. In the Milky Way, there is observational evidence that this transition may have been rapid (several hundreds of Myr, see Lehnert et al. 2014), posing a challenge for galaxy formation models.
Strong outflows with gas velocities larger than the escape velocity of their host halo have been observed for example by Heckman (2003) (see also the recent review on observations of galactic winds by Rupke 2018), with outflow velocities correlating with the surface density of star formation rate (SFR) SFR (Heckman & Borthakur 2016;Davies et al. 2019, and references therein). Interestingly, strong galactic outflows occur only above a certain threshold SFR . In local galaxies this was found to be SFR > 0.1M ⊙ yr −1 kpc −2 whereas at higher redshift this value was found to be larger, namely SFR > 0.5 − 1M ⊙ yr −1 kpc −2 (Newman et al. 2012;F orster Schreiber et al. 2019).
Such a phenomenon is found in theoretical studies using stratified box simulations (Walch et al. 2015;Girichidis et al. 2016)w h e r e strong outflows only occur when the number surface density of SN explosions in the disc exceeds a similar threshold -clustered SN events appear to be a key ingredient for powering strong outflows.
While strong feedback is a promising way to explain the low efficiency of galaxy and star formation, the same strong feedback can also destroy the thin galactic discs, leading to non-physical thick and turbulent discs at low redshifts (Governato et al. 2007;Agertz et al. 2011;Gibsonetal.2013;Stinson et al. 2013b;Roškar et al. 2014;Kretschmer & Teyssier 2020). Invoking strong feedback alone is therefore not sufficient. We need an additional mechanism to suppress feedback once large and thin discs have formed at low redshifts.
This paper aims at studying in detail the transition between highand low-redshift galaxy formation, i.e. between the early feedback dominated epoch to the present-day quiescent regime. We believe that this transition is the key to our understanding of the settling and existence of extended discs and the morphological differentiation between disc-dominated and bulge-dominated galaxies, as well as the observed epoch of 'disc settling' (z 1), where the neutral gas velocity dispersion in star-forming galaxies is found to decrease over time (Kassin et al. 2012).
In this work, we investigate this problem using cosmological zoom-in simulations of two galaxies forming in 10 12 M ⊙ haloes with almost identical mass growth histories. After a starburst episode triggered by the last major merger at z ∼ 1.5, the two galaxies similar evolutionary paths diverge, with one galaxy ending up as a compact bulge-dominated galaxy, whereas the other evolves to a disc-dominated galaxy analogue to our Milky Way. The reason for these different outcomes is the angular momentum (AM) of the accreted gas, and whether it adds destructively or constructively to the forming disc.
We have structured this paper as follows: In Section 2, we present our numerical methodology and summarize the properties of the subgrid models we have used. In Section 3, we present our results, focusing on the description of the last major starbursts around z ≃ 1.5, followed by the fast assembly of a large gaseous discs leading to the formation of a disc-dominated system analogue to our Milky Way. We describe how the formation of such a large disc naturally explains why SNe feedback stops to be efficient, leading to the maintenance of a quiescent regime. We systematically compare the properties of our disc-dominated system to the other simulation leading to a bugledominated galaxy. We finally discuss the implications of our work in Section 4 and conclude in Section 5.

Halo selection and initial conditions
Our study is based on two cosmological zoom-in simulations performed with the adaptive mesh refinement (AMR) code RAMSES (Teyssier 2002). Our methodology is as follows: We first ran an Nbody simulation with 512 3 dark matter particles in a periodic box of size 25 h −1 Mpc. From this volume we selected, at z = 0, two haloes with virial masses M vir = 1.1 × 10 12 M ⊙ and M vir = 1.4 × 10 12 M ⊙ in relative isolation; no halo more massive than half M vir within five virial radii exists. Note that the virial mass is calculated here using a spherical overdensity according to the definition of Bryan & Norman (1998).
We then generated higher resolution initial conditions around the selected haloes using the MUSIC code (Hahn & Abel 2011), with an initial hierarchy of concentric grids from ℓ min = 7, corresponding to a coarse grid resolution of 128 3 covering the entire periodic box, to ℓ max,ini = 11, corresponding to an effective initial resolution of 2048 3 . This yields a dark matter particle mass of m dm = 2.0 × 10 5 M ⊙ and an initial baryonic mass resolution of m bar = 2.9 × 10 4 M ⊙ .
We subsequently carried out our final simulations including gas and galaxy formation physics, as summarized below. The maximum resolution was set to ℓ max = 19 at z = 0, with refinement levels progressively released to enforce a quasi-constant physical resolution, where the smallest cells have sizes x min = 55 pc. The adopted refinement criterion is the traditional quasi-Lagrangian approach, namely cells are individually refined when more than 8 dark matter particles are present or when the total baryonic mass (gas and stars) exceeds 8 × m bar . Only the Lagrangian volume corresponding to twice the final virial radius of the halo was adaptively refined, the rest of the box being kept at a fixed, coarser resolution to provide the proper tidal field.
We intentionally picked haloes without any major mergers at z < 1 and with overall similar mass accretion histories. The haloes feature a strong merger-induced starburst at z = 1.4 and z = 2.0, respectively, followed by quiescent accretions history until z = 0. The final dark matter halo spin parameters, as defined in Bullock et al. (2001), are very similar for both haloes, namely λ = 0.015 and λ = 0.014. This fact is central to our study, as the resulting galaxies end up with very different bulge-to-disc ratios at z = 0, indicating that the dark matter halo spin is not a robust predictor for the final galaxy angular momentum and size (e.g. Jiang et al. 2019).

Galaxy formation physics
For our star formation recipe, we adopt a traditional Schmidt law, for which the SFR density is given bẏ where ρ is the density of the gas, t ff = √ 3π/(32Gρ) is its freefall time and ǫ ff is the star formation efficiency per free-fall time. Traditionally, ǫ ff is chosen to be a constant, close to 1 per cent, and star formation is allowed only in gas cells above a prescribed density threshold and below a prescribed temperature threshold. This is motivated by observations of inefficient star formation on galactic kpc scales (Bigiel et al. 2008) as well as in Milky Way giant molecular clouds (GMCs) (Krumholz & Tan 2007), albeit with a 1 dex spread (Murray 2011;Lee, Miville-Deschênes & Murray 2016;G r i s d a l e et al. 2019).
In this paper, we use a novel approach for which ǫ ff depends on the turbulent state of the gas following the so-called multifree-fall approach (Federrath & Klessen 2012;Semenov, Kravtsov & Gnedin 2016;Kretschmer & Teyssier 2020). Indeed, in a turbulent medium, like the interstellar medium (ISM), the gas density distribution is well described by a lognormal probability distribution function (PDF): where s = ln ρ/ρ 0 , σ s is the variance of s, ρ is the gas density, and ρ 0 the mean density. Gas with densities larger than a critical density s crit , corresponding to a collapse criterion, is converted into stars. The SFR per free-fall time ǫ ff is obtained by integrating the PDF weighted by the density and the inverse free-fall time above s crit (Hennebelle & Chabrier 2011;Federrath & Klessen 2012): We compute s crit using the model of Krumholz & McKee (2005)for which s crit is estimated by requiring the virial parameter of the gas to be less than 1. We slightly modify the original model to account for subsonic and transonic cases such that for M < 1 we require that the entire computational cell becomes gravitational unstable. With this, we can use the modified model for both regime M ≤ 1andM ≥ 1. The obtained critical density for star formation is where is the local virial number which can be interpreted as an estimator for the local stability, x is the cell size and σ is the turbulent 1D velocity dispersion which is computed by the subgrid scale (SGS) turbulent energy model (see Schmidt, Niemeyer & Hillebrandt 2006;Semenov et al. 2016;Kretschmer & Teyssier 2020, for more details). As a consequence, the efficiency is not a constant anymore but a function of the state of the gas in a computational cell which is characterized through the local virial parameter α vir and the turbulent Mach number M. This model therefore has two star formation channels, either α vir < 1 where the whole cell is collapsing under gravity or if M is large, such that large density fluctuations caused by supersonic turbulence occur. Stellar feedback is implemented as both thermal and momentum injection from SNe as described in Agertz et al. (2013)a n d Kretschmer & Teyssier (2020). A supernova remnant evolves from the energy conserving (adiabatic) phase to the momentum conserving phase. The transition happens when radiation losses become important, or when the shock radius reaches the so-called cooling radius R cool . Resolving this scale by at least three computational cells is crucial in order for the adiabatic phase to inject the correct amount of momentum into the surrounding gas (Kim & Ostriker 2015). If the cooling radius is unresolved, which can occurs at high gas densities, the thermal energy of the supernova is spuriously radiated away. In this case, in addition to the thermal energy injection, we also inject the predicted amount of terminal momentum (P SN ), as computed by Martizzi et al. (2015).
In this work, we use for the cooling radius and the terminal momentum the values computed by Martizzi et al. (2015)f o ra n inhomogeneous medium (as opposed to the values for an homogeneous medium). Specifically, we adopt where n H is the gas density of the cell where the star is exploding, and The injected thermal energy is set to E SN = 10 51 erg. We resolve individual SNe explosions in time, with each star particle triggering multiple SNe explosion, spread randomly between 3 and 20 Myr after the birth of the star particle according to the age distribution of massive stars (Kretschmer & Teyssier 2020). Each supernova injects a mass of metals set to 10 per cent of the progenitor mass. The metallicity is advected by the fluid as a passive scalar. Gas cooling and heating are implemented where equilibrium chemistry for hydrogen and helium is assumed (Katz, Weinberg & Hernquist 1996) as well as metal cooling at both low and high temperature. Additionally, heating by the UV background and its self-shielding is taken into account where the self-shielding density is assumed to be n H = 10 −2 Hcc −1 (Aubert & Teyssier 2010). Heating from the cosmic UV background is turned on at z reion = 10 (Haardt & Madau 1996).
Subgrid physics always remains highly tentative but we emphasize that the robustness of our models have been demonstrated in Kretschmer & Teyssier (2020) as well as in other independent studies (Semenov et al. 2016;Hopkins et al. 2018a;Gensior, Kruijssen & Keller 2020;N uñez-Castiñeyra et al. 2020). Even more important, gas accretion is a more robust feature of our model than the details of our subgrid physics. The behaviour of the CGM gas and its AM, is less dependent on the adopted subgrid physics. This is because gas at larger scales is more predictable than star formation and its feedback inside galaxies (apart from the CGM-wind interaction).

RESULTS
We now present our results for our two selected haloes. Recall that both haloes have a similar accretion history, with an early assembly and a quiet late evolution. The final total masses at z = 0( M vir = 1.1 × 10 12 M ⊙ and M vir = 1.4 × 10 12 M ⊙ ) are both similar to the estimated halo mass of the Milky Way (e.g. Callingham et al. 2019; Figure 1. Projected maps of the two galaxies at z = 0. Left two columns are the disc-dominated galaxy face-on and edge-on. Right two columns are the bulge-dominated galaxy face-on and edge-on. The top row shows the gas surface density and the bottom row shows true colour images with dust absorption allowing us to render both stars and gas. Cautun et al. 2020;W angetal.2020). Despite the fact that the two haloes also have identical halo spins, the two emerging galaxies have very different properties. From Fig. 1, it is immediately apparent that at z = 0 one simulation produced a large, extended thin gas disc. In contrast, the other galaxy contains a small nuclear gas disc within a spherically distributed stellar system. We classify the first galaxy as disc-dominated and the second as bulge-dominated. This visual classification will be confirmed later using a more detailed kinematic analysis (Section 3.4).

Star formation history and gas disc evolution
It is today established that strong feedback is required, especially at high redshift, to obtain a realistic amount of stars in presentday galaxies (Naab & Ostriker 2017). In our simulations, the early evolution is dominated by strong outflows driven by intense starbursts. Fig. 2 shows the evolution of the instantaneous SFR within 10 per cent R vir for our two galaxies. They both experience a sequence of starbursts with SFR ≃ 10 M ⊙ yr −1 until redshift 1.5. The dashed lines indicate the moment of the last starbursts in each galaxy (defined as the last peak in the SFR) which is accompanied by vigorous outflows with mass-loss ratesṀ out ∼ 15-25 M ⊙ yr −1 . It is after this event that the epoch of disc formation truly begins. Fig. 2 also shows the time evolution of the effective radius r e (halfmass radius) of the cold (T < 5 × 10 4 K) gas disc inside 0.1R vir .After the last starburst, both galaxies rapidly develop extended discs (see also Agertz et al. 2020;Dekel et al. 2020;Renaud et al. 2020a, b). More precisely, both galaxies feature r e ≃ 1 kpc before the last starburst, and end up with r e ∼ 5 kpc after only a few Gyr. In fact, the half-mass radius doubles in only 750 Myr. We call this event of spectacularly fast disc growth the Grand Twirl. In the case of the disc-dominated galaxy, r e continues to grow, reaching 10 kpc at z = 0, while the bulge-dominated galaxy shrinks back to r e ∼ 2 kpc.  We also show the evolution of the effective radius of the stars r e⋆ in 0.1R vir . For the disc-dominated system, r e⋆ growth after the Grant Twirl significantly from ∼0.5 to ∼3 kpc at z = 0 similar to the value of the Milky Way of 4-5 kpc (McMillan 2011;Cautun et al. 2020).
Interestingly, during this Grand Twirl, the instantaneous SFR drops significantly in both galaxies. This is due to quenching by gas depletion and outflows following the merger driven starburst (e.g. Renaud et al. 2014) which decreases the gas surface densities across the entire (rapidly expanding) discs. After this short period, gas re-accretion leads to the SFR increasing to 3-5 M ⊙ yr −1 in both galaxies, albeit with no more strong outflows in the disc-dominated case. Indeed, for this Milky Way-like galaxy, the Grand Twirl marks the transition between the high-redshift, outflow dominated regime and the low-redshift, quiescent, disc-dominated regime. At z = 0, the SFR of the disc-dominated galaxy is SFR = 1.3M ⊙ yr −1 which is in good agreement with values obtained for the MW (see e.g. Licquia & Newman 2015, where SFR ∼ 1.7M ⊙ yr −1 ). We note that the star formation history obtained from the simulated galaxies is to some extend determined by the chosen subgrid models. However, possible deviations from observed properties are not detrimental to our results since the evolution of the angular momentum of the gas is mostly dependent on physical processes outside the disc at larger scales. In the following sections we will outline why the galaxies diverge in terms of their overall properties and why outflows shut down in the disc-dominated galaxy.

Stellar disc evolution
In order to outline the differences between the two galaxies, we now analyse the evolution of the stellar discs. For this, we compute the specific angular momentum of all stars within 0.1R vir and compare its time evolution to observations from Fall & Romanowsky (2013) in Fig. 3. The evolution of the angular momenta for the two galaxies is very similar, with j ⋆ ≃ 100−200 km s −1 kpc up to the time of the last starburst, which is marked as a small outlined circle in the figure. During this early epoch, j ⋆ increases roughly in proportion to the stellar mass of the system.
After the last starburst and the Grand Twirl, on the other hand, the two evolutionary tracks diverge: For the disc-dominated galaxy j ⋆ increases by more than one order of magnitude and reaches the regime of observed spiral galaxies (for similar examples, see Agertz & Kravtsov 2016;Sokołowska et al. 2017). The specific angular momentum here grows as j ⋆ ∝ M 2 ⋆ (see also Pedrosa, Tissera & De Rossi 2014), which is steeper than the expected relation for dark matter haloes (j ∝ M 2/3 vir , Fall & Efstathiou 1980). We see that the specific angular momentum increased by a factor of ∼9 and the total stellar AM increases after the Grand Twirl at z ∼ 2 by a remarkable factor of ∼24 from 7.9 × 10 11 to 1.9 × 10 13 M ⊙ km s −1 kpc. The growth by such a remarkably large factor is in good agreement with recent studies (Swinbank et al. 2017;Marasco et al. 2019;Peng & Renzini 2020;Renzini 2020). In contrast, the specific angular momentum of the bulge-dominated galaxy decreases by one order of magnitude and ends up in the observational regime of elliptical galaxies. This classification in the j ⋆ versus M ⋆ plane is in clear agreement with the visual classification in Fig. 1.

Cold gas accretion and angular momentum build-up
To understand the origin of the diverging evolutionary tracks in the j ⋆ /M ⋆ plane, we next turn to the angular momentum of the cold gas (T < 5 × 10 4 K) within the virial radius, but outside 0.1R vir . Fig. 4 shows face-on maps of the mass-weighted distribution of J z at different epochs. The reference frame of the galaxy is computed using the total angular momentum of the cold gas inside 0.1R vir .
At high redshift (z ≃ 2.5), both galaxies show significant amounts of cold gas both counter (red)-and co-rotating (blue). However, going forward in time, it is apparent that for the disc-dominated galaxy, most of the cold gas is co-rotating. In particular, we see co-rotating gas streams, originating from the cosmic web, directly merging with the galactic disc. We also see the size of the galactic disc steadily increasing with time, with all of the disc's gas settled into co-rotation.
In contrast, the halo of the bulge-dominated galaxy is filled with mostly counter-rotating gas originating at large radii and merging with the disc. Quite strikingly, the Grand Twirl in this scenario, while producing an outer extended disc, is entirely counter-rotating in relation to the inner disc (see panels for z ≤ 1.6). This negative angular momentum gas is slowly diminishing the net angular momentum of the central galactic disc, until there is only a tiny cold gas disc left with a counter-rotating outer ring.
To quantify the above points, the distribution of the cold halo gas angular momentum is shown below each image in the bottom rows of Fig. 4 in units of j vir = V vir R vir . At early times (z ≃ 2.5), the distribution is dominated by almost radially in-falling gas. For the disc-dominated galaxy, after the Grand Twirl, the AM distribution shows a clear peak around j z /j vir ≃+ 0.1. For reference, recall that the spin parameter of the halo is more than five times smaller than this (see also Stewart et al. 2017, where a similar factor was found in a study of independent hydrodynamical simulations). This peak persists and signals a consistent feeding of the galactic disc with fresh gas and a constructive build-up of angular momentum. For the bulgedominated galaxy, however, we find a peak in the AM distribution around j z /j vir ≃− 0.2. This stream of gas contributes destructively Figure 4. Evolution of the mass-weighted cold gas angular momentum (AM) in the virial radius of the disc-dominated galaxy (a) and the bulge-dominated galaxy (b). Plots are done at z ∼ 2.5, 2.0, 1.5, 1.0, and 0.5. Top rows: Mass-weighted map of cold gas angular momentum inside the virial radius. Colours indicate the sign of the angular momentum where red corresponds to counter-rotating gas and blue to co-rotating. For the disc-dominated galaxy, most of the cold gas is joining the disc constructively. An extended disc of cold gas is formed. In contrast, for the bulge-dominated galaxy, most of the cold gas is counter-rotating with respect to the disc. This causes the galactic disc to shrink. Bottom rows: Mass-weighted distribution of specific angular momentum of infalling cold gas measured between 0.1 < r/R vir < 1.0 normalized by the specific virial AM of a given snapshot. j z is measured with respect to the disc AM inside 0.1R vir .Forthe disc-dominated galaxy, most infalling cold gas has positive j z /j vir resulting in constrictive AM. For the bulge-dominated galaxy, the distribution shows positive and negative contributions where sometimes the destructive part is dominating. Infalling cold gas therefore adds destructive AM to the disc.
to the AM of the central galaxy disc, explaining why the disc size shrinks at late time.
In summary, the disc morphology at z = 0 comes from cold gas being accreted constructively relative to the initial disc AM, while the bulge-dominated morphology (with a small surviving nuclear disc) originates from late accretion with a destructive impact on the disc AM. This gas is a mixture of pristine gas accreted from the cosmic web and metal-enriched gas recycled from previous galactic outflows.

Kinematics analysis of the stellar discs
We next analyse the angular momentum distributions of the stars in the two galaxies. For this we use the circularity parameter ǫ,defined as ǫ = j z /j c (E), where j z is the z-component of the specific angular momentum and j c (E) is the specific angular momentum for a pure circular orbit at the energy E of the corresponding star particle. For a prograde circular orbit in the mid-plane of the disc, ǫ =+1andǫ = −1 for a retrograde circular orbit. For inclined or eccentric orbits, ǫ will be −1 <ǫ<1. For a classical bulge, for example, we expect a symmetric distribution around ǫ = 0. Fig. 5 shows the evolution of the circularity distribution at different epochs, uniformly spread between the last starburst and the present day z = 0, in steps of ≃ 1.5 Gyr. The disc-dominated galaxy right after the last starburst shows a clear bulge component centred on ǫ = 0 and a small disc component.
Following the (ad hoc) convention described in the literature (e.g. Abadi et al. 2003b;El-Badry et al. 2018;Obrejaetal.2018;Parketal. 2019), we use a fixed cut-off at ǫ = 0.5 (shown as a vertical dotted line in Fig. 5) to separate bulge stars from disc stars. The time evolution in Fig. 5 reveals that the stellar disc component in the disc-dominated system steadily grows between the last starburst and the final epoch, with a three-fold increase in disc mass at an almost constant bulge mass. This level of disc growth, at constant bulge mass, is in excellent agreement with expected L ⋆ galaxy evolution since z ∼ 1−1.5 (van Dokkum et al. 2013). Furthermore, this confirms recent studies which highlighted that at high redshift, galaxies were bulge-dominated and only in the last ∼10 Gyr started to form strong disc components (Park et al. 2019;Bucketal.2020). During this last epoch from z = 2t oz = 0, the stellar mass of the system grew from 1.2 × 10 10 to 3.4 × 10 10 M ⊙ by a factor of 2.7, the radius from 0.5 to 3.1 kpc by a factor of 5.8 and the rotation velocity from 94 to 144 km s −1 by a factor 1.5. As a result the total stellar AM increased by a factor of ∼24 (Swinbank et al. 2017;Peng & Renzini 2020;Renzini2020). For the bulge-dominated galaxy, on the other hand, the disc component centred at ǫ =+ 1 does not grow significantly between the last starburst and the present epoch, while the bulge component sees its mass grow by a factor of 2. Furthermore, we measure a spectacular rise of a counter-rotating stellar disc centred at ǫ =−1, whose origin is directly related to the gaseous filament with negative contribution to the galaxy's angular momentum content identified in Fig. 4;new stars form in the counter-rotating disc that emerged from accretion of counter-rotating gas.
Another spectacular feature in Fig. 5 can be found for the last considered instance, where a dramatic evolution of the circularity distribution in the last 1.5 Gyr occurs due to a minor merger. Indeed, a satellite galaxy with stellar mass M ⋆ ≃ 10 9 M ⊙ on a nearly radial orbit hits the centre of our bulge-dominated galaxy, leading to a strong orbital reconfiguration of the disc stars into more eccentric orbits, ultimately increasing significantly the bulge mass (Naab 2013;Naab & Ostriker 2017).
In summary, the disc-dominated galaxy emerges as a large disc because of the consistent (over several Gyr) and smooth accretion of large angular momentum material. Using the above mentioned discbulge decomposition, the post-starburst disc-to-total ratio is D/T = 0.37 and evolves to D/T = 0.61 at z = 0. In contrast the bulgedominated galaxy suffers from multiple and combined destructive events, leading to the formation of a counter-rotating disc and a massive bulge. For this galaxy the post-starburst D/T = 0.58, but evolves to D/T = 0.19 at z = 0. We reiterate that the z = 0d a r k matter halo spin parameters are almost identical in the two galaxies, illustrating the fact that λ is a poor predictor of the final galaxy morphology (e.g. Jiang et al. 2019).
Note that our bulge-dominated galaxy, which features a gas-rich, star-forming nuclear disc, would be an ideal environment to grow a supermassive black hole and trigger strong AGN feedback. This would change the conditions in the galaxy and probably results in a completely dead, gas-poor spheroidal galaxy.

Evolution of the star formation intensity
It is still unclear why the disc-dominated galaxy, after an early phase of strong starbursts and outflows, managed to form such a quiescent, thin disc of gas and stars that survived until today. Other studies, based on different star formation and feedback implementations, have demonstrated that strong feedback, necessary at high redshift, leads to thick and disturbed discs at low redshift (Agertz et al. 2011;Stinson et al. 2013b;Roškar et al. 2014). Moreover, inspecting the gas density maps of Fig. 1 shows that the bulge-dominated galaxy has a significantly denser circumgalactic medium than the disc-dominated one, by almost a factor of 10.
To interpret these results, we propose the following scenario, based on the time evolution of the surface density of SFR, also referred to as the star formation intensity (SFI), defined as where r SFR/2 is the half-SFR radius. In the last row of Fig. 2,weshow the evolution of SFR for both galaxies as a function of redshift. We find that at high redshift both galaxies have high SFI with SFR ≃ 10 M ⊙ yr −1 kpc −2 maintained until z ≃ 1.5. After the Grand Twirl, however, the SFI in the disc-dominated galaxy drops very rapidly to a value smaller than 0.1M ⊙ yr −1 kpc −2 whereas the bulge-dominated galaxy maintains a high SFI until z = 0. The significant drop of the SFI in the disc-dominated case is naturally explained by (1) starburst quenching of star formation (reduced SFR), due to rapid gas depletion and outflows, and (2) the rapid size evolution and emergence of an extended gas disc (increased r SFR/2 ).
Since our disc-dominated galaxy's key properties (r e⋆ , SFR, M ⋆ and SFR ) as well as the resulting disc-to-total ratio at z = 0a r e similar to those obtained for the MW (e.g. Cautun et al. 2020, found D/T ∼ 0.6), we are interested in comparing the evolution of the SFI of our simulations with observational data for the MW. Although our disc-dominated galaxy might differ in certain quantities from the MW, such a comparison can shed light on the overall evolution of disc galaxies with respect to a possible Grand Twirl event. We compare our results to the SFI of the Milky Way using the analysis performed by Lehnert et al. (2014) based on stars in the solar neighbourhood (Haywood et al. 2013;Snaith et al. 2014), shown as a solid black line in Fig. 2. The agreement is striking, suggesting that a similar Grand Twirl phenomenon may have occurred in the Milky Way, around z ≃ 1.5. Also shown in the figure is a range of characteristic values of SFR that have been proposed as critical thresholds for driving large-scale outflows. This value, may range from 0.1M ⊙ yr −1 kpc −2 based on observations of nearby galaxies (Lehnert & Heckman 1996;Heckman 2003)to∼ 1M ⊙ yr −1 kpc −2 at high redshift (Newman et al. 2012;F orster Schreiber et al. 2019). Our simulations support the same idea: the disc-dominated galaxy shows very little outflows after the last starburst, while the bulge-dominated galaxy continues to eject significant quantities of gas in the halo. This explains the very different gas densities in the two galactic corona see in Fig. 1. We will discuss this point further in the next section.
The radial profiles of SFR , shown in Fig. 6, summarize and confirm this picture. First, we find the rapid growth of the gas disc after the Grand Twirl of the disc-dominated case, while the bulge-dominated galaxy maintains a compact star formation region. We note that the central gas-rich, star-forming nuclear disc of the bulge-dominated galaxy should probably host a supermassive black hole and trigger strong AGN feedback which are not included in our subgrid physics models. Strong AGN feedback would probably turn the galaxy into a completely quenched and gas-poor spheroidal galaxy. However, the main results of this paper concerning the importance of co-rotating accreted gas to form extended disc galaxies should not be affected by the lack of such a model. Secondly, we see that the disc-dominated galaxy SFI drops below the critical threshold for starburst-driven outflows everywhere at late times, even in the central region, while the bulge-dominated galaxy retains regions above the threshold, leading to a consistent launching of outflows, even at late times.

DISCUSSION
We now discuss in more details the rapid growth of the large disc we have observed in our simulations, an episode that we called the Grand Twirl. We speculate this also happened in the Milky Way, as suggested by the analysis of Lehnert et al. (2014). We also discuss how this rapid disc growth sets the conditions for the emergence of a quiescent disc, in which SNe-driven feedback becomes inefficient at launching outflows and at disturbing the thin disc. We note in passing that our analysis in based only on two simulations and cannot be used to draw any statistical conclusions on the frequency of large discs predicted by current galaxy formation simulations (for recent statistical studies, see e.g. Dubois et al. 2014;Genel et al. 2015;Zavala et al. 2016;DeFelippis et al. 2017;Lagos et al. 2017;Zjupa & Springel 2017;Parketal.2019).

The emergence of large discs
The traditional picture of disc assembly is based on tidal torque theory, for which haloes obtain their spin from their cosmic environment through tidal torques (Peebles 1969;Doroshkevich 1970;White 1984). Assuming that gas follows the dark matter dynamics on large scales, emerging star-forming discs are expected to have similar spins (Fall & Efstathiou 1980). However, simulations with gas cooling and star formation have demonstrated that this simple picture is invalid. Several studies have shown that the spin of baryons and dark matter inside the halo are barely correlated and that the two angular momentum distributions widely differ (van den Bosch et al. 2002;Genel et al. 2015;Teklu et al. 2015;Stewart et al. 2017;Jiang et al. 2019).
In both our simulations, we observe at high redshift (z ≃ 1.5−2) a spectacularly fast disc growth triggered by a combination of inflowing pristine gas through cold streams and recycled metal enriched outflows that combined and quickly fuelled the gas disc in the rapid event we called the Grand Twirl. This is in line with previous studies that highlighted the role played by cold stream accretion from the cosmic web of large-scale filaments in setting the spin of gaseous discs (Pichon et al. 2011;Dubois et al. 2014). Moreover, Stewart et al. (2013) showed that the high angular momentum of gaseous disc observed in simulations is a direct consequence of this filamentary accretion. This cold stream accretion is the main component of the new theory of disc evolution proposed by Danovich et al. (2015).
The second important ingredient for the fast build-up of a large disc is the re-accretion of outflow material. Indeed, feedback at these high redshifts is strong enough to remove low angular momentum gas from the disc and make it available in the halo to be efficiently torqued to large angular momentum by incoming cold streams. This combination of pristine and recycled material with large angular momentum is then accreted quickly into this Grand Twirl. The same process was described inÜbler et al. (2014), who argued that strong stellar feedback and re-accreted material can actually promote the formation of discs.
Comparing the later evolution of the discs in our two different haloes, it appears crucial for the disc survival to have a sustained accretion of positive angular momentum gas that is added constructively to the disc. For the disc-dominated galaxy, the incoming cold gas from the cosmic web is almost exclusively co-rotating leading to the steady build-up of an extended gas disc. On the other hand, for the bulge-dominated galaxy, the incoming material is almost purely counter-rotating leading to the steady build-up of a large bulge. The gas-rich and star-forming nuclear disc would be an ideal environment to grow a supermassive black hole and trigger strong AGN feedback.
The galaxy would probably end up as a passive bulge-dominated galaxy similar to those common at low redshift.
This picture is in agreement with previous studies that showed that the final galaxy morphology depends on the coherent alignment (or lack thereof) of the accreted angular momentum. In particular, Sales et al. (2012) compared the alignment of the angular momentum of baryons within a spherical shell measured at the turn-around radius with the total halo spin. The more consistently aligned the angular momentum in these shells is, the more likely a large disc will emerge at z = 0. Our simulations confirm this interpretation but instead of considering the turn-around radius and the entire halo, we used instead the alignment of angular momentum within R vir with the one of the central galaxy within 0.1R vir . This allows us to remove from the analysis satellites and gas streams within the haloes that never reach the central region.

The emergence of quiescent discs
A spectacular consequence of the Grand Twirl is the emergence of a very quiescent disc, in striking contrast with the earlier epoch that was dominated by strong starbursts and associated outflows, resulting in very disturbed environments. After the rapid formation of the extended gas disc, the star formation surface density drops abruptly, in agreement with the model of Lehnert et al. (2014)f o r the Milky Way. This is at odd with some earlier works for which strong feedback completely destroyed the final thin discs (Scannapieco et al. 2012;Stinson et al. 2013a;R ǒ skar et al. 2014), outlining an apparent contradiction between the strong feedback required at high redshift to regulate star formation and the weak feedback required at low redshift to allow the formation of quiescent and thin gas discs. One issue in these earlier studies is the unrealistic modelling of stellar feedback. Indeed, instead of resolving successive individual SNe (like we do here), these studies used individual star particle to represent entire star clusters and made them explode in one single energetic event.
In order to launch strong outflows, the SFI has to be high enough, so that overlapping SN explosions create a volume-filling hot phase such that gas is removed in a thermally driven outflow. In the Milky Way and at the present epoch, the SFI is too low and disfavour the build-up of overlapping hot bubbles and the formation of strong outflows (Gatto et al. 2015;L ie ta l .2015;Naab & Ostriker 2017;Steinwandel et al. 2020). This has been demonstrated for example byKim&Ostriker (2015), Walch et al. (2015), and Girichidis et al. (2016) using high-resolution simulations of stratified galactic discs with multiple SN explosions. These studies have also shown that the degree of clustering of SN explosions matters, with overlapping and percolating SNe driving much stronger outflows than individual, spatially separated events (Kim, Ostriker & Raileanu 2017;Fielding, Quataert & Martizzi 2018;Gentryetal.2019;Martizzi 2020). The Grand Twirl is the main mechanism causing this transition from highly clustered SN explosions driving strong outflows to a more uniform distribution of isolated events unable to perturb the disc anymore.
We believe these effects have been seen already in several past cosmological simulations of galaxy formation. For example, Muratov et al. (2015) have computed the mass loading factor in several zoomin simulations of disc galaxies, showing a striking and abrupt decline of the outflow rate at late time and for Milky Way sized haloes. Although this effect was attributed to the increased escape velocity in larger mass haloes, we speculate that the rapid formation of a large disc must have also played an important role. Agertz & Kravtsov (2015) modelled the formation of a Milky Way analogue and obtained also a rapid transition of the SFR, from very bursty at high redshift to very smooth at low redshift. The transition into an angular momentum rich disc occurred after the last major merger at z>1, and was likely due to a similar Grand Twirl episode.

CONCLUSIONS
In this paper, we have investigated the evolution of two galaxies using hydrodynamical simulations of two haloes that have masses similar to the estimated Milky Way halo mass at z = 0. We have seen the emergence of two very different galaxies albeit their very similar halo properties. Our main results and conclusions can be summarized as follows: (i) Our simulated galaxies undergo an epoch of rapid disc growth at z ≃ 1.5, where the gas half-mass radius doubles in less than a Gyr. We called this event the Grand Twirl. During this epoch, both pristine gas from cold streams and re-accreted outflow material combine to form an extended gas disc in a very short time.
(ii) For one galaxy, an extended and thin disc forms and persists until z = 0, as the result of the consistent accretion of cold gas with positive angular momentum, adding constructively material to the initial disc that formed after the Grand Twirl episode. For the other galaxy, accreted gas with negative angular momentum favours the growth of the bulge component.
(iii) The build-up of this extended gas disc explains the transition from a high-redshift, outflow-dominated era to a very quiescent epoch at low redshift, for which feedback is not capable of launching outflows anymore.
For our disc-dominated galaxy, the time evolution of the surface star formation density shows remarkable agreement with data obtained for the Milky Way. We speculate that the Milky Way may have undergone such a Grand Twirl episode at z ≃ 1.5. This echoes with recent work speculating a similar massive accretion event at the origin of a peculiar structure in the metal enrichment properties of stars in the Milky Way as seen by Gaia (Belokurov et al. 2020). Star formation in such a rapidly formed outer disc could favour the emergence of a characteristic chemically bimodality in stellar abundance ratios between α-elements and iron ([α/Fe]) over a wide range of [Fe/H]. This is the focus of a follow-up project (Agertz et al. 2020;Renaud et al. 2020a, b). In future work, we would also like to perform the same analysis presented here over a larger sample of haloes. This will allow us to determine whether such Grand Twirl episodes are a common feature in the past history of large galactic haloes.