On the origin of the peak of the stellar initial mass function: exploring the tidal screening theory

: Classical theories for the stellar initial mass function (IMF) predict a peak mass that scales with the properties of the molecular cloud. In this work, we explore a new theory proposed by Lee and Hennebelle. The idea is that the tidal field around first Larson cores prevents the formation of other collapsing clumps within a certain radius. The protostar can then freely accrete the gas within this radius. This leads to a peak mass of roughly 10 M 1LC , independent of the parent cloud properties. Using simple analytical arguments, we derive a collapse condition for clumps located close to a protostar. We then study the tidal field and the corresponding collapse condition using a series of hydrodynamic simulations with self-gravity. We find that the tidal field around protostars is indeed strong enough to prevent clumps from collapsing unless they have high enough densities. For each newly formed protostar, we determine the region in which tidal screening is dominant. We call this the tidal bubble. The mass within this bubble is our estimate for the final mass of the star. Using this formalism, we are able to construct a very good prediction for the final IMF in our simulations. Not only do we correctly predict the peak, but we are also able to reproduce the high- and low-mass ends. We conclude that tidal forces are important in determining the final mass of a star and might be the dominant effect in setting the peak mass of the IMF. ABSTRACT Classical theories for the stellar initial mass function (IMF) predict a peak mass that scales with the properties of the molecular cloud. In this work, we explore a new theory proposed by Lee and Hennebelle. The idea is that the tidal ﬁeld around ﬁrst Larson cores prevents the formation of other collapsing clumps within a certain radius. The protostar can then freely accrete the gas within this radius. This leads to a peak mass of roughly 10 M 1LC , independent of the parent cloud properties. Using simple analytical arguments, we derive a collapse condition for clumps located close to a protostar. We then study the tidal ﬁeld and the corresponding collapse condition using a series of hydrodynamic simulations with self-gravity. We ﬁnd that the tidal ﬁeld around protostars is indeed strong enough to prevent clumps from collapsing unless they have high enough densities. For each newly formed protostar, we determine the region in which tidal screening is dominant. We call this the tidal bubble. The mass within this bubble is our estimate for the ﬁnal mass of the star. Using this formalism, we are able to construct a very good prediction for the ﬁnal IMF in our simulations. Not only do we correctly predict the peak, but we are also able to reproduce the high- and low-mass ends. We conclude that tidal forces are important in determining the ﬁnal mass of a star and might be the dominant effect in setting the peak mass of the IMF.


INTRODUCTION
Star formation is a long-standing problem in astrophysics. While a general idea of how stars form is well established, see e.g. the reviews by McKee & Ostriker (2007) and Krumholz (2014), it is unclear which processes determine their final mass. A star's mass sets other important quantities, such as the luminosity and whether it will end its life with a supernova. Recent work has established that supernova feedback is important for the formation and evolution of galaxies. On smaller scales, the formation of planets is influenced by the luminosity of the parent star. So, knowing the distribution of stellar masses at birth, also called the stellar initial mass function or IMF, and understanding its origin are not only crucial for a complete theory of star formation but also have an impact on both planet and galaxy formation theories.
The IMF has been measured quite extensively inside the Milky Way and appears to be remarkably similar for different molecular clouds (Bastian, Covey & Meyer 2010;O f f n e re ta l .2014; Dib, Schmeja & Hony 2017;Hopkins 2018). This hints that star formation is governed mainly by processes that are independent of ⋆ E-mail: tinecolman@hotmail.com the global cloud properties. Recently, much effort is being put into measuring the IMF in other galaxies (e.g. van Dokkum & Conroy 2010;Gunawardhana et al. 2011) though this is not an easy task due to the limited resolution. It is still debated whether the results from these observations indicate an IMF that is universal across galaxies or not (Elmegreen 2009;Kroupa et al. 2013;Hopkins 2018).
Since Salpeter (1955) first introduced the concept, many different analytic forms have been proposed, defining the IMF ξ as a probability function with several free parameters, which are determined by fitting the observations. The most commonly used is the Chabrier IMF, consisting of a lognormal and a power-law tail (Chabrier 2003(Chabrier , 2005, but several others exist (e.g. Kroupa 2001;Parravano, McKee & Hollenbach 2011). These different forms are compared in Fig. 1.
While it is commonly accepted that the high-mass tail is a power-law with a slope close to the Salpeter value, there is a large uncertainty at the low-mass end. All of these forms, however, show a peak around 0.1-0.3 M ⊙ , which indicates a characteristic stellar mass. In this work, we focus on explaining this observed peak in the IMF. This paper is structured as follows: In Section 2, we explore some of the existing explanations for the IMF peak mass. Typical for these is that they predict a scaling with global cloud properties, something which is in contradiction with observations and recent simulations. However, such a scaling is not predicted by a recent theory proposed by Lee & Hennebelle (2018b), which states that the peak mass is set by tidal screening of the first hydrostatic core (also called the first Larson core). This tidal theory is explored further in Section 3, where we give a simple and intuitive analytical description. In Section 4, we describe the simulations used to investigate this theory. In Section 5, we present the results from the simulations. We study the tidal field around newly born stars and try to build a predicted IMF, which is compared to the IMF measured at the end of the simulation. In Section 6, we discuss what our results mean in the context of a general star formation theory, as well as the caveats associated with our approach. We conclude the paper in Section 7.

THEORIES FOR THE IMF'S PEAK
Over the years, many people have come up with theories to explain the origin of a typical mass scale for stars. In this section, we briefly review several of these with a focus on understanding why the peak of the IMF lies around 0.1 M ⊙ . Other interesting related questions are as follows: What is the smallest mass of a star? What is the largest possible mass? However, this is beyond the scope of this paper.

The Jeans mass at the opacity limit
A strong property of the isothermal gas that defines the interstellar medium (ISM) is the absence of a characteristic mass scale owing to the scale-free nature of gravity. From the collapse condition of an isothermal sphere, one can derive the Jeans mass as with c s being the sound speed, given by and the free-fall time t ff is defined here as One sees immediately that for an isothermally collapsing gas, we cannot choose a unique characteristic density that would result in a unique characteristic mass. The larger the density of the gas during collapse and fragmentation, the smaller the Jeans mass. At very large densities, though, owing to the increased dust absorption coefficient, the gas becomes opaque to its own radiation. This is sometimes referred to as the opacity limit (Low & Lynden-Bell 1976;Rees1976). At this point, the fragmentation is halted, setting a scale below which no fragments form. It is possible to estimate the characteristic density at which this process occurs by requiring the optical depth of one Jeans radius to be unity For typical molecular clouds conditions in the Milky Way with T = 10 K, μ = 2.2, and κ dust = 0.1 cm 2 g −1 , we find the critical density that defines the opacity limit At this critical density, we obtain a unique Jeans mass of the order of Note that the critical density can also be derived using slightly more complicated arguments, leading to a very similar value for typical Milky Way conditions (Krumholz 2017). This mass corresponds to the smallest gas clumps that can overcome the pressure gradients and collapse. At this characteristic density, the gas evolution transitions from isothermal to adiabatic, which prevents the collapse of smaller fragments, even at higher densities. The Jeans mass at the opacity limit is obviously much smaller than the observed characteristic mass of the IMF. It does not even correspond to the minimal mass of a star, as the smallest of these fragments will not be able to collapse enough to reach stellar densities in their centres.

The mass of the first Larson core
Once one of these small stable fragments forms, it can accrete mass from the surroundings and build up a hydrostatic core, a very important concept in the theory of star formation. Larson (1969) was the first to show that during the collapse of an isothermal cloud, a small core forms in hydrostatic equilibrium. He called this the first hydrostatic core, now also known as the first Larson core. Its formation is triggered when the density reaches the critical value defined above, changing the nature of the collapse from isothermal to adiabatic, with the formation of an accretion shock at the boundary of the core. Accretion continues and more and more mass accumulating in the hydrostatic envelope, compressing and heating the central core. At some point, the centre reaches a temperature high enough for molecular hydrogen to dissociate. Most of the energy now goes into dissociation and the collapse resumes, roughly isothermal again. This second collapse leads to the formation of a second hydrostatic core, which is the actual protostar. The remainder of the first core is then quickly absorbed by the protostar, which continues to accrete mass from the surrounding envelope.
Following this scenario, one computes the mass of the first Larson core by solving the Lane-Emden equation for a hydrostatic sphere with an equation of state (EOS) P ∝ ρ Ŵ . The properties of the resulting polytrope will depend on the following three parameters: the polytropic index Ŵ, the internal density in the centre, ρ int ,and the density at the edge of the core, ρ ext . Setting Ŵ = 5/3, ρ ext = ρ crit (corresponding to the transition from the isothermal to the adiabatic regime), and ρ int = 10 3 ρ crit fully specifies the mass of the core around M 1LC ≃ 0.02 M ⊙ , in agreement with the numerical value Larson obtained from his more complex calculations. Many authors have since repeated and improved Larson's numerical experiments, but the mass of the first Larson core M 1LC is still found to be of the order of 0.02 M ⊙ and has proven to be very robust against variations in the initial gas and dust temperatures and other properties of the surrounding envelop (e.g. Vaytet et al. 2013;Krumholz 2017;Bhandare et al. 2018). Choosing a different EOS and critical density will slightly change the mass of the hydrostatic core. A smaller Ŵ, corresponding to a shallower slope in the adiabatic part of the EOS, results in a smaller M 1LC , while lower ρ crit results in an increased M 1LC (Lee & Hennebelle 2018b). The EOS varies little between different environments, which makes M 1LC a universal quantity.
When the second Larson core forms, it accretes the leftovers of the first Larson core on a very short time-scale. All protostars will therefore have a mass of about M 1LC ≈ 0.02 M ⊙ shortly after their birth. This is why the mass of the first Larson core is probably a good estimate of the mass of the smallest stars, but it is still one order of magnitude lower than the observed peak of the IMF. This indicates that the majority of protostars continue to accrete material from a surrounding envelope after they have been formed and that the typical mass reservoir in the envelope should be about 0.1 M ⊙ .

The Bonnor-Ebert mass
On scales larger than the first Larson core, the flow is fully isothermal and scale invariant. Extracting a characteristic mass is therefore quite difficult. One strategy is to look for hydrostatic envelopes emerging from the parent molecular cloud. Hydrostatic spheres in this case are described using the isothermal Lane-Emden equation, leading to a family of solutions, with here again the internal and external densities ρ int and ρ ext as the two parameters. A particular solution, known as the Bonnor-Ebert sphere, corresponds to the largest stable isothermal self-gravitating sphere, with still one remaining free parameter ρ ext (Bonnor 1956;Ebert1957). Its mass, the Bonnor-Ebert mass, is given by and looks very similar to the Jeans mass, except that now the density is the density at the external boundary. We still have to find the value for ρ ext . One possibility is to use the typical density of molecular clouds. With n H ≃ 100 H cm −3 , we obtain M BE ≃ 21 M ⊙ , a value too large to explain the characteristic mass of the IMF. Moreover, molecular clouds are turbulent, strongly structured, and have different sizes, resulting in a wide range of densities and thus a wide range of possible Bonnor-Ebert masses, in contrast to the requirement of having a universal characteristic mass. A variation of this idea is based on the concept of turbulent Bonnor-Ebert spheres (Haugbølle, Padoan & Nordlund 2018). We know indeed that the ISM is turbulent and follows an observed universal scaling relation, known as Larson's relation (Larson 1981;Heyer et al. 2009), which related the velocity dispersion σ to the diameter of a cloud L where σ 0 = 1k ms −1 and L 0 = 1 pc are constant parameters valid throughout the Milky Way. This scaling relation is typical for supersonic (or Burgers) turbulence, in which strong isothermal shocks compress the gas to much higher densities than the initial cloud mean density. Hydrostatic envelopes will emerge from these compressed regions and deliver the envelopes to be accreted by the first Larson cores. Using the Rankine-Hugoniot relations for isothermal shocks, we can derive the external density of the Bonnor-Ebert sphere as where ρ is the mean cloud density and M is the Mach number of the turbulence, defined by The resulting Bonnor-Ebert mass still depends on the cloud density, but if one assumes that star-forming molecular clouds are in rough virial equilibrium, one has which results in a constant post-shock density with Injecting this in the Bonnor-Ebert mass formula finally gives for the turbulent Bonnor-Ebert mass (adding superscript 'turb') whose numerical value spectacularly fits the observed characteristic mass of the IMF (Haugbølle et al. 2018). One problem of this approach is that other galaxies have different global turbulent properties with different values for L 0 and σ 0 , resulting in a different IMF characteristic mass. A fair question to ask also is whether the formation of hydrostatic envelopes is at all possible within a selfgravitating, turbulent isothermal fluid. A more common situation in simulations of isothermal turbulence is a free-falling envelope around first Larson cores. Moreover, the numerical experiments we perform in this paper do not seem to support this scenario.

Turbulent fragmentation theory
The turbulent fragmentation scenario, introduced by Padoan, Nordlund & Jones (1997) and further developed by Hennebelle & Chabrier (2008) and Hopkins (2012), provides a more realistic framework for the origin of individual star-forming regions. This theory builds on the result that turbulence creates overdensities that follow a lognormal density PDF with δ = log(ρ/ρ)a n dδ =−σ 2 0 /2. The width of the distribution depends on the Mach number where the Mach number is defined at the scale of the entire molecular cloud and b is a parameter related to the nature of the turbulence forcing, either compressive with b = 1 or solenoidal with b = 1/3. The spectrum of turbulence is here again described using Burgers theory, with velocity fluctuations depending on each scale ℓ as where ℓ sonic is the sonic length,definedby corresponding to scales where turbulence becomes subsonic and density fluctuations become small. Using the analytical Press-Schechter theory, Hennebelle & Chabrier (2008) and Hopkins (2012) computed the mass fraction in collapsed structures, using a scale-dependent collapse criterion based on the virial parameter, which writes The excursion set formalism allowed them to compute the molecular cloud mass function using the first crossing of the collapse barrier, while the molecular core mass function was obtained using the last crossing of the collapse barrier. In other words, molecular clouds are the largest bound objects in a galaxy, while molecular cores are the smallest bound objects, at the bottom of the gravitational fragmentation hierarchy. Interestingly enough, this analytical formalism allowed these authors to compute the characteristic mass of the molecular cores as the mass of the smallest gravitational bound objects whose sizes are equal to the sonic length and the density is equal to the minimum density for gravitational collapse. Using equations (18) and (46) of Hennebelle & Chabrier (2008), we obtain This peak mass, under Milky Way conditions, is close to the value of the turbulent Bonnor-Ebert mass, although it has a slightly different dependence on the turbulent flow parameters. The numerical value we quote here was obtained using Larson's relation, b = 0.3 and M = 10.

Competitive accretion theory
A different theory was proposed by Bonnell et al. (2001a), who argued that observationally most young stars are seen in star clusters, sometimes still embedded in their parent gas cloud. The stellar density can be high, leading to dynamical interactions causing star ejections, mass segregation, and various scenarios of secular evolution. In this complex dynamical environment, young protostars, in particular the first Larson cores, interact with one another and compete for accreting gas. This scenario, called competitive accretion, is based on early first Larson core accreting gas at the core of their parent gas envelope, which is ejected after some typical interaction time by a more massive core, so that the accretion abruptly stops. As a result, smaller stars are the one ejected early, while more massive stars are the ones that managed to remain in their envelope until the exhaustion of the full gas reservoir. One can estimate the characteristic mass of the resulting IMF in the competitive accretion scenario, using whereṀ acc is a typical accretion rate and t dyn is the dynamical interaction time-scale. For the latter, one uses the local gas dynamical time-scale in the gas envelope t dyn = 3π 32Gρ env (22) while for the accretion rate, one can use the singular isothermal sphere model that giveṡ This leads to a characteristic mass equal within a factor of 2 to the Bonnor-Ebert mass. Here again, the difficulty is to adopt the right value for ρ env , using either the cloud mean density or the typical post-shock density emerging from the supersonic turbulence. For the latter, more realistic scenario, the resulting characteristic mass is half of the turbulent Bonnor-Ebert mass with exactly the same scaling than for the other previously discussed scenarios. Note that a more complicated derivation based on Bondi-Hoyle accretion was performed by Bonnell et al. (2001a), leading to a similar result, in very good agreement with the observed IMF.

Tidal screening theory
The previously discussed theories are based on supersonic selfgravitating isothermal turbulence, and always lead to a characteristic mass close to the turbulent Bonnor-Ebert mass. As a consequence, one expects the IMF to vary widely with the two parameters controlling the turbulence, namely the injection scale L 0 and the total velocity dispersion σ 0 , and a third parameter, namely the molecular gas sound speed c s . It is quite possible that the Universe conspires to keep the turbulent Bonnor-Ebert mass constant across very different environments. On the other hand, the mass of the first Larson core depends only on the thermodynamics of the gas and dust in the deepest and densest regions of molecular cores, which is expected to be universal for clouds with the same composition. This makes the first Larson core a very robust candidate for explaining why the characteristic mass of the IMF would remain constant in the Milky way.
In a recent series of papers, Lee & Hennebelle (2018b proposed to explain the characteristic mass of the IMF by tidal forces caused by the first Larson core and its envelope, preventing the formation of new cores in their surroundings. This allows the central core to gather all the mass available in this tidally screened reservoir, without sharing it with the other cores that would have collapsed otherwise. A rough estimate of the tidal screening process can be made using a simplified version of the virial theorem, applied to a small spherical region in the vicinity of an existing first Larson core. A more rigorous derivation follows in the next section.
We compare the tidal field of the first Larson core to the selfgravity of a collapsing gas clump by writing Assuming that gas around the newly formed core is distributed in an envelope described by a singular isothermal sphere, with we can deduce the radius within which the tidal field is strong enough to prevent the collapse as Integrating the total mass enclosed within this tidally screened region, we get independent on the envelope's properties (this is only valid when the exponent of the envelope profile is −2). Note that this very naive estimate is quite close to the numerical value found by Lee & Hennebelle (2018b), which was M tidal ≈ 10 M 1LC ≃ 0.2M ⊙ and also quite close to the characteristic mass of the observed IMF. Note that it is completely independent of any parent cloud physical parameters or turbulent galactic scaling relations. The goal of this paper is precisely to investigate this tidal screening theory with high-resolution simulations, using a complementary approach to Lee & Hennebelle (2018b).

TIDAL DISRUPTION THEORY
In this paper, we study specifically tidal screening of the first Larson core as the origin of the peak of the IMF. The mathematical aspects of the tidal screening theory are quite complex. For the original calculations, we refer to Lee & Hennebelle (2018b). Here, we adopt a slightly different point of view, exploiting directly our simulation data and highlighting the physical mechanism at work.

Model and assumptions
We consider here a basic configuration with a protostar (or a first Larson core) embedded in its gaseous envelope described by a power-law density profile In reality, star-forming cores are not spherically symmetric. Moreover, the parent clouds are highly turbulent and the envelopes are expected to contain filaments and voids. Nevertheless, we consider that a small, spherical clump with a uniform density ρ clump is forming at a distance r from the protostar in its spherically symmetric envelope. Certain conditions are required for this clump to be able to collapse and form a star. The fundamental idea is that due to the already existing protostar (and its power-law envelope), a tidal force is exerted on the clump, stretching it in the radial direction. If this stretching is too strong, the clump will be tidally stripped and will be prevented from forming a companion star. In this case, the gas in the envelope has no other option than to be accreted on to the central protostar. This process is called tidal screening and the gas mass in this tidally stripped region will determine the final mass of the star.

Tidal tensor
In general, the tidal tensor is defined as with φ the gravitational potential and g the gravitational acceleration. We have used the relation in the second equality. The derivatives commute so the tidal tensor is symmetric and can be diagonalized to the form with λ i the three eigenvalues. The sign of an eigenvalue indicates whether a local gas clump is expanding (λ i < 0) or contracting (λ i > 0) in the direction of the corresponding eigenvector. In the case of a spherically symmetric potential in spherical coordinates, the tidal tensor is given by

Collapse condition
Our clump will feel a total tidal force that is the sum of three contributions: the central protostar, the envelope, and the clump itself. We now derive a condition for the clump to collapse. We first need to compute the gravitational field of the protostar. It is simply that of a point source resulting in a tidal tensor of The tidal force is strongest in the direction towards the star. The clump is stretched in the radial direction and contracts in the two tangential directions, turning a spherical object into an elongated ellipsoid.
The gravity of the envelope with density profile from equation (30) is given by for α<3. Its tidal tensor thus becomes If α = 2, this generates a contraction in both the tangential directions, which is equally strong to the expansion in the radial direction.
The clump itself is considered to be homogeneous, with which results in a tidal tensor that is proportional to the clump density and contraction equally strong in all directions: Adding up these three contributions, the non-zero components of the total tidal tensor are which are also its eigenvalues. As noted before, a positive eigenvalue implies collapsing in the associated direction, and a negative one means expanding. To resist the effect of tidal stripping, and thus to be able to form a star, all three eigenvalues must be positive. The tangential eigenvalues cannot be negative for 0 ≤ α<3, so in order to be collapsing a clump must satisfy the condition T rr > 0or

Tidal radius and tidal mass
Once a first Larson core has formed inside its envelope, we know M * and ρ gas (r). For a given clump density, the collapse condition is a function of the distance to the central protostar. One can define the tidal disruption radius r tidal as the radius for which the inequality in equation (43) becomes a strict equality. In other words, if the clump with density ρ clump is located within r tidal , it will be tidally stripped. If it is outside this radius, it will resist the tidal forces. As a consequence, all the gas mass within the tidal radius will eventually end up in the central object. We call this mass the tidal mass, which gives us an estimate for the end mass of the protostar Contrary to the other theories discussed so far, it is difficult to estimate this tidal mass analytically. The typical clump density can be estimated using gravitational fragmentation theory, considering the smallest collapsing clumps, namely clumps with radius equal to the sonic length. Under typical Milky Way conditions, this gives ℓ sonic ≃ 0.04 pc, and α vir < 1 translates into ρ clump > 3.8 × 10 −18 gcm −3 . As discussed earlier, the first Larson core has a well-defined mass with M 1LC = 0.01 M ⊙ . The gaseous envelope properties, on the other hand, emerge from the supersonic turbulence in the cloud, and are not uniquely defined. We assume here that the envelope is a strict singular isothermal sphere, connected to the first Larson core with the Jeans density at the Jeans radius, both defined in Section 2.1, with α = 2, ρ 0 = ρ J = 5 × 10 −14 gcm −3 ,andr 0 = R J ≃ 12.5 au.
If we first ignore the tidal forces from the envelope, we can apply the collapse condition from equation (43), and find r tidal = 656 au for ρ clump = 10 −17 gcm −3 and r tidal = 305 au for ρ clump = 10 −16 gcm −3 . We see that denser clumps can survive tidal stripping closer to the first Larson core. If we consider now only the gaseous envelope, we find similar values, with r tidal = 1531 au for ρ clump = 10 −17 gcm −3 and r tidal = 484 au for ρ clump = 10 −16 gcm −3 . Overall, this gives a rather large range of values between 300 and 1500 au. The corresponding tidal masses range from 0.05 M ⊙ < M tidal < 0.25 M ⊙ , in good agreement with the characteristic mass of the IMF. A complete view of the parameter dependence of the tidal radius and tidal mass is shown in Fig. 2. For envelope profiles with a slope of α = 2, we recover the previous estimates, with tidal masses around 0.1 M ⊙ . For shallower profiles with the same normalization, it is much more difficult for clumps to collapse close to the centre. The tidal mass is thus larger. The opposite is true for steeper profiles, where the clump density quickly becomes much larger than that of the background envelope. In this case, the tidal mass typically becomes smaller than 0.1 M ⊙ .
Following Lee & Hennebelle (2018b), we can also consider that the clump density is proportional to the gaseous envelope, with Injecting this into equation (43), we see that a condition for the clump to collapse is A> 3(α−1) (3−α) . Only clumps significantly denser than the background gas envelope can collapse, while less dense clumps are always tidally stripped. If we consider again the singular isothermal sphere and choose A = 4 as a minimum clump overdensity, the collapse condition from equation (43) becomes identical to our simple derivation of Section 2.6 and the tidal mass becomes M tidal = 6M 1LC .
The arguments discussed here give us some analytical estimates for the tidal mass. Our next step is to use simulations to try and estimate the tidal mass directly from the actual gas distribution.

SIMULATION SET-UP
To explore the idea that tidal forces play an important role in setting the peak of the IMF further, we run a suite of turbulent box simulations. The goal is to study the tidal stripping in a realistic turbulent setting. The simulations are designed to form many stars while maintaining a resolution for which individual star can be distinguished. This allows us to investigate specific locations of star formation in detail, as well as providing good statistics on the cloud scale. The simulations are performed with the Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002) taking into account hydrodynamics, self-gravity, and a polytropic EOS to mimic the effect of radiation. Sink particles are used to represent stars. We choose a base refinement level of 8 covering a periodic box of size 0.25 pc and a maximum refinement level 12, which gives us a maximal spacial resolution of 12.6 au. This maximum resolution corresponds exactly to R J , the Jeans length at the opacity limit. This is also roughly the size of the first Larson core. In Appendix B, we show that this resolution is good enough for the purpose of this study and that our results are converged. We refine the grid based on a Jeans-length refinement criterion, namely that we always resolve the local Jeans length with at least four cells until we reach the maximum level of refinement.

Initial conditions
Our simulation set-up consists of a periodic box of size 0.25 pc containing 260 M ⊙ with fully developed turbulence with a Mach number M 1D = M 3D / √ 3 ≃ 7.3. In this simulation suite, we use five initial conditions, which only differ by the random seeds used to generate the Gaussian random field for the turbulent velocity. All the other average properties are strictly identical. We will refer to these as box 1 to box 5. Using different randomized turbulence gives us different large-scale structures in each box. This allows us to examine many configurations with different tidal fields and gives us an idea of how the direct environment of the protostar can influence its growth. For details about how the initial conditions are generated, we refer to Appendix A.
Note that for the sake of simplicity, we have no external turbulence driving force. During the run, the turbulence only evolves through the combined effects of hydrodynamics and selfgravity. Note also that our box is quite far from typical Milky Way conditions. The observed Larson's relation exhibits quite some scatter and inside a molecular cloud there are many subclumps and voids. This small 0.25 pc box should be considered as a small, particularly high-density region within a much bigger molecular cloud for which the 1D velocity dispersion approaches 1.5 km s −1 . High Mach numbers boost star formation, since more high density fluctuations are created. A small box size allows us to use high resolution. Using a periodic box without driving means that this piece of cloud is considered to be relatively decoupled from the parent cloud. This simplified set-up is well suited for the purpose of our experiment, which is to focus on the small mass end of the IMF and investigate the effect of local tidal forces around first Larson cores.

EOS
Initially, the temperature is set to T 0 = 10 K everywhere in the box. Since we do not include radiative transfer in our simulations, we adopt instead a simple recipe based on a polytropic EOS For densities below a certain threshold ρ knee , the gas is strictly isothermal. Above this density, T(ρ) follows a power law with a slope Ŵ − 1. This mimics the fact that, when the optical depth becomes higher than 1 and the density exceeds the critical density at the opacity limit, radiation gets trapped and the evolution becomes adiabatic. We use a value of ρ knee = 5 × 10 −14 gcm −3 or n knee = 2.5 × 10 10 Hcm −3 consistent with results found in RHD simulations (Vaytet et al. 2013). After reaching these densities, the adiabatic phase consists of two parts. At first, the gas evolves with a polytropic index of 5/3 until a temperature of about 100 K is reached and rotation modes of molecular hydrogen become excited. After this, it heats up further with a shallower index of 7/5 (Bhandare et al. 2018). Since we do not reach these very high densities and the index close to ρ knee is shown to be most important in setting the mass of the first Larson core (Lee & Hennebelle 2018b), we neglect the second phase and simply use Ŵ = 5/3. It has already been pointed out that the EOS is important in shaping the IMF (Kitsionas et al. 2007). Some form of polytropic EOS is certainly needed, since it was shown that pure isothermal simulations produce an IMF peak that directly depends on the resolution (Guszejnov et al. 2018). Without a transition to the adiabatic regime, there is no mechanism that allows for the formation of first Larson cores with a well-defined physical mass.

Sink particle formation
When a density peak reaches conditions for which it is expected to form a first Larson core that will trigger the second collapse towards the protostar, we introduce a sink particle (Bate, Bonnell & Price 1995;Bleuler & Teyssier 2014). Because our simulations are highly resolved, we believe sink particles can represent individual stars (see Appendix B). Resolving the second collapse is, however, quite demanding in terms of spatial resolution. It has been shown by Lee & Hennebelle (2018b) that the internal structure of the first Larson core allows one to adopt a relatively low density threshold, slightly larger than the critical density corresponding to the opacity limit, and still get the right statistics for the second collapse sites. The main requirement for us is thus to resolve the spatial extent of the first Larson cores, detect when the peak density of the first Larson core exceeds a density threshold ρ threshold ≥ ρ crit , and introduce a sink particle. These are the key elements that allow us to investigate the true IMF (instead of binary or cluster IMF).
The details of our sink particle formation criteria are as follows. Sinks are formed within density peaks, called here clumps, found by the native RAMSES clumpfinder PHEW (Bleuler et al. 2015). Individual cells are therefore not all eligible for sink formation, only cells located at density maxima. A clump must fulfil the following criteria to be eligible for sink formation: (i) it has to be dense enough: ρ peak >ρ threshold ; (ii) it has to be relevant enough: ρ peak > 2ρ saddle (see Bleuler et al. 2015, for the exact meaning of relevance); (iii) the sphere of radius R acc and centred on r peak has to contain a mass larger than M seed ; (iv) it is not kinetically supported: 2E kin < −E grav ; (v) it is not thermally supported: 3PdV < −E grav ; (vi) there is no existing sink nearby: |r peak − r sink | > 2R acc .
We remark that we compare the viral terms separately as first Larson cores are virialized: 2E kin + 3PdV + E grav ≈ 0. The density threshold and seed mass are chosen to match the conditions in the first Larson cores that trigger the second collapse. We adopt here and at our fiducial resolution ρ threshold = 5 × 10 −14 gcm −3 and M seed = 5 × 10 −3 M ⊙ . Once a sink is formed, it will accrete gas at the Bondi rate (Bleuler et al. 2015), with all Bondi flow parameters estimated within the accretion radius R acc = 50 au, or equivalently four cells at the maximum level of refinement. In practice, after the sink is created with an initial mass equal to M seed ,t h er e s t of the first Larson core mass is quickly accreted and the pressure support disappears, mimicking quite realistically the second collapse.
Another important parameter in our sink particle algorithm is the sink merging time-scale. New sinks are allowed to merge with another older sink only during a time-scale t merge after their birth (Bleuler et al. 2015). This time-scale is associated with the first Larson core lifetime, and is usually chosen around 1000 yr, which corresponds to the dynamical time of the first Larson core (Bhandare et al. 2018). During this time-scale, we expect the first Larson core to be still gaseous and a collision with another core will lead to one protostar rather than a binary. If the merging time-scale is chosen too small, one can get too many small-mass stars in multiple systems. If the chosen value is too high, not a single binary will form. In our simulation, this parameter is set to t merge = 1500 yr. Choosing any value within a factor of 3 from this gives similar results (see Appendix B).
The resolution of our fiducial set of simulations is fixed to 12.6 au, corresponding to AMR level 12. In Appendix B, we perform a resolution study to quantify the robustness of our results. In this study, we adopt the same polytropic EOS, so that the properties of our first Larson cores are preserved. We do, however, adjust the sink particle formation parameters, so that the density threshold for sink formation is larger and closer to the onset of the second collapse, and the seed mass at sink creation and the sink accretion region are smaller with increasing resolution. The conclusion of this resolution study is that our fiducial case (ℓ max = 12) is properly converged.

Evolution of the clouds
The evolution of the five clouds is visualized in Fig. 3. In an initial stage, large-scale overdense regions start to collapse under gravity, forming filaments, which are the locations of the formation of the first sinks. More and more stars form as the simulation continues and the mass in the box is slowly depleted. Fig. 4 shows the evolution of the total mass in stars compared to the initial cloud mass, also called the star forming efficiency (SFE) Once sink formation starts, the SFE initially evolves roughly parabolic and later transitions into a linear regime. What is striking is that the different SFEs vary by a factor of up to 3 for our different random seeds, even though the boxes have the same exact global mass, size, and Mach number. Box 5 is most notable, with an SFE much higher than all the other boxes at all times. As one can see in the density maps (Fig. 3), the primary locations of sink formation are the filaments formed from largescale overdensities. We see that box 5 has one massive filament in the centre of the box, while in box 1, for example, there are many small filaments that are spread out in a more uniform way over the entire volume. The presence of large-scale filaments thus has a strong effect on the overall star formation rate.
The origin of this difference in geometry might be that the turbulent energy is in different modes for the various initial conditions. It is indeed well known that star formation is more efficient in the presence of strong compressive modes, than for pure solenoidal turbulence (Federrath et al. 2010;Federrath & Klessen 2012;Orkisz et al. 2017) but this is not the topic of this paper, so we do not elaborate further. We only use this diversity of compressive and solenoidal modes to deliver a variety of possible environments in our star-forming regions. Fig. 5 shows the number of sinks as a function of mass for the five different boxes. Since in our simulations, a sink particle represents an individual star, this is equivalent to the IMF. Even though they have very different SFEs, for each box there is a clear peak around 0.1 M ⊙ . To improve the statistics, we can stack the data of the five boxes together and plot the resulting mass function, which is shown in Fig. 6. The stacking is done at different times, when the simulations all have the same SFE, to not let one simulation dominate the result. Also here, after stacking, we see a clear peak around 0.1 M ⊙ . We remark that the position of the peak remains roughly constant during the simulation. The exact moment we evaluate the IMF is thus not important. This also shows that the mechanism that sets the peak in our simulation does not depend on the evolution of the cloud.

The simulated IMF
Comparing our sink mass function with some of the commonly adopted analytic forms of the IMF (Fig. 6, bottom panel), we see that the simulated peak agrees quite well with the observed one, especially for the Kroupa and Chabrier IMFs. This indicates that even though our simulation only contains gravity and hydrodynamics, with a simplified EOS to account for radiation, we manage to capture a physical process that sets a clear peak in the IMF. Note that the agreement between our simulated results and the observed IMFs is not so good at the high-mass end. This might be due to our limited box size and total mass, giving us poor statistics at the high-mass end. It has also been shown that radiation and magnetic field are key ingredients for the formation of massive stars (Hennebelle & Commercon 2012;Tan et al. 2014). Since we do not include these in our simulations, it is quite possible that the most massive cores fragment too much and do not allow for the formation of massive enough massive stars. Another likely explanation comes from the recent work of Padoan et al. (2019). They claim that massive stars gather their mass through large-scale, converging, inertial flows. Since our set-up is a small periodic box, we do not have these largescale flows. Box 5 resembles closest this picture, as most of the gas happens to be in the centre of the box and a global collapse takes place. And indeed, as seen in Fig. 5, only in this box we find stars with a mass significantly higher than 1 M ⊙ . The fact that we do not reproduce the high-mass end of the IMF is not important here, as the goal of this paper is to focus on the characteristic mass of the IMF. On the low-mass end, however, we have a clear cut-off at a mass slightly lower than the first Larson core mass. This could be associated with our adopted numerical seed mass, which is, for our fiducial runs, only half the first Larson core mass. We show, however, in Appendix B that this cut-off persists even with a much lower seed mass (see Fig. B6).
Can any of the theories discussed in Section 2 explain the characteristic masses we observe in these simulations? First, we see very few sinks with masses lower than the first Larson core mass, around 0.01 M ⊙ , so this characteristic mass seems to corresponds to a lower limit for the mass of stars. The high-resolution run shows a few sinks with mass below the first Larson core mass, because our seed mass is significantly smaller (by a factor of 4) in this case. But these very low mass sinks are rapidly accreting their surrounding first Larson cores and quickly reach the low-mass end cut-off. Given our initial condition parameters, the Bonnor-Ebert mass is M BE = 0.26 M ⊙ , where we use the average box density M box /L 3 as the external density, namely ρ ext = 1.1 × 10 −18 gcm −3 . This is almost a factor of 3 too large to explain our peak value, which lies around 0.1 M ⊙ . As discussed in Section 2, a better approach in such turbulent environments is to use the compressed post-shock density as external density. This boils down to dividing the previous mass by the 1D Mach number of the turbulence M ≃ 7.3. We now obtain M turb BE = 0.03 M ⊙ , which is a factor of 3 too low. The reason for this disagreement comes from our adopted cloud properties that do not match the typical cloud conditions in the Milky Way. Our set-up is more compact, or for a fixed size, more turbulent. But this analysis also shows that the origin of the peak in the IMF in our simulations cannot be the turbulent Bonnor-Ebert mass. Most of the other theories discussed in Section 2 give similar characteristic masses, and therefore fail to explain our simulations. The only theory that could still work here is the tidal screening theory.

Studying the local tidal field
The tidal screening theory states that because of the presence of large tidal forces around young protostars, no other stars can form within a certain volume, called the tidally screened region. To study this in our simulations, we look at the tidal tensor and its eigenvalues in the envelope around newly formed sink particles.
In each AMR cell, we compute the three eigenvalues and the corresponding eigenvectors of the tidal tensor, using the computed gravitational force stored in the RAMSES outputs. If the eigenvalue is negative, we know that the gas at this location is currently being stretched in the direction of the corresponding eigenvector. If the eigenvalue is positive, the gas is collapsing along the direction of the corresponding eigenvector. To be able to collapse and form a first Larson core, a clump of gas needs to be collapsing in all three directions, otherwise it will simply turn into a long thin spaghetti, fall towards the centre of the envelope, and be accreted by the older first Larson core sitting there.
If we order the eigenvalues, it is thus enough to look at the smallest (or most negative) eigenvalue. If this one is positive, we know that the others will also be positive and the gas parcel will collapse. We can express the eigenvalue in units of density by dividing by 4π G/3. This is handy, because if we were to place a small homogeneous clump of certain density in the tidal field, we can immediately see whether it will be tidally stripped or not, since this simply adds a constant positive term to the most negative eigenvalue (see equation 41).
Using the definition of the tidal tensor (equation 31), the most direct approach is to calculate the spatial derivatives of the gravitational acceleration. Taylor expanding the acceleration around a fiducial point x 0 , we can write For each grid cell in the gas envelope around the first Larson core, we estimate the local tidal tensor's components T ij using a least-square approach on the closest neighbouring cells, using the relative positions and accelerations and fitting the above equation in each spatial direction. We then symmetrize the obtained tensor and calculate its eigenvalues using NUMPY's LINALG.EIGH. An example of this tidal tensor eigenvalue analysis for an isolated sink is shown in Fig. 7. The map shows the smallest (most negative) of the three eigenvalues in a 2D slice around the sink. The red colours are for a negative eigenvalue, corresponding to tidal stripping in the direction of the eigenvector, while the blue colours are for a positive eigenvalue, corresponding to fully collapsing fluid elements.
Looking at the map, we can see the presence of a strong tidal field in a region of several hundreds of au (up to 1000 au in diameter) around the sink particle, which is located in the centre of the image. The three grey contours mark constant value of −10 −16 , −10 −17 ,and−10 −18 gcm −3 , respectively. If we were to place a new, potentially collapsing clump within one of these contours, it would need to have a density higher than this value to be able to collapse and resist tidal stripping. Closer to the sink, the density has to be so high to overcome the tidal forces that it is virtually impossible for a potentially collapsing clump to form a star. These measurements confirm our analytical calculations in Section 3.4, where we have estimated the tidal radius and the corresponding tidally screened mass for typical potentially collapsing clump densities.

Computing the tidally screened mass
In the previous section, we have computed the tidal field around an isolated first Larson core, and showed that potentially collapsing clumps would be tidally stripped and will not collapse. In this section, we want to go one step further and compute the actual mass that is tidally screened around the same first Larson core. We will then apply the same method around all first Larson cores, right after they form during the simulation. Since this requires analysing all the simulation snapshots, we need a fast method to measure the tidal field around sink particles.

Spherical approximation of the tidal field
Instead of measuring the tidal force directly from the simulation like described in the previous section, we want to develop and test a simple model inherited from our spherical analysis of Section 3. The tidal field in our star-forming regions consists of the contribution of the sink and its surrounding turbulent envelope. We approximate the contribution of the gaseous envelope using a spherical model with a power-law density profile of the form ρ(r) = Ar −α .
To estimate the power-law slope and the normalization of the profile, we measure the enclosed mass at two different radii (R 1 = 300 au and R 2 = 3000 au). From this, we then find where M 1 and M 2 are the masses enclosed in radii R 1 and R 2 , respectively. Note that this spherical fit can be a poor approximation of the actual density distribution, if the sink is not in the centre of the envelope or the region is dominated by a filament. The corresponding gravitational acceleration is usually a much better fit, even in these unfavourable cases. The most negative eigenvalue of the tidal tensor is then The dashed circles in Fig. 7 compare this simplified model to the exact tidal field. We can see that it compares quite well with the exact measurement. More importantly, the volume enclosed by the isocontours matches quite closely to the actual tidally screened region, for a given clump density.

Applying the collapse condition on the actual density distribution
In our current simplified picture, we now want to estimate a local collapse criterion, using the gas density in each cell and comparing it to the eigenvalue of our spherically symmetric tidal tensor. We do not consider the density of a potentially collapsing clump virtually  . 2D slice of the gas density around an isolated sink, in the plane of the sink. The green circle in the centre represents the sink accretion radius. The tidal bubble is marked in black. The grey contours show the places where the exactly measured tidal eigenvalue tends towards collapsing (−10 −18 gcm −3 in Fig. 7). located at a given position like before, but we use instead the actual gas density at that position to compute the collapse condition described in Section 3.3. Usually, the dense filaments generated by the turbulence provide the overdensities within the envelope that will survive the tidal forces and collapse. The remainder of the gas in the envelope will be tidally stripped (negative collapse criterion) and will define the tidally screened region. Fig. 8 visualizes this collapse condition in a 2D slice around our isolated example. In this map, the gas that is expected to be tidally stripped is shown in light to dark orange, while the gas that can collapse is shown in light to dark green. When comparing with the density map in Fig. 9, we can see that collapsing regions closely follow the filaments surrounding the first Larson core. A large positive value for the collapse criterion (the units are again in density units) means that the gas will collapse in a short time-scale, while a large negative value corresponds to strong tidal stripping.

The tidal bubble and the corresponding tidal mass
Note that the previous analysis is performed right after the formation of the first Larson core. Our objective is therefore to predict at this very specific moment in time the mass that will be efficiently tidally stripped and end up in the central first Larson core, while the rest of the mass can eventually collapse and form other first Larson cores later on. We will call the region that is dominated by tidal stripping the tidal bubble. While a natural definition for the tidal bubble would be the region for which the collapse condition is negative, not all of this gas will make it to the centre. If the tidal field is only weakly stripping, a parcel of gas will not have the time to react to the effect before it is swept away by the turbulence or collapses on to another clump. In practice, we thus have a stricter criterion and a threshold value for the collapse criterion that is smaller than zero. In Fig. 8, we overplot several possible threshold values. We see that only for more negative thresholds, the tidal bubble is closed. The exact threshold value chosen is based on a time-scale argument, as we will discuss in Section 5.5.1. In Fig. 9,weshowa map of the gas density, with superimposed the region corresponding to the tidal bubble, shown as a black hatched area. The total gas mass within the tidal bubble is the tidal mass (including the initial mass of the first Larson core) and, in the tidal screening theory, corresponds to the final mass of the star. We have also indicated as a grey contour the collapsing regions identified using the exact tidal tensor analysis performed in the previous section. We see that they correspond nicely to all regions outside of the tidal bubble. As a caveat of our current approach, note that there is no guarantee that the gas within the tidal bubble at the formation time of the first Larson core will indeed end up in the final star, since the subsequent accretion will proceed in a highly dynamical environment.

Multiple systems and competitive accretion
In the previous sections, we have analysed the simple case of an isolated, newly born first Larson core. Many first Larson cores are, however, born in a much more complex environment, with multiple sinks forming a star cluster and sharing the same gas reservoir.
The first difficulty is to estimate the tidal field from the common gaseous envelope. A young core appearing in a crowded envelope will see its gaseous envelope severely depleted of its mass by the older sink nearby. The gas reservoir will therefore be significantly less massive compared to the first sink that formed in the same envelope. Its shape might also have changed. In order to estimate the envelope parameters, we apply the same technique as before, using as centre of the density profile the position of the newly born core.
The second difficulty is to compute the tidal field resulting from all the sink particles in the vicinity of the newly born first Larson core. Note that our goal is to assign to each new born core a tidally screened region, the tidal bubble. In case of a multiple system, we do so by assigning a gas cell to a tidal bubble depending on the relative strength of the different sink plus envelope tidal fields. Large sink particles with a strong tidal field will be assigned a larger tidal bubble.
In Fig. 10, we show an example of a binary system, with a younger core forming close to an older one. The young sink is shown as a green circle while the old one is shown as a red circle. The collapse criterion is shown on the top panel, and again, like for the isolated case, dense clumps and filaments are collapsing, while the rest of the envelope appears as tidally stripped. In the bottom panel, we show the map of the gas density in the same region, showing the tidal bubble as the black hatched area. We see that our analysis only assigns to the young sink the region that is dominated by the tidal field of the young sink. The rest is assigned to the old sink and will end up accreted by the old sink according to the theory.
The same caveats we discussed in the previous section apply here, probably even more. These often crowded and highly nonspherical environments are difficult to describe, so that complex numerical simulations are fully justified. The theory of competitive accretion, as discussed in Section 2.5, is probably very useful in such cases. Although we do not explore it in this paper, combining tidal screening theory and competitive accretion appears as an interesting idea. Overall, our current approach allows us to estimate roughly the final mass of the stars, by matching it with the tidal bubble at formation time. We will see in the next section whether this can indeed explain the numerical IMF we have obtained.

Predicting the IMF
Now that we are equipped with the previously described techniques to compute the tidally stripped mass in the region surrounding newly born first Larson cores (the region we called the tidal bubble), we can repeat the procedure for every new sink in the simulation.   Fig. 11 using different definitions for the edge of the tidal bubble. There is a bimodality for a collapse condition threshold close or equal to zero.

The edge of the tidal bubble
First, defining the edge of the tidal bubble as a collapse condition of exactly zero, corresponding to all the tidally stripped gas around the sink, we plot in Fig. 11 the relation that compares the tidal mass with the actual final mass for each sink for one of the simulations. We immediately see that there are two groups of sinks: the lowmass sinks, for which the prediction is very good, and the higher mass sinks, for which the mass is overestimated. Instead of taking a threshold of zero, we can use a stricter definition for the edge of the tidal bubble and choose a threshold for the collapse condition equal to some negative value. In Fig. 12, we show the corresponding mass functions and compare our predicted IMF when using different values for the tidal bubble threshold, corresponding to the contours in Fig. 8 and the top panel of Fig. 10. We see that there is a clear bimodality when using all the tidally stripped gas or a threshold Figure 13. The IMF of all sinks in the five simulations, categorized based on how many other sinks were close by at the moment of birth. One can see that sinks that will become low-mass stars formed in already existing clusters.
value that is close to zero to define the edge of the tidal bubble (light grey histograms). Interestingly, the minimum between the two peaks is located exactly at the position of the measured peak. When taking a stricter threshold for the tidal bubble edge, the highmass peak shifts to the left, until the bimodality disappears (dark grey histograms).
The presence of this bimodality indicates that there is a fundamental difference between the two groups of sinks. When we inspect several of these sinks more closely, we find that practically all sinks in the low mass group form in crowded regions deviating strongly from spherical symmetry. This situation is very dynamical, making it difficult for newly born sinks to accrete mass. The other members of the star cluster have also already partially depleted the gas in the shared envelope. In this case, the tidal mass is dominated by the effect of the neighbours. The fact that in this regime our prediction is quite good tells us that the procedure for dividing the mass between neighbouring sinks as discussed in Section 5.4.4 is reasonable. The sinks in the high-mass peak tend to form in relative isolation and can accrete undisturbed for a significant amount of time. The few sinks in between the peaks have configurations that are neither of the two extremes: isolated but with another star-forming region close by or forming in a binary or small star cluster. We can show this more quantitatively by dividing all sinks into groups based on their conditions of birth: sinks born in isolation, sinks born as part of a binary or triplet system, and sinks that form in clusters that already contain three or more members. In Fig. 13, we compare the sink mass function for each of these groups. This clearly shows that low-mass sinks formed in clusters. Note that isolated sinks do not always stay in isolation, as it is quite common for regions to merge during the course of the simulation.
A threshold for the collapse criterion close to zero clearly overestimates how much mass an isolated sink can accrete and thus how big the tidal bubble is. A threshold of −10 17 gcm −3 results in a peak in the same location as the measure one, but it also underestimates the number of high-mass sinks. As we can see in Fig. 11, high-mass sinks generally form early on in the simulation. They have more time to accrete mass and this needs to be taken into account. In order to estimate the mass that will be accreted by isolated first Larson cores, we need to apply this additional selection criterion within the tidally stripped region. We require the tidal stripping time-scale, evaluated at the birth epoch, to be shorter than the time during the initial and final ages of the sink particle. Converted in density units, this translates into This threshold defines a region in which gas is expected to be tidally stripped at a fast enough rate so that it can be accreted on to the sink before the end of the simulation. Gas experiencing weaker tidal forces is kept outside of the tidal bubble, as it will not react to the tidal field fast enough. Turbulent motions might have moved the gas parcel away from the sink or gravity might have made it fall on to another clump. We have adjusted the fudge factor to C = 0.5, which seems to give us the best match between the tidal mass and the final mass of the sink.

Final prediction for the IMF
We now finally analyse all the sinks in our five different simulations. In Fig. 14, we plot the comparison of the tidal mass as defined with the time-scale-dependent threshold, with the actual final mass for each sink. We see that both mass are well correlated, supporting the idea that the tidal mass, estimated at the birth epoch of the first Larson core, is a good proxy for the final mass of the star. The values of the Pearson correlation coefficient are listed in Table 1 and show that, although there is certainly some scatter, the correlation is very good. We also see in this figure, as pointed out before, the emerging of different families of sink particles forming at different epochs. There is a general trend that sinks that form early in the simulation end up having larger final masses, with values at or above the peak mass of the IMF. They are seen as blue circles in the top right part of each diagram. These stars are formed in isolation but typically end up being the most massive member of a relatively large star cluster. Sinks that form late in the simulation are shown as orange to dark red circles. They tend to have masses lower or equal to the peak mass of the IMF. These stars generally formed in relatively crowded environments or in already significantly depleted envelopes. Intermediate cases form the bulk of the stellar population in our simulations. They form in relative isolation, with typical zero or one neighbour. In Fig. 15, we plot the mass functions of our estimated tidal masses. They match quite nicely the measured final IMFs. Not only are the predicted peak positions quite close to the measured ones (except in box 5 where there is a larger offset), but also the highand low-mass ends are properly recovered. Since the individual box IMFs are quite noisy, in order to improve the statistics, we have again stacked them at the same SFE in Fig. 16. Here, we see that the agreement is quite good, with a small shift towards lower masses. This confirms our visual impression in Fig. 14 that our predicted masses using the tidal bubble are slightly (by 0.1 dex) underestimating the true final mass.

DISCUSSION
We now discuss our results in light of previous works on the same topics. We will also discuss caveats in our modelling approach, as well as remaining issues of the tidal screening theory.

Problems with established IMF theories
The fundamental issue in classical IMF theories is that the predicted peak mass depends on global cloud properties, e.g. the cloud average density and the turbulence Mach number. Using these theories, one would thus expect different values for the characteristic star mass if these global properties were to vary between different starforming clouds. Observations in the Milky Way show, however, a universal IMF for all clouds, even though they have a large variety of properties. As discussed in the beginning of this paper, the Larson relation comes to the rescue and imposes a magic cancellation of these possible variations. Whether the IMF is also universal in other galaxies is still unsure, since one expects the normalization of the Larson relation to differ widely in other galaxies. It remains unfortunately quite difficult to accurately measure the IMF, as well as the Larson relation, in distant galaxies (Pan et al. 2016;Conroy, van Dokkum & Villaume 2017;Gennaro et al. 2018). A recent theoretical study using cosmological galaxy formation simulations predicts extreme variations in these cloud properties (Guszejnov, Hopkins & Graus 2019). If the IMF truly depends on the global cloud properties as predicted by classical theories, the observed variations in the IMF would be much larger than what current observational evidence of IMF variations suggests (van Dokkum & Conroy 2012). Note, however, that small-scale predictions of current galaxy formation models must be taken with a grain of salt.
Unlike the classical theories, tidal screening theory does not depend on global cloud properties but only depends on the small-scale physics close to the opacity limit. The characteristic star mass is thus determined based on local gas conditions and not on global cloud conditions; in particular, the properties of the first Larson cores and of their gaseous envelopes are what matters. In this case, a universal IMF is a much more naturally outcome.
Many groups have also studied the IMF directly in numerical simulations. Early results from Bonnell & Bate (2006)a n dN o r dlund & Padoan (2003) confirmed both competitive accretion theory and turbulent fragmentation theory (the turbulent Bonnor-Ebert hydrostatic view). Note that the initial conditions in their simulated clouds were chosen to satisfy the Larson relation, and under these strict Milky Way-like conditions, the emerging IMFs were found to reproduce both the predicted IMF and the observed IMF. These early simulations were not as resolved as the more recent ones. Our resolution study in Appendix B, for example, strongly suggests that our IMF is fully converged. The work of Bate (2009) does have the necessary resolution, but it produces too many brown dwarfs. Otherwise, very similar to the set-up described in this paper, the polytropic index for the EOS used in Bate's work is 7/5, corresponding to a smaller M 1LC and thus an IMF peak at lower masses, in agreement with the tidal screening theory.
More recently, Lee & Hennebelle (2018a) performed simulations of clouds with very different properties, but with the same polytrope and first Larson core mass. They found that the resulting IMF characteristic mass was not evolving with the parent cloud properties, in contradiction with classical theories but in agreement with the tidal screening view.
Using a complementary approach, suites of simulations have been performed using different polytropic EOSs, corresponding to different first Larson core masses, but with identical initial turbulent cloud properties (Bertelli Motta et al. 2016;Lee & Hennebelle 2018b). According to classic theories, this should have led to an identical IMF characteristic mass, since the global cloud properties were unchanged. Surprisingly, they obtained instead a characteristic  These series of recent numerical results, including the one presented in this paper, show us that the characteristic mass of the IMF is still an open question. The tidal screening theory appears to alleviate a number of problems in the classic theories. It is, however, a relatively young theory and probably needs further numerical and analytical studies to fully understand what determines the peak mass of the IMF.

Open questions in tidal screening theory
Tidal screening theory appears as a nice alternative to turbulent fragmentation and competitive accretion. We believe that tidal stripping from the first Larson core and its surrounding envelope could be combined with these older theories to eventually deliver a fully consistent theory for the IMF.
In this paper, we have analysed in great details the tidal field around a newly formed first Larson core embedded in its gaseous, quasi-isothermal envelope. The magnitude of the tidal eigenvalue was used to estimate the region within which clumps of a given density are not able to collapse anymore. Using the lognormal distribution of the density fluctuations inside our turbulent box, we could estimate the probability of finding a collapsing clump of a certain density and at a certain distance from the central core. This is precisely what Lee & Hennebelle (2018b) and Hennebelle, Lee & Chabrier (2019) did in their earlier, albeit recent study using an analytical approach, providing also an analytical estimate of the tidal mass. In this paper, instead of adopting a probabilistic approach, we used the actual gas density field in the simulation to compute the tidal mass. We are, however, still missing a very important step, namely a complete statistical description of the IMF that accounts for tidal screening effects. We speculate that combining the Press and Schechter approach used in the theory of turbulent fragmentation of Hennebelle & Chabrier (2008) and Hopkins (2012) with tidal screening theory could be a promising path. The resulting collapse probability would depend on the distance of a protoclump to a nearby existing core through the two-point correlation function of the turbulent density fluctuations. We postpone the exploration of these ideas for future works.
A second issue that we have also addressed here still remains an open question: what is the real impact of the presence of two or more neighbouring first Larson cores on their respective tidal masses? It is indeed unclear how much mass will be assigned to each core in the final state. In Section 5.4.4, we have addressed this issue by comparing the relative strength of the tidal fields and defining unambiguously one tidal bubble per Larson core. This approach turns out to be a rather good predictor of the final mass for small-mass stars. We still believe, however, that this recipe is too simplistic. A better approach to deal with such cases is the theory of competitive accretion (Bonnell et al. 2001b;Bonnell & Figure 15. Predicted and measured IMFs of the sinks in Fig. 14.   Figure 16. Combining all the data from Fig. 15 for a comparison of the stacked measured IMF at 20 per cent SFE and predicted IMF. Bate 2006). This theory predicts that the accretion rate for stars in the centre of a star cluster is highest, and thus they become most massive. This is not unlike what we see in our tidal bubble scenario. Stars that form in already existing star clusters have small tidal masses because the tidal field of the older sinks and their envelope is quite strong. These stars indeed end up as small-mass stars in the majority of cases. In the tidal screening scenario, the cause of a small stellar mass is the presence of neighbours within a radius of a few hundreds of au (corresponding to the typical tidal radius for a tidal mass equal to the IMF peak mass). It is not entirely clear how this compares with observations, and the precise connection between tidal screening and competitive accretion remains to be explored in future studies. One option could be to explicitly follow the gas using tracer particles in our simulations.
Tidal screening theory also appears to have some difficulty explaining high-mass stars. Even with an age-dependent threshold for the collapse criterion, the final predicted IMF still slightly underestimates the number of stars with masses larger than the peak mass. As proposed by Padoan et al. (2019), massive stars might accrete mass from large-scale filaments that extend far beyond their tidal bubble.

Caveats in our numerical methodology
Modern simulations of turbulent, self-gravitating molecular clouds have come a long way since the early work of Bate (1998). Although resolution has increased significantly, thanks to better numerical techniques and more powerful computers, these simulations will never be a faithful representation of the reality. For example, as mentioned in Section 4, our computational boxes are periodic, without external driving to preserve the kinetic energy of turbulence on large scales. Since we aim at studying the effect of the local tidal field around first Larson cores, we believe our conclusions are not affected by this simplification.
A major caveat of this study is the absence of magnetic fields or radiation fields. While the effect of the latter is captured roughly by our polytropic EOS, the absence of magnetic field could be problematic in the vicinity of the first Larson core and could affect the dynamics of the surrounding envelope. Magnetic fields in the ideal MHD limit are known to prevent the formation of protostellar discs altogether (Hennebelle & Teyssier 2008), which may have an effect on the stellar masses. Fortunately, since what really matters is the geometry of the field lines during the collapse, rather than the strength of the field close to the first Larson core,  have shown that the effect of magnetism on the IMF is minor. Moreover, Masson et al. (2016) have shown that non-ideal effects play a very important role in the vicinity of the first Larson cores, and restore some of the properties (but not all) of the zero magnetic field limit. We therefore believe that it is probably more realistic to ignore magnetic fields in this study, than include them in the ideal MHD limit, and take the risk of grossly overestimating their effect on our results. A fully consistent, non-ideal MHD approach would ideally be needed, but it clearly lies beyond the scope of this study.
Feedback effects from the protostars are also completely ignored in our simulations. Both radiation from the massive end of the stellar population (Dale & Bonnell 2011;Gavagninetal.2017)and jets launched from magnetically supported protostellar discs (Frank et al. 2014;Bally 2016) could have a strong effect. In this paper, we use the simple approach of terminating the simulation early, when only 20 per cent of the mass has been converted into star, as a proxy for the regulating effect of these feedback processes. We also believe that the characteristic mass of the IMF is set up early and by only dynamical processes. This would ultimately need to be tested by implementing these feedback processes in similar simulations. Note, however, that, during all our simulation, the peak of the IMF remains roughly constant. This suggests that the characteristic mass of the IMF does not depend on the exact termination time of star formation and therefore our conclusions would not be affected by adding feedback processes.
While the missing physical effects described here might nevertheless alter the predicted IMF and bring it even closer to the observed one, we want to stress that tidal screening theory was used in this paper to explain the IMF that emerged from our simulations. The fact that it seems to match the observed IMF peak came as a bonus; any additional physics will thus not change our conclusions that tidal screening seems to be the main factor in determining the peak of the IMF in our current simulations. Note, however, that different physical effects play a role at different scales or evolutionary stages in the star formation process. Missing physics might be, for example, the reason why the high-mass tail in our IMF does not match the observed IMF. It has been pointed out before that high mass star formation does not necessarily follow the same mechanisms as low mass star formation (Hennebelle & Commercon 2012;Tan et al. 2014).
Finally, we discuss here possible caveats in our analysis of the simulation outputs and how it seems to support tidal screening theory. First, as presented in the previous sections, the tidal mass can predict the mass of individual stars quite well. The largest caveat here is that there is no guarantee that the mass inside the tidal bubble will end up fully in the sink particle. Star formation in a turbulent cloud is a complex dynamical process and the local environment of the first Larson core might change drastically from the moment of birth to the moment the star reaches its final mass. Individual cores within their gaseous envelopes are also not isolated and might interact with one another, which can strongly affect the local tidal field. Another caveat to discuss here concerns the tight correlation we have observed between the tidal mass and the final star mass. A correlation does not mean causation. While the good values for the Pearson correlation coefficient are surely encouraging to support tidal screening theory, there might be other confounding variables in the problem. More studies are therefore required to check for sure whether the tidal field is responsible for the characteristic mass of stars.

CONCLUSIONS
In conclusion, we summarize our findings as follows: (i) Observations show a universal IMF, something which is difficult to explain using classical IMF theories (turbulent Bonnor-Ebert spheres, turbulent fragmentation, and competitive accretion), which predict an IMF peak that scales with global cloud properties. Recent, high-resolution numerical experiments show that a polytropic EOS naturally introduces a characteristic stellar mass that is independent of the cloud properties and proportional to the first Larson core mass.
(ii) The recently proposed tidal screening theory states that strong tidal fields around first Larson cores prevent clumps from collapsing within a certain radius called the tidal radius. As a result, the mass within this radius will be accreted by the central first Larson core. This theory predicts an IMF peak around 10 × M 1LC and independent of the global cloud properties.
(iii) We explored this theory further and derived a collapse criterion for clumps that form close to an existing star (or first Larson core). Whether a clump can collapse or not depends on its density, the distance to the star, the mass of the star, and the profile of the envelope around the star.
(iv) We performed a series of numerical simulations of a small star-forming turbulent box representing a collapsing region inside a larger molecular cloud. Our simulations only include hydrodynamics, self-gravity, and a polytropic EOS. This way, we isolate the tidal effect. We measure the IMF and find it has a peak at about 0.1 M ⊙ throughout the simulations. Classical theories do not explain this result.
(v) We studied the tidal fields around first Larson cores (sink particles) and find they are strong enough to prevent potential clumps from collapsing unless they have high enough densities. Since these clumps are rare due to the form of the density PDF generated by turbulence, tidal forces effectively screen a large region around the first Larson core.
(vi) Using the actual density field of the simulation and a local collapse criterion, we construct a tidal bubble around every star at the moment of its birth, in which tidal forces are strong enough to strip anything within a time-scale shorter than the remaining simulation time. We then use the gas mass in this tidal bubble as a predictor for the final mass of the star.
(vii) We found a very good correlation between our predicted mass and the true final mass of the star in our simulations. Stars that formed initially in isolation generally have masses at or above the characteristic mass. Low-mass stars are formed in regions that already contain other stars within a radius of a few hundreds of au. This is something one would expect from the tidal screening theory.
(viii) Our predicted IMF, based on the tidal mass, matches the simulated one, based on the sink particle final mass, quite well, with a peak slightly shifted to lower masses. The predicted IMF also reproduces the simulated IMF at both the high-and low-mass ends. Note that both IMFs, predicted and simulated, also match quite well the observed IMF.
We conclude that tidal forces are an important ingredient in determining the final mass of a star. We speculate that a complete theory of the IMF would combine tidal screening theory with competitive accretion and turbulent fragmentation, possibly leading to more accurate predictions. It remains to confirm that the tidal screening formalism presented in this paper can still predict the IMF when additional physics is included. This will be the topics of a follow-up paper. Figure A1. Evolution of the volume-weighted density PDF in the initialization phase of box 1. The black curve is the PDF at the moment where the kinetic energy is half.  Figure A2. The volume-weighted density PDF at the end of the initialization for each box.
five boxes have different amplitudes for the velocity on large scales, they evolve quite differently. Table A2 lists some properties of the initial conditions at the end of our initialization phase. Note that the 1D Mach number is not necessarily equal in all directions. Fig. A2 compares the density PDF of these snapshots and one can see they are fairly similar.

APPENDIX B: RESOLUTION EFFECTS
We tested the convergence of our results by running the same simulation with three different maximum refinement levels: 10 (low), 12 (fiducial), and 14 (high), corresponding to a spatial resolution of 50.4, 12.6, and 3.15 au, respectively.

B1 Sink formation and resolving individual stars
The sink formation recipe has to be adapted for each resolution. We choose the sink formation threshold to correspond to the density at which a next refinement would have been triggered if it would have been allowed, The seed mass is chosen to be the isothermal Jeans mass at 10 K and the sink formation density, The exact values can be found in Table B1.
It is important to mention that the polytropic EOS is kept the same for all resolutions. The properties of the first Larson core depend on the EOS but not on the resolution. For the same reason, the merging time-scale for young sinks is kept constant for all resolutions. A list of the parameters that are kept fixed is given in Table B2. Fig. B1 shows a time sequence of the same star-forming region modelled at different resolutions. The top row corresponds to the first time in the sequence and shows the formation of the same sink in all cases. In the second row, the images correspond to a later time when a companion star is formed quite close to the first protostar. This second sink never formed in the low resolution case. This is a typical example that shows that the low-resolution simulation is not capable of resolving all binaries.   Figure B1. Example of an evolution sequence of a star-forming region. The maps show the maximum density along the line of sight. The green circles show the sink accretion radius with the colour indicating the sink's mass: the darker the circle, the more massive the sink. From left to right, the maximum refinement level is increased: levelmax 10, 12, and 14. This particular case has an unresolved binary in the low resolution case and merges with another region around 24 kyr. Top row: the primary sink is born. Second row: the companion sink's birth, which is not resolved in the low resolution case. Third row: the formationof another companion, this time present in all resolutions. Bottom row: the region merged with another star-forming region. Figure B2. The mass evolution of each sink in the star formation region displayed in Fig. B1. Each track represents a sink particle, while the coloured background shows the total mass in the region. The vertical lines mark the times indicated in Fig. B1.
In Fig. B2, we have represented a diagram that shows the mass evolution of the various sinks in this particular star-forming region. After 16 kyr, we see that the fiducial and high resolution runs do contain a second sink with a similar mass, while the low resolution has only one sink and hence failed to form a binary. Note, however, that the total mass in sinks at this time is very robust and is the same for all resolutions. Around time 22 kyr, a third member in the region appears for all resolutions. After time 25 kyr, another star-forming region merges with this one and the dynamical state of the resulting start cluster becomes a lot more complex. Note again, however, that the total mass in the region is still always the same for all resolutions.
In Fig. B3, we show four different star-forming regions at the same time but at different resolutions. They represent three typical configurations: an isolated sink, a small cluster with only two or three members and finally a large cluster with many sinks. As one might expect, the simpler the geometry, the faster the convergence and the more robust the numerical solution. Fig. B4 shows the sink mass evolution in these same four regions. For the isolated sink case, the agreement between the different resolutions is nearly perfect. For the other cases, we see that the convergence of the individual sink mass is slower, but the total mass in the region is very robust with increasing resolution.
In summary, our resolution study supports the claim that having a spatial resolution of 12.6 au is enough to get converged results for the number of sinks and for the value of their individual masses. This resolution corresponds to the size of the first Larson core.

B2 Influence on the IMF
Since we are mostly interested in the IMF and less in the exact mass of individual stars, we now study the effect of our different resolutions on the statistics of sinks. Fig. B5 shows the IMF of all initial conditions at 30 kyr and for different resolutions. In all cases, we see that the high-mass end agrees quite well. The peak and the low-mass end of the IMF for the fiducial and high resolution runs are also quite close, while the low resolution case is clearly not converged. In box 2 and box 5, we can also see a small excess of high-mass stars in the low resolution runs. These are unresolved star clusters with multiple members resolved only at higher resolution. Fig. B6 shows the stacked IMF (data from all initial conditions gathered together). Also, here we see that only the low resolution is incomplete, while the peak in the other two cases remains the same. From these results, and our findings in the previous section, we conclude that the IMFs in our fiducial runs are already numerically converged.

B3 Sink merging time-scale
To check how the results are affected by the chosen sink merging time-scale, we reran box 3 using a value different by a factor of 3. From Fig. B7, we can see that the results are not affected as long as the value chosen is of the order of the free-fall time at the critical density, corresponding to the time needed for the first Larson core to collapse entirely into the protostar. Figure B3. Examples of different star formation regions in our simulations (the layout of the maps is the same as in Fig. B1). From top to bottom, there is an increasing geometric complexity: an isolated sink, a binary, two small regions merging, and a complex star cluster. The more complex the star-formingregions, the more difficult it is to reproduce the results when changing the resolution.   This paper has been typeset from a T E X/L A T E X file prepared by the author.