Open Access
Issue
A&A
Volume 620, December 2018
Article Number A21
Number of page(s) 14
Section Interstellar and circumstellar matter
DOI https://doi.org/10.1051/0004-6361/201833161
Published online 23 November 2018

© ESO 2018

Licence Creative Commons
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1 Introduction

Star formation, a key phenomenon in our universe, is still not completely understood. Part of the reason for this is that many physical processes interact over a large range of temporal and spatial scales. Our understanding of star formation is directly linked to our comprehension of the cycle and the dynamics of the interstellar medium (ISM). Because of this large range of scales, it is not possible to simulate a galaxy (e.g. Tasker & Bryan 2006; Dubois & Teyssier 2008; Bournaud et al. 2010; Kim et al. 2011; Dobbs et al. 2011; Tasker 2011; Hopkins et al. 2011; Renaud et al. 2013; Goldbaum et al. 2016; Semenov et al. 2017) with a well-resolved interstellar medium everywhere. A resolution of ~2 pc is possible with adaptive mesh refinement, but only in very dense regions. One possible solution is to simulate a small part of a galaxy in order to have a better spatial resolution (e.g. Korpi et al. 1999; Slyz et al. 2005; de Avillez & Breitschwerdt 2005; Joung & Mac Low 2006; Kim et al. 2011, 2013; Hill et al. 2012; Gent et al. 2013; Hennebelle & Iffrig 2014; Gatto et al. 2015) at the cost of not solving the large galactic scales.

In this paper, we study the influence of stellar feedback (by supernovae and H II regions) and the galactic shear on the ISM and the star formation rate (SFR) by performing a series of simulations of a stratified shearing box, with various physical parameters such as the initial magnetic field or the velocity gradient in the box. The goal is to reproduce the observed Kennicutt law for the SFR (Kennicutt & Evans 2012). It is already known that magnetic fields and turbulence have a limiting effect on star formation (Ostriker & Shetty 2011; Hennebelle & Falgarone 2012). Feedback from stars can also greatly reduce the SFR (Hennebelle & Iffrig 2014; Kim & Ostriker 2015; Gatto et al. 2017). The values of the SFR obtained in Iffrig & Hennebelle (2017; with only supernovae feedback) were actually close to the observed ones. Supernovae inject momentum and energy into the gas around the dying star. In this way, they are capable of unbinding the gas in the host cloud of the massive star. The problem was that in these models, the supernovae directly occurred at the creation of the stars, and thus were not fully realistic.

Our new model for feedback introduces a delay for the supernovae, and H II regions. A H II region is composed of hydrogen that has been ionised by UV photons, and has a temperature of about 104 K. Because of the temperature difference with the surrounding neutral gas, a pressure difference is created across the ionisation front that triggers the expansion of the H II region. This expanding wave perturbs the interstellar medium, and therefore can prevent star formation. Theoretical features of H II regions have been previously described in many papers such as Kahn (1954), Franco et al. (1990), and Matzner (2002).

Galactic shear is another process that competes with gravity and slows down star formation. Previous work with simulations of shearing galactic gas disks has been done (e.g. Kim et al. 2002, 2013). Kim et al. (2013) find that with a gas column density Σgas ∝Ω (where Ω is the galactic angular velocity) so that the Toomre criterion Q is constant at ~ 2, the SFR follows the law . Kim & Ostriker (2017) find a numerically converged SFR ΣSFR ~ 5 × 10−3 M kpc−2 yr−1 for a gas surface density Σgas~ 10 M pc−2.

Section 2 describes the numerical and physical setup that we use, as well as the shearing box we have implemented. In Sect. 3, we give the results we obtained from the various simulations, in terms of SFR and structures properties. Section 4 discusses about the relevant shear value given the surface density, and about the physical limits of the model. Section 5 concludes this paper.

thumbnail Fig. 1

Sketch showing the galactic differential rotation and our rotating frame.

Open with DEXTER

2 Numerical methods

2.1 Numerical code

Our simulations are performed with a modified version of the RAMSES code (Teyssier 2002), a code using the Godunov method to solve the magnetohydrodynamics (MHD) equations (Fromang et al. 2006). We use a 2563 and 5123 grid with no adaptive mesh refinement. For a 1 kpc3 box, that corresponds to a physical resolution of ~4 or ~ 2 pc. The development of the shearing box also lead us to use a planar decomposition (along the z-axis) of the cubic domain, each partial domain being attributed to one core. We use sink particles (Bleuler & Teyssier 2014) to follow the evolution of the dense gas and model star formation.

2.2 Equations

We solve the ideal MHD equations with self-gravity, heating, and cooling processes, in a rotating frame (at Ω0 = Ω(R0)), where R0 is the distance from the centre of the box to the centre of the galaxy. In the galaxy, there is differential rotation, that is, Ω = Ω(R) (see Fig. 1). In the rotating frame, this implies for the velocity along the x-axis (the minus sign comes from the fact that here the x-axis is in the opposite direction to the usual polar axis) (1)

after a Taylor expansion (the box is assumed to be far enough from the galactic centre). In this approximation, there is a constant gradient of velocity along y. We define the velocity gradient (2)

with L being the size of the box.

Following Kim et al. (2002), assuming the angular velocity is a power law Ω (R) ∝ Rq, the effective tidal potential is (3)

We take the value q = 1 (flat rotation curve, Kim et al. 2002). Since Ω(R0) ∝ 1∕R0, we also notice that Vshear = −Ω′(R0)R0 = Ω(R0). From Binney & Tremaine (1987), Ω(R0) = 28 kms−1 kpc−1 corresponds to the value in the solar neighbourhood.

To this effective potential, we must add the potential from self-gravity. We also add an analytical gravity profile accounting for the distribution of stars and dark matter (Kuijken & Gilmore 1989): (4)

with K = 1.42 × 10−3 kpc Myr−2, F = 2.75 × 10−4 Myr−2, and z0 = 180 pc (Joung & Mac Low 2006). Through a Poisson equation, the potential is associated with a density, (5)

Adding the Coriolis force, the MHD equations that are solved are

with ρ, v, P, B, Φ, and E, respectively, being the mass density, velocity, pressure, magnetic field, total gravitational potential, and total (kinetic plus thermal plus magnetic) energy.

The total gravitational potential is written Φ = ϕ + ϕext, with ϕext being the analytical potential as explained above. For non-ionised gas, the loss function is such that , where n is the particle density, Γ represents constant and uniform UV heating, and Λ(T) is the cooling function, similar to the one used by Joung & Mac Low (2006). For photoionised gas (H II regions), radiative cooling and heating is performed by the radiative transfer module RAMSES-RT (Rosdahl et al. 2013). We use the same model as in Geen et al. (2016), where they add radiative heating and cooling functions for metals.

thumbnail Fig. 2

Drawing of the shearing box used at the y boundaries.

Open with DEXTER

2.3 Boundary conditions and shearing box

A schematic of the simulated volume in the context of the galactic disk, with axis directions, is shown in Fig. 1. We use periodic boundary conditions along the x-axis, and zero-gradient (except for Φ) boundary conditions along the z-axis (outside thegalactic plane). The most difficult part concerns the boundary conditions along the y-axis. Simple periodic conditions cannot be used, as the two boundaries have opposite and non-zero velocities in our rotating frame (see Sect. 2.2). A shearing box is therefore implemented (Hawley et al. 1995; see Fig. 2).

In each ghost cell in the upper y boundary, the values of ρ, vy, vz, and T are equal to those of the corresponding cell in the lower domain, shifted of − Vshear × t in the x direction. For the x-velocity vx, it is not simply the value from the lower cell, but the necessary differential of velocity that is added: vx + Vshear is applied in the ghost cell.

The treatment of the magnetic field is a bit different because zero-divergence must be assured in the cells, as well as the continuity of the normal field at the boundary (in RAMSES, each cell has six values for the magnetic field, one on each face). We do the following: first the values of the tangential fields (Bx and Bz) are applied in the ghost cell. They are read as usual in the corresponding cell of the domain. The continuity of B is then assured by taking the value of the normal By of the cell just below. Finally, the remaining unknown By in the ghost cell is specifically chosen in order to ensure the zero-divergence of the magnetic field. The same treatment is applied in the z-boundaries, adapted to zero-gradient conditions. Then, the energy per unit of mass must be defined in the ghost cell: the internal energy from the cell of the domain is taken, but the kinetic and magnetic energies have to be re-calculated.

The last physicial parameter that remains unknown in the ghost cell is the gravitational potential Φ. The same treatment cannot be applied, since the Poisson solver is an iterative solver (for a single timestep), but RAMSES does not update the boundary conditions at each iteration of the Poisson solver. Therefore, an analytical boundary condition for Φ has to be taken (at the y-boundaries and z-boundaries only, the x-boundaries being periodic and correctly solved), which does not need to be updated at each of these iterations (but that is updated at each timestep).

We proceed as follows: we notice that the external potential ϕext (Eq. (4)) is far stronger than the self-gravitating potential ϕ. Hence, in order to find the analytical boundary condition, an approximation is made that the self-gravitating potential in the box is of the form (11)

with α and z0 to be determined at each timestep. This approximate ϕapp will be used in the boundary condition (12)

Via the Poisson equation, this potential is associated with an approximate density (13)

By integrating this density, α is deduced from the total mass of gas M in our box, (14)

To estimate the scale height, z0, we simply find the altitude at which the density is given by (15)

Finally, the boundary condition used in the y and z directions is (16)

With this boundary condition, the gas is stable on large scales while otherwise a global collapse occurs if the usual shearing box treatment is applied. We have performed a test for the accuracy of these gravitational boundaries (see Appendix A) and the error is around 3040%, which is lower than the incertitude coming from our numerical setup (see Sect. 3.3).

Section 3 of Shen et al. (2006) describes hydrodynamic tests for shearing boxes. They have been performed in our case, and are in good agreement with analytical results.

2.4 Initial conditions

The size of our cubic box is 1 kpc. We impose initially a Gaussian density profile (for hydrogen atoms) along the z-axis, (17)

with n0 = 1.5 cm−3 and z0 = 150 pc. The initial column density is then , where mp = 1.4 × 1.66  × 10−24 g is the mean mass per hydrogen atom. The initial temperature is chosen to be 8000 K. This is a typical temperature of the warm neutral medium (WNM) phase of the ISM. We also add an initial turbulent velocity field with a root mean square (RMS) dispersion of 5 km s−1 and a Kolmogorov power spectrum with random phase (Kolmogorov 1941).

Finally, we add a Gaussian magnetic field, oriented along the x-axis, (18)

with B0 ≃ 0 or 4 μG depending on the simulation (magnetohydrodynamic or just hydrodynamic).

Our different runs are detailed in Table 1.

The Toomre criterion for our shearing disk of gas is (Binney & Tremaine 1987) (19)

where cs is the sound speed of the gas. This disk is stable to axisymmetric perturbations if (20)

With a speed of sound of 8 km s−1, the values for our simulations are listed in Table 2. These values are very close to the stability criterion, hence the consistency of our parameters.

Table 1

Different runs.

2.5 Sink particles

The star-forming regions are represented by Lagrangian sink particles (Krumholz et al. 2004; Bleuler & Teyssier 2014): the resolution being too low to describe each star individually, they are described as a group by subgrid models. The sink particles are created when a density threshold is reached (usually 103   cm−3; see Sect. 3.3), and when virial criteria are satisfied (see Bleuler & Teyssier 2014). As the particles represent a group of stars and because the resolution is coarse, it is unclear that virial criteria should be used. Therefore, we have performed simulations with and without the criteria. In the following, unless explicitly mentioned, the virial criteria are being tested. The sink particles gain mass by accreting gas from their surroundings: when the density at the position of the sink is higher than this same threshold, the mass excess is accreted.

Table 2

Values of the Toomre criteria (Eq. (19)) for different values of Vshear.

2.6 Supernova and H II feedback

The supernova feedback implemented in RAMSES is based on scheme C of Hennebelle & Iffrig (2014): a massive star is formed each time a sink particle accretes more than 120 M. Its mass is determined randomly assuming a Salpeter initial mass function (Salpeter 1955; Chabrier 2003). This massive star is placed randomly within a 10 pc radius around the sink. Its lifetime is calculated from its mass, and the supernova is triggered at the end of its life at its current position: we assume the star has moved with a velocity of vdisp = 1 km s−1 in the sink frame throughout its life. Simulations with a lower or higher density threshold, or with vdisp = 0 have also been performed, and are discussed in Sect. 3.3. The supernova then consists of an injection of 4 × 1043 g cm s−1 (Iffrig & Hennebelle 2015) of radial momentum around the place of the star.

We manually limit the temperature to Tsat = 106 K and the velocity provided by the supernova to Vsat = 500 km s−1 in order to prevent too high sound speeds and too small timesteps. We also performed some simulations with greater values (see Table 1): Tsat = 107 K and Vsat = 1000 km s−1, to analyse their influence.

The maximum of velocity Vsat is set higher than in the previous simulations because now, as previously explained, the supernova is not set immediately but rather at the death of the star: therefore, the explosion may happen in diffuse gas and the maximum velocity would be too easily reached, thus causing a loss of precision.

Feedback by H II regions is also implemented in the code. The energy and momentum injected by the H II region (via the pressure gradient) in the interstellar medium is highly dependent on the flux of ionising photons S* emitted by the star. Using the data of Vacca et al. (1996), we predicate the following relation between S* and the mass of the star: (21)

We then have the asymptotic relations (22)

The results of this empirical model are given in Table 3 and Fig. 3. The evolution of H II regions is then computed with the radiative transfer module RAMSES-RT (Rosdahl et al. 2013). Comparisons between analytical models and simulations results for the expansion of H II regions with RAMSES-RT have been done by Geen et al. (2015). RAMSES-RT has been validated by the STARBENCH comparison project (Bisbas et al. 2015), which tested the D-type expansion of an H II region with different codes. Other works have recently studied the role of H II regions in star formation (Butler et al. 2017; Peters et al. 2017).

Table 3

Parameters of Eq. (21).

thumbnail Fig. 3

Fitting of Eq. (21) with the data of Vacca et al. (1996).

Open with DEXTER

3 Results

3.1 Global considerations

Figure 4 shows the column density in two of the MHD simulations (at a resolution of 2563), with feedback by supernovae and H II regions, at 80 Myrs. The coloured points are sink particles, as explained in Sect. 2.5. We can see the effects of the shearing box in these images: the y-boundaries are now well defined, with gas and sink particles getting through (with a corresponding new position and velocity). Figure 5 shows the same physical simulations at the same time, at a higher resolution. The global patterns are thus broadly the same, but the small scales and the small high-density regions become more visible. The impact of the shear is that the gas is more scattered at higher shear, and more condensed at lower shear. We notice that the galactic shear tends to stretch the large structures of gas: it competes with gravity and opposes the collapse of the molecular clouds. Hence, the galactic shear is expected to reduce the star formation rate. We will quantify this effect in Sects. 3.2 and 3.5.

Supernovae and H II regions are also implemented, as explained in Sect. 2.6. The temperature and density maps of Fig. 6 show the effects of supernovae on the interstellar gas: there are peaks of temperature and minima of density at the supernovae locations. The surrounding gas is heated and dispersed (see Sect. 3.2 for the effect on star formation). Some of the explosions are at the sink particle location, but some others are lightly shifted from the star cluster. This is the consequence of the scheme described in Sect. 2.6: the star is moving, and the supernovae is triggered at the end of its life.

3.2 Influence of the physical parameters on the star formation rate (SFR)

Figure 7 displays the total mass of stars (in units of M) formed over time in our 1 kpc3 box, for different physical simulations. The slope of the dashed line is the observed SFR given our column density, which we want to reproduce. This value ΣSFR ≃ 7 × 10−3 M yr−1 kpc−2 comes from Fig. 11 of Kennicutt & Evans (2012).

Figure 7a shows the influence of the feedback (by supernovae and H II regions) on the total mass of stars. Without any kind of feedback (black curve), the mass of stars (hence the SFR) is very high. When feedback by H II regions is added (blue curve), it decreases by20%. The feedback by supernovae is, however, dominant: with supernovae and no H II regions (green curve), the asymptotic SFR decreases by a factor approximately three. The red curve is obtained with both kinds of feedback: it is now far closer to the observed SFR, and it shows the huge impact of the feedback on the ISM. However, in these simulations, the impact of the supernovae is lesser and the SFR is greater than in Iffrig & Hennebelle (2017), because the massive stars now explode with some delay (see Sect. 2.6). The explosions can now happen outside the clouds, in more diffuse gas, where they are less effective at suppressing star formation.

There isa non-linear coupling between the supernovae and the H II regions that can be seen in the gas distributions of these four simulations, shown in Fig. 8. When one kind of feedback is absent, there is a drop in the distribution function, between 1 and 2 cm−3. Both kinds of feedback must be present for the hollow to disappear (red). The effect of both supernovae and H II regions (on the SFR or on the PDF) is not the sum of the effects of the individual processes. Such non-linear coupling was also noticed in the simulations of Bournaud (2016) where the outflow rate when both processes are present was above the sum of the outflow rates in the individual cases.

Figure 7b shows the influence of the magnetic field: without magnetic field, the total stellar mass increases by a factor of two. The magnetic field still strongly suppresses the SFR, and has a strong impact on the shape of the ISM: in the presence of magnetic field, the clumps are more filamentary than they are in pure hydrodynamic cases (Hennebelle 2013). In MHD simulations, the matter is gathered in more stable filaments and that contributes to reduce the SFR.

Figure 7c and d: these plots show how the total mass of stars depends on the shear. The velocity gradient Vshear takes the values: 0 (black), 14 (blue), 28 (red), and 56 (green) km s−1 kpc−1. For the MHD cases, the higher the gradient, the lower are the mass and the SFR. This is exactly the behaviour that was expected, as stated in Sect. 3.1.

This is not so true for hydrodynamic simulations. This comes from the fact that in the hydrodynamic case, the filaments are very unstable, and thus the shear is less effective.

In MHD, the case with Vshear = 14 km s−1 kpc−1 (blue) is very similar to the one without shear (black). This is probably because the Toomre stability criterion Eq. (20) is not reached at this Vshear.

For the MHD simulations, the asymptotic behaviour of the green curve is close to the observed SFR (too high by a factor of ~ 1.5). Nevertheless, this is for a gradient Vshear = 56 km s−1 kpc−1, while the value in the Sun’s neighbourhood is 28 km s−1 kpc−1 (see Sect. 2.2). As for the SFR for a velocity gradient Vshear = 28 km s−1 kpc−1 (red curve), it is still too high by a factor of approximately four in the simulations compared to observations. However, as described in Sect. 4.1, 56 km s−1 kpc−1 may be the right value to use.

thumbnail Fig. 4

Column density for our MHD simulations (2563), with feedback by supernovae and H II regions, at 80 Myrs. For each plane, there are two rows: the first one does not display sink particles, but the second one does, each one representing a star cluster (the colours are random). Left panel: Vshear = 28 km s−1 kpc−1. Right panel: Vshear = 56 km s−1 kpc−1.

Open with DEXTER
thumbnail Fig. 5

Column density for our MHD simulations (5123), with feedback by supernovae and H II regions, at 80 Myrs. Left panel: Vshear = 28 km s−1 kpc−1. Right panel: Vshear = 56 km s−1 kpc−1.

Open with DEXTER
thumbnail Fig. 6

Density and temperature maps of our SH28_mhd simulation (2563) at 85 Myrs. The arrows on the density map represent the velocity of the gas. The regions above 105 K are due to the supernovae.

Open with DEXTER
thumbnail Fig. 7

Total mass of stars formed over time for different simulations with different physical parameters. The slope of the dashed line is the observed SFR for our Σgas (Kennicutt & Evans 2012). Top-left panel: importance of the two types of feedback. Top-right panel: importance of magnetic field. Bottom panels: illustrate the importance of the shear. Altogether we see that there is no dominant effect but instead all physical processes, namely, feedback by supernovae, ionising radiation, magnetic field, and shear, all contribute appreciably to reduce the SFR.

Open with DEXTER

3.3 Numerical parameters and numerical convergence

We now turn to Fig. 9 to study the influence of several numerical parameters on the star formation process.

Figure 9a: as explained in Sect. 2.6, simulations with higher values for Tsat and Vsat have been performed. In these simulations, the total mass of stars formed is lower because the supernovae are more effective at injecting momentum into the ISM. However, in the end, the asymptotic slopes (and then the asymptotic SFR) are quite similar. Thus, at these values, these two parameters do not have a strong impact on the asymptotic SFR. This justifies the choice of limiting the temperature and the velocity of the supernovae explosions that reduce the timestep.

Figure 9b: as stated in Sect. 2.6, a simulation with vdisp = 0 has also been performed, meaning that the massive stars leading to supernoave now stand motionless in the sink particle’s frame. In this simulation, the total mass of stars and the SFR is roughly the same as in the simulation with vdisp = 1 km s−1. Thus, the value of this parameter has a negligible influence on the SFR. This is significantly different from the conclusion of Hennebelle & Iffrig (2014; see also Gatto et al. 2015) and could stem from the fact that support from galactic shear and H II regions is now considered, or that the supernovae can happen in more diffuse gas and be less effective because of the delay of the explosion.

Figure 9c shows that at a higher resolution, the dependence of the SFR on the shear is still the same: when the velocity gradient increases, the SFR decreases. This result is consistent with the distribution functions of the gas for both simulations, which are simultaneously displayed in Fig. 10. The function is higher for greater shear simply because, as previously stated, the gravity is less efficient and less gas is accreted into the sink particles. However, the pattern is the same for both simulations.

Going back to the SFR plots, while at 2563 the final slope/SFR of the green curve (SH56) matches the observed one (see the end of Sect. 3.2), it is no longer the case at 5123. The final SFR is too high (by a factor of about two) at higher resolution for the column density used in this paper. In order to explain this, we have performed complementary runs.

Figure 9d shows the mass of stars formed for two identical simulations (Vshear = 28 km s−1 kpc−1), except for the resolution (2563 vs. 5123, blue and red): the higher the resolution, the higher the mass, and the higher the SFR. The problem is that the asymptotic behaviours do not match, the final SFR depends on the resolution (as stated above for the case Vshear = 56 km s−1 kpc−1): therefore, it seems numerical convergence for the SFR has not been attained.

However, it turns out that the numerical creation and accretion threshold of the sink particles (see Sect. 2.5) can explain why we do not have numerical convergence of the SFR. Indeed, this parameter has been chosen until now to be constant and independent of the resolution (103   cm−3). This threshold value was not well defined. Kim & Ostriker (2017) argue that it should vary as (23)

where Δx is the size of one cell. The quantity Δx2 is then proportional to the surface of the cell. Figure 9d displays a new simulation (2563, green, with nsink four times lower (i.e. 250 cm−3) than the threshold of the high resolution simulation (red). This way, the previous proportionality rule is respected. With this prescription, the SFR of the low resolution case is now very close to that of the high resolution case, and thus numerical convergence is reached.

Figure 9e shows the total mass of stars for the same resolution simulations but with three different sink thresholds: 250 cm−3, 103   cm−3, and 4 × 103 cm−3. The one with the higher nsink (purple) has some bursts and more mass, but the asymptotic SFR (from 100 Myrs) is actually lower, as one would expect. It is, however, still higher than the observed one. Between the two extreme cases, there is a factor of two in the asymptotic SFR. The lower threshold case has an asymptotic SFR too high by a factor approximately six whereas the higher threshold case has an asymptotic SFR too high by a factor of approximately three compared to the observed one. This leaves an open question: even if numerical convergence is reached, what proportionality coefficient (or what reference value) must be taken in the numerical formula (23)? Without answering that question, there is still an incertitude of a factor of two in the SFR. This uncertainty should be understood as a limit to the present models.

thumbnail Fig. 8

Distribution functions of the density of the gas for SH28_mhd simulations, with or without feedback, at 80 Myrs.

Open with DEXTER

3.4 Star formation in simulations without virial criteria for the sink particles

As explained in Sect. 2.5, we have performed runs without the virial criteria for the sink particles. Figure 11 shows their star formation.

Figure 11a shows the total mass of stars for two MHD runs, with both kinds of feedback, with Vshear = 28 or 56 km s−1 kpc−1. Compared to the cases where the virial criteria are being tested (Fig. 7e), the total masses of stars are here lower by a factor of approximately two. This is because without the virial tests, the sink particles are created sooner, and thus stellar feedback takes effect sooner. With the virial criteria, more mass of gas is accreted under gravity before the particles areformed, causing the many bursts we see in Fig. 7 but not in Fig. 11.

The total masses are here lower, however the asymptotic SFRs are quite similar to the ones in the runs with the virial criteria. When Vshear = 28 km s−1 kpc−1 (red), the SFR is too high by a factor of approximately four compared to the observed one, and when Vshear = 56 km s−1 kpc−1, it is too high by a factor of ~1.5 as stated in Sect. 3.2.

Figure 11b shows the dependence of star formation on the threshold of the sink particles threshold nsink. The three curves are associated with the three previous values of nsink: 250 cm−3, 103 cm−3, and 4 × 103 cm−3. Contrary to the cases where the virial criteria must be satisfied (Fig. 9e), the total mass of stars and the SFR do not depend significantly on the value of this numerical parameter. This is likely because when virial criteria are considered, enough gas must have piled up to satisfy them. Here, the sink particles are being introduced almost immediately and this makes the feedback less dependent on the system history.

3.5 Structure properties

We now look at the properties of the dense clouds. The Larson relations (Larson 1981) constitute observables that are worth reproducing. More specifically, the relation between the velocity dispersion σ in the clumps and their size R is expected to be (24)

From the data and the results of Falgarone et al. (2009), typical values of the parameters are σ0 = 1.1 km s−1 and α = 0.5. These values are not well established in the literature (e.g. Heyer et al. 2009; Miville-Deschênes et al. 2017) infer a relation that depends upon the column density of the cloud: ). However, a slightly different σ0 or power index α would not change the following results.

In order to identify structures in our simulations, we proceed as in Iffrig & Hennebelle (2017) and use a friends to friends algorithm with a density threshold of 50 cm−3. For each structure, the velocity dispersion σ and the size R are then computed like this: (25)

where the summations are done on the cells that define the clump. The λi are the eigenvalues of the inertia matrix defined in the centre of mass of the clump, (26)

Figure 12 displays the velocity dispersions of the clumps versus their sizes, for a velocity gradient Vshear = 28 or 56 km s−1 kpc−1. The red straight line is the observed Larson relation (Eq. (24)). The structures of our simulations follow a similar power law, with a slope of the same order (α = 0.5). However, the velocity dispersions are well below the observed ones, by a factor of approximately three, as already noted in Iffrig & Hennebelle (2017). This shows that the shear does not modify this conclusion and this may indicate that other energy sources are present and drive the turbulence, like, for example, the large scale gravitational instabilities that are discussed in Krumholz & Burkhart (2016) and Krumholz et al. (2018). They show that for massive galaxies with high gas column density and for rapidly star-forming galaxies, gravity-driven turbulence is dominant over feedback-driven turbulence. However, such large scale gravitational instability is absent in the model of this paper.

Another possible reason for these low-velocity dispersions is that for smaller structures (with a size of a few cells), the dispersions are damped by numerical dissipation. However, even the structures larger than approximately ten cells are below the expected velocity dispersion. Moreover, Iffrig & Hennebelle (2017) performed runs at various resolutions and did not find significantly different dispersions between their more or less resolved runs (see their Fig. 11 with runs B1 and B2L). If the low-velocity dispersions were caused by numerical diffusion, the high-resolution simulations would have less dissipation for a given physical size and the velocity dispersions would be higher, but this is not the case. Therefore, it seems more likely that this is actually a consequence of a missing injection of energy.

For small sizes (log(R) ≲ 0.8), the simulation with the highest shear (right) has less clumps and they are less massive than in the simulation with the lowest shear (left). For large sizes (log(R) ≳ 0.8), the simulation with the highest shear (right) has more clumps and they are more massive than that of the lowest shear. This confirms what is stated in Sect. 3.1: the galactic shear stretches the clouds (hence more large-sized structures and less small-sized structures at higher shear) and competes with gravity (hence less gas accreted into the stars).

thumbnail Fig. 9

Same as Fig. 7 but showing the influence of numerical parameters. Panels a, b: SFR does not really depend on the limitations on velocity and temperature of the gas and on the velocity of the massive stars in the sink particles. Panel c: influence of the shear on simulations at higher resolution. Panels d, e: star formation depends on the density threshold of the sink particles, giving a numerical incertitude factor of approximately two on the SFR.

Open with DEXTER
thumbnail Fig. 10

Distribution functions of the density of the gas for MHD simulations (at 5123), at 80 Myrs. Red: Vshear = 28 km s−1 kpc−1. Green: Vshear = 56 km s−1 kpc−1.

Open with DEXTER

4 Discussion

4.1 Shear and column density

We have seen in Sects. 3.2 and 3.3 that the galactic shear competes with gravity and greatly reduces the SFR. However, at the solar value Vshear = 28 km s−1 kpc−1, the SFR was too high compared to the observed Kennicutt value. The reference value Vshear in the solar neighbourhood may not be the correct choice to reproduce the Kennicutt relation given our initial column density Σgas = 19.1 M pc−2.

Star-forming clouds are generally in spiral arms, where the shear is usually high (Roberts 1969). In order to quantify this effect, we use simulations of a galaxy that is similar to the Milky Way, with the same radial density profile of the gas, the same stellar mass distribution, and the same rotation curve. These simulations are presented in Renaud et al. (2013) and Kraljic et al. (2014). Figure 13 of Renaud et al. (2013) shows the high shear in spiral arms: there can be a velocity gradient of 40 km s−1 over 200 pc. From these simulations, 1 kpc boxes are isolated in the outer disk, between 5 and 10 kpc from the galactic centre. This corresponds to regions located within or between spiral arms. From each one of these boxes, we compute the mass-weighted average of Σgas, the usual surface-weighted average Σgas,mean (which respects mass conservation), and the mass-weighted mean velocity gradient Vshear. The surface-weighted Σgas are used to define some bins, and the final results are presented in Table 4. Figure 13 is the corresponding scatter plot of Vshear versus log(Σgas,mean).

The higher the column density, the higher is the velocity gradient. Then, for the column density of this paper, Σgas,mean = 19.1 M pc−2, the relevant value should not be Vshear = 28 km s−1 kpc−1 but about Vshear ~ 63 km s−1 kpc−1. As previously stated (Figs. 7d and 9c), the SFR is very close to the observed one for Vshear = 56 km s−1 kpc−1. That shows that density and shear must not be taken independently, and that with a relevant value of Vshear we are able to reproduce Kennicutt’s law (within our intrinsic uncertainty of a factor approximately two).

We will explore the full parameter space of Σgas,mean and Vshear in a future work.

thumbnail Fig. 11

Same as Fig. 7 but for simulations without the virial criteria for the sink particles. Left panel: total masses of stars are now lower but that the SFRs are similar to the ones measured in the runs with the virial tests. Right panel:without the virial criteria, star formation no longer depends on the numerical density threshold of the sink particles.

Open with DEXTER
thumbnail Fig. 12

Left panel: velocity dispersion–size relation of the clumps at 80 Myr for the simulation SH28_mhd_512. Right panel: for the simulation SH56_mhd_512. Each square is a two-dimensional bin, and the colour depends on the total mass of the clumps associated with this bin. The red line is the observed Larson relation (Eq. (24)), from Falgarone et al. (2009).

Open with DEXTER

4.2 Physical limits of the model

The star formation efficiency (SFE) could be another key to restrain star formation: all the gas accreted onto the clusters is not necessarily converted into stars. A large fraction can be re-injected into the interstellar medium through stellar feedback operating at the cluster scale. The SFE is defined for an embedded cluster as Mgas∕(Mgas + Mstars) and its value is usually inferred to be around 20% (Lada & Lada 2003). However, as shown in Lee & Hennebelle (2016) this value is not well defined and should be regarded with great care. In the simulations of this paper, the SFE is simply equal to 100%. All the gas accreted onto the sink particles is considered as stars. Having a model with a lower efficiency would obviously reduce the SFR.

The model of this paper also does not take stellar winds into account. Unlike supernovae, winds are not delayed, and Gatto et al. (2017) find that they regulate the growth of young star clusters by suppressing the gas accretion on them shortly after the first massive star is born. In their model, when winds and supernovae are present, the SFR is close to the one with stellar winds only because the gas surrounding the cluster is already unbound by the winds and supernovae have little additional effect. They also find that adding stellar wind feedback to supernovae decreases ΣSFR by a factor of approximately two, compared to the case with supernovae feedback only. Thus, stellar winds could be another missing physical process that would reduce the SFR even more in our model.

Table 4

Mean gradient of velocities for different bins of log Σgas.

thumbnail Fig. 13

Scatter plot of Table 4. The mean value of the shear in the galaxy simulation of Renaud et al. (2013) is given as a function of Σgas,mean.

Open with DEXTER

5 Conclusions

We have performed a set of simulations of 1 kpc3 galactic regions, with differential rotation, magnetic field, and stellar feedback. We have confirmed that the feedback from the stars has a huge limiting role on the star formation. Our simulations also point out that with both H II regions and supernovae, there is a more effective suppression of star formation than the linear combination of both effects would suggest.

As expected, the value of the galactic shear has a great impact: the higher the velocity gradient, the lower the SFR. Using a galactic simulation, it seems this value actually depends on the column density of gas Σgas. With a gradient Vshear = 56 km s−1 kpc−1 and a column density Σgas = 19.1 M pc−2, we get an asymptotic SFR very close to that of the observed Kennicutt law.

We manage to get numerical convergence of the SFR by making the threshold parameter of the sinks dependent on the spatial resolution: nsink ∝ 1∕Δx2. However, the proportionality coefficient is not well defined, and gives an incertitude of a factor of approximately two on the asymptotic SFR.

Finally, the analysis of the structures of the properties confirmed that the galactic shear stretches the gas clouds. We also find that the velocity dispersion is below the observed Larson law, by a factor of approximately three. This indicates that we may have missed someenergy injection from the large galactic scales that drives the turbulence. We would expect such an injection to reduce the SFR even more.

Acknowledgements

This work was granted access to HPC resources of CINES and CCRT under the allocation A0010407023 made by GENCI (Grand Equipement National de Calcul Intensif). This research has received funding from the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 306483). Wethank the referee for his helpful report.

Appendix A Relevance of our gravitational boundary conditions

thumbnail Fig. A.1

Star formation in the simulation (Vshear = 0) with “shearing box” boundary conditions for gravity (Eq. (16)) compared with the same simulation with true periodic boundary conditions.

Open with DEXTER

The boundary conditions used for gravity in Sect. 2.3 are still an approximation.For instance, as the boundary potential only depends upon z, there is no force in the x and in the y directions. This can be troublesome, as gas and clouds are supposed to cross the boundary. To test the influence of this approximation for the boundary conditions, we have performed a comparison between a fully periodic box and one with these conditions, for a run without shear (Vshear = 0). These runs were performed without the virial criteria of the sink particles (see Sect. 2.5). The corresponding star formation curves are displayed in Fig. A.1.

The simulation of the fully periodic box (blue) has a jump of star formation at around 110 Myr, making it difficult to define a proper SFR. However, the total mass of stars is about only 30% higher in the true periodic case at the end of the runs. Then, estimating the slope from the end of the jump (at 115 Myr), the case with the approximated boundary conditions has a SFR higher by ~ 40% compared tothe fully periodic box. These differences are lower than the numerical incertitude factor of two coming from the sink particles, so we neglected it in our discussion. Moreover, due to the large fluctuations that the SFR is experiencing, a finer quantification would require the realisation of many runs. One improvement of our numerical setup would be to define better boundary conditions for gravity when differential rotation is present.

References

All Tables

Table 1

Different runs.

Table 2

Values of the Toomre criteria (Eq. (19)) for different values of Vshear.

Table 3

Parameters of Eq. (21).

Table 4

Mean gradient of velocities for different bins of log Σgas.

All Figures

thumbnail Fig. 1

Sketch showing the galactic differential rotation and our rotating frame.

Open with DEXTER
In the text
thumbnail Fig. 2

Drawing of the shearing box used at the y boundaries.

Open with DEXTER
In the text
thumbnail Fig. 3

Fitting of Eq. (21) with the data of Vacca et al. (1996).

Open with DEXTER
In the text
thumbnail Fig. 4

Column density for our MHD simulations (2563), with feedback by supernovae and H II regions, at 80 Myrs. For each plane, there are two rows: the first one does not display sink particles, but the second one does, each one representing a star cluster (the colours are random). Left panel: Vshear = 28 km s−1 kpc−1. Right panel: Vshear = 56 km s−1 kpc−1.

Open with DEXTER
In the text
thumbnail Fig. 5

Column density for our MHD simulations (5123), with feedback by supernovae and H II regions, at 80 Myrs. Left panel: Vshear = 28 km s−1 kpc−1. Right panel: Vshear = 56 km s−1 kpc−1.

Open with DEXTER
In the text
thumbnail Fig. 6

Density and temperature maps of our SH28_mhd simulation (2563) at 85 Myrs. The arrows on the density map represent the velocity of the gas. The regions above 105 K are due to the supernovae.

Open with DEXTER
In the text
thumbnail Fig. 7

Total mass of stars formed over time for different simulations with different physical parameters. The slope of the dashed line is the observed SFR for our Σgas (Kennicutt & Evans 2012). Top-left panel: importance of the two types of feedback. Top-right panel: importance of magnetic field. Bottom panels: illustrate the importance of the shear. Altogether we see that there is no dominant effect but instead all physical processes, namely, feedback by supernovae, ionising radiation, magnetic field, and shear, all contribute appreciably to reduce the SFR.

Open with DEXTER
In the text
thumbnail Fig. 8

Distribution functions of the density of the gas for SH28_mhd simulations, with or without feedback, at 80 Myrs.

Open with DEXTER
In the text
thumbnail Fig. 9

Same as Fig. 7 but showing the influence of numerical parameters. Panels a, b: SFR does not really depend on the limitations on velocity and temperature of the gas and on the velocity of the massive stars in the sink particles. Panel c: influence of the shear on simulations at higher resolution. Panels d, e: star formation depends on the density threshold of the sink particles, giving a numerical incertitude factor of approximately two on the SFR.

Open with DEXTER
In the text
thumbnail Fig. 10

Distribution functions of the density of the gas for MHD simulations (at 5123), at 80 Myrs. Red: Vshear = 28 km s−1 kpc−1. Green: Vshear = 56 km s−1 kpc−1.

Open with DEXTER
In the text
thumbnail Fig. 11

Same as Fig. 7 but for simulations without the virial criteria for the sink particles. Left panel: total masses of stars are now lower but that the SFRs are similar to the ones measured in the runs with the virial tests. Right panel:without the virial criteria, star formation no longer depends on the numerical density threshold of the sink particles.

Open with DEXTER
In the text
thumbnail Fig. 12

Left panel: velocity dispersion–size relation of the clumps at 80 Myr for the simulation SH28_mhd_512. Right panel: for the simulation SH56_mhd_512. Each square is a two-dimensional bin, and the colour depends on the total mass of the clumps associated with this bin. The red line is the observed Larson relation (Eq. (24)), from Falgarone et al. (2009).

Open with DEXTER
In the text
thumbnail Fig. 13

Scatter plot of Table 4. The mean value of the shear in the galaxy simulation of Renaud et al. (2013) is given as a function of Σgas,mean.

Open with DEXTER
In the text
thumbnail Fig. A.1

Star formation in the simulation (Vshear = 0) with “shearing box” boundary conditions for gravity (Eq. (16)) compared with the same simulation with true periodic boundary conditions.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.