Converging on the Initial Mass Function of Stars

Understanding the origin of stellar masses—the initial mass function (IMF)— remains one of the most challenging problems in astrophysics. The IMF is a key ingredient for simulations of galaxy formation and evolution, and is used to calibrate star formation relations in extra-galactic observations. Modeling the IMF directly in hydrodynamical simulations has been attempted in several previous studies, but the most important processes that control the IMF remain poorly understood. This is because predicting the IMF from direct hydrodynamical simulations involves complex physics such as turbulence, magnetic fields, radiation feedback and mechanical feedback, all of which are difficult to model and the methods used have limitations in terms of accuracy and computational efficiency. Moreover, a physical interpretation of the simulated IMFs requires a numerically converged solution at high resolution, which has so far not been convincingly demonstrated. Here we present a resolution study of star cluster formation aimed at producing a converged IMF. We compare a set of magnetohydrodynamical (MHD) adaptive-mesh-refinement simulations with three different implementations of the thermodynamics of the gas: 1) with an isothermal equation of state (EOS), 2) with a polytropic EOS, and 3) with a simple stellar heating feedback model. We show that in the simulations with an isothermal or polytropic EOS, the number of stars and their mass distributions depend on the numerical resolution. By contrast, the simulations that employ the simple radiative feedback module demonstrate convergence in the number of stars formed and in their IMFs.


Introduction
The IMF is the distribution of stellar masses when they are born. We know from observational surveys that most stars have masses of about half the mass of our Sun (M ). Stars with smaller masses are rarer and there is a natural cutoff at at the brown dwarf mass of 0.08 M , because gaseous, self-gravitating objects with masses M 0.08 M cannot fuse hydrogen to helium. Stars more massive than the Sun also become rarer with increasing mass. The highmass tail of the IMF is indeed a steeply decreasing power-law function with the number of stars N (M ) ∝ M −1. 35 [1,2,3,4,5,6,7,8]. Understanding this seemingly universal power-law tail and the turnover at around 0.5 M is one of the most challenging open problems in astrophysics. The IMF has far-reaching consequences and applications, e.g., for calibrating extra-galactic star formation relations used to understand galaxy formation and evolution.
A number of numerical simulations exist in the literature that aim to model and understand the IMF from direct hydrodynamical simulations [9,10,11,12,13,14,15,16]. A critical problem in modeling the IMF is that one must resolve the fragmentation of magnetized, turbulent gas, including radiation feedback [17,18,19,20,21,22,23] and mechanical feedback by jets and outflows [24,25,26,27]. Modeling all of the required physical mechanisms is difficult, and achieving sufficient resolution at the same time is even harder. Furthermore, forming a sufficiently large number of stars in a simulation requires one to either run many smallscale clouds with different realizations of the turbulence or one large-scale, high-mass cloud that spawns many stars to obtain a statistically representative and converged mass function. Currently, there is no set of simulations that fulfills all of these requirements and produces a converged IMF.
Here we present a simple model to include stellar heating feedback in MHD simulations of star cluster formation. We compare this stellar heating model to simulations that use an isothermal or polytropic EOS. Running low-resolution, mid-resolution and high-resolution simulations for all of these models, we show that only the simulations that include the stellar heating feedback converge on the number of stars formed and their cumulative mass distributions.

Equations of magnetohydrodynamics
We use a modified version of the adaptive mesh refinement (AMR) [28] code FLASH [29,30] (v4) to solve the three-dimensional, ideal MHD equations including gravity, where the gravitational acceleration of the gas g, is the sum of the self-gravity of the gas and the contribution of sink particles (see Section 2.3 below): In these equations, ρ, v, P tot = P + 1/(8π) |B| 2 , B, and E = ρ int + (ρ/2) |v| 2 + 1/(8π) |B| 2 denote the gas density, velocity, pressure (thermal+magnetic), magnetic field, and total energy density (internal+kinetic+magnetic), respectively. We use the positive-definite HLL3R approximate Riemann scheme [31] to solve the MHD equations (1). The HLL3R solver guarantees stability and accuracy of the numerical solution of the hydrodynamical equations [32,33,34,31].

Turbulence driving
All our simulations include the turbulence driving module in the FLASH code [35]. Turbulence in the interstellar medium is likely driven by a combination of supernova explosions, stellar feedback in the form of jets, outflows and winds, gravitational contraction and accretion, as well as magneto-rotational instability and galactic spiral-arm/disk dynamics [36,37,38].
Turbulence is indeed a critical ingredient for star formation [39,36,40,41,42,43,44,45,46,47,48,49]. Since turbulence decays quickly [50,51], but is observed on virtually all spatial scales in the ISM, the turbulence needs to be driven. Here we use a turbulence driving module that produces turbulence similar to the observed turbulence in real molecular clouds, i.e., driving on the largest scales [52,53] and with a velocity power spectrum ∼ k −2 , consistent with supersonic, compressible turbulence [54,55,56]. This type of velocity spectrum is further consistent with simulations of supersonic turbulence [57,58,35,59,60]. Our turbulence driving module [35] applies a stochastic Ornstein-Uhlenbeck process [61,62] to construct an acceleration field F stir , which serves as a momentum and energy source term in the MHD equations. Following observations, F stir only drives large-scale modes, 1 < |k| L/2π < 3, where most of the power is injected at the k inj = 2 mode in Fourier space, i.e., on half of the size of the computational domain. The turbulence on smaller scales is not directly affected by the driving and develops self-consistently from the cascade of energy starting on the large scales. In this study we use mixed driving of the turbulence with a turbulence driving parameter b = 0.5, typical for clouds in the Milky Way disk [63]. However, cloud-to-cloud variations in the turbulence driving parameter b, which is related to the mixture of driving modes [64,35,65,66,67,68], suggest variations from b ∼ 1/3 (purely solenoidal driving, ∇ · F stir = 0) to b ∼ 1 (purely compressive driving, ∇×F stir = 0) in observed clouds in the Milky Way [69,70,71,72,73,63]. Investigating potential variations of the IMF between purely solenoidal driving and purely compressive driving will be the subject of a follow-up study.

Sink particles and AMR
In order to model star formation and accretion of gas, we use the sink particle method in the FLASH code [74]. Sink particles form dynamically in our simulations when a local region in the simulation undergoes gravitational collapse and forms stars. This is technically achieved by first flagging each computational cell that exceeds the Jeans resolution density, where c s is the sound speed, G is the gravitational constant, and λ J = [πc 2 s /(Gρ)] 1/2 is the local Jeans length. This criterion defines a sink particle radius, which we set to 2.5 grid cell lengths at the highest level of AMR, to avoid artificial fragmentation [75]. If the gas density in a cell exceeds ρ sink , a spherical control volume with radius r sink is gathered around the cell and it is checked that all the gas within this control volume is Jeansunstable, gravitationally bound and converging towards the central cell of the control volume. A sink particle is formed there, only if all of these checks are passed. This avoids spurious creation of sink particles and guarantees that only bound and collapsing gas forms stars in our simulations [74]. On all AMR levels lower than the maximum AMR level (where sink particles form), we use an AMR criterion based on the local Jeans length, such that λ J is always resolved with at least 16 grid cell lengths. This resolution criterion is conservative and computationally costly, but allows us to resolve turbulence on the Jeans scale [76], potential dynamo amplification of the magnetic field in dense cores [77], and allows us to capture the basic structure of accretion discs forming on the smallest scales resolved in the simulations [27]. Accretion of gas occurs if a cell within the accretion radius of an existing sink particle exceeds ρ sink , is gravitationally bound to the sink particle and converging towards it. If these accretion criteria are met, the excess mass above ρ sink is accreted onto the sink particle, with conservation of mass, center-of-mass, momentum, and angular momentum enforced by the algorithm. Self-gravity of the gas is computed with a multigrid solver [78], while all gravitational interactions of the sink particles with the gas and with one another are computed by direct summation over all sink particles and grid cells. This procedure avoids inaccuracies in the particle accelerations if particle-mesh and mesh-particle mapping methods were used instead. A second-order leapfrog integrator is used to advance the sink particles with a timestep that allows us to resolve highly eccentric orbits of stars [74].

Equation of state (EOS)
To model the thermal evolution of the gas, it is common practice in star formation simulations to use an isothermal or polytropic equation of state for the gas pressure, 1 which, using the ideal gas EOS, defines the temperature due to the EOS, where the polytropic exponent Γ may be set to Equations (5-7) approximate previous detailed radiation-hydrodynamic simulations of dense star-forming gas [82]. The polytropic constant K in Equation (5) is adjusted such that K = c 2 s , i.e, the square of the sound speed. When the polytropic exponent Γ changes according to the density regimes given by Equation (7), the temperature and sound speed change, and K is adjusted, such that the pressure and temperature are continuous functions of density. In the isothermal regime (Γ = 1), which serves as the normalization, the sound speed c s = 0.2 km s −1 and the temperature T = 11 K for gas with a typical molecular weight of 2.3 m H (with m H being the mass of a hydrogen atom). The initial isothermal evolution for ρ ≤ 2.5 × 10 −16 g cm −3 is a reasonable approximation for dense, molecular gas of solar metallicity, over a wide range of densities [83,84,85,86,87,80,88,89,90]. We note that for the numerical resolutions and respective sink particle density thresholds used in the simulations presented below, only the first two density regimes in Equation (7) are relevant.
Equations (5)(6)(7) are an approximation to hydrodynamic calculations that take radiative transfer into account [91,92,82,17,93,18,94,20,21,14], covering the phases of isothermal contraction, adiabatic heating during the formation of the first and second core and the effects of H 2 dissociation in the second collapse. However, it does not include the radiative heating feedback from the accreting protostar once it has formed. This additional stellar heating feedback, however, is crucial, because it dominates the heating in the direct vicinity of the protostar, out to distances of the order of ∼ 1000-10000 AU from the star. Our main goal here is to account for this additional heating by the protostars.

A simple model for stellar heating
We have implemented a new stellar feedback module into the FLASH MHD code, which accounts for the effect of gas heating around newborn protostars until they reach the main sequence. Stellar heating is a direct feedback mechanism (with a spatial correlation based on the position of each star) and must be distinguished from the EOS (which relates pressure, density and Since there is no spatial information in the EOS, feedback cannot be modeled with a standard EOS. Here we provide a simple algorithm to add stellar heating to the simulations, which allows us to approximate the direct radiative heading feedback from stars.

Protostellar and stellar luminosities
The first step in accounting for the direct heating feedback from stars and protostars is to calculate the star's luminosity. Here we use the stellar evolution model developed by Offner et al. [93], which is based on one-zone protostellar evolution models [95,96,97] calibrated to match detailed numerical calculations of protostellar evolution [98,99]. This module provides stellar luminosities from the birth of a star, while it is still accreting, to the main sequence, depending on the evolutionary stage. The following evolutionary phases are distinguished: 'pre-collapse', 'no burning', 'Deuterium burning at fixed core temperature', 'Deuterium burning at variable core temperature', 'Deuterium burning in the shell', and 'main sequence'. The 'pre-collapse' state applies to a newly formed sink particle with M sink < 0.01 M and zero luminosity. During the 'no burning' stage, the star accretes, but the luminosity is determined solely based on the energy released by gravitational contraction of the protostar. As the radius and Deuterium mass of the star evolve further, Deuterium burning first starts in the core and eventually in the shell. Finally, the star reaches the zero-age main sequence, and after that stellar evolution continues on the main sequence.

Calculating the stellar heating temperature
Based on the stellar luminosity (L ) from the stellar evolution module discussed in the previous subsection, which includes accretion luminosity, we can now estimate the heating of the gas around each sink particle. An accurate calculation of the temperature around each stellar radiation source would require solving the radiative transfer equation [93]. Here we use a simple model for the gas heating that assumes spherical symmetry and a homogenous gas distribution around the star, which allows us to approximate the heating of the gas without solving the radiative transfer equation. We estimate the radial dependence of the temperature due to stellar heating (T heat ) by where σ SB is the Stefan-Boltzmann constant and r is the distance from the star. This allows us to calculate T heat (r) based on L for each star in the computational domain. For practical reasons and to speed up the computation, we only calculate T heat for r ≤ 10 4 AU, because the contribution of stellar heating at greater distances from the star is essentially negligible [93]. We further apply an inner heating radius of 100 AU, below which T heat = T heat (100 AU) = const, to approximate the radial dependence of T heat obtained in full radiation hydrodynamical simulations [93]. If the heating radii of multiple stars overlap, we add T 4 heat (r) from each star locally in each cell where the overlap occurs.

Coupling to the gas EOS
Finally, we add the stellar heating T heat (r) from Equation (8) to the existing gas temperature (from the EOS, Eq. 6) at every grid point r and for every stellar source. This is similar to the procedure in recent semi-analytic models of the IMF [100] and gives the final gas temperature,  (FB). Column 3: maximum effective grid resolution (refinement based on Jeans length, guaranteeing a minimum of 16 cells and a maximum of 32 cells per Jeans length). The base grid is always resolved with at least 256 3 cells. Columns 4-6: minimum cell size, sink particle radius, and sink density threshold (see Eq. 3 and 4). Columns 7-9: number of sink particles formed, mean, and median sink particle mass when SFE = 10%, i.e., when a fraction of 10% of the original cloud mass has formed stars. or in terms of the total gas pressure (directly relevant for integrating the hydrodynamical equations; see Eq. 1), We note that adding the fourth power of the EOS and heating temperatures is a simplification that may lead to overestimating the heating in regions where T EOS ∼ T heat . This will be addressed in a forthcoming paper. Here we use the simple summation given by Equation (9), suggested in previous work [100].

Simulation parameters and initial conditions
All our simulations were carried out in a three-dimensional triple-periodic computational domain with side length L = 2 pc. The initial gas density is uniform with ρ = 6.56 × 10 −21 g cm −3 , giving a total cloud mass of M = 775 M and a mean freefall time of t ff = 0.82 Myr. The initial velocities are zero and the turbulence driving (see §2.2) excites turbulent motions. A fullydeveloped turbulent state is reached after 2 crossing times, L/(Mc s ) [101,35]. The amplitude of the driving was adjusted to drive turbulence with a steady-state sonic Mach number of M = σ v /c s = 5.0, where the velocity dispersion σ v = 1.0 km s −1 and the initially isothermal sound speed c s = 0.2 km s −1 , appropriate for the initial density in the molecular phase of the ISM. The mass and velocity dispersion define the initial virial parameter, the ratio of twice the kinetic energy to the gravitational energy, α vir = 2 E kin /E grav = 0.5, which is in the range of observed virial parameters [102,103,104]. We add an initially uniform magnetic field along the zaxis of the computational domain with B = 10 −5 G, which is subsequently compressed, tangled and twisted by the turbulence, approximating the magnetic field structure in real molecular clouds [105]. This results in an initial plasma beta (ratio of thermal to magnetic pressure)  [102,106,107,108].
We run a total of 9 simulations (see Table 1): 3 different numerical resolutions (low-resolution, mid-resolution, and high-resolution) times 3 different methods of treating the thermodynamics (isothermal EOS, polytropic EOS, and polytropic EOS + stellar heating feedback). A detailed list of the different simulation parameters (EOS and numerical resolution) and main outcomes (number of sink particles formed and their mean and median mass) is provided in Table 1.

Results
We now compare our 3 physical models with different EOS (isothermal, polytropic, and stellar heating) for different numerical resolutions. The main aim is to determine which model may achieve convergence with increasing numerical resolution. Figure 1 shows the column density along the z-axis (the gas density integrated along the lineof-sight) of simulations 1, 2, and 3 (top panels) and simulations 7, 8, and 9 (bottom panels) from Table 1 turbulence is fully developed, and self-gravity and sink particle formation had been switched on. At this time, all simulations have reached a star formation efficiency (SFE) of about 10% (i.e., a fraction of 10% of the original cloud mass has formed stars). However, the number of sink particles formed varies significantly, depending on the EOS and the numerical resolution. The isothermal simulations form 37 and 84 sink particles in the low-and high-resolution runs, respectively, and the polytopic EOS runs from 33 and 62, respectively. This means that these simulations clearly lack convergence in the number of stars formed. Only the simulations that include stellar heating feedback show convergence, with 28 and 29 sink particles formed in the low-and high-resolution simulations, respectively.

Time evolution of the number sink particles and cumulative mass distribution functions
We now investigate the numerical convergence in the number of sink particles as a function of time. The left-hand panels of Figure 2 show the number of sink particles in all our simulations as a function of SFE. We see excellent convergence in our stellar heating runs. Finally, we show the cumulative mass functions in all our simulations in the right-hand panels of Figure 2. We see that in the isothermal and polytropic simulations, the mass functions depend on the numerical resolution. By contrast, the stellar mass functions are well converged for the simulations with stellar heating feedback.

Summary and conclusions
We implemented a new simple stellar heating module in the FLASH MHD code ( §3). This heating feedback module takes the luminosity of a protostar or evolved star and distributes the radiated energy in a spherically symmetric approximation around each star. This allows us to approximate the effects of full radiation transfer of the stellar radiation field in a computationally fast and simple implementation.
In a series of simulations that employ an isothermal EOS (Eq. 5 with Γ = 1), a polytropic EOS (with Γ given by Eq. 7), and our new stellar heating module, we investigate the spatial distribution of formed stars, the time evolution of the number of stars, and the cumulative stellar mass functions. We compare low-resolution, mid-resolution, and high-resolution simulations and demonstrate that only the simulations with stellar heating feedback achieve numerical convergence, while simulations with a purely isothermal or polytropic EOS fail to converge on the number of stars and stellar masses formed. We conclude that the direct heating feedback from stars is essential to understanding the IMF of stars and needs to be included in simulations to achieve numerical convergence.
Finally, we caution that given our current resolutions of ∼ 100 AU, these simulations may not be able to accurately capture physical fragmentation on the disk scale. Moreover, the radiation feedback model presented here assumes spherical symmetry, which breaks down on disk scales. Future improvements of the model should take the asymmetry of the accretion disks and the resulting shielding of radiation into account.