Home Search Collections Journals About Contact us My IOPscience Converging on the Initial Mass Function of Stars This content has been downloaded from IOPscience. Please scroll down to see the full text. 2017 J. Phys.: Conf. Ser. 837 012007 (http://iopscience.iop.org/1742-6596/837/1/012007) View the table of contents for this issue, or go to the journal homepage for more Download details: IP Address: 131.215.225.187 This content was downloaded on 31/05/2017 at 17:41 Please note that terms and conditions apply. You may also be interested in: STARS AS THEY LOOK AND AS THEY ARE P. W. Merrill ROTATIONAL VELOCITIES OF SOME B STARS Hugo Levato and Stella Malaroda CORPUSCULAR RADIATION FROM STARS Su-Shu Huang ANALYTICAL THEORY FOR THE INITIAL MASS FUNCTION. II. Patrick Hennebelle and Gilles Chabrier TECHNIQUES FOR SURFACE IMAGING OF STARS N. E. Piskunov and J. B. Rice NOTE ON THE RADIATION FROM STARS W. W. Coblentz Thermal Conduction in Galaxy Clusters K. Dolag, M. Jubelgas, V. Springel et al. Molecular gas in active environments S Mühle, C Henkel, T de Maio et al. Molecular Hydrogen and Star Formation Brant E. Robertson and Andrey V. Kravtsov ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 Converging on the Initial Mass Function of Stars Christoph Federrath1, Mark Krumholz1, Philip F. Hopkins2 1Research School of Astronomy and Astrophysics, The Australian National University, Canberra, ACT 2611, Australia 2TAPIR, Mailcode 350-17, California Institute of Technology, Pasadena, CA 91125, USA E-mail: [email protected] Abstract. 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 simulationsofgalaxyformationandevolution,andisusedtocalibratestarformationrelations 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 IMFremainpoorlyunderstood. ThisisbecausepredictingtheIMFfromdirecthydrodynamical simulationsinvolvescomplexphysicssuchasturbulence,magneticfields,radiationfeedbackand 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 thermodynamicsofthegas: 1)withanisothermalequationofstate(EOS),2)withapolytropic EOS, and 3) with a simple stellar heating feedback model. We show that in the simulations withanisothermalorpolytropicEOS,thenumberofstarsandtheirmassdistributionsdepend 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. 1. 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 (cid:12) smaller masses are rarer and there is a natural cutoff at at the brown dwarf mass of 0.08M , (cid:12) because gaseous, self-gravitating objects with masses M (cid:46) 0.08M cannot fuse hydrogen to (cid:12) helium. Stars more massive than the Sun also become rarer with increasing mass. The high- mass 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 andtheturnoverataround0.5M isoneofthemostchallengingopenproblemsinastrophysics. (cid:12) 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 probleminmodelingtheIMFisthatonemustresolvethefragmentationofmagnetized,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, ContentfromthisworkmaybeusedunderthetermsoftheCreativeCommonsAttribution3.0licence.Anyfurtherdistribution ofthisworkmustmaintainattributiontotheauthor(s)andthetitleofthework,journalcitationandDOI. PublishedunderlicencebyIOPPublishingLtd 1 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 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 small- scale 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. HerewepresentasimplemodeltoincludestellarheatingfeedbackinMHDsimulationsofstar 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. 2. Numerical methods 2.1. 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, ∂ρ +∇·(ρv) = 0, ∂t (cid:18) (cid:19) ∂ (B·∇)B ρ +v·∇ v = −∇P +ρg, tot ∂t 4π (cid:20) (cid:21) ∂E (B·v)B +∇· (E +P )v− = ρv·g, tot ∂t 4π ∂B = ∇×(v×B) , ∇·B = 0, (1) ∂t 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): g = −∇Φ +g , gas sinks ∇2Φ = 4πGρ. (2) gas In these equations, ρ, v, P = P +1/(8π)|B|2, B, and E = ρ(cid:15) +(ρ/2)|v|2 +1/(8π)|B|2 tot int 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]. 2.2. 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 2 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 consistent with simulations of supersonic turbulence [57, 58, 35, 59, 60]. Our turbulence driving module[35]appliesastochasticOrnstein-Uhlenbeckprocess[61,62]toconstructanacceleration fieldF ,whichservesasamomentumandenergysourcetermintheMHDequations. Following stir observations, F only drives large-scale modes, 1 < |k|L/2π < 3, where most of the power stir is injected at the k = 2 mode in Fourier space, i.e., on half of the size of the computational inj 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 mixeddrivingoftheturbulencewithaturbulencedrivingparameterb = 0.5,typicalforcloudsin theMilkyWaydisk[63]. However, cloud-to-cloudvariationsintheturbulencedrivingparameter 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 = 0) to b ∼ 1 (purely compressive driving, stir ∇×F = 0)inobservedcloudsintheMilkyWay[69,70,71,72,73,63]. Investigatingpotential stir variations of the IMF between purely solenoidal driving and purely compressive driving will be the subject of a follow-up study. 2.3. 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, πc2 πc2 ρ = s = s , (3) sink Gλ2 4Gr2 J sink where c is the sound speed, G is the gravitational constant, and λ = [πc2/(Gρ)]1/2 is the local s J s Jeans length. This criterion defines a sink particle radius, r = λ /2, (4) sink J whichwesetto2.5gridcelllengthsatthehighestlevelofAMR,toavoidartificialfragmentation [75]. If the gas density in a cell exceeds ρ , a spherical control volume with radius r is sink sink gathered around the cell and it is checked that all the gas within this control volume is Jeans- unstable, 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]. OnallAMRlevelslowerthanthemaximumAMRlevel(wheresinkparticlesform),weusean AMR criterion based on the local Jeans length, such that λ is always resolved with at least 16 J 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 scalesresolved inthe simulations [27]. Accretion of gasoccurs if a cell within the accretion radius of an existing sink particle exceeds ρ , is gravitationally bound to the sink sink particle and converging towards it. If these accretion criteria are met, the excess mass above ρ is accreted onto the sink particle, with conservation of mass, center-of-mass, momentum, sink 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]. 3 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 2.4. 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 P = KρΓ, (5) EOS which, using the ideal gas EOS, defines the temperature due to the EOS, µm µm T = HP = HKρΓ−1, (6) EOS EOS k ρ k B B where the polytropic exponent Γ may be set to   1 for ρ ≤ ρ ≡ 2.50×10−16gcm−3,  1    1.1 for ρ < ρ ≤ ρ ≡ 3.84×10−13gcm−3,  1 2  Γ = 1.4 for ρ < ρ ≤ ρ ≡ 3.84×10−8 gcm−3, (7) 2 3    1.1 for ρ < ρ ≤ ρ ≡ 3.84×10−3 gcm−3,  3 4     5/3 for ρ > ρ . 4 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 = c2, 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 = 0.2kms−1 s and the temperature T = 11K for gas with a typical molecular weight of 2.3m (with m being H H the mass of a hydrogen atom). The initial isothermal evolution for ρ ≤ 2.5×10−16gcm−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–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 dissociation in the second collapse. However, it does not include the radiative heating 2 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–10000AU from the star. Our main goal here is to account for this additional heating by the protostars. 3. A simple model for stellar heating WehaveimplementedanewstellarfeedbackmoduleintotheFLASHMHDcode,whichaccounts 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 1 Another common alternative is to use an ideal gas equation of state and to include heating and cooling terms, eitherbyusingacoolingcurve[79],followingthechemicalevolutionincludingheatingandcoolingterms[80,81], or by directly solving the radiative transfer equation. 4 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 temperature of the gas); see §2.4. 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 heatingtothesimulations, whichallowsustoapproximatethedirectradiativeheadingfeedback from stars. 3.1. 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], whichisbased onone-zone protostellarevolutionmodels [95,96,97] calibratedtomatch 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’stateappliestoanewlyformedsinkparticlewithM < 0.01M andzero sink (cid:12) 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. 3.2. Calculating the stellar heating temperature Based on the stellar luminosity (L ) from the stellar evolution module discussed in the previous (cid:63) 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 ) by heat L T4 (r) = (cid:63) , (8) heat 4πσ r2 SB where σ is the Stefan-Boltzmann constant and r is the distance from the star. This allows SB us to calculate T (r) based on L for each star in the computational domain. For practical heat (cid:63) reasons and to speed up the computation, we only calculate T for r ≤ 104AU, because the heat contribution of stellar heating at greater distances from the star is essentially negligible [93]. We further apply an inner heating radius of 100AU, below which T = T (100AU) = const, heat heat to approximate the radial dependence of T obtained in full radiation hydrodynamical heat simulations [93]. If the heating radii of multiple stars overlap, we add T4 (r) from each star heat locally in each cell where the overlap occurs. 3.3. Coupling to the gas EOS Finally, we add the stellar heating T (r) from Equation (8) to the existing gas temperature heat (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, T = (cid:0)T4 +T4 (cid:1)1/4, (9) EOS heat 5 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 Table 1. Key simulation parameters and results. Simulation EOS N3 ∆x r ρ N Mean Median eff sink sink sinks [AU] [AU] [gcm−3] M M sink sink (1) (2) (3) (4) (5) (6) (7) (8) (9) 1. LowResIso Isothermal 20483 200 500 8.3×10−17 35 2.2M 1.1M (cid:12) (cid:12) 2. LowResPoly Polytropic 20483 200 500 8.3×10−17 31 2.5M 2.0M (cid:12) (cid:12) 3. LowResHeat Heating FB 20483 200 500 8.3×10−17 29 2.7M 1.4M (cid:12) (cid:12) 4. MidResIso Isothermal 40963 100 250 3.8×10−16 44 1.7M 1.1M (cid:12) (cid:12) 5. MidResPoly Polytropic 40963 100 250 3.8×10−16 47 1.7M 0.9M (cid:12) (cid:12) 6. MidResHeat Heating FB 40963 100 250 3.8×10−16 29 2.6M 1.4M (cid:12) (cid:12) 7. HighResIso Isothermal 81923 50 125 1.8×10−15 88 0.9M 0.4M (cid:12) (cid:12) 8. HighResPoly Polytropic 81923 50 125 1.8×10−15 62 1.3M 0.4M (cid:12) (cid:12) 9. HighResHeat Heating FB 81923 50 125 1.8×10−15 30 2.6M 1.8M (cid:12) (cid:12) Notes. Column 1: simulation model. Column 2: equation of state: isothermal, polytropic, or polytropic plus stellar heating feedback (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 2563 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), (cid:34) (cid:18) k ρ (cid:19)4 (cid:35)1/4 P = P4 + B T4 . (10) EOS µm heat H 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 ∼ T . This will be EOS heat addressed in a forthcoming paper. Here we use the simple summation given by Equation (9), suggested in previous work [100]. 4. Simulation parameters and initial conditions Alloursimulationswerecarriedoutinathree-dimensionaltriple-periodiccomputationaldomain with side length L = 2pc. The initial gas density is uniform with ρ = 6.56 × 10−21gcm−3, giving a total cloud mass of M = 775M and a mean freefall time of t = 0.82Myr. The initial (cid:12) ff velocities are zero and the turbulence driving (see §2.2) excites turbulent motions. A fully- developed turbulent state is reached after 2 crossing times, L/(Mc ) [101, 35]. The amplitude s of the driving was adjusted to drive turbulence with a steady-state sonic Mach number of M = σ /c = 5.0, where the velocity dispersion σ = 1.0kms−1 and the initially isothermal v s v sound speed c = 0.2kms−1, appropriate for the initial density in the molecular phase of the s ISM. The mass and velocity dispersion define the initial virial parameter, the ratio of twice the kinetic energy to the gravitational energy, α = 2E /E = 0.5, which is in the range of vir kin grav observedvirialparameters[102,103,104]. Weaddaninitiallyuniformmagneticfieldalongthez- axis of the computational domain with B = 10−5G, 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) 6 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 Figure 1. Each image shows a (2pc)2 projection of gas density along the z-axis. Each column shows simulations with an isothermal EOS, Eq. (5) with Γ = 1 (left-hand panels), a polytropic EOS with Γ given by Eq. (7) (middle panels), and the polytropic EOS combined with our new stellar heating feedback introduced in §3 (right-hand panels). The top panels show the low-resolution runs with an effective maximum resolution of 20483 cells and the bottom panels show the high-resolution runs with 81923 effective resolution (see Tab. 1 for details). Only the simulations that include stellar heating feedback show convergence in the number of sink particles formed. of β = 0.66 and an Alfv´en Mach number (ratio of velocity dispersion to Alfv´en speed) of M = 2.9. The magnetic field strengths in the simulations are consistent with observed values A [102, 106, 107, 108]. Werunatotalof9simulations(seeTable1): 3differentnumericalresolutions(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. 5. 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. 5.1. Spatial distribution of gas and stars Figure 1 shows the column density along the z-axis (the gas density integrated along the line- of-sight) of simulations 1, 2, and 3 (top panels) and simulations 7, 8, and 9 (bottom panels) from Table 1. All panels show each simulation at the same time, 0.5Myr (or 0.55t ) after the ff 7 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 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. 5.2. 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, weshowthecumulativemassfunctionsinalloursimulationsintheright-handpanels ofFigure2. Weseethatintheisothermalandpolytropicsimulations,themassfunctionsdepend on the numerical resolution. By contrast, the stellar mass functions are well converged for the simulations with stellar heating feedback. 6. 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 approximatetheeffectsoffullradiationtransferofthestellarradiationfieldinacomputationally 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 distributionofformedstars,thetimeevolutionofthenumberofstars,andthecumulativestellar 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, wecautionthatgivenourcurrentresolutionsof∼ 100AU, thesesimulationsmaynot 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. Acknowledgements We thank the anonymous referee for their useful comments on the manuscript. C.F. gratefully acknowledges funding provided by the Australian Research Council’s Discovery Projects (grants DP150104329 and DP170100603). The simulations presented in this work used high performance computing resources provided by the Leibniz Rechenzentrum and the Gauss Centre for Supercomputing (grants pr32lo, pr48pi and GCS Large-scale project 10391), the Partnership for Advanced Computing in Europe (PRACE grant pr89mu), the Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre with funding from the Australian Government and the Government of Western Australia, in the 8 ASTRONUM 2016 IOP Publishing IOP Conf. Series: Journal of Physics: Conf. Series 1823374 (526071879) 0012007 doi :10.1088/1742-6596/837/1/012007 Figure 2. Left-hand panels: the number of sink particles as a function of SFE in runs with the isothermal EOS (top panels), the polytropic EOS (middle panels), and with stellar heating feedback (bottom panels). Right-hand panels: same as left-hand panels, but showing the cumulative stellar mass functions for all our simulations at SFE=10%. Each panel includes runs with three different numerical resolutions: high-resolution (solid), mid-resolution (dashed), and low-resolution (dotted). Only the simulations with stellar heating show numerical convergence in the number of sink particles and their stellar mass functions. framework of the National Computational Merit Allocation Scheme and the ANU Allocation Scheme. The simulation software FLASH was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago. References [1] Salpeter E E 1955 ApJ 121 161 [2] Kroupa P 2001 MNRAS 322 231–246 (Preprint arXiv:astro-ph/0009005) [3] MuenchAA,LadaEA,LadaCJandAlvesJ2002ApJ 573366–393(Preprint arXiv:astro-ph/0203122) [4] Chabrier G 2003 Publ. Astron. Soc. Pacific 115 763–795 (Preprint arXiv:astro-ph/0304382) [5] Alves J, Lombardi M and Lada C J 2007 A&A 462 L17–L21 (Preprint arXiv:astro-ph/0612126) 9