© ESO, 2014

1. Introduction

Lambda-CDM cosmological simulations tend to show that a major merger is at work shaping galaxy properties at high redshifts (Stewart et al. 2009). Although it is often set as a competitor of the smooth cold gas accretion along cosmic filaments, which is believed to be very efficient at feeding star formation (Dekel et al. 2009a; Kereš et al. 2009b), mergers still contribute to around a third of the baryonic mass assembly history (Brooks et al. 2009; Dekel et al. 2009b). The pioneering work of Toomre & Toomre (1972) first highlighted that disk galaxy mergers are able to drive large amounts of baryons in tidal tails. Mihos & Hernquist (1994) showed that the redistribution of gas can fuel star formation enhancement in the core of the remnants. Furthermore, stars during a merger event are gravitationally heated and can form spheroids (Barnes & Hernquist 1996; Mihos & Hernquist 1996), underlining a convincing link between the late-type and early-type galaxies of the Hubble sequence. Boxy slowly rotating ellipticals, however, probably formed at much higher redshifts through multiple minor mergers or in-situ star formation (Oser et al. 2010; Feldmann et al. 2010; Johansson et al. 2012).

The study of stars and gas kinematics is a good way to detect signatures of merger in the recent history of galaxies (Barnes 2002; Arribas & Colina 2003; Bois et al. 2011), and allows constraining the role of mergers in galaxy mass assembly. The past decade has seen the first resolved observations of galaxies in the redshift range 0.5 < z < 3 using integral fields unit spectrographs (IFU; Yang et al. 2008; Förster Schreiber et al. 2009; Law et al. 2009; Gnerucci et al. 2011; Contini et al. 2012), where a peak is observed in the cosmic star formation history. This peak located around z ~ 2 (Hopkins & Beacom 2006; Yang et al. 2008) could arise from intense merger activity, since it is an efficient mechanism for producing starbursts in the local Universe. Shapiro et al. (2008) did a kinematical analysis of the high-z IFU SINS sample to determine the fraction of mergers. To calibrate this analysis, a set of local observations (Chemin et al. 2006; Daigle et al. 2006; Hernandez et al. 2005), hydrodynamical cosmological simulations (Naab et al. 2007), and toy models (Förster Schreiber et al. 2006) were used. The halos in the Naab et al. (2007) cosmological simulations were selected to host a merger around z = 2. Although the baryon accretion history makes these simulations credible in terms of cosmological mass assembly, the resulting low number of halos could be considered as insufficient at statistically detecting various merger signatures.

The GalMer database (Chilingarian et al. 2010) favors a statistical approach with hundreds of idealized merger simulations, which are probing the orbital configurations. GalMer is relevant for studying such merger signatures at low redshift (Di Matteo et al. 2008); however, the low gas fractions makes the comparison with high redshift galaxies impossible. Additionally, simulating the interstellar medium (ISM) of high redshift galaxies requires correctly resolving the high redshift disk scaleheights, which can otherwise artificially prevent the expected Jeans instabilities. Indeed, it is now commonly accepted that high redshift disks are naturally subject to such instabilities (Elmegreen et al. 2009, 2007). The high gas fractions at z > 1 (Daddi et al. 2010a; Tacconi et al. 2010) are strongly suspected of driving violent instabilities that fragment the disks into large star-forming clumps (Bournaud et al. 2008) and generate turbulent velocity dispersions (e.g., Epinat et al. 2012; Tacconi et al. 2008). Therefore, the canonical image of smooth extended tidal tails falling onto the merger remnant cannot be valid in the context of gas-rich interactions (Bournaud et al. 2011).

The ability to form such clumps is important to understand the complex behavior of high redshift galaxies, but it is also essential to prevent the overconsumption of gas expected at these very high gas densities from the classical Schmidt law. To match the Kennicutt-Schmidt (KS) relation (Kennicutt 1998) and to have acceptable gas consumption timescales, an efficient stellar feedback is required to deplete the gas reservoir of the star-forming regions. Indeed, cosmological simulations with weak or no feedback models produce galaxies with too many baryons in the galactic plane (Kereš et al. 2009a) when compared to the abundance-matching techniques (Guo et al. 2010). The constraints on the intergalactic medium (IGM) metal enrichment also imply that baryons entered galaxies at some points and underwent star formation (Aguirre et al. 2001). It has been demonstrated that scaling supernovae stellar winds in cosmological simulations to the inverse of the mass of the host galaxy produces models in reasonable agreement with the local mass function (Oppenheimer et al. 2010). It is therefore essential to constrain the parameters controlling the stellar feedback processes in order to better understand the scenarios of galaxy evolution.

To get insight into the various processes of galaxy mass assembly, such as mergers, the Mass Assembly Survey with SINFONI in VVDS (MASSIV, Contini et al. 2012) aims to probe the kinematical and chemical properties of a significant and representative sample of high redshift (0.9 < z < 1.8), star-forming galaxies. Observed with the SINFONI integral-field spectrograph at the VLT and built upon a simple selection function, the MASSIV sample provides a set of 84 representatives of normal star-forming galaxies with star formation rates ranging from 5 to 400 M_⊙ yr^-1 in the stellar mass regime 10⁹ − 10¹¹ M_⊙. Compared to other existing high-z IFU surveys, the main advantages of the MASSIV sample are its representativeness since it is flux-selected from the magnitude-limited VVDS survey (Le Fèvre et al. 2005) and its size, which allows probing different mass and star formation rate (SFR) ranges, while keeping enough statistics in each category. Together with the size of the sample, the spatially-resolved data therefore allows galaxy kinematics and chemical properties to be discussed across the full mass and SFR ranges of the survey to derive robust conclusions for galaxy mass assembly on cosmological timescales. By studying strong kinematic signatures of merging and detecting pairs in the first-epoch MASSIV, Epinat et al. (2012) have shown that the fraction of interacting galaxies is up to at least one third of the sample and that more than a third of the galaxies are non-rotating objects. In addition, there are more non-rotating objects in mergers than in isolated galaxies. This suggests that a significant number of isolated non-rotating objects could be mergers in a transient state in which the gas is not dynamically stable. Furthermore, based on the whole MASSIV sample, López-Sanjuan et al. (2013) find a gas-rich major merger fraction of ~20% in the redshift range 1 < z < 1.8 and a gas-rich major merger rate of ~0.12. Quantification of the kinematical signatures of interacting galaxies and mergers and the understanding of the high fraction of non-rotating systems, the existence of inverse metallicity gradient in some disks (Queyrel et al. 2012), and more generally, a comprehensive view of the process of formation of turbulent and clumpy gaseous galaxy disks, has motivated building a set of simulations of merging galaxies in the redshift range probed with MASSIV, i.e. the Merging and Isolated high redshift Adaptive mesh refinement Galaxies (MIRAGE) simulations.

We describe, a set of 20 idealized galaxy mergers and three isolated disks using adaptive mesh refinement (AMR) simulations with a physically motivated implementation of stellar feedback¹. This paper focuses on presenting the MIRAGE sample, the numerical technique employed, and the physical properties deduced. The analysis is extended in a companion paper Bournaud et al. (2014) that presents a study of the clumps properties in the three isolated disk simulations of MIRAGE sample. The paper is organized as follows. In Sect. 2, we describe the numerical technique used to build our simulation sample. In Sect. 3, we specifically describe the idealized initial conditions generation. For this purpose we introduce the new public code DICE and summarize the different numerical techniques used to generate stable galaxies models. Section 4 reviews the MIRAGE sample definition of galactic models and orbital parameters. Section 5 describes the global properties of the sample. For each simulation, we present the star formation histories, the disk scalelengths evolution, and their position on the KS relation.

2. Simulations

We ran a set of idealized AMR high redshift galaxy simulations. The sample encompasses 20 major and minor galaxy mergers and three isolated disks, with a high gas fraction (>50%) typical of 1 < z < 2 galaxies (Daddi et al. 2010a), evolved over 800 Myr. In this work, we choose to balance the available computational time between high-resolution and statistical sampling of the orbital parameters to provide new insight into the galactic mass assembly paradigm.

2.1. Numerical technique

To build our numerical merger sample, we use the AMR code RAMSES (Teyssier 2002). The time integration of the dark matter and the stellar component is performed using a particle-mesh (PM) solver, while the gas component evolution is insured by a second-order Godunov integration scheme. The code has proven its ability to model the complexity of interstellar gas on various galaxies simulations (e.g., Dubois & Teyssier 2008; Teyssier et al. 2013). The computational domain of our simulations is a cube with a side l_box = 240 kpc, and the coarsest level of the AMR grid is ℓ_min = 7, which corresponds to a Cartesian grid with (2⁷)³ elements and with a cell size of Δx = 1.88 kpc. The finest AMR cells reach the level ℓ_max = 15, where the cell size corresponds to Δx = 7.3 pc. The grid resolution is adapted at each coarse time step between the low refinement (ℓ_min = 7) and the high refinement levels (ℓ_max = 15). Each AMR cell is divided into eight new cells if at least one of the following assertions is true: (i) it contains a gas mass greater than 1.5 × 10⁴ M_⊙; (ii) it contains more than 25 particles (dark matter or stars); or (iii) the local Jeans length is less than four times the current cell size. This quasi-Lagrangian refinement scheme is comparable to the one introduced in Teyssier et al. (2010) and Bournaud et al. (2010).

The star formation is modeled with a Schmidt law triggered when the density ρ_gas overcomes the threshold ρ₀ = 100 cm^-3, with an efficiency ϵ_⋆ = 1%: $\mathit{\rho ̇}\mathit{\star }\mathrm{=}\{\begin{array}{c}\\ \mathrm{0}& \mathrm{if}& {\mathit{\rho }}_{\mathrm{gas}}\mathit{<}{\mathit{\rho }}_{\mathrm{0}}\\ \mathrm{0}& \mathrm{if}& \mathit{T}\mathit{>}\mathrm{2}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\mathrm{K}\\ {\mathit{ϵ}}_{\mathit{\star }}{\mathit{\rho }}_{\mathrm{gas}}\mathit{/}{\mathit{t}}_{\mathrm{ff}}& \mathrm{else}\mathit{,}\end{array}$(1)where is the local SFR, ${\mathit{t}}_{\mathrm{ff}}\mathrm{=}\sqrt{\mathrm{3}\mathit{\pi }\mathit{/}\mathrm{\left(}\mathrm{32}\mathit{G}{\mathit{\rho }}_{\mathrm{gas}}\mathrm{\right)}}$ is the free-fall time computed at the gas density ρ_gas, and T is the temperature of the cell considered. AMR cells with temperature greater than 2 × 10⁵ K are not allowed to form stars.

The gravitational potential is computed using a PM scheme with a maximum level ℓ_max,part = 13 for the grid, which ensures gravitational softening of at least 29 pc for Lagrangian particles. This choice prevents a low number of dark matter particles per cell, often synonymous with N-body relaxation. We use a thermodynamical model modeling gas cooling provided by the detailed balance between atomic fine structure cooling and UV radiation heating from a standard cosmic radiation background by using tabulated cooling and heating rates from Courty & Alimi (2004). In this model, the gas metallicity acts like a scale factor on the cooling rate.

The gas is forced to stay within a specific area in the density-temperature diagram to prevent multiple numerical artifacts (see Fig. 1):

• In the low IGM density regime (ρ < 10^-3 cm^-3), we ensure a gravo-thermal equilibrium for the gas by introducing a temperature floor in the halo following a gamma polytrope at the virial temperature T_min(ρ) = 4 × 10⁶(ρ/10^-3)^2/3 K, as in Bournaud et al. (2010).

• For densities between 10^-3 cm^-3 < ρ < 0.3 cm^-3, the temperature floor is isothermal and set to T_min(ρ) = T_floor. At full resolution (i.e. ℓ_max = 15), we have T_floor = 300 K. The densest IGM can reach the ρ = 10^-3cm^-3 limit and can condense on the gaseous disk.

• For densities above 0.3 cm^-3, we use the temperature floor T_min(ρ) = 300 × (ρ/0.3)^−1/2 K; this choice allows us to have a dynamical range in the thermal treatment of the gas up to 30 times colder than the slope of the thermodynamical model used in Teyssier et al. (2010) and Bournaud et al. (2010).

• A density-dependent pressure floor is implemented to ensure that the local Jeans length is resolved at least by n_Jeans = 6 cells in order to avoid numerical fragmentation, as initially proposed by Truelove et al. (1997). This Jeans polytrope acts like a temperature floor for the very dense gas: it overcomes the cooling regime of the temperature floor starting from ρ = 2.6 cm^-3 when the resolution is maximum, i.e. cells with a size of 7.3 pc. The Jeans polytrope is described by the equation ${\mathit{T}}_{\min }\mathrm{\left(}\mathit{\rho }\mathrm{\right)}\mathrm{=}\mathit{\rho G}{\mathit{m}}_{\mathrm{H}}\mathrm{(}{\mathit{l}}_{\mathrm{box}}{\mathit{n}}_{\mathrm{Jeans}}\mathit{/}{\mathrm{2}}^{{\mathit{\ell }}_{\max }}{\mathrm{)}}^{\mathrm{2}}\mathit{/}\mathrm{(}\mathit{\pi }{\mathit{k}}_{\mathrm{B}}\sqrt{\mathrm{32}}\mathrm{)}$, with m_H the proton mass and k_B the Boltzmann constant.

• We impose a maximum temperature for the gas T_max = 10⁷ K. Indeed, the clumps generated by Jeans instabilities typical of gas-rich disks (Bournaud et al. 2008) lead to regions of low density inside the disk, where supernovae can explode. This thermal explosion is thus able to produce sound speed greater than 1000 km s^-1, which affects the time step in the Godunov solver. Setting an upper limit to the temperature is not fully conservative in terms of energy, but our choice of T_max ensures a viable time step and a reasonably low energy loss. This issue typical of grid codes is handled in the same way in the recent work of Hayward et al. (2013).

Owing to non-periodic boundary conditions applied to the AMR box, we impose a zero density gradient in the hydrodynamical solver at the boundaries. To avoid galaxies passing close to the edges of the box, which could induce numerical artifacts, this gradient is required to have sufficiently large AMR volume to encompass the whole trajectories of both galaxies through the simulation duration.

2.2. Feedback models

Because we do not resolve individual stars, each stellar particle models a population that contains massive OB-type stars with masses M > 4 M_⊙ (Povich 2012) and which is responsible for injecting energy into the surrounding ISM. Assuming a Salpeter (1955) initial mass function, we consider that a fraction η = 20% of the mass of a stellar particle contributes to stellar feedback, which is effective during 10 Myr after the star particle is spawned. We use the Renaud et al. (2013) physically-motivated model implementation for the OB-type stars feedback, and summarize the three main recipes below.

• Photo-ionization: OB-type stars produce highly energeticphotons capable of ionizing the surrounding ISM. Using a simplemodel for the luminosity of the star, the radius of theStrömgren (1939) sphere is computed accordingto the mean density of electrons n_e. The gas temperature inside the HII regions is replaced by an isothermal branch at T_HII = 10⁴ K. The radius of the HII sphere is computed via the equation: ${\mathit{R}}_{\mathrm{HII}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}\mathit{\pi }}\frac{{\mathit{L}}_{\mathrm{\ast }}}{{\mathit{n}}_{\mathrm{e}}^{\mathrm{2}}{\mathit{\alpha }}_{\mathrm{r}}}\mathit{,}$(2)where L_∗ is the time-dependent luminosity of the star in terms of ionizing photons, and α_r is the effective recombination rate.

• Radiative pressure: Inside each HII bubble, a kinetic momentum Δv is distributed as a radial velocity kick over the time interval Δt, matching the time step of the coarsest level of the simulation. This velocity kick is computed using ionizing photons momentum, which is considered to be transferred to the gas being ionized, i.e. the gas within the radius R_HII of the Strömgren sphere: $\mathrm{\Delta }\mathit{v}\mathrm{=}\mathit{k}\frac{{\mathit{L}}_{\mathrm{\ast }}\mathit{h\nu }}{{\mathit{M}}_{\mathrm{HII}}{\mathit{v}}_{\mathrm{c}}}\mathrm{\Delta }\mathit{t,}$(3)with h the Planck constant, v_c the speed of light, M_HII the gas mass of the bubble affected by the kick, and ν the frequency of the flux representative of the most energetic part of the spectrum of the ionizing source. In this model, it is considered that Lyman-α photons dominate this spectrum, implying ν = 2.45 × 10¹⁵ s^-1. The distribution of the momentum carried by ionizing photons is modeled by the trapping parameter k = 5, which basically counts the number of diffusion per ionizing photon and energy loss. This value may appear to be rather high compared to recent works (e.g., Krumholz & Thompson 2012), but is more acceptable once considered that we miss other sources of momentum such as protostellar jets and stellar winds (Dekel & Krumholz 2013).

• Supernova explosions: we follow the implementation of supernova feedback of Dubois & Teyssier (2008): the OB-type star population that reaches 10 Myr transforms into supernovae (SNe) and releases energy, mass, and metals into the nearest gas cell. The gas that surrounds the supernovae receives a fraction η = 20% of the stellar particle mass, as well as a specific energy E_SN = 2 × 10⁵¹ ergs/10 M_⊙, which is the product of the thermo-nuclear reactions. The energy injected by each SN is higher by a factor of two compared to some other works (e.g., Teyssier et al. 2013; Dubois et al. 2012), but simulations of individual type II SN releasing such energy could be frequent in the early Universe (Joggerst et al. 2010). Moreover, the use of an IMF with a lower statistical contribution of low mass stars would imply higher values for η (e.g. η ≃ 35% for Kroupa 2001 IMF), which balances our choice of a high value for E_SN. Each supernovae event also releases into the surrounding ISM metals derived form the nucleosynthesis following the equation $\mathit{Z}\mathrm{=}\mathit{y}\mathrm{+}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{y}\mathrm{)}{\mathit{Z}}_{\mathrm{ini}}\mathit{,}$(4)with Z the mass fraction of metals in the gas, Z_ini the initial metal fraction of the supernova host, and y the yield that is set to y = 0.1, as in Dubois et al. (2012). To account for non-thermal processes due to gas turbulence on subparsec scales, we follow the revised feedback prescriptions of Teyssier et al. (2013). The numerical implementation is similar to introducing a delayed cooling in the Sedov blast wave solution. At each coarse time step, the fraction of gas released by SNe in AMR cells is evaluated in a passively advected scalar; the gas metal line cooling is switched off as long as $\frac{{\mathit{m}}_{\mathrm{ejecta}}}{\mathit{m}}\mathit{>}{\mathrm{10}}^{-3}\mathit{,}$(5)with m the gas mass of the cell, and m_ejecta the total mass of the gas ejected by SNe in the same cell. To model the turbulence dissipation, the mass of gas contributing to the Sedov blast wave is lowered by a factor γ at each coarse time step following $\mathit{\gamma }\mathrm{=}\exp \left(\mathrm{-}\frac{\mathrm{d}{\mathit{t}}_{\mathrm{cool}}}{{\mathit{t}}_{\mathrm{dissip}}}\right)\mathit{,}$(6)with dt_cool the cooling time step, and t_dissip a typical timescale for the turbulence induced by the detonation. The dissipation timescale for the unresolved subgrid turbulent structures is the crossing time (Mac Low 1999), i.e the ratio of the numerical resolution over the velocity dispersion. Since our simulations are able to resolve structures down to 7 pc, we presume that the non-thermal velocity dispersion in the smallest AMR cells is close to 5 km s^-1, which is a typical supersonic speed in regions of star formation with gas temperatures below 10³ K (Hennebelle & Falgarone 2012). Under these assumptions, we set t_dissip = 2 Myr.

Our feedback model does not assume the systematic destruction of the clumps by the star formation bursts following their formation, unlike what is done in some other works (e.g., Hopkins et al. 2013; Genel et al. 2012). Smaller clumps are subject to disruption, but larger clumps may survive such thermal energy injection. This model clearly favors the scenario of long-lived star-forming clumps, which we aim to address in this study.

3. DICE: a new environment for building disk initial conditions

The initial conditions of the MIRAGE sample are constructed using software developed for the purpose of the task, named disk initial conditions environment (DICE). DICE is an implementation of the numerical methods described in Springel et al. (2005a). It is able to set up multiple idealized galaxies in a user friendly context. The software is open source and available online².

3.1. Density distributions

DICE initial conditions are generated using Lagragian particles whose distributions are built using a Metropolis-Hasting Monte-Carlo Markov Chain algorithm (Metropolis et al. 1953). The strength of this algorithm lies in its ability to build a distribution for a sample of Lagrangian particles having only the knowledge of the probability distribution function. After having initialized the first Lagrangian particles of each component (disk, bulge, gas, halo, etc.) to a probable location, the algorithm loops over the desired number of Lagrangian particles and iteratively produces a candidate position for each of them. The probability of setting a Lagrangian particle to the randomly picked candidate Cartesian position x′ depends on the Cartesian position x of the previous particle in the loop, and is written as $\mathit{\alpha }\mathrm{\left(}{x}\mathit{,}{{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\begin{array}{c}\min \end{array}\left(\mathrm{1}\mathit{,}\frac{\mathit{\rho }\mathrm{\left(}{x}\mathrm{\right)}\mathit{P}\mathrm{\left(}{x}\mathit{,}{{x}}^{\mathrm{\prime }}\mathrm{\right)}}{\mathit{\rho }\mathrm{\left(}{{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathit{P}\mathrm{\left(}{{x}}^{\mathrm{\prime }}\mathit{,}{x}\mathrm{\right)}}\right)\mathit{,}$(7)with ρ(x) the density function of the considered component at the position x, and P(x,x′) the probability of placing the particle at x′ considering the position x of the previous particle. Indeed, our implementation uses a Gaussian walk, meaning that the candidate coordinates are generated using the rule ${{x}}^{\mathrm{\prime }}\mathrm{=}{x}\mathrm{+}\mathit{\sigma W,}$(8)where W is a standard Gaussian random variable, and σ a dispersion factor tuned to a fixed fraction of the targeted scalelength of the component to build, ensuring satisfying convergence. For each particle, a uniform random value τ ∈  [0,1] is picked, and the position of the Lagrangian particle is set to x′ if τ ≤ α and to x, otherwise. The first 5% of the iterations to build the distribution are not taken into account because they are considered as a “burning period” to account for any eventual poor choice of initial values.

To fit the system in the finite AMR domain, we cut the density profiles of all the components. We apply these cuts using an exponential truncation profile at the edges of each component, in order to prevent strong discontinuities nearly the cut region, which would make the numerical differentiation quite unstable. The scalelength of the exponential truncation profile is set to be one percent of the gas disk scaleheight.

3.2. Gravitational potential

To set up the velocities in our initial conditions, we compute the gravitational potential using a PM technique. We first interpolate the densities of all the components onto a Cartesian grid using a cloud-in-cell scheme. We compute the gravitational potential Φ by solving the Poisson equation: $\mathrm{\Phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\mathrm{\int }\mathrm{𝒢}\mathrm{\left(}{x}\mathit{,}{{x}}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{4}\mathit{\pi \rho }\mathrm{\left(}{{x}}^{\mathrm{\prime }}\mathrm{\right)}{\mathrm{d}}^{\mathrm{3}}{{x}}^{\mathrm{\prime }}\mathit{,}$(9)where is the Green function, and ρ is the density function of all the mass components interpolated on the Cartesian grid. We compute this integral by performing a simple product on the Fourier plane, which is equivalent to a convolution in the real plane. We eliminate the periodicity associated to the fast Fourier transform algorithm using the zero-padding technique described in Hockney & Eastwood (1988).

3.3. Velocities

To fully describe our system, we assume that the mean radial and vertical velocities ⟨v_r⟩ and ⟨v_z⟩ are equal to zero. The velocities of each Lagrangian particle are determined by integrating the Jeans equations (Binney et al. 2009), assuming that the velocity distribution is shaped as a tri-axial Gaussian. For the dark matter halo and the stellar bulge, we numerically solve the equations $\begin{array}{ccc}\mathrm{⟨}{\mathit{v}}_{\mathit{z}}^{\mathrm{2}}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{\mathit{\rho }}{\mathrm{\int }}_{\mathit{z}}^{\mathrm{\infty }}\mathit{\rho }\mathrm{\left(}\mathit{r,}{\mathit{z}}^{\mathrm{\prime }}\mathrm{\right)}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial }{\mathit{z}}^{\mathrm{\prime }}}\mathrm{d}{\mathit{z}}^{\mathrm{\prime }}\mathit{,}\\ \mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}& \mathrm{=}& \mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}\mathrm{⟩}\mathrm{+}\frac{\mathit{r}}{\mathit{\rho }}\frac{\mathit{\partial }\mathrm{(}\mathit{\rho }\mathrm{⟨}{\mathit{v}}_{\mathit{r}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{)}}}{\mathit{\partial r}}\mathrm{+}\mathit{r}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\mathrm{·}\end{array}$The velocity dispersion can thus be computed using the relation ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{=}\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{⟩}\mathrm{-}\mathrm{⟨}{\mathit{v}}_{\mathit{\phi }}{\mathrm{⟩}}^{\mathrm{2}}\mathit{.}$(12)The dark matter halo is generally described with an angular momentum that is not specified by the Jeans equations. The streaming component is set to be a low fraction f_s of the circular velocity: i.e., ⟨v_φ⟩  = f_sv_c. The fraction f_s depends on the halo spin parameter λ and the halo concentration parameter c (Springel & White 1999), which are used as input parameters in our implementation.

For the stellar disk, we choose to use the axisymmetric drift approximation (Binney et al. 2009), which allows fast computation, although we caution against the risk of using this approximation with thick and dispersion supported disks³. This approximation relates the radial Gaussian dispersion to the azimuthal one: ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathit{r}}^{\mathrm{2}}}{{\mathit{\eta }}^{\mathrm{2}}}\mathit{,}$(13)with ${\mathit{\eta }}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}}{\mathit{r}}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}{\left(\frac{\mathrm{3}}{\mathit{r}}\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\mathrm{+}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathrm{\Phi }}{\mathit{\partial }{\mathit{r}}^{\mathrm{2}}}\right)}^{-1}\mathrm{·}$(14)The Toomre parameter for the stellar disk is written as ${\mathit{Q}}_{\mathrm{stars}}\mathrm{=}\frac{{\mathit{\sigma }}_{\mathit{z}}\mathit{\kappa }}{\mathrm{3.36}\mathit{G}{\mathrm{\Sigma }}_{\mathrm{stars}}}\mathit{,}$(15)where κ is the so-called epicyclic frequency, and Σ_stars is the surface density of the stellar disk. It is used to control the stability of the stellar disk by setting a minimum value for the velocity dispersion σ_z which prevents the local Toomre parameter from going below a given limit of 1.5 in the initial conditions of our simulations, although this parametrization cannot prevent the natural fragmentation of the gaseous disk at later stages.

The only component to specify for the gas is the azimuthal streaming velocity, derived from the Euler equation: $\mathrm{⟨}{\mathit{v}}_{\mathit{\phi ,}\mathrm{gas}}\mathrm{⟩}\mathrm{=}\mathit{r}\left(\frac{\mathit{\partial }\mathrm{\Phi }}{\mathit{\partial r}}\mathrm{+}\frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{gas}}}\frac{\mathit{\partial P}}{\mathit{\partial r}}\right)\mathit{,}$(16)where P is the gas pressure.

3.4. Keplerian trajectories

DICE is also able to set up the Keplerian trajectories of two galaxies involved in an encounter. Using the reduced particle approach, we can setup the position of the two galaxies with only three input parameters: (i) the initial distance between the two galaxies r_ini; (ii) the pericentral distance, r_peri i.e. the distance between the two galaxies when they reach the periapsis of the Keplerian orbit; (iii) the eccentricity of the trajectories, which are equal for both of the galaxies. The position of the barycenter of each galaxy in the orbital plane can be expressed in Cartesian coordinates as $\begin{array}{c}\\ {\mathit{x}}_{\mathrm{1}}\mathrm{=}{\mathit{r}}_{\mathrm{1}}\cos \mathrm{\left(}{\mathit{\psi }}_{\mathrm{1}}\mathrm{\right)}\mathit{,}& {\mathit{y}}_{\mathrm{1}}\mathrm{=}{\mathit{r}}_{\mathrm{1}}\sin \mathrm{\left(}{\mathit{\psi }}_{\mathrm{1}}\mathrm{\right)}\mathit{,}\\ \\ {\mathit{x}}_{\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{r}}_{\mathrm{2}}\cos \mathrm{\left(}{\mathit{\psi }}_{\mathrm{2}}\mathrm{\right)}\mathit{,}& {\mathit{y}}_{\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{r}}_{\mathrm{2}}\sin \mathrm{\left(}{\mathit{\psi }}_{\mathrm{2}}\mathrm{\right)}\mathit{,}\end{array}$(17)with ψ₁ and ψ₂ the true anomaly of the first and second galaxy, respectively. The Cartesian velocities v_x, v_y of the two galaxies in the orbital plane are computed using $\begin{array}{c}\\ {\mathit{v}}_{\mathit{x,}\mathrm{1}}\mathrm{=}{\mathit{k}}_{\mathrm{1}}\sqrt{\frac{\mathit{\gamma }}{\mathrm{ℒ}}}\sin \mathrm{\left(}{\mathit{\psi }}_{\mathrm{1}}\mathrm{\right)}\mathit{,}& {\mathit{v}}_{\mathit{y,}\mathrm{1}}\mathrm{=}\mathrm{-}{\mathit{k}}_{\mathrm{1}}\sqrt{\frac{\mathit{\gamma }}{\mathrm{ℒ}}}\left[\mathit{e}\mathrm{+}\cos \mathrm{\left(}{\mathit{\psi }}_{\mathrm{1}}\mathrm{\right)}\right]\mathit{,}\\ \\ {\mathit{v}}_{\mathit{x,}\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{k}}_{\mathrm{2}}\sqrt{\frac{\mathit{\gamma }}{\mathrm{ℒ}}}\sin \mathrm{\left(}{\mathit{\psi }}_{\mathrm{2}}\mathrm{\right)}\mathit{,}& {\mathit{v}}_{\mathit{y,}\mathrm{2}}\mathrm{=}{\mathit{k}}_{\mathrm{2}}\sqrt{\frac{\mathit{\gamma }}{\mathrm{ℒ}}}\left[\mathit{e}\mathrm{+}\cos \mathrm{\left(}{\mathit{\psi }}_{\mathrm{2}}\mathrm{\right)}\right]\mathit{,}\end{array}$(18)with k_i the mass fraction of the i-galaxy compared to the total mass of the system, γ the standard gravitational parameter, ℒ the semi-latus rectum of the reduced particle of the system, and e the eccentricity of the orbits. With these definitions, it is possible to set trajectories for any eccentricity. This parametrization holds for point-mass particles, while galaxies are extended objects that undergo dynamical friction. The galaxies quickly deviates from their initial trajectories because of the transfer of orbital energy towards the energy of each galaxy, which can lead to coalescence.

4. Sample definition

4.1. Galaxy models

The different parameters of our disk initial conditions are summarized in Table 1. We set up three idealized galaxy models based on the MASSIV sample stellar mass histogram (Contini et al. 2012). The choice of the initial stellar masses of our simulations was made in order to sample this histogram with all the available snapshots, i.e. in the redshift range 1 < z < 2. We chose to build our sample out of three disk models with the respective stellar masses: log (M_⋆/M_⊙) = 9.8 for our low mass disk, log (M_⋆/M_⊙) = 10.2 for our intermediate mass disk, and log (M_⋆/M_⊙) = 10.6 for our high mass disk. All of our models have a stellar disk and a gaseous disk with an initial gas fraction f_g = 65%. The stellar density profile is written as ${\mathit{\rho }}_{\mathrm{stars}}\mathrm{\left(}\mathit{r,z}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{stars}}}{\mathrm{2}\mathit{\pi }{\mathit{h}}_{\mathrm{stars}}^{\mathrm{2}}}\exp \left(\mathrm{-}\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{stars}}}\right)\exp \left(\mathrm{-}\frac{\mathit{z}}{{\mathit{h}}_{\mathrm{stars}}}\right)\mathit{,}$(19)with r_stars the scalelength of the stellar disk, h_stars the scaleheight of the stellar disk, and M_stars is the uncut stellar disk mass. We use the exact same exponential profile to set up the gaseous disk, with scalelengths 1.68 times shorter than the stellar counterpart as measured in the MASSIV sample data (Vergani et al. 2012). We initialize the metallicity in the gas cells modeling the ISM of the disks following an exponential profile to be consistent with the previous prescriptions: $\mathit{Z}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{Z}}_{\mathrm{core}}\exp \left(\mathrm{-}\frac{\mathit{r}}{{\mathit{r}}_{\mathrm{metal}}}\right)\mathrm{·}$(20)We choose to have negative initial metallicity gradients, with values of r_metal equal to the gaseous disk scalelength. The fraction of metals in the center Z(r = 0) = Z_core of each model is chosen to follow the mass-metallicity relation at z = 2 found in Erb et al. (2006). Such a choice combined with the exponential profile provides integrated metallicities that are 50 percent lower than the mass-metallicity relation at z = 2 for starburst galaxies, but this choice is consistent with our aim of modeling normal star-forming galaxies. The numerical implementation of metallicity treatment of the stellar particles ignores the stars present in the initial conditions. It is therefore not required to set a metallicity profile for these stars.

Dark matter halos were modeled using a Hernquist (1990) profile, with a spin parameter set close to the conservative value with λ = 0.05 (Warren et al. 1992; Mo et al. 1998): $\begin{array}{ccc}& & {\mathit{\rho }}_{\mathrm{halo}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{{\mathit{M}}_{\mathrm{halo}}}{\mathrm{2}\mathit{\pi }}\frac{\mathit{a}}{\mathit{r}\mathrm{(}\mathit{r}\mathrm{+}\mathit{a}{\mathrm{)}}^{\mathrm{3}}}\mathit{,}\\ & & \mathit{a}\mathrm{=}{\mathit{r}}_{\mathrm{halo}}\sqrt{\mathrm{2}\left(\ln \mathrm{(}\mathrm{1}\mathrm{+}\mathit{c}\mathrm{)}\mathrm{-}\frac{\mathit{c}}{\mathrm{1}\mathrm{+}\mathit{c}}\right)}\mathit{,}\end{array}$where M_halo is the total dark matter mass, a is the halo scalelength and r_halo the scalelength for an equivalent Navarro et al. (1997) halo with the same dark matter mass within r₂₀₀ (Springel et al. 2005b). We can therefore define our halo with the frequently used concentration parameter c, which is set to a value c = 5 as measured at z ~ 2 in N-body cosmological simulations (Bullock et al. 2001). We do not consider the mass dependence of the halo concentration function to ensure that our simulations are comparable in terms of disk instability between each other.

Finally, a bulge enclosing 8% of the total initial stellar mass is modeled again using a Hernquist profile, with a scalelength set to be equal to 20% of the stellar disk scalelength.

4.2. Orbital parameters

The MIRAGE sample is designed to constrain the kinematical signatures induced by a galaxy merger on rotating gas-rich disks. To this purpose, we built a sample that explores probable disks orientations that are likely to produce a wide range of merger kinematical signatures. It has been statistically demonstrated using dark matter cosmological simulations that the spin vectors of the dark matter halos are not correlated one to the other when considering two progenitors as Keplerian particles (Khochfar & Burkert 2006). We use this result to assume that no spin orientation configuration is statistically favored. Our galaxy models are placed on Keplerian orbits using θ₁ the angle between the spin vector of the first galaxy and the orbital plane, θ₂ the angle between the spin vector of the second galaxy and the orbital plane, and κ the angle between the spin vector of the first galaxy and the second one (see Fig. 2). If these angles are uncorrelated, the normalized spin vectors are distributed uniformly over the surface of a sphere. Consequently, all the spin orientations are equally probable. If one considers a random sampling of these disk orientations using a small finite solid angle, having the spin vector coplanar to the orbital plane produces the most configurations. Therefore, we favor configurations where we have at least one spin vector in the orbital plane, i.e. θ₁ = 90° in all the configurations. We specifically avoid configurations where both of the disks are in the orbital plane because they are highly unlikely and are subject to strong resonances that are not statistically relevant. We assume that the fourth angle ω that orients the first galaxy with respect to its line of node (Toomre & Toomre 1972) might not affect the kinematics and the shape of the merger remnant since this parameter does not affect the total angular momentum of the system. Consequently, we arbitrarily chose to have the spin vector of the first galaxy always collinear to its Keplerian particle velocity vector. We defined each orbit name with an identifier referring to the angles θ₁, θ₂, and κ (see Table 2).

The choice of studying a wide range of spin vector orientations was motivated by the requirement of detecting extreme signatures and binding the kinematical and morphological parameters of the merger remnants. However, we introduced a random angle δ when setting up the spin vector of our galaxies, picked using a uniform distribution introducing a ± 5° uncertainty. This method was implemented to prevent alignment with the AMR grid, which could produce spurious effects. The slight misalignment also increases the numerical diffusion typical of grid codes, which in our case can help relax our initial conditions.

The pericenter distance, i.e. the distance between the two galactic centers at the time of the closest approach along the Keplerian trajectory, is chosen to be r_peri = r_1,cut,gas + r_2,cut,gas, where r_1,cut,gas and r_2,cut,gas are the cut radii of the first and the second galaxies, respectively. This parametrization ensures that disks do not collide on the first pericentral passage, even if the disks are already fragmented at the pericentral time. This choice is supported by the argument that low pericentral distances are not statistically relevant because the collision cross section is directly proportional to the value of this parameter; i.e., low pericentral distance is less probable. By specifying one value for the specific orbital energy of the system, one can compute the eccentricities of the orbits. To better understand the effects of the interaction parameters on the kinematics of our merger remnants, we set this parameter to a fixed negative value: ${\mathit{E}}^{\mathrm{\ast }}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{ini}}^{\mathrm{2}}}{\mathrm{2}}\mathrm{-}\frac{\mathit{G}\mathrm{(}{\mathit{m}}_{\mathrm{1}}\mathrm{+}{\mathit{m}}_{\mathrm{2}}\mathrm{)}}{{\mathit{r}}_{\mathrm{ini}}}\mathrm{=}\mathrm{-}\mathrm{2.85}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{k}{\mathrm{m}}^{\mathrm{2}}\hspace{0.17em}{\mathrm{s}}^{-2}\mathit{,}$(23)where E^∗ is the specific orbital energy, v_ini the initial relative velocity of the galaxies, and r_ini the initial distance between the galaxies. This negative specific orbital energy means that all of our trajectories are elliptic (e < 1). The parameters of the Keplerian orbits are listed in Table 3. We acknowledge that such low eccentricities might not be statistically relevant (Khochfar & Burkert 2006), but we do save computational time.

Finally, we define the initial distance between galaxies with a conservative expression through different merger masses, using the pericenter time, i.e., the time for the galaxies to reach the pericenter with t_peri = 250 Myr. The choice of the specific orbital energy is achieved in order to be able to set t_peri to 250 Myr with elliptic orbits e < 1. Because the dynamical times of all the models are close, this formulation ensures that the models relax synchronously before the start of the interaction (see Sect. 5.1). Our sample encompasses 20 merger configurations (four sets of orbital angles, five sets of orbital parameters due to different galaxy masses), to which we add the three isolated disk models to have a reference for secular evolution. We have excluded the G3_G3 interaction to save computation time, since the relative resolution on the merger remnant is coarser than any other cases.

4.3. Environment

We aim to simulate the accretion from an idealized hot gaseous halo surrounding the galactic disks. To this purpose, we model the intergalactic medium (IGM) by setting an initial minimum gas density ρ_IGM = 2.3 × 10^-4  cm^-3 within the AMR box. The gas present in the IGM is initialized with no velocity, so that it collapses towards the central potential well at the free-fall velocity. After a dynamical time, the gas halo reaches a state close to a spherical hydrostatic equilibrium where the densest regions are allowed to cool down. The zero gradient condition imposed in the grid boundaries implies a continuous injection of pristine gas on the boundaries of the AMR box.

5. Global evolution of physical properties

Figure 3 shows the morphology of the gas and the stars after 400 Myr of evolution along two orthogonal line-of-sight (LOS) for the simulation G2. With the first LOS, we see the disk edge-on, while the second LOS provides a face-on view. For each LOS, the gas density, temperature, the morphology of the stellar component through a rest-frame SDSS mock composite (ugr bands), and the stellar mass maps are displayed. We used a pixel of 0.396′′, and we projected our simulations to a luminous distance of 45 Mpc, which gives a pixel size of 0.12 kpc assuming WMAP9 cosmological parameters values. The physical quantities computed for the gas are all mass-weighted averages along the LOS. The stellar emission is computed using the STARBURST99 model (Leitherer et al. 1999) given the age and the mass of each particle. Unlike Hopkins et al. (2013), we chose to neglect the dust absorption in the building of the SDSS mock images to emphasize the stellar light distribution. Projections misaligned with the AMR grid are always difficult to build. To palliate this common issue, we used multiple convolutions with smoothing kernel sizes adapted to the cell sizes.

The projections in Fig. 3 show a disk with clumps lying in a turbulent medium. The most massive clumps reach masses of ~ 10⁹ M_⊙ (Bournaud et al. 2014). We observe there a gaseous disk thickened by stellar feedback. The edge-on velocity field nevertheless shows clear ordered rotation. The clumps concentrate most of the stellar emission due to young stars, since they host most of the star formation. Figure 4 further emphasizes this highly complex behavior of the gas with substantial turbulence and disk instabilities by showing the mass-weighted average density of one of the most massive merger simulation in our sample (G1_G1_90_90_0). We observe star-forming clumps wandering in a very turbulent ISM where the spiral structures are continuously destroyed by the cooling induced fragmentation and the thermal energy injection from stellar feedback. The edge-on view displays a disk thickened by the tidal torque induced by the recent merger. In Appendix A, projections similar to Fig. 3 are given for three simulations of the MIRAGE sample, covering the evolution up to 800 Myr, displayed in 11 time steps. The whole MIRAGE sample maps, containing 26 figures, are available in Appendix A.

5.1. Initial conditions relaxation

The relaxation of the disk plays a fundamental role at the beginning of the simulation. The low halo concentration when compared to lower redshift, combined with a high gas fraction drives the gas disk towards an unstable state with Q < 1,even though we start our simulation with the requirement Q > 1.5 everywhere in the stellar disks. The high cooling rates of the gas in the initial disk allow very fast dissipation of internal energy. To prevent uncontrollable relaxation that occurs rapidly, we start our simulations with a maximum resolution of 59 pc (ℓ_max = 12) and with a temperature floor for the gas of T = 10⁴  K (see Table 4). To establish the turbulence smoothly afterwards, we progressively increase the resolution every 25 Myr starting from 85 Myr, until we reach a maximum resolution of 7.3 pc, and a temperature floor for the gas of 300 K. This allows the disks spiral features supported by the thermal floor to form quickly during the first time steps. Once the resolution is increased, Jeans instabilities arise and give birth to clumps of a few 10⁸ solar masses, which can quickly merge to form more massive ones. We observe a rapid contraction of the disks, reducing their radial size by ~ 20% during the first 80 Myr (about a third of the dynamical time), owing to the dissipation of energy by the gas component. This ad-hoc relaxation strategy insures that the internal energy of the gas disk dissipates gradually through cooling over 130 Myr, and also helps us to save computational resources. The refinement down to the level ℓ = 14 at t = 105 Myr allows reaching densities ρ > ρ₀, thus enabling the formation of stars and all the associated feedback of newly formed stars.

5.2. Gas accretion from hot halo

As mentioned in Sect. 4.3, the AMR box is continuously replenished with metal-free gas. The very low-density component is constrained by a gamma-polytrope, which ensures the formation of a hot stabilized halo. The central part of the halo reaches densities above 10^-3 cm^-3, where the pressure support from the gamma-polytrope ends. Thanks to metal lines cooling, the central part of the gaseous halo can cool down and condense on the galactic disk. We measure the accretion rates in a spherical shell of 20 kpc (typical value of the halo scalelength in the most massive galaxy model), by detecting the metal poor gas (Z < 10^-3) able to enter the sphere within a time step of 5 Myr. In Table 5 we display the mean values of these accretion rates, as well as the mean SFR for the different masses configurations of the sample. We compare these values to theoretical predictions from the baryonic growth rate formula found in Dekel et al. (2009a). The theoretical values are obtained using total halo masses (without considering radial cuts), and we assume that only two thirds of this accretion rate can be associated to smooth gas accretion, the remaining third being associated to mergers as observed in Dekel et al. (2009a).

Furthermore, at z = 2, Agertz et al. (2009) find an accretion rate of hot gas for a galaxy with a baryonic mass M_b ~ 10¹¹ M_⊙ in good agreement with our model having a close mass (namely the G1 model). They also show that the cold gas accretion that prevails at z ≳ 2 becomes dominated by hot gas accretion at lower redshifts, which makes our implementation agree with this statement. This scenario is also supported by recent work that uses moving-mesh codes, which find a substantially lower cold gas accretion rate than in comparable SPH simulations (Nelson et al. 2013). The gas accretion rate in the MIRAGE sample slightly increases with time (see Sect. 4.3), implying that simulating more than 1 Gyr of evolution would lead to unrealistic high accretion rates. The mean accretion rates measured over 800 Myr in our simulations remain consistent with theory and cosmological simulations (see Table 5).

5.3. Mass–size evolution

To follow the evolution of the mass–size relation in our simulation sample, we proceed to a centering and spin alignment with the z-axis of the AMR Cartesian grid. As the stars age, they experience a progressive gravitational heating that redistributes the oldest stars into a diffuse halo component with a smoother gravitational potential. The center of the each simulation is therefore found by using the peak of the mass-weighted histogram of the positions of the oldest stars, i.e. those stars only present in the initial conditions. This peak in the distribution of the old stars is hereafter associated to the center of the bulge. We recover the galactic disk orientation using the spin vector of the stars younger than 50 Myr for a given snapshot. They tend to still be located within the gas disk, which cannot be used to perform such a computation because of the turbulence and the outflows carrying consequent momentum. Once the orientation of the galactic disk is correctly recovered, we compute the stellar surface density profile. The distribution is decomposed into a bulge and a disk by performing a linear regression on the surface density profile in the interval [r_cut,bulge,r_cut,stars] to extract the disk profile (see Table 1). To each surface density measurement we associate a relative error proportional to the square root of the number of particles found within the radial bin.

The evolution of the disk scalelengths is displayed in Fig. 5. It confirms that both mergers and isolated disks can produce an inside-out growth (Naab & Ostriker 2006) regardless of the orbital configurations, and despite the proven ability of gas-rich mergers to produce compact systems (Bournaud et al. 2011). For each simulation, we estimate the growth time, which we define as the time needed for the stellar disk to double its size measured right after the coalescence (or at t = 400 Myr for the isolated systems). A mean growth time of 3.9 Gyr is measured for the mergers, with the fastest systems reaching growth time close to 2 Gyr. It appears that the less massive systems, therefore the less clumpy, are less efficient at driving stellar disk growth. This inside-out growth is taking place in an idealized framework, although the galaxies are accreting gas from the halo at a rate comparable to cosmological simulations (see Sect. 5.2). This continuous gas accretion fuels secular evolution processes that are able to drive such growth by performing a mass redistribution. Consequently, our results suggest that other mechanisms than late infall of cold gas from the cosmic web (Pichon et al. 2011) may alternatively build up high redshift disks inside-out.

The stellar mass–size relation for the MIRAGE sample is shown in Fig. 6. We plot the mass–size relation found in Dutton et al. (2011) and shifted at different redshifts. Our choice of stellar sizes in the initial conditions makes the simulations of the MIRAGE sample lie in the dispersion range computed for the z = 1.5 mass–size relation. Nonetheless, the size evolution is fast in the MIRAGE sample, although one can expect this rapid growth to stop once the clumpy phase ends. Indeed, the size growth is linked to the gas-rich clump interaction that is able to redistribute significant amount of stellar mass toward the outskirts of the disk. We overplot on the simulations data the values for the MASSIV sample, and the error bars show the 1σ standard deviation computed using the errors on the stellar mass and size (Vergani et al. 2012). We use the classification of López-Sanjuan et al. (2013) to differentiate isolated galaxies from minor and major mergers on the plot using different symbols. We observe that the majority of galaxies classified as major merger lie above or close to the z = 1.5 mass–size relation, which is straightforward once one considers that the size measurement is done on a extended system where the two disks are not yet mixed well. This gives credit for the major merger classification performed by López-Sanjuan et al. (2013). Overall, the bulk of the MASSIV sample ranges within the dispersion fork of the z = 1.5 relation, which makes our simulations consistent with observations. A fraction of the isolated and minor merger systems are more than 1σ below the z = 1.5 relation, suggesting a population of compact galaxies.

5.4. Star formation

Figure 7 presents the star formation histories of the MIRAGE sample for the different masses and merger orbital configurations. For each simulation, we indicate the coalescence time t_c expressed in Myr in the legend, as well as the pericentral time t_peri. The star formation histories exhibit stochastic behavior due to the clump interactions and the cycling energy injection by stellar feedback maintaining the gas turbulence. The mean ratio of the SFR dispersion over the average SFR (σ_SFR/ ⟨SFR⟩) for the whole MIRAGE sample is roughly equal to 30%. Quite surprisingly at first glance, we do not observe any SFR enhancement due to the galaxy merger. Neither orientations nor mass configurations appear to produce enhanced SFR.

Figure 8 shows the histogram of the normalized quantity in the interval t_c ± 100 Myr: $\begin{array}{ccc}\mathrm{\left(}\mathit{SFR}\mathrm{\right(}\mathit{t}\mathrm{)}\mathrm{-}{\mathit{SFR}}_{\mathrm{iso}}\mathrm{(}\mathit{t}\mathrm{\left)}\mathrm{\right)}\mathit{/}\mathrm{⟨}{\mathit{SFR}}_{\mathrm{iso}}\mathrm{⟩}& & \end{array}$with SFR_iso the summed SFR of the fiducial simulations evolved in isolation, and ⟨SFR_iso⟩ the mean value of SFR_iso. For each histogram, we display the value of barycenter of the distribution β in the legend, which allows estimating how much the interaction enhances the star formation in the time interval defined previously. We observe no trend to SFR enhancement due to merger (β ≤ 0), even if in the case G1_G3 two of the merger produces somewhat more stars than the summed fiducial isolated models (0.1 ≤ β ≤ 0.2). However, this value is too low to be considered as a starburst. Generally, the mergers are even less effective at producing stars compared to isolated simulations. This result contradicts other works (e.g. Bournaud et al. 2011; Teyssier et al. 2010; Cox et al. 2008; Powell et al. 2013). Since this paper does not intend to perform a full study of the starburst efficiency in high redshift galaxies, we list subsequently and briefly discuss the possible reasons for the suppression of starburst in our simulation sample:

• Our choice of elliptic Keplerian trajectories might affect the star formation efficiencies of our merger simulations. However, many works have demonstrated that the starburst efficiency of equal mass galaxy mergers is insensitive to the orbits, the disk orientations, and the physical properties of these galaxies (e.g., Mihos & Hernquist 1996; Springel 2000; Cox et al. 2004). Consequently, the initial configuration of the major merger simulations G1_G1 and G2_G2 should not be considered as responsible for the absence of starburst. However, longer interactions would lead to coalescence of more concentrated systems because of the clumps migration, which could enhance a nuclear starburst. That higher mass ratios simulations (G1_G2, G1_G3, G2_G3) that explore more elongated orbits (see Table 3) do not exhibit star formation enhancement might suggests that the orbits and the disk orientations are generally not to blame for this lack of star formation overactivity.

• As highlighted by Moster et al. (2011), the hot gaseous halos implied in a galaxy merger are likely to be heated by shocks, together with an acquisition of specific angular momentum increasing the centrifugal barrier. Both of these processes can push toward a lower starburst efficiency because isolated disks are more effective at accreting gas from the hot halo.

• A complex treatment of the ISM favors the production of hot gas, which systematically lowers SF, as Cox et al. (2006) point out. The simulations performed in Bournaud et al. (2011) constitute a good dataset for direct comparison, owing to the initial conditions definitions very close to our G1 model. The presence of starbursts in comparable simulations when the gas obeys to a 1D equation-of-state suggests a change in the gas response to a galactic interaction.

• Teyssier et al. (2010) demonstrate that the starburst in a low redshift major merger is mostly driven by the enhancement of gas turbulence and fragmentation as long as the numerical resolution allow it to be resolved. It may be more difficult to increase this turbulence and fragmentation at high redshift because both are already high in our isolated gas-rich disks. The isolated disks simulations are indeed able to maintain this high level of turbulence and fragmentation thanks to continuous gas refilling by the hot halo accretion and an efficient stellar feedback. This scenario would suggest that star formation can saturate and prevent starbursts in galaxy mergers of very turbulent and clumpy gas-rich disks.

• High gas fractions (>50%) are maintained throughout the duration of the mergers. These high gas fractions may prevent the formation of a stellar bar in the remnant, which would drag a large amount of gas toward the nucleus to fuel a starburst (Hopkins et al. 2009). The large fraction of cold fragmented gas prevents the formation of a bar in the stellar component. Additionally, the stellar feedback removes gas from the star-forming regions continuously and may also act against the formation of a large stellar disk by lowering the SFR Moster et al. (2011).

• The feedback model adopted in this study might not be energetic enough to succeed in ejecting important quantities of gas on very large scales, especially because of the isotropic hot gas accretion that systematically curbs the outflowing material. The adopted feedback model may be efficient enough to saturate the star formation during the premerger regime, but is not strong enough to deplete the disk from large quantities of gas, which would then be accreted again later, feeding a star formation burst.

Numerous processes can explain the starburst removal in very gas-rich clumpy and turbulent galaxy mergers. The star formation histories of the MIRAGE sample remain difficult to interpret without a complete study in a full cosmological environment to weight each configuration according to its occurrence probability. Generally, the link between mergers and starburst may be more fuzzy at high redshift than at lower redshift.

Figure 9 displays the SFR as a function of the stellar mass. We compared the MASSIV “first-epoch” data with the MIRAGE sample, for which the SFR has been estimated from the integrated Hα luminosity, and stellar mass within the optical radius. The MASSIV error bars were computed using the errors on the Hα flux measurement found in Queyrel et al. (2012). The scatter observed for a given simulation stellar and gas mass underlines the stochastic nature of the star formation in gas-rich clumpy disks. This scatter is nevertheless still lower than the one observed in the MASSIV data, which encompasses very varied gas fractions.

Figure 10 displays the position of the MIRAGE sample on the KS diagram between 200 Myr and 800 Myr for the isolated disks and between the coalescence and 800 Myr for the merger simulations. We computed the gas surface density Σ_gas and the star formation surface density Σ_SFR quantities within r_stars, the stellar disk scalelength estimated with the method described in Sect. 5.3. We also rejected all the gas in the IGM by only considering cells with densities greater than ρ = 2 × 10^-3 cm^-3, which typically corresponds to the frontier between ISM and IGM in all of the simulations. The quantities Σ_gas and Σ_SFR are measured on face-on projections, after having centered our referential on the peak of the old stars probability distribution function and aligned the spin of the young stars disks with our LOS. We note that our sample lies close to the relation found in Daddi et al. (2010b), with a slight shift towards lower star formation efficiencies (within the 1σ dispersion), which can be attributed to the shutdown of star formation at high gas temperatures. The MIRAGE sample does not show any bimodality, as expected from the star formation histories displayed in Fig. 7. One should take into account that, by construction, our sample does not provide the statistical cosmological weight of a volume-limited sample of the 1 < z < 2 galaxy population.

We also overplot the position of the MASSIV sample on this diagram for comparative purposes. The MASSIV error bars on the quantity Σ_SFR are computed again using the uncertainties on the Hα flux. We also take an error proportional to the spatial sampling of the SINFONI data into account, which we propagate to the measurement of the radius of the ionized emission region. We do not have a strong observational constraints for the amount of gas in the MASSIV galaxies. Nevertheless, we mark out the gas mass for each galaxy assuming a mean gas fraction f_g = 45%, which is the mean value obtained on the dynamical/stellar mass diagram of the MASSIV sample (Vergani et al. 2012). Using the relation M_gas = f_g M_stars/(1 − f_g), we can overplot the MASSIV data on the KS diagram. We compute the errors bars of the Σ_gas quantity assuming a minimal gas fraction of f_g,min = 25% and f_g,max = 65%, a range where we can expect the MASSIV sample to lie. We then propagate the errors on the stellar mass using f_g,min and f_g,max. Therefore, the distribution of the MASSIV data on the KS relation is close to the “normal” regime of star formation, considering our assumptions on the gas fractions. Our merger simulations match the area covered by both the isolated and merging galaxies of the MASSIV sample on the KS diagram.

6. Summary and prospects

In this paper, we introduce a new sample of idealized AMR simulations of high redshift (1 < z < 2) mergers and isolated disks referred to as MIRAGE. The sample is originally designed to study the impact of galaxy merger on the gas kinematics in a clumpy turbulent medium. We focus on presenting the methods used to build the MIRAGE sample and on the first results obtained for the evolution of the masses, sizes, and star formation rates.

The key points of the goals and methods used in this paper can be summarized as follows:

• We presented the MIRAGE sample, a series of mergers andisolated simulations using the AMR technique in an idealizedframework that compares extreme signatures in terms of gaskinematics. The MIRAGE sample initial conditions probe fourdisk orientations (with κ ranging from 0° to 180°), five total baryonicmerger masses (ranging from 4.9 to 17.5 × 10¹⁰ M_⊙) and three galaxy mass ra-tios (1:1, 1:2.5, 1:6.3) among 20 merger simulations designedfrom three disk models (with baryonic masses of 1.4, 3.5, and 8.8 × 10¹⁰ M_⊙).The case of low gas fractions has been extensively studied in theliterature, so we choose here to only study gas-rich galaxies (f_g ~ 60%) tostudy the impact of the presence of giant star-forming clumps onmerging turbulent disks.

• We introduced DICE, a new public code designed to build idealized initial conditions. The initialization method is similar to what has been done in Springel et al. (2005b). The use of MCMC algorithm to build a statistical distribution requiring only the 3D-density function as input allows us to consider building components in future developments with more complex density functions compared to the canonical ones used in this paper.

• We used a new implementation of stellar feedback from the young, massive part of the IMF (Renaud et al. 2013), coupled to a supernova feedback with non-thermal processes modeled using a cooling switch (Teyssier et al. 2013). The new physically-motivated implementation of young-star feedback allowed us to track the formation of Strömgren spheres where the energy from the massive young stars is deposited, allowing future comparisons with simulations using feedback recipes parametrized with wind mass-loading factors.

The key results of this paper can be summarized as follows:

• Star formation in disks – We find that the star formation history ofisolated disk galaxies fluctuates strongly throughout the durationof simulations, with a SFR dispersion close to 30% around itsmean value. This star formation proceeds mostly in giant clumpsof gas and stars and naturally gets a stochastic behavior. The smallstar formation bursts may account for the intrinsic scatter of the“main sequence” of star forming galaxies at z = 1–2 (Daddiet al. 2010b).

• Star formation in mergers – The minor and major gas-rich mergers of our sample do not induce major bursts of star formation significantly greater than the intrinsic fluctuations of the star formation activity. The mechanisms for triggering active starburst at high redshift could be different from the ones at low redshift due to large differences in the amount of gas available for accretion in the circumgalactic medium around stellar disks. This suggests that a complex modeling of the gas capturing a high level of fragmentation and turbulence maintained by stellar feedback and gas accretion may offer a mechanism of saturation for the star formation activity in high-z galaxies. The remarkable homogeneity of the observed specific SFR in high redshift galaxies (Elbaz et al. 2007, 2011; Nordon et al. 2012), coupled to the prediction of a high occurrence of minor mergers in this redshift and mass range (Dekel et al. 2009a; Brooks et al. 2009), may support the scarcity of star formation bursts triggered by very gas-rich mergers.

• Star formation scaling laws – Overall, our sample of disks and mergers is compatible with the evolution of the mass-SFR relation observed for a complete sample of star-forming galaxies in the same mass and redshift range (namely the sample MASSIV, Contini et al. 2012), and independently of the assumed disk and merger fraction in the sample. In a Kennicutt-Schmidt diagnostic, the majority of mergers are close to the “normal” regime of disk-like star formation as defined by Daddi et al. (2010b), Genzel et al. (2010), with a slight deviation towards lower star formation efficiencies.

• Size evolution – A stellar mass–size relation in accordance with Dutton et al. (2011) has been obtained by our models, and the evolution with redshift of this relation was also reproduced. In particular, inside-out growth can be obtained as a natural outcome of internal dynamical processes redistributing angular momentum mostly through clumps interactions: these processes can naturally make disks become larger over time, for any given stellar mass, even if mergers are expected to produce more compact systems. Our simulations only include infall of low-angular momentum material through hot gas accretion, suggesting that the radial inside-out growth at the observed rate might not need to be achieved through a cold mode in the context of our idealized modeling.

The significant fraction of active galactic nucleus (AGN) in the redshift range 1 < z < 2 means that the inclusion of black-hole particles and the associated AGN feedback should be addressed specifically. Nevertheless, the recent work of Gabor & Bournaud (2013) shows that AGN feedback is unable to disrupt the clumps of 10⁸–10⁹M_⊙ formed in-situ in comparable idealized gas-rich disk simulations. Dubois et al. (2013) also show that massive clumps may survive AGN feedback during their migration towards the bulge in a fully cosmological context. Since the clumps drive most of the SF, we do not expect major changes on short term star formation histories by including AGN feedback. However, at later stages of evolution, the strong heating of the gaseous halo driven by shocks due to AGN feedback should lower the accretion rate and lead to lower gas fractions in the merger remnants. Finally, the results obtained in the MIRAGE sample call for further investigations to assess the effect of AGN feedback in this kind of simulations.

With the MIRAGE sample we have the opportunity to investigate further questions of galaxy evolution. The combination of statistical probing of the orbital parameters, the controlled input parameters due to idealized framework, the parsec scale resolution, and the explicit physically motived implementation of stellar feedback make it a useful database to use for studying (i) the properties and lifetime of the giant (10⁸–10⁹ M_⊙) star-forming clumps (already presented in Bournaud et al. 2014); (ii) the impact of the migration and interaction of the clumps on the galaxies properties; (iii) the metallicity evolution in mergers and isolated disks; (iv) the classification of velocity fields of high-z galaxies based on a large set of mock observations derived from the MIRAGE sample, among other studies.

Online material

Appendix A: Simulations maps

---

Acknowledgments

We thank the anonymous referee for the useful comments that greatly improved this paper. We thank the whole MASSIV team for encouraging this project with stimulating discussions. The simulations presented in this work were performed at the Très Grand Centre de Calcul of CEA under GENCI allocations 2012-GEN 2192 and 2013-GEN2192, and at the LRZ SuperMUC facility under PRACE allocation number 50816. F.R. and F.B. acknowledge funding from the EC through grant ERC-StG-257720. We thank J. Billing for distributing of the starscream code under the GPL license. Snapshots format conversions were managed using the UNSIO library (http://projets.lam.fr/projects/unsio).

References

All Tables

All Figures