Galactic discs mass transport.

Multiscale mass transport in z6 galactic discs: fueling black holes.


By using AMR cosmological hydrodynamic N-body zoom-in simulations, with the RAMSES code, we studied the mass transport processes onto galactic nuclei from high redshift up to . Due to the large dynamical range of the simulations we were able to study the mass accretion process on scales from to few . We studied the BH growth on to the galactic center in relation with the mass transport processes associated to both the Reynolds stress and the gravitational stress on the disc. Such methodology allowed us to identify the main mass transport process as a function of the scales of the problem. We found that in simulations that include radiative cooling and SNe feedback, the SMBH grows at the Eddington limit for some periods of time presenting throughout its evolution. The parameter is dominated by the Reynolds term, , with . The gravitational part of the parameter, , has an increasing trend toward the galactic center at higher redshifts, with values at radii few contributing to the BH fueling. In terms of torques, we also found that gravity has an increasing contribution toward the galactic center at earlier epochs with a mixed contribution above . This complementary work between pressure gradients and gravitational potential gradients allows an efficient mass transport on the disc with average mass accretion rates of the order few . These level of SMBH accretion rates found in our cosmological simulations are needed in all models of SMBH growth that attempt to explain the formation of redshift quasars.

galaxies: formation — large-scale structure of the universe — stars: formation — turbulence.

1 Introduction

The mass transport (MT) process in astrophysical environments has relevance for different phenomena in our Universe. It is important for planet formation in proto-planetary discs, it triggers the AGN activity associated to super massive black hole (SMBH) accretion at high redshift and it is responsible for mass accretion from the filamentary structures around dark matter (DM) haloes into the central regions of the first galaxies. A full understanding of this phenomenon is very important in the construction of a galaxy formation theory.

In particular, the MT phenomenon has a crucial relevance in models of black hole (BH) formation and their growth in the early stages of our Universe. The observation of very bright quasars at redshift with luminosities implies the existence of BHs with masses of the order when our Universe was about Gyr old (Fan et al., 2001), i.e. SMBH should be formed very early in the history of our Universe and they should grow very fast in order to reach such high masses in the first of our Universe. To understand such a rapid early evolution is one of the main challenges of current galaxy formation theories. For a more extended discussion on massive BH formation at high redshift see Volonteri (2010) and Haiman (2013).

There are three main scenarios for the formation of SMBH seeds:

  • Seeds from the first generation of stars and their subsequent accretion to form the SMBH: In this scenario a top heavy initial mass function (IMF) associated with population III (pop III) stars (Abel et al., 2002; Bromm & Larson, 2004) can leave BHs with masses of the order of M (Heger & Woosley, 2002). We note that more recent studies have shown that pop III stars could be formed in clusters of low mass stars, e.g. Stacy et al. (2010); Greif et al. (2011); Clark et al. (2011).

  • Massive seeds formed by the direct collapse of warm ( K) neutral hydrogen inside atomic cooling haloes: In this scenario the initial BH seeds are formed by the direct collapse of warm gas triggered by dynamical instabilities inside DM haloes of mass at high redshift, (e.g. Oh & Haiman, 2002; Lodato & Natarajan, 2006; Begelman et al., 2006). If there is no molecular hydrogen to avoid fragmentation (e.g. Agarwal et al., 2012; Latif et al., 2013, 2014) and there is an efficient outward AM transport (Choi et al., 2015) such a scenario could favor the formation of a BH seed of , e.g. Begelman et al. (2008).

  • BH seeds formed by dynamical effects of dense stellar systems like star clusters or galactic nuclei, e.g. Schneider et al. (2006): In this scenario the formed BH seed can be of mass (Devecchi & Volonteri, 2009)

Studies related to the initial mass for SMBH formation favor massive M seeds inside primordial atomic cooling haloes (Volonteri et al., 2008; Tanaka & Haiman, 2009; Lodato & Natarajan, 2006). Despite that, due to the lack of observational evidence it is not clear yet if one of these scenarios is preferred by nature or all of them are working at the same time in different haloes.

There is dynamical evidence for the existence of SMBH in the center of nearby galaxies (Ferrarese & Ford, 2005) with masses in the range suggesting that the BHs formed in the first evolutionary stages of our Universe are now living in the galactic centers around us, including our galaxy (Ghez et al., 2005). Besides their ubiquitous nature there is evidence of scaling relations connecting the BH mass with its host galaxy properties, namely the galactic bulge - BH mass relation (e.g. Gültekin et al., 2009) and the bulge stars velocity dispersion - BH mass relation (e.g. Tremaine et al., 2002; Ferrarese & Merritt, 2000). Such relations suggest a co-evolution between the BH and its host galaxy.

In a cosmological context, motivated by the theoretical study of Pichon & Bernardeau (1999) a series of recent simulations Pichon et al. (2011) and Codis et al. (2012) have argued that the galactic spin may be generated entirely in the baryonic component due to the growth of eddies in the turbulence field generated by large-scale ( few Mpc) mass in-fall relating the large scale angular momentum (AM) acquisition with mass transport phenomena inside the virial radius.

Danovich et al. (2015) studied the AM acquisition process in galaxies at redshift identifying 4 phases for AM acquisition: i) the initial spin is acquired following the Tidal Torque Theory (Peebles, 1969; Doroshkevich, 1970). ii) In this phase both the mass and AM are transported to the halo outer region in the virialization process following the filamentary structure around them. iii) In the disc vicinity the gas forms a non-uniform ring which partially suffers the effect of the galactic disc torques producing an alignment with the inner disc AM. iv) Finally, outflows reduce gas with low specific AM, increasing its global value at the central region and violent dynamical instabilities (VDI) associated to clump-clump interactions and clump-merged DM halo interactions remove AM, allowing a more centrally concentrated gas.

At smaller scales ( few ), Prieto et al. (2015) studied the MT and AM acquisition process in four DM haloes of similar mass and very different spin parameter (Bullock et al., 2001a) at redshift . The main result of this work is the anti-correlation between the DM halo spin parameter and the number of filaments converging on it: the larger the number of filaments the lower the spin parameter. Such a result suggests that DM haloes associated to isolated knots of the cosmic web could favor the formation of SMBH because the inflowing material would have to cross a lower centrifugal barrier to reach the central galactic region.

In a non-cosmological context, Escala (2006) and Escala (2007) has shown that the interplay and competition between BH feeding and SF can naturally explain the relation. Using idealized isolated galaxy evolution simulations Bournaud et al. (2007) showed that galaxies are able to form massive clumps due to gravitational instabilities (Toomre, 1964) triggered by its high gas mass fraction. Such clumpy high redshift galaxies evolve due to VDI to form spiral galaxies characterized by a bulge and an exponential disc. Similar results have been found in cosmological contexts by Mandelker et al. (2014). They show that the formation of massive clumps is a common feature of galaxies. Due to the fast formation and interaction between them the VDI dominate the disc evolution. A similar clump migration has been observed at higher redshift in a DM halo at in Dubois et al. (2012) and Dubois et al. (2013). In these works the migration has been triggered by DM merger induced torques. In contrast to the scenario presented above there are studies supporting the idea that clump interaction in high redshift discs are not the main source to build up the galactic bulges (e.g. Hopkins et al., 2012; Fiacconi et al., 2015; Tamburello et al., 2015; Behrendt et al., 2016; Oklopcic et al., 2016).

Inspired by the parametrization in the seminal paper of Shakura & Sunyaev (1973), Gammie (2001) studied the gravitational stability in cool thin discs. In his work Gammie (2001) quantified the rate of angular momentum flux in terms of the Reynolds and gravitational stress. In this work we will use a similar -formalism to study the MT process on galactic discs at redshift performing N-body and hydrodynamic numerical simulations from cosmological initial conditions. We will study a halo of few , a mass value not studied already in this context. It is the first time that such an approach is being used to study the MT process on galaxies at high redshift. Furthermore we will compute directly the torques working on the simulated structures from scales associated to the cosmic web around the central DM halo to few pc scales associated to the galactic disc. Such an approach will allow us to have clues about the main source of mass transport on these objects and then to have some insights about the SMBH growth mechanisms at high redshift.

The paper is organized as follows. Section §2 contains the numerical details of our simulations. Here we describe the halo selection procedure, our refinement strategy and the gas physics included in our calculations. In section §3 we show our results. Here we present radial profiles of our systems, star formation properties of our galaxies, a gravitational stability analysis and show a mass transport analysis based on both the formalism and the torques analysis on small and large scales. In section §4 we discuss our results and present our main conclusions.

2 Methodology and Numerical Simulation Details

2.1 RAMSES code

The simulations presented in this work were performed with the cosmological N-body hydrodynamical code RAMSES (Teyssier, 2002). This code has been written to study hydrodynamic and cosmological structure formation with high spatial resolution using the Adaptive Mesh Refinement (AMR) technique, with a tree-based data structure. The code solves the Euler equations with a gravitational term in an expanding universe using the second-order Godunov method (Piecewise Linear Method).

2.2 Cosmological parameters

Cosmological initial conditions were generated with the mpgrafic code (Prunet et al., 2008) inside a side box. Cosmological parameters where taken from Planck Collaboration (2013): , , , , and .

2.3 Halo selection

Using the parameters mentioned above, we ran a number of DM-only simulations with particles starting at . We selected one DM halo of mass at redshift z=6. We gave preference to DM haloes without major mergers through its final evolution in order to have a more clean and not perturbed system to analyze.

After the selection process we re-simulated the halo including gas physics. For these simulations we re-centered the box on the DM halo position at redshift . We set a coarse level of (level 7) particles and allowed for further DM particle mass refinements until level 10 inside a variable volume called mask (as we wil explain below). In this way we were able to reach a DM resolution equivalent to a particle grid inside the central region of the box, which corresponds to a particle mass resolution , in other words we resolved the high redshift halo with particles and our final halo with particles.

2.4 Refinement strategy

In order to resolve all the interesting regions we allowed refinements inside the Lagrangian volume associated to a sphere of radius 1 around the selected DM halo at . Such a Lagrangian volume is tracked back in time until the initial redshift of the simulation, =100. In this way we ensure that the simulation is resolving all the interesting volume of matter throughout the experiment, i.e. all the material ending inside the at the end of the simulation2. The Lagrangian volume (the mask) is defined by an additional passive scalar advected by the flow throughout the simulation. At the beginning the passive scalar has a value equal to 1 inside the mask and it is 0 outside. We apply our refinement criteria in regions where this passive scalar is larger than .

In our simulations a cell is refined if it is in regions where the mask passive scalar is larger than and if one of the following conditions is fulfilled:

  • it contains more than 8 DM particles,

  • its baryonic content is 8 times higher than the average in the whole box,

  • the local Jeans length is resolved by less than 4 cells (Truelove et al., 1997), and

  • if the relative pressure variation between cells is larger than 2 (suitable for shocks associated to the virialization process and SNe explosions).

Following those criteria the maximum level of refinement was set at , corresponding to a co-moving maximum spatial resolution of and a proper spatial resolution of at redshift . With this resolution we were able to resolve the inner DM region with computational cells.

2.5 Gas physics

Our simulations include optically thin gas cooling. The gas is able to cool due to H, He and metals following the Sutherland & Dopita (1993) model until it reaches a temperture of . Below this temperature the gas can cool until due to metal lines cooling. We note that the cooling functions assume collisional ionization equilibrium. The metals are modeled as passive scalars advected by the gas flow. In order to mimic the effect of H cooling in primordial environments all our simulations started with an initial metallicity (Powell et al., 2011). Furthermore, a uniform UV background is activated at , following Haardt & Madau (1996).

Star formation

The numerical experiments include a density threshold Schmidt law for star formation: above a given number density, set to in our case, the gas is converted into stars at a rate density, , given by (e.g. Rasera & Teyssier, 2006; Dubois & Teyssier, 2008):


where is the local gas density, is the constant star formation efficiency and is the density dependent local free fall time of the gas. The number density for star formation, , corresponds to the value at which the local Jeans length is resolved by 4 cells with a temperature 3

When a cell reaches the conditions to form stars we create star particles following a Poisson distribution


with the number of formed stars and


where is the cell grid side, is the mass of the stars, is the time step integration and is the star formation time scale. This process ends up with a population of stars inside the corresponding cell. In order to ensure numerical stability we do not allow conversion of more than of the gas into stars inside a cell.

SNe feedback

After 10 Myr the most massive stars explode as SNe. In this process a mass fraction (consistent with a Salpeter initial mass function truncated between and ) of the stellar populations is converted into SNe ejecta:


In this case is not the single stellar particle mass but the total stellar mass created inside a cell. Furthermore, each SNe explosion a releases specific energy into the gas inside a sphere of :


As mentioned above, metals are included as passive scalars after each SNe explosion and then they are advected by the gas flows. This means that after each SNe explosion a metallicity


is included as metals in the gas in the simulation. Such an amount of metals is consistent with the yield of a type II SNe from Woosley & Weaver (1995).

In this work we used the delayed cooling implementation of the SNe feedback (discussed in Teyssier et al., 2013; Dubois et al., 2015). This means that in places where SNe explode, if the gas internal energy is above an energy threshold , the gas cooling is turned off for a time in order to take into account the unresolved chaotic turbulent energy source of the explosions.

As written in Dubois et al. (2015) the non-thermal energy evolution associated to the SNe explosions can be expressed as


In an equilibrium state it is possible to write


If we assume that non thermal energy is associated to a turbulent motion with a velocity dispersion and that this energy will be dissipated in a time scale of order the crossing time scale associated to the local jeans length then , and it is possible to write


Then, expressing the local Jeans length as , with the proper cell side at the highest level of refinement, it is possible to write the dissipation time scale as:

Given our parameters, we set the non-thermal energy dissipation time scale as .

In this model the gas cooling is switched off when the non-thermal velocity dispersion is higher than a given threshold:

which for us is .

2.6 Sink particles and black hole accretion

In order to follow the evolution of a black hole (BH) in the simulations we introduced a sink particle (Bleuler & Teyssier, 2014) at the gas density peak inside a DM halo of at redshift . The BH seed mass is 10, roughly following the relation of McConell et al. (2011). Such a black hole mass is in the range of masses associated to direct collapse of warm gas inside atomic cooling haloes at high redshift (e.g. Oh & Haiman, 2002; Lodato & Natarajan, 2006; Begelman et al., 2006, 2008; Agarwal et al., 2012; Latif et al., 2013, 2014; Choi et al., 2015). We did not allow more BH formation after the formation of the first one. In order to compute the mass accretion rate onto the BH we use the modified Bondi-Hoyle accretion rate described below.

Modified Bondi accretion

Bleuler & Teyssier (2014) implemented the modified Bondi mass accretion rate onto sink particles in the RAMSES code. They use the expression presented in Krumholz et al. (2004) based on the Bondi, Hoyle and Lyttleton theory (Bondi & Lyttleton, 1939; Bondi, 1952). There the Bondi radius


defines the sphere of influence of the central massive object of mass and its corresponding accretion rate is given by


where is the gravitational constant, is the average sound speed and is the average gas velocity relative to the sink velocity, is an equation of state dependent variable and it is in the isothermal case. The density is the gas density far from the central mass and it is given by


where is the solution for the density profile in the Bondi model (Bondi, 1952). The variable is the dimensionless radius and the corresponding density. In our case , with the minimum cell size in the simulation.

The modified Bondi accretion rate is limited by the Eddington accretion rate. In other words, the sink particle can not accrete at a rate larger than the Eddington rate, given by:


where is the proton mass, is the Thomson scattering cross section, is the speed of light and is the fraction of accreted mass converted into energy.

3 Results

In this work we will analyze three simulations:

  • NoSNe simulation: Includes star formation and modified Bondi-Hoyle-Lyttleton (BHL) accretion rate onto sinks,

  • SNe0.5 simulation: Includes star formation, BHL accretion rate onto sinks and SNe feedback with a consistent delayed cooling , and

  • SNe5.0 simulation: Similar to SNe0.5 but with an out of model

We are not including AGN feedback in these experiments. This important ingredient has been left for a future study to be presented in an upcoming publication.

Figure 1 shows a gas number density projection of our systems at redshift at two different scales. The top rows show the large scale () view of the systems and the bottom rows show a zoom-in of the central () region.

In the top panels it is possible to recognize a filamentary structure converging on to the central galaxy position. Such filaments work as pipes channeling cold baryonic matter onto the converging region: the place for galaxy formation, i.e. knots of the cosmic web. Aside from the accreted low density gas it is possible to recognize a number of over-densities associated to small DM haloes merging with the central dominant halo, a common feature of the hierarchical structure formation. Such mini haloes certainly perturb the galactic disc environment as we will see later.

Whereas at large scale we can see that the difference between runs are the small scale features associated to the shocks produced by SNe explosions, the bottom panels show a clear difference between simulations at the end of the experiments. The NoSNe run developed a concentrated gas rich spiral galaxy whereas both SNe runs have less concentrated gas and much more chaotic matter distribution. Such a difference certainly is a consequence of SNe explosions: the energy injected into the environment is able to spread the gas out of the central region and then to decrease the average density due to effect of the expanding SNe bubbles. Such a phenomenon is able to destroy the galactic disc as can be seen in the central and right-bottom panel where the proto-galaxy is reduced to a number of filaments and gas clumps.

Figure 2 shows the rest frame face on (top panels) and edge-on (bottom panels) stellar populations associated to our systems in three combined filters: i, u and v. The images were made using a simplified version of the STARDUST code (Devriendt et al., 1999)4.

3.1 Radial profiles

Before computing any radial average from our AMR 3D data we aligned the gas spin vector with the Cartesian direction. After this procedure, we performed a mass weighted average in the direction in order to have all the interesting physical quantities associated to the disc surface:


with and the grid size.

In order to get our radial quantities we have performed a mass weighted averaging in the cylindrical direction of our disc surface data:


where and are the radial and azimuthal cylindrical coordinates, with and .

Figure 3 shows the gas surface density5 (SD) (top panels), the stars SD (central panel) and the SD star formation rate (bottom panel) for different redshifts as a function of radius. (The line style-redshift relation will be kept for the following plots). All our simulations show a number of peaks in the gas SD profiles at almost all the sampled redshifts, a proof of the irregular and clumpy structure of the gas.

The second row of the figure shows more clear differences between our runs. Apart from an almost monotonic increment on the stars SD for each simulation it is also possible to see lower central peaks from left to right. Such a trend is a consequence of the different strength of the feedback which also increases from left to right.

The bottom panels show the SD star formation rate (SFR). In order to compute this quantity we took into account all the stars with an age formed in the time elapsed between consecutive outputs (which is in the range ). Inside one tenth of the virial radius we compute the height and the radius where of the gas and stars are enclosed. We averaged these scales for stars and gas to define a cylinder where we compute the SFR (a similar procedure will be used in the next sub-section to compute the Kennicutt-Schmidt law). We can see that our NoSNe run shows higher SD SFR peaks above and the SNe5.0 run has the lower SD SFR. Such a fact will be confirmed after computing the global SD against the global SD SFR in the next section.

Figure 4 shows different gas velocities associated to our systems at different redshift as a function of radius. It plots the radial velocity (in solid black line), azimuthal velocity (in long-dashed blue line) and the spherical circular velocity of the disc (in dot-dashed cyan line) defined as:


where is Newton’s constant and is the total (gas, stars and DM) mass inside the radius .

In our simulations the radial velocity fluctuates from negative to positive values. Such a feature is a proof of a non stationary disc where at some radii there is inflowing material whereas at other radii there are gas outflows. Such features can be produced by virialization shocks, DM halo mergers or SNe explosions. It is worth noticing that at kpc, which roughly correspond to the outer edge of the galactic disc, the gas is inflowing in most of the cases. Such a feature is a remarkable signal of radial gas inflows at distance from the center of the system. This fast radial material comes from larger scales channeled by the filamentary structure shown in the top panels of figure 1 and, as mentioned above, they supply the central DM halo region with cold gas at rates as high as as we will see in the following sections.

The orbital velocity tends to be roughly similar to the spherical circular velocity at large radii in most of the cases but in general the circular velocity does not follow the spherical circular orbit. Such deviations can be explained due to the shocked gas inflows, the mergers suffered by the central halo and due to SNe explosions which enhance the pressure support against gravity. We emphasize that these kinds of interactions have a gravitational effect due to tidal forces (mergers and clump-clump interaction) on the disc and also have a hydrodynamical effect (shocks). In our SNe runs it is clear that the spherical circular curves are lower than the NoSNe curve. In other words the enclosed mass inside is lower in the SNe runs. That is because the SNe explosions spread the gas out of the central region. Actually from the shocked gas features at the top right panel of figure 1 it is possible to see that the outflows can reach regions at from the central galaxy, i.e. at .

3.2 Star formation

Despite the lack of observations of the Kennicutt-Schmidt (KS) law at high redshift () it is interesting to compare the KS law from our simulations with its currently accepted functional form from Kennicutt (1998) (hereafter K98). Furthermore, it is also interesting to compare our data with more recent literature from the Daddi et al. (2010) (hereafter D10) results for normal and star burst galaxies.

Figure 5 shows the KS law for our runs. Each point marks the SD SFR as a function of the total gas SD. The SD SFR was computed following a similar procedure as in the previous sub-section.

There is a correlation between the lower gas SD and the level of feedback in our results: the higher the feedback the lower the gas SD, which is a natural consequence of the gas heating due to SNe events. Whereas the NoSNe run shows points covering decade at high SD with a large scatter in SD SFR from below the D10 normal galaxies sequence to above the D10 star burst sequence the SNe runs cover a larger SD range. Both SNe runs are in agreement with the star burst sequence of D10 and the SNe0.5 simulation shows a lower scatter in the points. Such a behavior could be due to the number of mergers suffered by these kind of haloes at high redshift.

Figure 6 shows the stellar mass normalized by , where is the universal baryonic fraction. At the end of the simulation our SNe0.5 galaxy has a stellar metallicity and our SNe5.0 galaxy a metallicity of . It is clear from the figure that the NoSNe run is producing much more stars than our SNe runs and that due to the extreme feedback our SNe5.0 simulation form less stars than our SNe0.5 simulation. When we compare our results with the one shown in Kimm et al. (2015) we can see that our results are in the range of their MFB and MFBm simulations at similar . Despite the uncertainties and the lack of robust observational constrains, such values are not far (a factor of few for SNe5.0 and just in the limit of the order of magnitude for SNe0.5) from the prediction from Behroozi et al. (2013) (hereafter B13) where the stellar to halo mass ratio is of the order of few at the same mass range and high redshift.

Figure 7 shows the SFR for our runs as a function of redshift. When we compare the NoSNe run with our SNe runs the main difference arises in the continuity of the SFR history. Whereas the NoSNe run shows a continuous line the SNe runs present periods of almost zero SFR. Such periods last few and are more frequent in our SNe5.0 run due to the stronger feedback. Despite the large fluctuations in the SFR data our SNe runs tend to be in the range . Such numbers are in line with the one found by Watson et al. (2015) (here after W15) and references therein for high redshift galaxies. Taking into account the uncertainties of the predictions, if we compare our SNe run results with B13 they tend to be below or similar to their data for halo at .

Figure 1: Mass weighted projection of the gas number density for our simulations: NoSNe left column, SNe0.5 central column and SNe5.0 right column. The top row is a large scale ( ckpc square side) view of our systems and the bottom row is a zoom-in of the central region of the system ( ckpc square side). From the top panels it is possible to identify the filamentary structure converging at the central region of the system: the galaxy position. Such filaments channel and feed the galactic structure. At large scales it is possible to recognize shock waves associated to the SNe explosions of our SNe runs. Beside the low density gas, there are a number of over-densities associated to small DM haloes about to merge with the central structure. The bottom panels show a dramatic difference between our simulations: a compact gas rich spiral galaxy for the NoSne experiment, a rough spiral galaxy disturbed by SNe feedback in our SNe0.5 run and a group of clumps in our SNe5.0 simulation.

3.3 Disc stability

High redshift galactic environments have a high gas fraction (e.g. Mannucci et al., 2009; Tacconi et al., 2010). Figure 8 shows the gas fraction of our systems as a function of redshift. Here we define the gas fraction as the ratio between the galactic gas mass and the mass of the gas plus the stars in the galaxy: . All our systems shows a high gas fraction with values fluctuating around . In fact the average values for our simulations are , and below . If we average below redshift we find , and . It is interesting to compare such numbers with those found by W15. In this work the authors describe the properties of a galaxy. The galaxy at this redshift has a gas fraction , in other words our values of are inside the errors associated to their observations as we can see in figure 8 with the SNe runs closer to the observational expectations.

The non-stationary and highly dynamic nature of the high gas fraction systems makes them susceptible to gravitational instabilities. In order to analyze the disc stability throughout its evolution we will use the Toomre parameter, , stability criterion (Toomre, 1964):


A convenient modification of the Toomre parameter to take into account the turbulent velocity dispersion of the fluid has the form . Despite the ad hoc modification of the parameter it is not straightforward to interpret the turbulent velocity dispersion of the gas as a source of pressure counteracting the gravity (Elegreen & Scalo, 2004). This comes from the fact that this pressure term could only be defined in the case where the dominant turbulent scale is much smaller than the region under consideration, which is in fact not the case of the ISM. Rigorous analysis indeed shows that turbulence can be represented as a pressure only if the turbulence is produced at scales smaller than the Jeans length (micro-turbulence in Bonazzola et al., 1992). Therefore the gravitational instability analysis is not strictly applicable with a turbulent pressure term that could stabilize and dampen all the substructure below the unstable scale associated to .

The left column of figure 9 shows the Toomre parameter for our three runs at different redshifts. The gray dashed horizontal line marks the state. For completeness, the right column of figure 9 shows the Toomre parameter associated to the turbulent velocity dispersion. Due to the high Mach numbers (see appendix Appendix 3) of these systems it is 1 order of magnitude above the thermal Toomre parameter.

Our NoSNe run tends to have lower values with a smaller dispersion compared with our SNe runs. In the case of no feedback the Toomre parameter fluctuates around 1 above showing an unstable disc at high redshift. At the parameter is of order inside and above this radius it increases due to the combined effect of low density and higher sound speed (high temperature) stabilizing the system at these radii.

Due to the higher temperature associated to SNe explosions the Toomre parameter tends to be larger in our SNe runs showing a more stable system in these cases. Despite that it is also possible to find regions with in our feedback runs. We have to take into account that after each SNe explosion a given amount of metals is released into the gas. Such a new component allows the gas to reach lower temperatures creating unstable regions.

We applied the clump finder algorithm of Padoan et al. (2007) to our galactic disc inside a box. The clump finder algorithm scans regions of density above , with the average density inside the analyzed box which is of order . In practice it means that we look for gas clumps at densities above . The scan is performed increasing the density by a fraction until the maximum box density is reached. For each step the algorithm selects the over-densities with masses above the Bonnor-Ebert mass in order to define a gravitationally unstable gas clump. This algorithm gave us clump masses in the range few . Figure 10 shows the clump mass function found in each of our simulations at different redshifts. In order to complement this analysis we have computed the mass associated to the maximum unstable scale length of a rotating disc (Escala & Larson, 2008)


The vertical lines of figure 10 mark the average at each sampled redshift. Our NoSNe run formed the bigger gas clumps. In this case due to the lack of feedback the most massive objects () can survive at different redshifts. Such mass is of the order of the expected .

The SNe runs could form objects as big as . These masses are below the few . The SNe0.5 simulation forms much more massive objects compared with the SNe5.0 run. Due to the extreme feedback of the SNe5.0 experiment it is not easy for the clump to survive in such a violent environment. This is why the SNe5.0 simulation forms less clumps throughout its evolution.

All the clumps formed in our simulations have sizes in the range of few to few (note that as this size is associated to all the cells above the threshold , then it is a minimum size because it could increase if we reduce ). These sizes are below the unstable length scale (averaged on the inner ) associated to the maximum clump mass: few .

3.4 Mass transport on the disc

It is well known that in a cosmological context the large scale () gas cooling flows associated to DM filaments converging onto DM haloes have influence on the small scales () galactic AM (e.g. Powell et al., 2011; Prieto et al., 2015; Danovich et al., 2015). Such an interplay between large and small scales suggests that the mass/AM transport analysis should be performed taking into account both regimes.

Figure 2: Combined rest frame stars visualization for our three runs using SDSS , and filters in blue green and red colors, respectively. The images correspond to the end of our simulations and there is no dust extinction. The face on view of the NoSNe system shows a smoother star distribution compared with our SNe runs where feedback is able to create a non-homogeneous star distribution characterized by green-blue star clumps around the center of the system.

Stresses on the disc

The MT on the galactic disc can be studied based on the momentum conservation equation. Written in its conservative form this equation tell us that the local variation of momentum is due to the rate of momentum fluxes:


where is the gas density, are the Cartesian coordinates and are the Cartesian components of the gas velocity. All the terms inside the divergence are related with the rates of momentum flux and they can be written as follow:


is the term associated to the Reynolds (or hydrodynamic) stress. It is the momentum flux term associated to the total fluid movement. Instead of being a momentum flux source it quantifies the transported momentum due to the addition of different phenomena on the disc, namely gravitational stresses, magnetic stresses, viscous stresses or pressure stresses.


This is the pressure term, where is the Kronecker delta symbol, is the gas pressure and its gradient will be a source of torque as we will show in the following lines.


with the gravitational potential and Newton’s constant. is the term associated to the gravitational stress and it is related with the movements of the fluid due to the gravitational acceleration. This term also will be a source of torques acting on the fluid as we will show later.

Because we are not including magnetic fields we have neglected the term associated to it. Furthermore, the dissipative-viscous term is negligible in this context and will not be taken into account in the following discussion (e.g. Balbus, 2003).

In the disc MT context it is useful to quantify the momentum transport in the direction due to processes in the direction where and are the radial and the azimuthal cylindrical coordinates, respectively. If is the rate of momentum flux in the direction due to the processes in the direction associated to any of the stresses mentioned above, in general we can write (see appendix Appendix 1)


After some algebra it is possible to write the momentum fluxes for each of our sources as follow (e.g. Balbus, 2003; Fromang et al., 2004, and references therein):


It is worth noticing that in the case of symmetry the gravitational term vanishes. In other words, any density perturbation in the direction, e.g. an asymmetric density distribution of gas clumps in the disc, will cause a momentum flux in the direction. This will be the term associated to the VDI as we will show later.

The terms associated to the Reynolds and the gravitational stress as defined in the above expressions are averaged in space in order to quantify the radial momentum flux associated to perturbations in the azimuthal direction (Hawley, 2000) The Reynolds and the gravitational stress are defined as follow6:


where , the average circular velocity of the fluid and the averages are computed as


In this context it is useful to define an parameter for each of our stresses. For a given rate of momentum flux, following Gammie (2001), we define:


Each parameter is interpreted as the rate of momentum flux associated with a given process normalized by the gas pressure. Because the gas pressure is it can be interpreted as a “thermal momentum” advected at the sound speed or as a “thermal rate of momentum flux”. In this sense an is a sign of super-sonic movements in the fluid. This parameter is for ionized and magnetized discs (Fromang et al., 2004; Nelsen & Papaloizou, 2003) and observations of proto-stellar accretion discs (Hartmann et al., 1998) and optical variability of AGN (Starling et al., 2004) give an alpha parameter . Due to the turbulent (i.e. high Mach number, see appendix Appendix 3) nature of the environments studied here, the alphas will typically be higher than 1. In fact, with the gas Mach number.

Figure 11 shows the radial values of in the top row and in the bottom row for our three simulations in different columns.

The first thing that we should notice from this figure is that the Reynolds parameters are not constant neither in time nor in space and furthermore they reach values well above unity. In other words, our high redshift galactic discs are not in a steady state. Such a dynamical condition does not allow use of the Shakura & Sunyaev (1973) mass accretion rate expression as a function of the computed parameter. (See appendix Appendix 2 for a more detailed discussion.) Instead of that we must compute a mass accretion rate directly from our data.

From figure 11 it is clear that the Reynolds stress tends to be much larger than the gravitational stress and then it dominates the MT process in most of the cases (note that both top and bottom panels are not in the same range). In other words, the rate of momentum flux associated to the gravitational potential gradients is lower than the rate of momentum flux associated to local turbulent motions of the gas in most of the cases. Such high values of are associated to high velocity dispersions which can be an order of magnitude above the sound speed. We note that our two SNe runs have lower due to the higher sound speed in their environment.

Here we emphasize that the Reynolds tensor is not a source of momentum flux, in the sense that if we start the disc evolution from a spherical circular rotation state, i.e. without a radial velocity component, with null viscosity and one small gravitational potential perturbation in the direction the variation in momentum will be associated to the gravitational stress and the appearance of the Reynolds stress will be a consequence of this process.

It is interesting to note that the parameter in our NoSNe and SNe0.5 has a decreasing trend with the galactic radius at some redshifts: the smaller the radius the larger the gravitational stresses. If we take into account that the accreted material tends to concentrate in the inner part of the galaxy then it is reasonable that the larger gravitational stresses act at small radii. In the NoSNe run it is of the order of the pressure at the galactic center at all redshifts, whereas in the SNe0.5 run it is comparable to the pressure at high z. Due to the high feedback the SNe5.0 run is dominated by the Reynolds stress in all the sampled redshifts.

Figure 3: Radial profiles of the gas SD (top panels), stars SD (central panel) and SD SFR (bottom panel) for the NoSne run in the left column, SNe0.5 run in the central column and SNe5.0 run in the right column. All the quantities are plotted for different redshifts (from 6 in blue to 10 in green): (solid line), (long-dashed line), (short-dashed line), (dotted line) and (dot-dashed line). The vertical lines mark the at each following the same line style as the profiles. The top and central panels show density fluctuations associated to the non-homogeneous nature of the galactic disc at high redshift. SNe feedback clearly decreases the amount of stars formed in our galaxies.

Torques on the disc

Figure 4: The figure shows the gas radial velocity (solid black line), orbital velocity (long-dashed blue line) and the spherical circular velocity (dot-dashed cyan line) as a function of radius at different redshifts for the NoSNe run (left columns), the SNe0.5 run (central column) and the SNe5.0 run (right column). The vertical lines mark the position.

After observing that the Reynolds stress associated to the gas turbulent motions dominates the rate of momentum flux in the disc and that the gravitational tends to reach its maximum at the central galactic region, it is relevant for the MT study to analyze the torques acting in the disc associated to forces in the direction. In order to do that we compute the torques associated to both the gravity and the gas pressure for our systems. We define these two quantities as:


which actually are specific torques, i.e. torques per unit gas mass. These two terms will act as a source of AM transport in the galactic disc and will give us some clues about the MT process in high redshift galactic discs. In order to compute this we have defined the radial origin to be in the cell where the sink particle is set.

Figure 5: KS law for each of our snapshots from to . Each point marks the KS relation for our galaxies computed as mentioned in the text at different redshifts in different colors. From top to bottom: NoSNe, SNe0.5 and SNe5.0. The solid black line marks the K98 fit. The thick long-dashed green line marks the D10 fit for normal galaxies and the short-dashed blue line marks the D10 fit for star burst galaxies. Our SNe0.5 run has the lowest scatter, and most closely follows the D10 sequence of star burst galaxies.

Figure 12 show the ratio between and , with . The NoSNe run shows a decreasing trend with radius, like in the profile. The pressure gradients tend to dominate above and the gravity force dominates in the innermost region. As already shown in the alpha profiles, in the SNe0.5 run the gravity dominates the central part of the system at high and at lower redshifts the pressure torques are the source of mass transport. And finally, due to the high feedback which is able to create strong shocks and destroy gas clumps, the SNe5.0 simulation tends to be dominated by torques associated to pressure gradients in line with the previous alpha results.

As a complement to our findings it is useful to take a look at torques at large scales. Figure 13 shows the ratio of the total torques (not only the component in the disc) for our three runs at two different redshifts, at in the top row and at in the bottom row.

The maps take into account the gas with density above , where is the critical density of the Universe. Such a cut in density was set by inspection in order to have a clear view of the filaments around the central DM halo.

It is interesting to note that in our three simulations the border of the filaments is clearly dominated by the pressure torque: material from voids falls onto the filamentary DM structure creating large pressure gradients (Pichon et al., 2011; Danovich et al., 2015). There it loses part of its angular momentum and flows onto the DM halo.

At high redshift it is possible to see that inside the filaments the gravitational torque is of the pressure torque. The picture changes when we look at the bottom panels, there the shocked filaments have a ratio .

It is possible to find regions with a ratio around gas over-densities and near the main central halo. All the gas over-densities, in general associated to DM haloes at these scales, have a higher gravitational to pressure torque ratio. In particular at we can see that the central galactic region for the NoSNe simulation is dominated by the gravitational torque. This is not the case for both our SNe runs where at low redshift the pressure torque is dominating the AM re-distribution. Such behavior confirm the radial profile results of figure 12 and the parameters of figure 11.

At the edge of the galactic disc the pressure torque associated to the in-falling shocked material tends to dominate AM variations whereas at the central region the gravitational potential gradient is the main source of torque in the NoSNe simulation. In the SNe runs the energy injection spread out the high density material and there is a more flat potential at the center of the galaxy implying a non clear gravity domination there. Such behavior is more evident in our SNe5.0 run where it is possible to see a dark region in the center of the system.

Figure 6: Halo mass - stellar mass (normalized by the halo baryonic content assuming that it has exactly the universal baryonic fraction ) relation as a function of redshift in different colors for our runs. From top to bottom: NoSNe, SNe0.5 and SNe5.0. The filled circles mark the mass in stars inside the virial radius and the empty triangles mark the mass in stars inside one tenth of the virial radius. Our SNe5.0 run is in closest agreement with B13 at while for high galaxies our SNe0.5 experiment is still in the limit of the order of magnitude predicted by B13.

Having clarified that the source of pressure gradients are the shocks associated to both the filamentary incoming material from the cosmic web and the SNe explosions, it is interesting to elucidate the origin of the gravitational torque acting mainly in the central region of the galactic disc. In order to do that it is useful to study the density distribution in the disc. In particular, it is worth computing the Fourier modes associated to the gas mass surface density:


Figure 14 shows the square of each Fourier mode (from to ) normalized by for five different redshifts for our three runs. It is clear from the figure that the mode has the highest power in the spectrum for all the shown redshifts. Despite that, it is also possible to see that the difference in power between the first and the second mode is not too much for all the sampled redshifts, i.e. . In this sense it is not possible to say that the first mode is the only contribution to the surface density spectrum because the second mode (and even the third one) is also important. Furthermore, it is worth noticing that the powers with are also there and they have values below . It is interesting to compare our result with the one shown in Krumholz et al. (2007) where they found that the source of the torques on a proto-stellar disc is associated to the domination of the mode due to the SLING instability (Adams et al., 1989; Shu et al., 1990). In their case the first power is at least one order of magnitude higher than the mode with an increasing domination of mode with time. They argue that the spiral mode produces global torques which are able to efficiently transport AM. In our case, the global perturbation will be associated to a more complex disc structure. The reason for this difference will be clear after looking at the surface gas density projections.

Figure 15 reflects the fact that the power spectrum of the gas surface density shows power for different modes . There we can see a complex spiral-clumpy structure defining the galactic disc. Such density field features create an inhomogeneous gravitational potential field which will exert torques on the surrounding media. In particular, the clumps formed on the disc by gravitational instabilities interact between themselves migrating to the central galactic region: the VDI acts on these high redshift clumpy galactic discs (Bournaud et al., 2007).

Figure 7: SFR as a function of redshift for our experiments: NoSNe (black solid line), SNe0.5 (dashed blue line) and SNe5.0 (dot-dashed line). It is worth to notice that above redshift the NoSne simulation presents a higher SFR compared with the SNe simulations: the NoSne experiment can form stars continuously without feedback. The SNe runs shows a “bursty” nature with peaks of SF in each few .

It is worth noticing that due to the SNe energy injection in the SNe runs the disc takes a longer time to appear comparable to the NoSne experiment. Whereas the NoSne simulation develops a disc that is progressively disturbed in time by no more than mergers, the SNe runs show a disturbed clumpy environment characteristic of turbulent gas where the SNe explosions disrupt the galaxy with a strong effect on the central BH accretion rate as we will see in the next sub section. Furthermore, due to the metal release in the SNe runs the gas can cool more efficiently than in the NoSNe simulation. Such an important difference allows the gas to form more self gravitating over-densities and produce the clumpy galaxies shown at redshift in the second and third columns of figure 15.

Mass accretion and BH growth

High redshift galaxies are far from isolated systems. In fact, as has been shown above, they are very dynamic, in the sense that they are being built up in environments disturbed by filamentary accretion from the cosmic web, mergers and SNe explosions which affect the transport of AM. In the context of BH evolution at high redshift it is relevant to study and quantify the mass accretion rate in the galactic disc due to the processes described in the previous sub section, and more relevant yet is the quantification of the mass accretion rate onto the central BH and the relation of its mass accretion with the large scale filamentary inflows.

Figure 16 shows the radial gas mass accretion rate on the disc as a function of radius inside at different redshifts for our three runs. We defined the mass accretion on the disc as:


The radial coordinate is defined in the disc plane and the gas SD and the radial velocity are cylindrical shell averages in the direction.

We note that the NoSNe simulation shows continious lines at almost all redshifts and radii, but our SNe simulations present non continuous lines due to gas outflows. Such features are another proof of the highly dynamic environment where the first galactic discs are formed. The mass accretion rate fluctuates roughly between and in the range of radius shown in the figure. As mentioned above such a huge dispersion reflects the fluctuating conditions of the galactic disc environment at high , where due to the continuous gas injection through filamentary accretion the evolution is far from the secular type we see in low redshift galaxies.

At the NoSNe run has an accretion rate of on the disc. This value is higher than the in the case of SNe0.5 and for the SNe5.0 simulation. This trend can be related to the SNe feedback strength. As we will see below, the SNe feedback will have important effects on the BH growth and on the mass accretion at larger scales also.

Figure 8: Gas mass fraction as a function of redshift for our three runs: NoSNe (solid black line), SNe0.5 (dashed blue line) and SNe5.0 (dot-dashed cyan line). The dotted thick green line at marks the W15 observed value. The dotted thin green line are the errors associated to that observation. In general our SNe runs are well inside the error bars of W15. Our SNe0.5 run has closer to the observed value at and it is just in the limit of when we average below .
Figure 9: Toomre parameter as a function of radius for our three simulations. The thermal Toomre parameter in the left column and the turbulent Toomre parameter in the right column. From top to bottom: NoSNe, SNe0.5 and SNe5.0. From the thermal Toomre panels we see that our NoSNe run presents lower values compared with our SNe experiments. Despite that the SNe runs do have some regions of the Toomre parameter due to the effect of metal line cooling. Our turbulent Toomre parameter shows much higher values. This is due to the high Mach number of our runs. The difference is more dramatic in our NoSNe run due to the low gas temperatures reached without SNe heating. See the text for a discussion about this parameter.

Figure 17 shows the ratio of the BH mass accretion rate to the Eddington mass accretion rate . The NoSNe run BH accretes matter at the Eddington limit until . At this redshift the central galactic region undergoes several mergers losing a lot of gas, leaving the BH with almost no material to consume. After this event the accretion rate fluctuates until the end of the simulation. Our SNe run has an Eddington limited BH accretion rate throughout its evolution and a below .

From figure 17 the effect of SNe feedback on the BH growth is clear. SNe feedback perturbs the BH accretion rate from the beginning of its evolution decreasing its value until in the SNe0.5 run and reaching even lower values in the SNe5.0 simulation. Such a difference in the BH accretion rate is translated into a for SNe0.5 throughout its evolution and a for the SNe5.0 experiment. These values do not change significantly when we average below . The mass accretion rate onto the BH at the end of the simulations is , and for the NoSNe, SNe0.5 and SNe5.0 run, respectively (see appendix Appendix 4).

Figure 10: The gas clumps mass function. From top to bottom: NoSNe, SNe0.5 and SNe5.0. The vertical lines are the average maximum unstable mass for rotating disc . This scale mass has values for all our runs. We divided the mass range in the bins . Our simulations form clumps in the range few to few The NoSNe simulation reaches the higher mass clump with due to its null feedback. Our SNe runs reached lower maximum clump mass due to the SNe heating and our SNe0.5 experiment form more massive clumps compared with our SN5.0 run.

Figure 18 shows the BH mass evolution as a function of redshift. From this figure it is possible to see the effect of the different behavior in the accretion rate, as shown in figure 17. Whereas the NoSne sink has an approximately exponential evolution ending with a mass , due to the SNe feedback our SNe runs show episodes of no growth at some redshift. Such a feature is much clearer in our SNe5.0 run (see between or for example). The final mass in these two runs was for SNe0.5 and for SNe5.0.

Figure 11: The parameters in the disc as a function of radius for different redshifts. The gray dashed line marks the position. From left to right: NoSNe, SNe0.5 and SNe5.0. The top row shows the hydrodynamic and the bottom row the gravitational . The are much larger than one in most of our cases. Due to the high Mach number of the NoSNe run it shows higher values. Due to the high temperatures reached after SNe explosions our SNe runs have lower Reynolds alphas. The gravitational is much lower than the Reynolds stress. It reaches higher values at the central galactic region for NoSNe and SNe0.5 (). The SNe5.0 run shows lower values due to the destruction of dense gas features and higher gas temperatures.

3.5 Mass transport on larger scales

At high redshift we cannot study the small-scale galactic phenomena without taking into account the effects of the large scale structure in a cosmological context. Here we study the behavior of the mass accretion rate above the kpc scales, i.e. beyond the galactic disc edge.

Figure 19 shows the mass accretion rate out to . The mass accretion has been computed taking into account all the mass crossing a spherical shell at a given radius centered at the sink cell position:


The left column of figure 19 shows the total mass accretion rate for our three simulations. In the right column we have plotted the mass accretion rate associated to gas densities below , with the critical density of the Universe. The vertical lines mark the DM virial radius at each sampled redshift.

The right column of the figure can be interpreted as smooth accretion associated to non collapsed objects. The NoSNe panel shows a smooth decreasing behavior almost independent of redshift above . The smooth mass accretion rate presents roughly constant values above the virial radius, with accretion rates of the order . Such a value is consistent with the one found by e.g. Dekel et al. (2009) and Kimm et al. (2015) for a halo at high redshift (see also Neistein et al., 2006). On the other hand, the SNe run panels have a more irregular decreasing behavior with a notable dependence on redshift due to SNe explosions. In these runs the SNe feedback is able to heat up the gas and create hot low density gas outflows almost depleting the system of low density gas. In particular, for a number of redshifts it is possible to see that the smooth accretion is practically erased at radii , a clear signal of low density gas evaporation. In other words, due to SNe explosions only the dense gas is able to flow into the inner region of the galaxy.

Figure 12: Gravitational torque to pressure gradient torque ratio as a function of radius for different redshifts. From top to bottom: NoSNe, SNe0.5 and SNe5.0. The gray dashed line marks the state. The NoSNe run MT tends to be dominated by gravitational torque in the inner galactic region at all redshifts. Beyond that radius the pressure gradients and gravity work together to re-distribute AM. The SNe0.5 run tends to be dominated by pressure with a central gravity domination at high , inside . The pressure gradient domination is more clear in our SNe5.0 due to the extreme SNe feedback.
Figure 13: Modulus of the mass weighted gravitational to pressure gradients torque ratio at in the top row and at in the bottom row. From left to right: NoSNe, SNe0.5 and SNe5.0. It is interesting that the pressure torque dominates over the gravitational torque in most of the mapped filamentary dense regions. The gravitational torque increases its influence at the central region of filaments and around gas over-densities. Such a fact confirms our previous finding based on the torques ratio radial profiles: gravitational torque increases its influence in the central galactic region. Such a behavior is not true in our SNe runs where the SNe feedback creates a region dominated by pressure gradients at the galactic center.

3.6 Gas-stars-DM spin alignment

As a complementary analysis it is interesting to study the alignment between the AM of the different components of the system, namely DM, gas and stars. Figure 20 shows the alignment between the AM of the different components of our systems. The misalignment angle between the gas AM and the component of the system was computed as:


The rotational center to compute the gas, DM and stars AM was set at the sink position, . This point coincides with the gravitational potential minimum cell within around the sink particle.

Figure 14: The Fourier modes associated to the mass surface density spectrum on the disc. From top to bottom: NoSNe, SNe0.5 and SNe5.0. Despite of the mode has the higher power, the mode is also important being a fraction of the at all redshifts. Furthermore, all the other modes have a contribution of roughly similar order between them. In other words, the disc have developed a complex azimuthal structure allowing gravitational torques on the galaxy.

We computed the AM of the different components as


where the sum is calculated inside for each component7. is the the center of the cell where the sink particle is located, is the average gas velocity of all cells inside a radius of around the sink position and is the mass of our different quantities: gas, stars and DM.

From the figure we can see that the gas and the DM spins are far from aligned. The misalignment angle between them fluctuates from a parallel alignment to an almost anti-parallel configuration in our NoSNe experiment. In our SNe runs the fluctuations are more dramatic due to SNe explosions. Such a non-correlation between the AM vector of these two components has been studied before in, e.g. Prieto et al. (2015). In their work the authors noticed that after the cosmological turn around the gas can decouple from the DM due to its collisional nature: while the DM can feel only the gravity the gas can also feel the gas pressure gradients. Such pressure gradients are responsible for an extra torque on the baryonic component and its AM vector deviates from the DM AM orientation. As already shown in figure 12 the pressure gradients are not negligible inside the virial radius of our haloes. Such torques are able to change the orientation of the gas AM and then create a misalignment between gas and DM AM vectors.

Figure 15: Gas surface density projections for our runs at different . From left to right: NoSNe, SNe0.5 and SNe5.0. The evolution of the density maps show that the galaxy develops a complex spiral clumpy structure supporting the existence of high powers in the Fourier analysis of figure 14. Due to SNe explosion the spiral shape of the object appears only at in the SNe0.5 run. Below this redshift the galaxy is successively destroyed by SNes and re-built by gravity. In our SNe5.0 run it is almost impossible to see a spiral shape due to the extreme feedback.

The alignment between gas and stars has a different behavior in our runs. For the NoSNe case it is possible to see that at high redshift, between and , the stars and the gas had a very different spin orientation. This is because at this stage the galaxy is starting to be built by non spherical accretion and mergers, conditions which do not ensure an aligned configuration. After the gas and stars reach a rather similar spin orientation with a misalignment angle fluctuating around the value of . There the proto-galaxy can not be perturbed easily by minor mergers and acquires a defined spiral shape allowing the gas-star alignment. Such an alignment is perturbed at redshift . At this redshift the main DM halo suffers a number of minor mergers which can explain the spin angle perturbation. After that the gas and stars again reach an aligned configuration which will be perturbed by mergers again at lower redshifts.

The SNe runs show a much more perturbed gas-stars AM evolution. In this case, as well as the merger perturbations the systems also feel the SNe explosions which continuously inject energy. The strong shocks associated to this phenomenon are able to decouple the gas AM from the stellar AM as we can see from the blue solid line. Such perturbations are more common in our SNe5.0 run compared with our SNe0.5 run due to the stronger feedback as we can see from figure 15 where the SNe5.0 simulation shows a number of clumps instead of a defined spiral shape.

Figure 16: Mass accretion rate radial profiles for our three runs. The vertical lines mark at each redshift. In all simulations the accretion rate has huge fluctuations between few and few , another proof of the highly dynamic nature of the system. The SNe runs show a less continuous accretion with lower values. In particular, our SNe simulations have a few at the end of the simulation.

4 Discussion and Conclusions

By using cosmological hydrodynamic zoom-in simulations we have studied the MT process from few to few scales on a DM halo of at redshift . We have studied the evolution of the system without SNe feedback (NoSNe run) with the delayed cooling model for SNe feedback (SNe0.5 run) and with an extreme case of delayed cooling SNe feedback (SNe5.0 run).

We found that the SNe0.5 run is the best match with the D10 star burst galaxy sequence. It covers about two decades in SD with the lowest scatter among our simulations. When we look at the stellar mass of the systems our SNe5.0 run shows a stellar mass close to the expected value from B13, . Looking at this quantity our SNe0.5 run is still in the order of magnitude compared with B13 for a DM halo at high redshift. Such an offset can be related to the “bursty” nature of high redshift galaxies. In terms of the SFR, due to the extreme feedback, our SN5.0 run has the lowest values with a SFR at . At the same range our SNe0.5 run has a SFR in agreement with results from the W15 high galaxies. Despite this both SNe runs present low () episodic SFR values due to the SNe heating.

Our SNe experiments show lower gas fractions among our three simulations. They have values below . If we look at the gas fraction below our SNe0.5 run has the lowest value with which is just within the upper limit for the SFR of the galaxy found by W15.

Following Gammie (2001) we have computed the parameters associated to both the Reynolds and the gravitational stresses. In other words, we have computed both the Reynolds and the gravitational rate of momentum fluxes on the disc normalized by the gas pressure. Gammie (2001) showed that the parameters associated to radial mass transport are of order , reasonable values for a subsonic stationary accretion disc. In our case the parameters reach values above unity, meaning that the rate of momentum flux has values higher than the gas pressure . Such high values are characteristic of a turbulent super-sonic environment associated to dynamical systems like the ones in our simulations. The highly non-stationary gas behavior is confirmed also by the highly fluctuating values of at all redshift.

We found that the Reynolds stress dominates over the gravitational one in most of the analyzed redshifts. Here it is worth noting that the Reynolds stress tensor is a measurement of the turbulent motions in the gas. In these systems the gas falls from large scales, channeled by filaments almost freely onto the DM halo central region, gaining super-sonic velocities. Through the virialization process strong shocks are created developing a turbulent environment which is enhanced due to SNe explosions. Under such conditions the rate of momentum flux associated to this term normalized by the gas pressure will be much higher than 1 if the rms gas velocity in the and directions are super-sonic.

We emphasize that the Reynolds stress is not a source of mass transport but it is a measurement of the local rate of momentum transport triggered by other processes, namely pressure gradients, gravitational forces, magnetic fields or viscosity. In this sense, its high value simply tells us that throughout galaxy evolution there exists processes capable of transporting mass from large scales to small scales very efficiently. In fact, in our systems, gravity triggers the mass flows through the DM filamentary structure around the central halo and then a combined effect of gravity and pressure gradients allows the MT in the disc. The Reynolds term tends to be higher for our NoSNe run where Mach numbers are higher due to the null SNe heating.

The gravitational parameter has a different behavior in our three experiments. The NoSNe run shows a clear decreasing trend in radius until for all sampled redshifts. Beyond that fluctuates between values lower than 1. It reaches values at the central region. Such a behavior is telling us that the gravitational term is more important at the central galactic region where matter is more concentrated.

Our SNe0.5 run has a peak above unity in the central region at high redshift decreasing until . Beyond that radius it has a similar behavior compared with our NoSNe simulation. Below redshift the gravitational alpha parameter reduces its value to around with a lot of dispersion but always below . In this case the SNe feedback is able to deplete the central galactic region of gas after reducing the stresses associated to the gravitational gradients.

The SNe5.0 simulation has no peak at the galactic center. It seems to have the lowest values at the central regions. Due to the extreme feedback adopted in this simulation it is much more difficult for gas to create dense structure producing important gravitational forces. Furthermore, in this case the gas maintains higher temperatures implying higher pressures counteracting the gravitational effect.

Figure 17: BH mass accretion rate normalized by the Eddington accretion, as a solid black line for the NoSNe run, as a dashed blue line for the SNe0.5 run and as a dot-dashed line for the SNe5.0 run. From the figure the increasing feedback perturbation on the BH growth is clear. It has an average value of in the NoSne case and in our SNe runs throughout the BH evolution.

The torques acting in the disc show the sources of angular momentum variations triggering the MT process in these galaxies. In our systems the sources of torques are the pressure gradient due to the shocks created through the virialization process and SNe explosions, and the gravitational forces associated to gas inhomogeneities.

As in the analysis when we compute the gravitational to pressure gradient ratio our NoSNe run shows a decreasing trend in radius. Gravity dominates over pressure and has a maximum at the central galactic regions reaching values at radius . Beyond this radius the pressure gradients tend to dominate the AM re-distribution. Without SNe feedback the pressure domination at large radius is associated to shocks created by the large scale in-falling material to the central region of the host DM halo. Despite the domination of the pressure gradients in the outer regions the system shows a number of regions where gravity acts showing a mixed contribution for the MT process.

In our SNe runs the domination of gravity at the central regions is not as clear as in the NoSNe run. In the SNe0.5 simulation the gravitational gradients dominate above inside . At lower redshifts the pressure gradients clearly dominate the torques at the inner . Beyond that radius again it is possible to see a mixed torque contribution to the MT. A similar scenario is shown in our SNe5.0 simulation. In this case the gravity can dominate the very central regions ( few ) at high redshift and the pressure gradients have a more clear domination at larger radii, but there is still a mixed contribution to the AM re-distribution.

When we look at the large scales related with the filamentary structure around the central DM halo it is possible to see that pressure torques dominate over gravitational torque in filaments. The central region of the filaments show an enhanced gravitational contribution, but it is not enough to be dominant. These results are consistent with the picture in which the material filling the voids falls onto the filamentary over-densities where it is channeled to the central DM halo region (Pichon et al., 2011; Danovich et al., 2015) by gravity. Once the gas reaches the filaments it feels the pressure gradient on the edge of the filaments and it loses part of its AM. Then gravity acts and transport the mass almost radially inside the cold filaments to the central region of the DM halo. Such a process allows the gas to reach the galactic edge almost at free-fall. There the gas pressure acts reducing its initially high radial velocity and at the same time exerting torques allowing the MT. Throughout this process the gravitational torques also work in the galactic gas helping the MT process in the disc.

Figure 18: BH mass evolution for our three simulations: NoSne (solid black line), SNe0.5 (dashed blue line) and SNe5.0 (dot-dashed cyan line). The NoSNe BH reaches a mass of at the end of the simulation. Such a high mass was reached because most of the time the BH was accreting at the Eddington limit. In the SNe0.5 run the sink particle reaches a final mass of and our SNe5.0 BH mass reaches due to the extreme feedback.

A Fourier analysis of the disc gas surface density field for our runs shows that the density power spectrum has a number of excited modes. Despite the and modes dominating the power spectrum, the other modes do exist and have roughly comparable values between them. Such features tell us that the gas SD develops a complex structure throughout its evolution. The information given by the Fourier analysis is confirmed by visual inspection. The galactic discs develop spiral arms and gas clumps which interact between them by gravity. The gas clumps are formed from the cold gas flowing from the cosmic web onto the central DM halo region. The high gas fraction (which is ) and cold environment is a perfect place to produce a clumpy galactic disc. The interaction between gas clumps, spiral arms and merged DM haloes exert gravitational torques which are capable of transporting mass onto the galactic center in times comparable to the dynamical time of the system: this is the so-called VDI (Mandelker et al., 2014; Bournaud et al., 2007).

Due to the process described above, i.e. large scale gravitational collapse inducing filamentary accretion onto the DM central region and both gravitational and pressure torques acting in the galaxy, the mass can flow through the galactic disc and reach the galactic center. The radial mass accretion rate inside has huge fluctuations with values in the range for our SNe runs, a clear proof of a non-stationary and highly dynamic environment.

Figure 19: Left column: Same as figure 16 but for larger radii taking into account material till around the central halo. From top to bottom: NoSNe, SNe0.5 and SNe5.0. Right column: Same as left column but for the smooth accretion, i.e. for gas with a density below the collapse density . Beyond the virial radius the total accretion has a floor similar to the smooth accretion. Inside the virial radius the accretion rate is dominated by dense gas. The SNe explosions have a clear effect on the smooth accretion. At kpc scales the smooth accretion is practically erased due to the SNe heating.

The high mass accretion rate in the high gas fraction disc allows the central BH to grow at the Eddington limit most of the time for the NoSNe run whereas in the SNe runs it is clearly affected by the SNe explosions showing an intermittent Eddington-limited accretion rate. Despite this it can increase its mass substantially throughout the simulation. The violent events, namely mergers (which can also trigger mass accretion torquing the gas in the disc) and SNe explosions are not enough to stop the BH growth. The BH seed can evolve until in our NoSNe experiment, in our SNe0.5 and in our SNe5.0.

When we look at the mass transport beyond the virial radius we find that the large scale mass accretion rate has a floor of the order with peaks associated to gas inside DM haloes of few in all our runs (consistent with Dekel et al., 2009; Kimm et al., 2015). Inside the virial radius the smooth accretion decreases monotonically reaching values in the galactic outer regions, i.e. . These values change dramatically when we look at our feedback simulations. In these cases the SNe feedback practically depletes the galactic central region of low density gas. Due to the strong feedback effect only dense gas is able to reach the outer regions of the central galaxy. The mass accretion rate associated to dense gas is of the order in our NoSNe systems and it is almost devoid of discontinuities. On the other hand, despite the SNe runs reaching similar accretion rates in the disc, they do have discontinuities, i.e regions of zero accretion rate, affecting the amount of gas reaching the outer galactic region.

At the end of the simulation our NoSNe run shows an accretion rate which is similar to the total mass accretion in the disc. Contrarily, our SNe0.5 run ends with and our SNe5.0 run reaches at the end of the experiment showing how important are the SNe explosions to the BH accretion rate.

The gas AM vector orientation fluctuates a lot with respect to the DM spin vector through out the system evolution. The gas and DM start their evolution with spin vectors roughly aligned but once the pressure gradients increase due to virialization shocks, mergers (e.g. Prieto et al., 2015) and SNe explosions they decouple reaching an almost anti-parallel orientation at some stages. The alignment between these two components is more clear in our NoSNe run where the angle between them is below . The picture changes when we look at our SNe simulations, there the effect of SN feedback is capable of changing the alignment from to in few . Such an effect is stronger in our SNe5.0 where the big angle fluctuations are present throughout the entire system evolution.

The inclusion of AGN feedback in our simulations certainly could change both the galaxy and the BH evolution. The strong energy release in the gas can increase the gas temperature and may suppress the star formation changing SFR properties of those objects. Furthermore, due to the outflows associated to BH feedback the gas may not reach the central galactic region as easily as in the simulations presented here. A more detailed study of mass transport on high redshift galaxies with AGN feedback is left for a future study in preparation.

Figure 20: The misalignment angle between the gas AM and stellar AM (solid blue line), and the gas AM and DM (short-dashed cyan line). From top to bottom: NoSNe, SNe0.5 and SNe5.0. The gas and DM show a fluctuating misalignment angle with a high value at the end of the simulation. Due to the collisional nature of the gas it decouple from the DM once the pressure torques start to work on it. For the same reason the SNe simulations have a larger misalignment through out the simulation.

To summarize: In a cosmological context galaxies are formed inside knots of the cosmic web surrounded by filaments. The gas flows from voids to the DM filaments from all directions. There the gas piles up in the filamentary structure and its pressure gradient cancels part of its angular momentum. The pressure torques dominate the filamentary structure whereas the gravitational torques have a non-dominant enhancement at the center of filaments. Part of the material inside the filaments formed by dense cold gas, flows into the DM central halo in almost free-fall due to the host DM halo gravitational attraction. The transported cold gas reaches the DM halo with a high radial component of the velocity producing strong pressure gradients at the edge of the galactic disc. The constant inflowing of cold gas creates a high gas fraction and cold environment at the inner . Such conditions promote a low Toomre parameter in the disc and so it becomes gravitationally unstable, forming very efficiently gas clumps of masses in the range which interact between them, with the galactic spiral arms and with merged DM haloes. Such a clumpy environment produces regions in the disc dominated by gravity torques, in other words the gravitational torque due to the VDI acts as a source of MT in high redshift galaxies. The other, dominant, source of torques in our system is the pressure gradients. It is produced by SNe explosions and virialization shocks complementing the gravitational MT effect on this galaxy. The mass accretion rate in the disc triggered by pressure gradients and gravity can reach peaks of few and average values of few allowing an efficient BH mass growth.


J.P. and A.E. acknowledges the anonymous referee for the invaluable comments to improve this work. J.P. acknowledges the support from proyecto anillo de ciencia y tecnologia ACT1101. A.E. acknowledges partial support from the Center of Excellence in Astrophysics and Associated Technologies (PFB06), FONDECYT Regular Grant 1130458. Powered@NLHPC: This research was partially supported by the supercomputing infrastructure of the NLHPC (ECM-02). The Geryon cluster at the Centro de AstroIngenieria UC was extensively used for the analysis calculations performed in this paper. The Anillo ACT-86, FONDEQUIP AIC-57 and QUIMAL 130008 provided funding for several improvements to the Geryon cluster. J.P. acknowledges the valuable comments and discussion from Yohan Dubois and Muhammad Latif. J.P. and A.E. acknowledge to Marta Volonteri for her enlightening comments on this work.

Appendix 1

In order to compute the momentum flux in the direction due to processes in the direction we projected the tensor in the and then in direction: