Impact of galactic shear and stellar feedback on star formation

A numerical shearing box is used to perform three-dimensional simulations of a 1 kpc stratified cubic box of turbulent and self-gravitating interstellar medium (in a rotating frame) with supernovae and HII feedback. We vary the value of the velocity gradient induced by the shear and the initial value of the galactic magnetic field. Finally the different star formation rates and the properties of the structures associated with this set of simulations are computed. We first confirm that the feedback has a strong limiting effect on star formation. The galactic shear has also a great influence: the higher the shear, the lower the SFR. Taking the value of the velocity gradient in the solar neighbourhood, the SFR is too high compared to the observed Kennicutt law, by a factor approximately three to six. This discrepancy can be solved by arguing that the relevant value of the shear is not the one in the solar neighbourhood, and that in reality the star formation efficiency within clusters is not 100%. Taking into account the fact that star-forming clouds generally lie in spiral arms where the shear can be substantially higher (as probed by galaxy-scale simulations), the SFR is now close to the observed one. Different numerical recipes have been tested for the sink particles, giving a numerical incertitude of a factor of about two on the SFR. Finally we have also estimated the velocity dispersions in our dense clouds and found that they lie below the observed Larson law by a factor of about two. Conclusions. In our simulations, magnetic field, shear, HII regions, and supernovae all contribute significantly to reduce the SFR. In this numerical setup with feedback from supernovae and HII regions and a relevant value of galactic shear, the SFRs are compatible with those observed, with a numerical incertitude factor of about two.


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 wellresolved 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. 2011Kim et al. , 2013Hill 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 10 4 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). A&A 620, A21 (2018) is the distance from the centre of the box to the centre 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. 2002Kim et al. , 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 Σ SFR ∝ Σ 2 gas . 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.

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 256 3 and 512 3 grid with no adaptive mesh refinement. For a 1 kpc 3 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.

Equations
We solve the ideal MHD equations with self-gravity, heating, and cooling processes, in a rotating frame (at Ω 0 = Ω(R 0 )), where R 0 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) 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 with L being the size of the box. Following Kim et al. (2002), assuming the angular velocity is a power law Ω(R) ∝ R −q , the effective tidal potential is We take the value q = 1 (flat rotation curve, Kim et al. 2002).
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): with K = 1.42 × 10 −3 kpc Myr −2 , F = 2.75 × 10 −4 Myr −2 , and z 0 = 180 pc (Joung & Mac Low 2006). Through a Poisson equation, the potential is associated with a density, Adding the Coriolis force, the MHD equations that are solved are with ρ, u, 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 L is such that n 2 L = n 2 Λ(T ) − nΓ, 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.

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 the galactic 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 ρ, v y , v z , and T are equal to those of the corresponding cell in the lower domain, shifted of −V shear × t in the x direction. For the x-velocity v x , it is not simply the value from the lower cell, but the necessary differential of velocity that is added: v x + V shear 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 (B x and B z ) 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 B y of the cell just below. Finally, the remaining unknown B y in the ghost cell is specifically chosen in order to ensure the zerodivergence 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 yboundaries 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 with α and z 0 to be determined at each timestep. This approximate φ app will be used in the boundary condition Via the Poisson equation, this potential is associated with an approximate density By integrating this density, α is deduced from the total mass of gas M in our box, To estimate the scale height, z 0 , we simply find the altitude at which the density is given by Finally, the boundary condition used in the y and z directions is 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 30-40%, 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.

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, with n 0 = 1.5 cm −3 and z 0 = 150 pc. The initial column density is then where m p = 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,  with B 0 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 where c s is the sound speed of the gas. This disk is stable to axisymmetric perturbations if 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.

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 10 3 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.  (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 v disp = 1 km s −1 in the sink frame throughout its life. Simulations with a lower or higher density threshold, or with v disp = 0 have also been performed, and are discussed in Sect. 3.3. The supernova then consists of an injection of 4 × 10 43 g cm s −1 (Iffrig & Hennebelle 2015) of radial momentum around the place of the star. We manually limit the temperature to T sat = 10 6 K and the velocity provided by the supernova to V sat = 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): T sat = 10 7 K and V sat = 1000 km s −1 , to analyse their influence.   (21) with the data of Vacca et al. (1996).
The maximum of velocity V sat 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: We then have the asymptotic relations 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). Figure 4 shows the column density in two of the MHD simulations (at a resolution of 256 3 ), 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. 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 by 20%. 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.

Global considerations
There is a 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)    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. 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.
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 V shear = 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 V shear .
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 V shear = 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 V shear = 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.

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 T sat and V sat 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. 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 v disp = 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 v disp = 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 256 3 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 512 3 . 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 (V shear = 28 km s −1 kpc −1 ), except for the resolution (256 3 vs. 512 3 , 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 V shear = 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 (10 3 cm −3 ). This threshold value was not well defined. Kim & Ostriker (2017) argue that it should vary as where ∆x is the size of one cell. The quantity ∆x 2 is then proportional to the surface of the cell. Figure 9d displays a new simulation (256 3 , green, with n sink 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 , 10 3 cm −3 , and 4 × 10 3 cm −3 . The one with the higher n sink (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.

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 V shear = 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 are formed, 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 V shear = 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 V shear = 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 n sink . The three curves are associated with the three previous values of n sink : 250 cm −3 , 10 3 cm −3 , and 4 × 10 3 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.

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 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: 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, I 11 = (y 2 + z 2 ) ρdV, I 22 = (x 2 + z 2 ) ρdV, I 33 = (x 2 + y 2 ) ρdV, I 12 = I 21 = − xy ρdV, Figure 12 displays the velocity dispersions of the clumps versus their sizes, for a velocity gradient V shear = 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).

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 V shear = 28 km s −1 kpc −1 , the SFR was too high compared to the observed Kennicutt value. The reference value V shear 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 massweighted average of Σ gas , the usual surface-weighted average Σ gas,mean (which respects mass conservation), and the massweighted mean velocity gradient V shear . 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 V shear 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 V shear = 28 km s −1 kpc −1 but about V shear ∼ 63 km s −1 kpc −1 . As previously stated (Figs. 7d and 9c), the SFR is very close to the observed one for V shear = 56 km s −1 kpc −1 . That shows that density and shear must not be taken independently, and that with a  relevant value of V shear 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 V shear in a future work.

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 M gas /(M gas + M stars ) 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. Notes. The values of Σ gas used to define these bins are mass-weighted averages whereas Σ gas,mean is the usual surface-weighted average. The mean V shear is mass-weighted. As can be seen, higher V shear are bound towards higher Σ gas . Section 4 gives more details.  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 .
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.

Conclusions
We have performed a set of simulations of 1 kpc 3 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 V shear = 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: n sink ∝ 1/∆x 2 . 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 some energy injection from the large galactic scales that drives the turbulence. We would expect such an injection to reduce the SFR even more.