Clumps with and without Radiative Feedback

# Giant Clumps in Simulated High-z Galaxies: Properties, Evolution and Dependence on Feedback

## Abstract

We study the evolution and properties of giant clumps in high- disc galaxies using AMR cosmological simulations at redshifts . Our sample consists of 34 galaxies, of halo masses at , run with and without radiation pressure (RP) feedback from young stars. While RP has little effect on the sizes and global stability of discs, it reduces the amount of star-forming gas by a factor of , leading to a similar decrease in stellar mass by . Both samples undergo extended periods of violent disc instability (VDI) continuously forming giant clumps of masses at a similar rate, though RP significantly reduces the number of long-lived clumps (LLCs). When RP is (not) included, clumps with circular velocity , baryonic surface density and baryonic mass are short-lived, disrupted in a few free-fall times. More massive and dense clumps survive and migrate toward the disc centre over a few disc orbital times. In the RP simulations, the distribution of clump masses and star-formation rates (SFRs) normalized to their host disc is similar at all redshifts, exhibiting a truncated power-law with a slope slightly shallower than . The specific SFR (sSFR) of the LLCs declines with age as they migrate towards the disc centre, producing gradients in mass, stellar age, gas fraction, sSFR and metallicity that distinguish them from the short-lived clumps which tend to populate the outer disc. Ex situ mergers comprise of the mass in clumps and of the SFR. They are more massive and with older stellar ages than the in situ clumps, especially near the disc edge. Roughly half the galaxies at redshifts are clumpy, with of their SFR and of their stellar mass in clumps.

###### keywords:
cosmology — galaxies: evolution — galaxies: formation — galaxies: kinematics and dynamics — stars: formation
12

## 1 Introduction

The peak phase of galaxy formation is at redshifts , when star-formation is at its peak and most of the mass is assembled into galaxies (Madau, Pozzetti & Dickinson, 1998; Hopkins & Beacom, 2006; Madau & Dickinson, 2014). Observations of massive star-forming galaxies (SFGs) of in baryons at this epoch reveal high star-formation rates (SFRs) of order (Genzel et al., 2006; Förster Schreiber et al., 2006; Elmegreen et al., 2007; Genzel et al., 2008; Stark et al., 2008). A large fraction of these galaxies have been spectroscopically confirmed to be rotating discs (Genzel et al., 2006; Shapiro et al., 2008; Förster Schreiber et al., 2009; Wisnioski et al., 2015). They are perturbed, thick and turbulent, with high velocity dispersions of , and low rotation to dispersion ratios of , as opposed to in today’s spiral galaxies (Elmegreen & Elmegreen, 2005; Genzel et al., 2006; Förster Schreiber et al., 2006, 2009; Cresci et al., 2009). Estimates of their gas fractions from CO measurements are in the range (Daddi et al., 2010; Tacconi et al., 2010, 2013; Genzel et al., 2015), much higher than the fractions of in today’s discs (Saintonge et al., 2011). Many of these galaxies exhibit irregular morphologies, in both rest-frame UV and rest-frame optical emission (Genzel et al., 2008; Förster Schreiber et al., 2009, 2011), with much of the UV light concentrated in a few large “clumps”, each a few percent of the disc mass and on the order of a in size, much larger than the star-forming complexes observed in local galaxies. Recent observations of over 3000 galaxies in the Cosmic Assembly Near-IR Deep Extragalactic Legacy Survey (CANDELS, Grogin et al., 2011; Koekemoer et al., 2011) reveal that roughly of the SFGs at are clumpy, with clumps accounting for of their UV light and of their SFR (Guo et al., 2015).

According to our developing theoretical framework of high redshift galaxy formation, these SFGs are in a perpetual state of violent-disc-instability (VDI, Dekel, Sari & Ceverino, 2009). Intense inflow of cold gas in narrow streams along the filaments of the cosmic web (Birnboim & Dekel, 2003; Kereš et al., 2005; Dekel & Birnboim, 2006; Ocvirk, Pichon & Teyssier, 2008; Dekel et al., 2009) maintains high gas fractions throughout the disc, repleneshing losses to star-formation, outflows, and compaction events which drive large amounts of gas into the disc centre (Bournaud et al., 2011; Forbes et al., 2014; Dekel & Burkert, 2014; Zolotov et al., 2015; Tacchella et al., 2016a, b). The high gas fraction together with the high density in the early Universe, leads to a violent gravitational disc instability (Toomre, 1964), which involves giant clumps and operates on short, orbital timescales (Dekel, Sari & Ceverino, 2009; Ceverino, Dekel & Bournaud, 2010), as opposed to the slow, secular instability in today’s discs.

The basic idea, summarized in Dekel, Sari & Ceverino (2009), is that during VDI the high surface density of gas and “cold” young stars, , drives the Toomre parameter below unity, , where is the one-dimensional velocity dispersion and is the angular frequency, a proxy to the epicyclic frequency , which is related to the potential well (Toomre, 1964). It has been established that under such conditions the disc will fragment and produce large star-forming clumps. This has been shown using idealized simulations of isolated galaxies (Noguchi, 1999; Gammie, 2001; Immeli et al., 2004a, b; Bournaud, Elmegreen & Elmegreen, 2007; Bournaud & Elmegreen, 2009; Elmegreen, Bournaud & Elmegreen, 2008; Hopkins et al., 2012), as well as cosmological simulations (Agertz, Teyssier & Moore, 2009; Ceverino, Dekel & Bournaud, 2010; Ceverino et al., 2012; Genel et al., 2012; Mandelker et al., 2014; Oklopcic et al., 2016). The ratio of clump mass to the mass of the cold disc scales as , where is the ratio of the cold disc mass to the total mass within the disc radius, which includes the bulge and dark matter halo (e.g. Dekel, Sari & Ceverino, 2009). This leads to much larger clumps at than the low-redshift giant molecular clouds (GMCs). Gravitational interactions in the perturbed disc drive turbulence causing the disc to self-regulate in a marginally stable state with (Dekel, Sari & Ceverino, 2009; Ceverino, Dekel & Bournaud, 2010; Krumholz & Burkert, 2010; Cacciato, Dekel & Genel, 2012; Forbes, Krumholz & Burkert, 2012; Forbes et al., 2014) that can last for more than a Gyr so long as the accretion is not interrupted. Some recent works have called into question the validity of linear Toomre analysis in the context of these highly non-linear galaxies (Behrendt, Burkert & Schartmann, 2015; Tamburello et al., 2015; Inoue et al., 2016) and others have suggested alternate fragmentation mechanisms related to turbulence (e.g. Hopkins, 2013). However, since clump formation is largely determined by the balance between self-gravity, turbulent pressure and the centrifugal force, the largest clumps are always roughly at the Toomre scale. Larger clumps would be disrupted due to the shear and/or tidal forces within the disc, or would not collapse in the first place due to the centrifugal force. Therefore, regardless of the full validity of linear Toomre analysis, it is plausible that the Toomre parameter can serve as a crude criterion for instability, possibly with a critical value that is larger than unity.

If clumps survive stellar feedback, their large masses cause them to migrate to the galactic centre on short timescales of orbital times at the disc edge (Dekel, Sari & Ceverino, 2009; Ceverino et al., 2012; Bournaud et al., 2014; Mandelker et al., 2014). This was proposed as a mechanism for the formation of galactic spheroids, in parallel with the traditional merger scenario (Bournaud, Elmegreen & Elmegreen, 2007; Elmegreen, Bournaud & Elmegreen, 2008; Genzel et al., 2008; Dekel, Sari & Ceverino, 2009; Ceverino et al., 2015) and the fueling of a central AGN (Bournaud et al., 2011, 2012). However, simulations show that most of the gas inflowing to the galactic centre is in fact inter-clump gas (e.g. Hopkins et al., 2012; Zolotov et al., 2015). Indeed, wet inflows to the galactic centre induced by gravitational torques are a generic feature of VDI, not necessarily associated with clump migration (Krumholz & Burkert, 2010; Forbes, Krumholz & Burkert, 2012; Forbes et al., 2014; Dekel & Burkert, 2014). While clump migration can contribute to the dry growth of the bulge, cosmological simulations indicate that this contribution is limited to of the bulge mass at (Zolotov et al., 2015). Nevertheless, the giant clumps themselves, their properties, and their evolution, are all major observational indicators of VDI. Clump migration is related to the overall inflow, and if signatures of this migration are observed, they would indicate the validity of VDI and VDI-driven gas inflow.

There is much debate regarding the ability of these giant clumps to survive feedback. Based on rather crude arguments, massive clumps of are not expected to be disrupted by supernova feedback alone (Dekel, Sari & Ceverino, 2009; Ceverino, Dekel & Bournaud, 2010; Ceverino et al., 2012; Bournaud et al., 2014; Mandelker et al., 2014). However, winds are observed from giant high- clumps (Genzel et al., 2011; Newman et al., 2012), with mass loading factors of order unity. If these outflows are very intense on timescales of free fall times, they could lead to significant mass loss and maybe even disruption of the clump (Genel et al., 2012; Hopkins et al., 2012). Murray, Quataert & Thompson (2010) argued that the high-z clumps are likely to be disrupted by momentum driven feedback, similar to GMCs at low redshift, but Krumholz & Dekel (2010) showed that this was unlikely, unless the star-formation efficiency per free fall time in the clumps is much higher than observed locally, which does not appear to be the case (Freundlich et al., 2013). Furthermore, Dekel & Krumholz (2013) argued that a steady wind generated by radiation pressure was not expected to unbind the clump before it had turned most of its mass into stars. Their argument was based on simulations by Krumholz & Thompson (2012, 2013) which showed that radiation trapping within the giant clumps is negligible because it destabilizes the wind. Thus, the total momentum injection from radiation pressure, stellar winds and supernovae is limited to , significantly lower than the values of advocated by Hopkins et al. (2012).

In previous work, we studied the formation and evolution of giant clumps at redshifts in the weak feedback limit (Mandelker et al., 2014, hereafter M14), using a large suite of 29 cosmological zoom-in simulations of haloes at , run with adaptive mesh refinement (AMR) with a maximal resolution of . Using a three-dimensional clump finder on the gas, we identified clumps in the simulations, distinguishing in situ clumps formed during VDI from ex situ clumps that joined the disc as minor mergers. The simulations included only thermal feedback from supernovae and stellar winds, under which the massive clumps survived and migrated towards the disc centre. Clump evolution during migration generated gradients in clump properties across the disc, so that clumps closer to the centre tended to be more massive, with older stellar populations, lower specific SFR (sSFR), lower gas fractions and higher metalicities. Observational indications for similar gradients (Förster Schreiber et al., 2011; Guo et al., 2012) provided support for the survival of massive clumps.

In a preliminary attempt to study the effect of radiation pressure (RP) feedback on clump survival, we examined a small subset of the simulations used in this work (Moody et al., 2014), consisting of 8 pairs of cosmological zoom-in simulations of haloes at , with a maximal AMR resolution of . Each pair consisted of two simulations run from the same initial conditions, one implementing only thermal feedback from supernovae and stellar winds and the second adding a model for RP feedback from UV photons, with no IR trapping (Ceverino et al., 2014). Using a two-dimensional clump finder on stellar mass maps, we found that the number of clumps with stellar mass less than of the disc stellar mass was drastically reduced with the inclusion of RP, but that more massive clumps were largely unaffected. We also found that clump counts in mock HST observations of the RP simulations roughly matched observational estimates. However, this study had several limitations. Firstly, we did not account for the significant gas content of clumps, studying only the stellar clumps. Second, while identifying clumps in two-dimensional stellar mass maps is closer to what is done observationally, it limits the analysis of the three-dimensional physical properties of the clumps and the role played by RP in their formation, evolution and disruption.

In this work, we expand on M14 and Moody et al. (2014). We use a large suite of 34 cosmological zoom-in simulations of galaxies, most evolved to , whose halo masses are in the range at . These 34 simulations include the same model for RP feedback, in addition to supernova feedback and stellar winds, as in Ceverino et al. (2014) and Moody et al. (2014), and 29 of them have counterparts run from the same initial conditions without the inclusion of RP. We identify clumps in three-dimensions, in both gas and stars, and study their physical properties. By tracking clumps through time using their stellar particles, we study the effect of RP feedback on the formation and evolution of giant clumps and find the conditions leading to clump disruption or survival. Our two main goals are to asses theoretically the conditions for clump survival when RP is included, and to find observational tests that can distinguish short-lived clumps from long-lived clumps and place constraints on models of stellar feedback and clump evolution.

This paper is organized as follows: In §2 we provide an overview of the simulations, including their limitations. In §3 we describe our analysis method and how we define galaxy and clump properties. In §4 we compare global properties of the two galaxy samples, with and without radiative feedback. In §5 we address the effect of radiative feedback on clumps. In §6 we address properties of clumps, observational signatures of clump migration and disc clumpiness. We discuss our results in §7 and summarize our conclusions in §8.

## 2 Simulations

We use zoom-in hydro-cosmological simulations of 34 moderately massive galaxies that comprise the VELA simulation suite (Ceverino et al., 2014; Zolotov et al., 2015). The fiducial sample includes an implementation of RP feedback from young stars (RP simulations, hereafter RPs), and 29 of the 34 galaxies were resimulated without RP (NoRP simulations, hereafter NoRPs), but with all other sub-grid physics unchanged. Several recent works have used the VELA suite to study a variety of issues relating to high redshift galaxy formation. These include the stellar to halo mass relation (Moody et al., 2014); galaxy clumpiness and morphological evolution (Moody et al., 2014; Snyder et al., 2015); the evolution of galaxy shapes (Ceverino, Primack & Dekel, 2015; Tomassetti et al., 2016); the link between metal inhomogeneities and gas inflows (Ceverino et al., 2016a); galaxy outflows (Ceverino et al., 2016b); compaction, quenching and the formation of blue and red nuggets (Zolotov et al., 2015) and their relation to the evolution of density profiles (Tacchella et al., 2016a) and the confinement of the main sequence of star formation (Tacchella et al., 2016b); and the relation between clump formation and Toomre instability (Inoue et al., 2016). In this section, we give an overview of the key aspects of the simulations. Further details regarding the numerical methods can be found in Ceverino et al. (2014).

### 2.1 Simulation Method and Sub-Grid Physics

The VELA simulations utilize the Adaptive Refinement Tree (ART) code (Kravtsov, Klypin & Khokhlov, 1997; Kravtsov, 2003; Ceverino & Klypin, 2009), which accurately follows the evolution of a gravitating -body system and the Eulerian gas dynamics, with an AMR maximum resolution of in physical units at all times3. The dark matter particle mass is and the stellar particles have a minimum mass of , similar to the stellar mass of an Orion-like star cluster. Each AMR cell is split into 8 cells once it contains a mass in stars and dark-matter higher than , equivalent to 3 dark-matter particles, or a gas mass higher than . This quasi-Lagrangian strategy ends at the highest level of refinement that marks the minimum cell size at each redshift. We often refine based on stars and dark-matter particles rather than gas, so within the central halo and the star-forming disc the highest refinement level is reached for gas densities between , and occasionally for densities as low as . In the outer circumgalactic medium, near the halo virial radius, the median resolution is .

Beyond gravity and hydrodynamics, the code incorporates the physics of gas and metal cooling, UV-background photoionization, stochastic star formation, gas recycling, stellar winds and metal enrichment, and thermal feedback from supernovae (Ceverino, Dekel & Bournaud, 2010; Ceverino et al., 2012), plus an implementation of feedback from radiation pressure (Ceverino et al., 2014). A subsample of the full VELA suite, 29 of the 34 galaxies, was run using the same initial conditions and sub-grid physics, but without the implementation of radiative feedback.

We use the CLOUDY code (Ferland et al., 1998) to calculate the cooling and heating rates for a given gas density, temperature, metallicity, and UV background, assuming a slab of thickness 1 kpc. We assume a uniform UV background, following the redshift-dependent Haardt & Madau (1996) model, except for gas densities higher than where we use a substantially suppressed UV background () in order to mimic the partial self-shielding of dense gas. This allows dense gas to cool down to temperatures of . The equation of state is assumed to be that of an ideal monoatomic gas. Artificial fragmentation on the cell size is prevented by introducing a pressure floor, ensuring that the Jeans scale is resolved by at least cells (Ceverino, Dekel & Bournaud, 2010). The pressure floor is given by

 Pfloor=Gρ2N2Δ2πγ (1)

where is the gas density, is the cell size, and is the adiabatic index of the gas.

Star formation is allowed to occur at densities above a threshold of and at temperatures below . Most stars () form at temperatures well below , and more than half of them form at in cells where the gas density is higher than . New stellar particles are generated with a timestep of . We implement a stochastic model, where the probability to form a stellar particle in a given timestep is

 P=min(0.2,√ρgas1000cm−3) (2)

In the formation of a single stellar particle, its mass is equal to

 m∗=mgasdtSFτ∼0.42mgas (3)

where is the mass of gas in the cell where the particle is being formed, and is a parameter of the simulations which was calibrated to match the empirical Kennicutt-Schmidt law (Kennicutt, 1998). We assume a Chabrier 2003 stellar initial mass function. Further details can be found in Ceverino et al. (2014).

The thermal stellar feedback model releases energy from stellar winds and supernova explosions as a constant heating rate over following star formation. The heating rate due to feedback may or may not overcome the cooling rate, depending on the gas conditions in the star-forming regions (Dekel & Silk, 1986; Ceverino & Klypin, 2009), as we do not explicitly switch off cooling in these regions. The effect of runaway stars is included by applying a velocity kick of to of the newly formed stellar particles. The code also includes the later effects of Type Ia supernova and stellar mass loss, and it follows the metal enrichment of the ISM.

In the fiducial (RP) models, radiation pressure is incorporated through the addition of a non-thermal pressure term to the total gas pressure in regions where ionizing photons from massive stars are produced and may be trapped. This ionizing radiation injects momentum in the cells neighbouring massive star particles younger than whose column density exceeds , isotropically pressurizing the star-forming regions. The expression for the radiation pressure is

where is set to half the cell size for the cell hosting a stellar mass , and to the cell size for its closest neighbours. The value of is taken from the stellar population synthesis code STARBURST99 (Leitherer et al., 1999). We use a value of which corresponds to the time averaged luminosity per unit mass of the ionizing radiation during the first of the evolution of a single stellar population. After , the number high mass stars and ionizing photons declines signifcantly. Since the significance of radiation pressure also depends on the optical depth of the gas within a cell, we use a hydrogen column density threshold of , above which ionizing radiation is effectively trapped and radiation pressure is added to the total gas pressure. See the “RadPre” model of Ceverino et al. (2014) for further details.

The initial conditions for the simulations are based on dark-matter haloes that were drawn from dissipationless N-body simulations at lower resolution in three comoving cosmological boxes (box-sizes of 10, 20, and 40 Mpc/h). We assume the standard CDM cosmology with the WMAP5 cosmological parameters, namely , , , and (Komatsu et al., 2009). Each halo was selected to have a given virial mass at and no ongoing major merger at that time. This latter criterion eliminates less than of the haloes, which tend to be in dense proto-cluster environments at . The target virial masses at were selected in the range , with a median of . If left in isolation, the median mass at would be . In practice, the actual mass range is broader, with some of the haloes merging into more massive haloes that host groups at .

### 2.2 Galaxy Samples

More than half the sample was evolved with RP to and all but 5 were evolved to . The NoRP runs were often stopped at higher redshifts, due to an overproduction of stars which slowed the simulations down considerably. The simulation outputs were stored and analyzed at fixed intervals in the cosmic expansion factor , , which at corresponds to about . Table 1 lists the final available snapshot for each of the 34 RP runs and the 29 corresponding NoRP runs in terms of expansion factor, , and redshift, . We detect the main progenitor in the final available output by using the AdaptaHOP group finder (Aubert, Pichon & Colombi, 2004; Tweed et al., 2009) on the stellar particles. It is then traced back in time until it contains fewer than 100 stellar particles, typically between ().

In §4 and §5 we directly compare the RP and NoRP samples using only snapshots available in both sets of simulations. In these cases, visual inspection of each snapshot is performed to ensure that the same galaxy is being traced in both samples. Occasionally, when there is a roughly equal mass merger, the automated algorithm may trace different progenitors in the RP and NoRP samples. Whenever this happens, we change the NoRP version to track the same galaxy as in the RP sample. In §6 we analyze the RPs independently of their NoRP counterparts, making use of the full sample.

### 2.3 Limitations of the Current Simulations

The cosmological simulations used in this paper are state-of-the-art in terms of high-resolution AMR hydrodynamics and the treatment of key physical processes at the subgrid level, highlighted above. However, like other simulations, they are not yet perfect in their treatment of the star formation and feedback processes. While the SFR recipe was calibrated to reproduce the KS relation and a realistic SFR efficiency per free fall time (Ceverino et al., 2014), the code does not yet follow in detail the formation of molecules and the effect of metallicity on SFR (Krumholz & Dekel, 2012). Additionally, the resolution does not allow the capture of Sedov-Taylor adiabatic phase of supernova feedback. The radiative stellar feedback assumed no infrared trapping, in the spirit of low trapping advocated by Dekel & Krumholz (2013) based on Krumholz & Thompson (2012, 2013). Other works assume more significant trapping (Murray, Quataert & Thompson, 2010; Hopkins et al., 2012; Hopkins, Quataert & Murray, 2012), which makes the strength of the radiative stellar feedback here lower than in other simulations4. Finally, AGN feedback and feedback associated with cosmic rays and magnetic fields are not yet incorporated. Nevertheless, the star formation rates, gas fractions, and stellar to halo mass ratios are all in the ballpark of the estimates deduced from abundance matching and provide a better match to observations than earlier simulations (Ceverino et al., 2014; Moody et al., 2014 and Fig. 3).

The uncertainties, and any possible remaining mismatches by a factor of order 2, are comparable to the observational uncertainties. For example, the stellar-to-halo mass fraction is not well constrained observationally at . Recent estimates by Burkert et al. (2015) (see their Fig. 5) based on the observed kinematics of SFGs reveal significantly larger ratios than the estimates based on abundance matching (Conroy & Wechsler, 2009; Moster et al., 2010; Moster, Naab & White, 2013; Behroozi, Conroy & Wechsler, 2010; Behroozi, Wechsler & Conroy, 2013) at . In §4, we present a detailed comparison of the stellar-to-halo mass relation in our simulations and the observational data (Fig. 3). We conclude that our simulations produce values that are in the ballpark of the observational estimates, and within the observational uncertainties.

If the SFR at very high redshifts is still overestimated, then this may lead to more clumps forming at earlier times when the gas fractions are still high, and less efficient clump formation at when gas fractions are underestimated. However, we find that the fraction of clumpy galaxies as a function of redshift roughly matches observations 6). Additional sources of feedback, from infrared trapping and AGN, may affect the lifetimes and properties of clumps. However, our simulations with intermediate feedback have the advantage of forming both long-lived bound clumps and short-lived disrupting clumps in the same galactic environment. By studying these two populations, we deduce observable differences between them (§6), which will place tighter constraints on models for clump evolution and for feedback. Furthermore, Ceverino et al. (2014) showed that photon trapping with boost factors of above the UV radiation, as advocated by Dekel & Krumholz (2013), changes the final stellar mass by and has only a minor effect on the density distribution of ISM gas. We therefore suspect that our current simulations capture the majority of the effect.

## 3 Analysis

In this section we describe how physical properties of the virial haloes, the galactic discs and the clumps are defined. We also describe the clump finding algorithm, the method for tracking clump histories and the distinction between in situ, ex situ and bulge clumps.

### 3.1 Physical Quantities

The virial mass, , is the total mass within a sphere of radius that encompasses an overdensity of , where and are the cosmological parameters at (Bryan & Norman, 1998). The virial properties for the 34 RP galaxies are listed in Table 1, together with the stellar mass, , gas mass, , and star-formation rate, SFR, within . These are quoted at except for the 5 galaxies that were stopped at higher redshift, marked by , for which we quote the properties at .

The stellar mass, , is the instantaneous mass in stars, after accounting for stellar mass loss. The simulation calculates stellar mass loss using an analytic fitting formula, where , and of the initial mass of a stellar particle is lost after , and respectively. We refer to the birth mass of stars, without accounting for stellar mass loss, as .

The SFR is obtained by , where is the mass at birth in stars younger than . The average is obtained for in the interval in steps of 0.2 Myr in order to reduce fluctuations due to the discreteness in stellar birth times in the simulation. The in this range are long enough to ensure good statistics. This represents the SFR on timescales, which is a crude proxy for H  based SFR measurements, while UV based measurements are sensitive to stars younger than .

We define the specific SFR as 5 and the gas fraction as . The stellar age of a system is defined as the instantaneous mass-weighted mean age of its stellar particles.

### 3.2 The Galactic Disc

We define a disc in each snapshot using gas with temperatures and stars younger than . Varying these thresholds in the range and has no significant impact on any of our results. We hereafter refer to this as the “cold” component of the system. The disc plane is determined by the angular momentum vector of the cold component. This vector, along with the disc radius, , and height (half thickness), , are computed iteratively until all converge to within , following the prescription in Appendix B of M14. Briefly, the disc radius, , contains of the cold mass within a cylinder of radius and height . The disc thickness, , contains of the cold mass within a cylinder of radius and height . The disc angular momentum is that of the cold material within a cylinder of radius and height . We varied in the range , in the range as well as the value of from the previous iteration, and the cold mass fraction in the range , and found no significant effect on our main results. See M14 for further details.

All gas within the final cylinder is assigned to the disc. However, a stellar particle is assigned to the disc only if it meets an additional kinematic criterion, whereby the -component of its specific angular momentum (parallel to the disc angular momentum), , is higher than a fraction of the maximum specific angular momentum for the same orbital energy, . Here, is the magnitude of the velocity of the stellar particle and is its cylindrical distance from the galactic centre. We adopt (Ceverino, Dekel & Bournaud, 2010, M14). Whenever there is a well defined disc, the vast majority of the cold mass obeys this criterion (see Ceverino, Primack & Dekel, 2015 and Tomassetti et al., 2016 for the typical stellar mass and redshift when galaxies develop well defined discs).

### 3.3 Clump Analysis

We identify clumps in 3D using a method similar to that of M14. The main difference is that in that work we only identified clumps in gas, while here we identify clumps in both the cold component and the stellar component. As was highlighted in M14 and in Moody et al. (2014), there is typically not a one-to-one correspondence between gas and stellar clumps, making it important to study both populations. We briefly review the clump finding method here and refer the reader to M14 for further details.

#### Clump finder

We detect clumps within a cube of side centred on the main galaxy. Using a cloud-in-cell (CiC) interpolation, we deposit mass in a uniform grid with a cell size of , so between 2 and 4 times the maximal AMR resolution, whose z axis is alligned with the normal to the disc plane. We then smooth this grid with a spherical Gaussian filter whose full width at half maximum (FWHM) is , and calculate the density residual at each point , where and are the unsmoothed and smoothed density fields respectively. We do this once for the cold mass and once for the stellar mass, and at each point take the maximum of the two residual values. We define clumps as connected regions containing at least 8 grid cells, corresponding to between 64 and 512 cells at the maximal AMR refinement, above a residual threshold of . This technique selects regions that are at least 10 times denser than their local surroundings in cold and/or stellar mass. A discussion on the sensitivity of our results to variations in the parameters of the algorithm, , and , is provided in §A.

Gas, stars and dark matter particles are deposited in the clump cells using CiC interpolation, and no attempt was made to remove unbound material from the clumps. We remove from the sample all clumps with baryonic mass . Hereafter, when we discuss clump mass we refer to the baryonic mass unless explicitly stated otherwise.

We define the clump centre as the baryonic density peak and the clump radius, , as the radius of a sphere with the same volume as the clump: , where is the total number of grid cells within the clump. This translates to a minimum clump diameter of , which is on the order of 7 cells at the maximal AMR resolution of . Since the pressure floor ensures the Jeans length is always resolved by at least 7 cells, this is the minimal clump size where we could in principle have a complete sample.

We assign a shape parameter to each clump, , using its inertia tensor, . The sum is over all cells in the clump, is the baryonic mass in the -th cell, is the displacement vector of the -th cell relative to the clump centre, , and is the Kronecker Delta. We diagonalize the tensor and find the three eigenvalues, . Finally, we define . For a perfect sphere, . For a flattened, oblate clump, . For a very prolate or filamentary clump, .

The clump finder is illustrated in Fig. 1, for one of the galaxies in the RP sample, V13 at (). On the left we plot the surface density of the cold component in a face on projection and on the right we plot the stellar surface density in the same projection. Besides being a large, clumpy, star-forming disc, this example highlights the importance of detecting clumps in multiple tracers, since there is not a one-to-one correspondence between stellar and cold clumps.

The disc radius, , is marked in both panels, as are all identified clumps with baryonic volume densities and masses . These thresholds were introduced here for clarity, because many low-mass, low-density clumps are not visible in 2D projection plots (recall that clumps are identified in 3D), and in any case for most of our analysis in §5 and §6 we will focus on clumps above these thresholds. There are two ex situ clumps (defined below) in this galaxy, which are marked with squares. The rest of the off-centre clumps are all in situ and are marked with circles, as is the central bulge clump (defined below). The size of the marker corresponds to the clump radius, , while the color represents its mass. Clumps with are marked in grey, clumps with are marked in red, clumps with are marked in blue, and clumps with are marked in black. The disc mass is and the gas fraction is .

#### Tracking individual clumps through time

The time between consecutive outputs is typically about half an orbital time at the disc edge. Since the expected migration time for clumps is roughly 2 disc orbital times, we are occasionally able to track individual clumps across several snapshots, which allows us to directly study the evolution and lifetimes of clumps. As the simulations were run with an Eulerian, grid based, code we cannot trace gas elements through time, and therefore track clumps using their stellar particles. We only attempt to track clumps that contain at least 10 stellar particles. For each such clump in a given snapshot, we search the preceding snapshot for all “progenitor clumps”, defined as clumps that contributed at least of their stellar particles to the clump in question6. If a given clump has more than one progenitor clump in the same snapshot, we rank them by baryonic mass, consider the most massive one the main progenitor and the others as having merged, thus creating a clump merger tree. If a clump has no progenitors in the preceding snapshot, we search the previous snapshots until we either find a progenitor or reach the initial timestep of the simulation.

We define the clump time, , as the time since clump formation. In the first snapshot when the clump is identified, we set to the stellar age of the clump, provided this is less than the time since the previous snapshot, when the clump was not identified. Otherwise, we set . In later snapshots when the clump is identified, we add to the physical time that has elapsed since its initial identification. Note that clump time differs from clump age, which is given by the mean stellar age and is thus affected by ongoing star-formation and tidal stripping of older stars (see Fig. 12 and Bournaud et al., 2014). The clump lifetime, , is the clump time at the final snapshot where the clump is identified. We refer to clumps that were identified in only one snapshot and have as zero-lifetime-clumps, or ZLCs. In §5 we will show that most of these are false positives, characterised by low densities and masses.

#### In situ, ex situ and bulge clumps

Following M14 we divide clumps into central “bulge” clumps, and off-centre clumps. A bulge clump forms the nucleus of the galactic bulge and is defined as a clump whose centre is within 2 grid cells of the galaxy centre. If there is more than one such clump, the most massive is considered the bulge clump and the rest are considered off-centre clumps. We discuss the prevalence of bulge clumps in the simulations in §5.

The off-centre clumps are divided into “in situ” clumps which formed within the disc through VDI, and “ex situ” clumps which are external mergers. As in M14 we distinguish in situ from ex situ clumps by their dark matter content. We compare the mean dark matter density within the clump, , to the mean dark matter density in a spherical shell (concentric with the host halo) thick around on the clump, , and define the dark matter overdensity . Figure 2 shows the probability density7, where , for the off-centre clumps in the NoRPs (blue line) and the RPs (red line). In both samples, the in situ clumps exhibit a roughly log-normal distribution centred on with a dispersion of , while the ex situ clumps form an extended tail to high values. We define ex situ clumps as clumps with , where is the maximal value of during the clump lifetime, and mark this threshold in Fig. 2. We discuss the contribution of ex situ clumps to the total clump population in §5.

An alternate definition of ex situ clumps could have been based on the birth place of the stellar particles within the clump. Among clumps with and , where our sample is complete and ZLCs are rare (see §5), these two definitions agree well. In of ex situ clumps selected by dark matter in the RPs (NoRPs), more than half the baryonic mass consists of stars formed outside the disc, as defined at the snapshot closest to the star particle’s birth. On the other hand, in both the RPs and the NoRPs, only of in situ clumps selected by dark matter have more than half their baryonic mass in ex situ stars.

## 4 Global Properties of the Galaxy Samples

In this section we examine the global effect of RP feedback on the galaxies in our sample, before studying its effect on giant clumps in §5. We limit our analysis here to snapshots that have both a RP and a NoRP version, which leaves 798 snapshots: 255 at , 131 at , 168 at , 244 at .

### 4.1 Stellar Mass

Figure 3 shows the stellar-to-halo mass ratio as a function of halo mass for the two galaxy samples at four redshifts, . 8 Stellar mass is defined here as the total within . For comparison, we show the abundance matching predictions of Behroozi, Wechsler & Conroy (2013), and recent estimates from (Burkert et al., 2015) which are based on observed kinematics of several hundred SFGs at redshifts from KMOS, with similar masses to our simulated galaxies. RP feedback typically lowers the stellar mass by a factor of (see also Moody et al., 2014). At low halo masses at high redshifts, the stellar mass can be decreased by up to a factor of , suggesting that radiative feedback is more efficient in low mass galaxies.

At redshifts , the RPs are still a factor of above the Behroozi, Wechsler & Conroy (2013) curve, though they appear consistent with the estimates of Burkert et al. (2015), at least for halos with , where most of our data lies. This reflects a factor of uncertainty in the observational estimates of the stellar-to-halo mass ratio for these halo masses and redshifts, and our simulations produce values within the observational uncertainties. At , most of our sample lies at low halo masses where no direct comparison with existing predictions can be made. By extrapolating the Behroozi, Wechsler & Conroy (2013) relation towards lower halo masses, the off-set between our sample and the prediction is comparable to that at . We conclude that the RPs produce stellar masses that are in the ballpark of realistic values, within the observational uncertainties of a factor .

### 4.2 Disc Properties and Toomre Q

In Fig. 4 we compare the galaxy samples in terms of global disc structure and stability. We examine as a function of time, represented by the expansion factor , the disc radius, thickness, mass fraction with and Toomre parameter, as described below. For each of these properties we show the median and scatter among all snapshots at a given expansion factor.

In the top left panel we examine the disc radius normalized to the virial radius, . On average, the disc sizes are very similar in the two galaxy samples. The NoRPs have a very slight tendency to produce smaller discs than the RPs, by at all times. However, given the scatter in the data this does not appear significant. Most of the scatter is caused by small differences in the times of major mergers and compaction events between the two samples, which mostly occur after (Zolotov et al., 2015).

The top right panel examines the disc thickness, normalized to its radius . To first order, this represents the velocity dispersion within the disc, . As and were both calculated using the cold component, this refers to the gas velocity dispersion, rather than the stellar value. The stronger feedback in the RPs results in thicker and more turbulent discs, especially at later times. The median ratio between individual snapshots paired in both sets of simulations increases monotonically from at to at . However, the scatter is large, at and at . It is unclear how much of this trend is the result of instantaneous feedback driving additional turbulence, and how much is the result of decreased gas fractions in NoRPs due to excessive star-formation at early times. Surely the latter contributes to the increased discrepancy at later times. However, this distinction is beyond the scope of the current paper.

The bottom left panel examines the ratio of the disc mass to the total mass within the disc radius, including the spheroid and dark matter halo components, . This is an important indicator of disc instability (e.g. Dekel, Sari & Ceverino, 2009). At early times, the NoRPs have a more prominent disc, which results from earlier settling of the galaxies into oblate structures (Tomassetti et al., 2016). Detailed analysis of the shapes of galaxies in the RPs indicate that galaxies tend to maintain a prolate shape at early times, developing a disc only after the first compaction event, when the baryons dominate the mass and potential at the inner radii (Ceverino, Primack & Dekel, 2015). This happens earlier in the NoRP galaxies (Tomassetti et al., 2016). At late times, the two sets of simulations have comparable . It is interesting to characterise the cold component in the disc by defining instead . At , the NoRPs and the RPs have and respectively. At , both have comparable values of . By , the NoRPs have saturated at , while the RPs have saturated at . This indicates that the gas fractions in the RP discs are higher than in the NoRP discs at these redshifts.

The bottom right panel shows the ratio of the global Toomre parameter in the RPs to in the NoRPs, using the crude approximation (e.g. Dekel, Sari & Ceverino, 2009):

 Q∝σΩΣ∝σVMtotMd∝δ−1HdRd (5)

We do this once for the cold component only, using , and once for the baryonic disc mass, using (which still only includes stars that are co-rotating with the gas disc). This crudely approximates for a one-component model, assuming the gas and stars have comparable velocity dispersions. Since the discs contain dynamically hot stars as well, a two component Toomre analysis is required to fully address the issue of disc instability. However, eq. (5) is not applicable for a two component system, and a more detailed analysis is beyond the scope of this paper. Therefore, rather than focus on the absolute values of , we only address the relative values in the two simulation suites, in order to compare the predicted strength of VDI. We note that in the RPs, the two component is not very different from computed using the cold component only (Inoue et al., 2016).

The NoRPs have slightly lower values of , indicating that they may be slightly more prone to instability and clump formation than the RPs. However, this difference is not large, when considering all baryons in the disc and less than that when considering only the cold component. We conclude that RP does not fundamentally stabilize the discs against VDI, so clump formation is not expected to be greatly suppressed.

### 4.3 ISM Structure

Figure 5 shows the density probability distribution function (PDF) of the ISM gas in the NoRPs and RPs, stacked over all common snapshots. We consider all gas with temperatures located within the disc cylinder at each snapshot9. Hot gas, with , contains only of the gas mass in the ISM and is not important for disc instability or for clump detection. We show both the mass-weighted PDF, (left), and the volume weighted PDF, (right). These are the mass fraction and volume fraction per logarithmic density interval, respectively.

At low to intermediate densities, the mass-weighted PDF of the NoRPs and the RPs have a very similar shape. Both are well fit by a log-normal distribution which has a mean density of and a standard deviation of . The mean density decreases from at to at , in agreement with the cosmological scaling, while the standard deviation decreases from to in the same redshift range. We note that the adopted density threshold for star formation, , is at the peak of the mass weighted PDF. A log-normal shape for the density PDF is characteristic of isothermal, supersonic turbulence (Vazquez-Semadeni, 1994; Padoan, Nordlund & Jones, 1997; Scalo et al., 1998; Federrath, Klessen & Schmidt, 2008; Price, Federrath & Brunt, 2011; Hopkins et al., 2012) with the width of the distribution weakly depending on the Mach number as (Padoan, Nordlund & Jones, 1997; Federrath, Klessen & Schmidt, 2008).

At high densities, , the distributions deviate from log-normality and develop power law tails, with slopes of and in the RPs and NoRPs respectively. The transition density increases with redshift, maintaining a factor above the mean. The development of a power law tail at high densities is typical when self-gravity becomes important and the gas can no longer be held up by turbulence (Vázquez-Semadeni et al., 2008; Elmegreen, 2011; Hopkins et al., 2012). Radiation pressure regulates the distribution of very dense, star-forming gas. As a result, the mass fraction of dense, self-gravitating gas in the power law tail is reduced by a factor of in the RPs, consistent with the typical decrease in the amount of stars formed in the simulations (Fig. 3). The mass in the hot component, with (not shown), is roughly twice as high as in the RPs compared to the NoRPs. Similar results have been obtained in studies of isolated galaxy simulations (Hopkins et al., 2012; Rosdahl et al., 2015) and in one of the cosmological simulations used here (Ceverino et al., 2014).

Both distributions are attenuated at the highest densities of , with a slight tendency for a higher threshold in the RPs vs the NoRPs, and at lower redshifts. This is likely due to the pressure floor (eq. (1)), to which we can associate an effective velocity

 σfloor=√γPfloorρ≃120n1/24Δ25kms−1 (6)

where and .10 In the RPs, the radial velocity dispersion, , is typically before compaction and after compaction (Zolotov et al., 2015, figure 11). So at densities of , the pressure floor becomes comparable to the turbulent pressure, preventing cells from reaching higher densities.

The volume-weighted PDF is well fit by a log-normal in both NoRPs and RPs, with peak densities of and respectively, and standard deviations of and . Dense gas, with comprises less than of the total volume. The volume filling factor of star-forming gas, with , is reduced by a factor of in the RPs, which in turn reduces the star formation in the RPs, as seen in Fig. 3. The hot component, with (not shown) occupies a comparable volume to the cooler component, despite containing only of the mass. The distribution of densities is also log-normal, with mean densities of and in the NoRPs and RPs respectively, and corresponding standard deviations of and . The volume fraction of this hot gas is roughly 10 times higher in the RPs than in the NoRPs, showing that RP heats the gas and drives it to lower densities.

To summarise, we expect bound star-forming clumps to form at densities , while most of the volume in the ISM is at densities of . Our threshold density residual for clump detection, 3.3), is thus well suited to detect the vast majority of clump candidates. Naturally, given the wide range of low densities, we also expect to detect many unbound, low density, transient structures, which form the bulk of the ZLCs discussed earlier.

## 5 The Effect of Radiative Feedback on clumps

In this section we examine the effect of RP feedback on the formation, lifetime, and properties of giant clumps in high- disc galaxies. As in the previous section, we limit our analysis to the snapshots with both RP and NoRP versions, and discuss separately the effects on in situ, ex situ and bulge clumps.

### 5.1 In Situ Clumps

#### Mass-size plane and sample completeness

Figure 6 shows the distribution of in situ clumps in the mass-radius plane, for the NoRPs (left) and the RPs (right). In all panels, the color represents the median value within each pixel of a respective property. The top row shows the clump shape, , and the bottom row shows the clump time, . Cases where were artificially assigned a value of so that they appear on the plot in dark blue. The red and magenta lines represent constant baryonic volume densities of as marked. RP feedback reduces the typical density of in situ clumps by roughly , similar to the factor by which the typical ISM density (by volume) is reduced (Fig. 5, right). This is expected since we define clumps based on local overdensity.

It is clear from the distribution of pixels in each diagram that our sample of clumps is incomplete at masses , due to our minimal clump volume, . At a given mass, , we can only resolve clumps with densities . Due to the maximal ISM density we resolve, there is a maximal clump density, , so that we are complete for all masses . For , and , we can resolve all clumps with , and respectively. We estimate the incompleteness at these three masses by calculating the fraction of clumps with that have densities above these thresholds. The corresponding fractions are , and , and are very similar in both sets of simulations. We conclude that our sample is reasonably complete for clumps with . Hereafter, when we discuss the sample of clumps, we refer only to this mass complete sample.

From Fig. 6, we learn that clumps with are rather round and oblate, , and dense, . This corresponds to the scale where the ISM becomes self-gravitating, allowing clumps to collapse to bound, oblate, long-lived structures. At lower densities, turbulence prevents the clumps from collapsing, and they remain extended, typically prolate structures. Such structures are rapidly mixed back into the ISM through a combination of turbulence and shear, or are torn apart by tidal forces. The oldest clumps are typically old, comparable to the expected migration time for clumps in high redshift discs (Dekel, Sari & Ceverino, 2009). This will be discussed further in §6.2.

Figure 7 shows the distribution of clump lifetimes in the NoRPs (blue) and the RPs (red). The axis shows the average number of clumps per snapshot per bin of , shown on the axis. The clump lifetime has been normalized by the mass-weighted average free-fall time during its lifetime

 tff=√3π32Gρ≃0.54√Gρ (7)

There is a bi-modal distribution with a minimum near . This is in the ball park of the output time step in the simulations, and it can therefore separate clumps that were identified in one or multiple snapshots. We also note that is the maximal clump lifetime observed in isolated galaxy simulations with much stronger radiative feedback (Hopkins et al., 2012). We therefore use this threshold to distinguish short-lived clumps (SLCs) from long-lived clumps (LLCs). In both sets of simulations, only of clumps that can be tracked for more than one timestep fall into the SLC category, and these make up less than of the SLCs themselves.

Ignoring ZLCs, the RPs have as many clumps as the NoRPs, and as many LLCs. In the NoRPs, of clumps are LLCs compared to only in the RPs. We conclude that RP has a minor effect on clump formation, but it dramatically reduces the number of long-lived clumps.

#### Clump structural and virial properties

In Fig. 8, we compare structural and virial properties of clumps in the RPs and NoRPs, to try and discern what causes clump disruption. In each panel we show the distribution of a different property, represented on the axis, where the axis shows the average number of clumps per snapshot per bin of , as in Fig. 7. In clockwise order from the top left the properties are clump shape, ; mass, ; circular velocity, ; surface density, ; volume density, ; and virial parameter , i.e. twice the ratio of kinetic to potential energies within the clump. For each property, we show the distributions including ZLCs, the distributions for non-ZLCs, and the distributions for LLCs only, as marked in the legend.

The virial parameter is measured in a sphere of radius centered on the peak baryonic density within the clump. This is a valid approximation for nearly spherical clumps and can serve as an order of magnitude estimate for elongated clumps as well. In the rest-frame defined by the center of mass velocity of baryons within the sphere, we measure the total baryonic kinetic energy, which gives us the mass-weighted 3D velocity dispersion,