The influence of the cloud virial parameter on the initial mass function

ABSTRACT

Crucial for star formation is the interplay between gravity and turbulence. The observed cloud virial parameter, |$\alpha _{\mathrm{vir}}$|⁠, which is the ratio of twice the turbulent kinetic energy to the gravitational energy, is found to vary significantly in different environments, where the scatter among individual star-forming clouds can exceed an order of magnitude. Therefore, a strong dependence of the initial mass function (IMF) on |$\alpha _{\mathrm{vir}}$| may challenge the notion of a universal IMF. To determine the role of |$\alpha _{\mathrm{vir}}$| on the IMF, we compare the star-particle mass functions obtained in high-resolution magnetohydrodynamical simulations including jet and heating feedback, with |$\alpha _{\mathrm{vir}}=0.0625$|⁠, 0.125, and 0.5. We find that varying |$\alpha _{\mathrm{vir}}$| from |$\alpha _{\mathrm{vir}}\sim 0.5$| to |$\alpha _{\mathrm{vir}}< 0.1$| shifts the peak of the IMF to lower masses by a factor of |$\sim 2$| and increases the star formation rate by a similar factor. The dependence of the IMF and star formation rate on |$\alpha _{\mathrm{vir}}$| is non-linear, with the dependence subsiding at |$\alpha _{\mathrm{vir}}< 0.1$|⁠. Our study shows a systematic dependence of the IMF on |$\alpha _{\mathrm{vir}}$|⁠. Yet, it may not be measurable easily in observations, considering the uncertainties, and the relatively weak dependence found in this study.

1 INTRODUCTION

One of the fundamental questions in star formation is the initial mass function (IMF), which is the stellar mass distribution in a newly formed cluster. The IMF is considered nearly universal based on observations (see reviews by Bastian, Covey & Meyer 2010; Kroupa et al. 2013; Krumholz 2014; Offner et al. 2014; Hopkins 2018; Lee et al. 2020; Hennebelle & Grudić 2024). Despite some variations in the functional form, most observational IMF models identify the existence of a peak at around |$0.3\, \mathrm{M_\odot }$| and a power-law tail at high masses given by |$\mathrm{ d}N \propto M^{-\Gamma }\, \mathrm{ d}\mathrm{log}M$|⁠, where N is the number of stars and M is the stellar mass (Salpeter 1955; Miller & Scalo 1979; Kroupa 2001; Chabrier 2005; Parravano, McKee & Hollenbach 2011). The slope (power-law on a logarithmic scale) of the high-mass end measured in the Milky Way and nearby galaxies is generally consistent with the Salpeter (1955) slope (Kroupa 2002; Bastian et al. 2010; Kroupa et al. 2013; Hopkins 2018), although cluster-to-cluster variations exist (Dib 2014; Dib, Schmeja & Hony 2017; Schneider et al. 2018; Hennebelle & Grudić 2024; Kroupa et al.2024).

Recent observations challenge the notion of a universal IMF. There is compelling evidence for a top-heavy IMF near the Galactic centre, i.e. having a higher fraction of high-mass stars compared to the IMF in the solar neighbourhood (e.g. Figer et al. 1999; Kim et al. 2006; Lu et al. 2013; Hosek Matthew et al. 2019). On the other hand, massive elliptical early-type galaxies (ETGs) are discovered to have a bottom-heavy IMF, i.e. a higher fraction of low-mass stars (e.g. Cenarro et al. 2003; van Dokkum & Conroy 2010; Treu et al. 2010; Cappellari et al. 2012; Conroy et al. 2013). Further, the slope of the IMF seems to be dependent on the star formation rate (SFR) of the galaxy, where the power-law slope becomes shallower with the increase in SFR (Gunawardhana et al. 2011). In addition, Kroupa et al. (2013) suggest that the underlying formation mechanism of brown dwarfs (BDs) and stars are likely different based on the observed BD binary statistics and propose separate mass functions for BDs and stellar objects. The need to treat the mass functions for BDs and stars separately would put further constraints on the modelling of the IMF parameters and interpreting the IMF shape in the subsolar range. An accurate theoretical model of the IMF must reproduce the typical IMF characteristics such as the presence of a peak in the low-mass end, the existence of a power-law tail at supersolar masses with a slope close to the Salpeter (1955) estimate, and at the same time address the observed outlier IMFs, apparent SFR dependence, and the BD-star issue discussed above. A pre-requisite to achieving such a model is understanding the dependence of the IMF on environmental conditions.

There have been significant contributions from theoretical and numerical works towards understanding the IMF. Some of the proposed candidates responsible for setting the peak mass include the Jeans mass at the density of the thermal coupling of gas to the dust (Li, Klessen & Mac Low 2003; Jappsen et al. 2005; Larson 2005; Elmegreen, Klessen & Wilson 2008; Grudić & Hopkins 2023), the Jeans mass associated with protostellar heating (Bate 2009b; Krumholz, Klein & McKee 2011; Guszejnov, Krumholz & Hopkins 2016; Federrath, Krumholz & Hopkins 2017; Hennebelle et al. 2020; Mathew & Federrath 2020), and the formation of the first hydrostatic core and tidal screening of the gas in the immediate environment by the core (Hennebelle, Lee & Chabrier 2019; Colman & Teyssier 2020). Another important mechanism extensively studied in the context of the IMF is the role of magnetic fields, which is found to reduce fragmentation by providing additional pressure support (Price, Bate & Dobbs 2008; Federrath 2015), although the extent of influence is debatable in the presence of protostellar feedback (Myers et al. 2014; Krumholz & Federrath 2019). Protostellar jets/outflows, which are driven by magnetic fields, induce fragmentation and reduce the median mass by a factor of |$\sim 2$|–3 (Li & Nakamura 2006; Wang et al. 2010; Federrath et al. 2014; Frank et al. 2014; Offner & Chaban 2017; Guszejnov et al. 2021; Mathew & Federrath 2021; Lebreuilly et al. 2024). Further, the impact of magnetic fields may differ for super-Alfv|$\acute{\text{e}}$|nic and trans-Alfv|$\acute{\text{e}}$|nic conditions, and also in a multiphase medium (Mac Low et al. 1998; Krumholz & Federrath 2019; Beattie et al. 2021; Seta & Federrath 2022). Assuming a correlation between the core mass function (CMF) and the IMF, supported by some observations (Motte, Andre & Neri 1998; Testi & Sargent 1998; Johnstone et al. 2000; Stanke et al. 2006; Alves, Lombardi & Lada 2007; Könyves et al. 2015), analytical theories of the CMF/IMF (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012) based on gravo-turbulent fragmentation suggest that the IMF is sensitive to the turbulent properties of the cloud, such as the Mach number and the mode of turbulence driving. Meanwhile, star formation and IMF theories based on a global hierarchical collapse (e.g. Ballesteros-Paredes et al. 2007; Vázquez-Semadeni, González-Samaniego & Colín 2017; Vázquez-Semadeni et al. 2019) and stochastic, competitive accretion are also prominent (e.g. Zinnecker 1982; Bonnell et al. 2001; Bate & Bonnell 2005; Basu, Gil & Auddy 2015). Recently, there have been numerous studies on the IMF in the early universe, with a focus on the influence of varying metallicities (e.g. Sharda & Krumholz 2022; Bate 2023; Li et al. 2023; Menon et al. 2024; Tanvir & Krumholz 2024; Yan et al. 2024).

While the above findings put crucial constraints, we are yet to reach a clear consensus on the origin of the IMF. A primary reason why it is difficult to isolate the mechanism(s) responsible for the observed IMF characteristics is because the effects of the mechanisms involved in star formation are inter-related. Thus, it is important to include all these mechanisms simultaneously and analyse the impact of the relative parameters associated with the cloud. The current analytical theories of the CMF/IMF are fundamentally based on the interaction between gravity and turbulence (Mac Low & Klessen 2004). The principal quantity that defines the relative importance of gravity versus turbulence is the cloud virial parameter |$\alpha _{\mathrm{vir}}$|⁠, defined as twice the ratio of turbulent kinetic to gravitational energy. It is typically found to be around unity, but the scatter in the measurements is large and the difference in |$\alpha _{\mathrm{vir}}$| between individual clouds can be more than an order of magnitude (Kauffmann, Pillai & Goldsmith 2013; Polzin et al. 2024). A strong dependence of the IMF on |$\alpha _{\mathrm{vir}}$| would place the concept of IMF universality under scrutiny. Therefore, we carry out a comprehensive study of the influence of the virial parameter on the star formation process, in particular the IMF, using simulations of star cluster formation. The simulations include gravity, turbulence, magnetic fields, protostellar heating, and mechanical feedback in the form of jets and outflows.

In Section 2, we outline the numerical methodology and the initial conditions of the simulations. In Section 3, we study the impact of the cloud virial parameter on the SFR and the IMF. We also compare the IMF and the multiplicity associated with simulations of different virial parameters with observations. In Section 4, we discuss some of the previous numerical studies on the effect of the cloud virial parameter on the IMF. The main conclusions are presented in Section 5.

2 METHODOLOGY

2.1 Numerical setup

Here, we provide a brief summary of the numerical methods used in this study, which build upon the simulation framework in Mathew & Federrath (2021). A more detailed description of the general setup can be found in section 2 of Mathew & Federrath (2021). We model star cluster formation by solving the magnetohydrodynamical (MHD) equations in the presence of gravity utilizing the flash (version 4) code (Fryxell et al. 2000; Dubey et al. 2008) with in-house modifications (Waagan, Federrath & Klingenberg 2011). In addition to gravity and magnetic fields, which are taken into account through the MHD equations, we also implement other important mechanisms such as turbulence, protostellar heating and jets/outflows in our simulations.

2.1.1 Turbulence driving

We employ a stochastic Ornstein–Uhlenbeck process (Eswaran & Pope 1988) to drive turbulence in the simulations, which generates an acceleration field that is included as a source term in the momentum equation of MHD (Federrath et al. 2010a, 2022). The turbulent energy is injected only on the largest scales, but it cascades down to smaller scales, producing a velocity power spectrum |$\sim k^{-2}$| or equivalently a velocity dispersion – size relation of |$\sigma _v \propto \ell ^{1/2}$| (Larson 1981; Ossenkopf & Mac Low 2002; Heyer & Brunt 2004; Roman-Duval et al. 2011; Federrath 2013; Federrath et al. 2021). In this study, we use a natural mixture of turbulence driving modes, i.e. equal power in compressive and solenoidal modes, which corresponds to a turbulence driving parameter of |$b \sim 0.4$| (see Federrath, Klessen & Schmidt 2008; Federrath et al. 2010a, 2022 for details on the turbulence driving method adopted here).

2.1.2 Sink particles and protostellar feedback

We use sink particles to replace the gas within the innermost parts of the collapsing regions in the simulations, which is a common practice in numerical works (Bate, Bonnell & Price 1995; Krumholz, McKee & Klein 2004; Federrath et al. 2010b). The introduction of sink particles is to prevent the simulations from reaching extremely small time-steps and stalling as a consequence of the rapid increase in the density within the bound regions. We set the sink particle radius |$r_{\mathrm{sink}} = 250\, \mathrm{au}$|⁠, i.e. sink particles will replace spherical volumes of gas with radius |$r_{\mathrm{sink}} = 250\, \mathrm{au}$| that satisfy several gravitational collapse and sink particle formation criteria (Federrath et al. 2010b), to avoid spurious formation of star particles. The sink particles inherit the position, linear momentum, and angular momentum of the gas enclosed within the respective spherical volume. The sink particles in our simulations therefore represent star + disc systems. To model protostellar heating and outflow feedback, we use subresolution models, which are approximations based on previous theoretical, numerical and observational studies. We encourage the reader to refer to Federrath et al. (2014) and Mathew & Federrath (2020, 2021) for the full description of the outflow and heating feedback models, respectively.

2.1.3 Limitations

We note that the minimum attainable grid cell size in our simulations is |$100\, \mathrm{au}$| (see Section 2.2) and the sink particles correspond to star–disc systems as mentioned above. While fragmentation in the extended discs (Bate, Bonnell & Bromm 2002; Goodwin & Whitworth 2007; Stamatellos, Hubber & Whitworth 2007; Stamatellos, Whitworth & Hubber 2011; Rogers & Wadsley 2012; Thies et al. 2015) can occur in our simulations, fragmentation on typical disc scales (see reviews by Kratter & Lodato 2016; Lee et al. 2020, and the references therein) is not fully resolved. However, we aim to analyse the relative differences in the IMFs on changing the virial parameter, which does not require resolving fragmentation entirely on all scales. Moreover, our comparison analysis with observations is aimed at the system IMFs (unresolved close binaries) rather than the canonical (individual star IMF).

2.2 Initial conditions

The computational domain of the simulations is a three-dimensional box with side length |$L = 2\, \mathrm{pc}$| and periodic boundary conditions. The grid structure embodies an adaptive mesh refinement framework with a maximum effective resolution of |$4096^3$| cells or equivalently a minimum cell size of |$100\,$| au. The initial gas density |$\rho _0$| and the magnetic field strength |$B_0$| in the simulations are uniform. The magnetic field is directed only along the z-axis initially. The initial cloud temperature is isothermal at |$10\, \mathrm{K}$|⁠, but we employ a polytropic equation of state (EOS) such that the gas pressure (and the temperature via the ideal gas equation) is dependent on the local density (see Mathew & Federrath 2021 for a detailed description of the polytropic EOS implemented here). The polytropic EOS is based on previous radiation-hydrodynamical simulations and theoretical works (Larson 1969; Yorke, Bodenheimer & Laughlin 1993; Masunaga & Inutsuka 2000; Offner et al. 2009). The steady state sonic Mach number is set as |$\mathcal {M} = \sigma _{\mathrm{v}}/c_{\mathrm{s}} = 5$|⁠, where |$\sigma _{\mathrm{v}} = 1.0\, \mathrm{km\, s^{-1}}$| is the velocity dispersion on the driving scale and |$c_{\mathrm{s}} = 0.2\, \mathrm{km\, s^{-1}}$| is the isothermal sound speed. Once turbulence driving starts, filamentary and clump-like overdense regions start to emerge via turbulent shocks. The stirring of the gas also makes the magnetic field morphology inhomogeneous due to the tangling, compression and elongation of magnetic field lines by the turbulent flow (Seta et al. 2020; Seta & Federrath 2021), producing a field structure comparable to that in actual molecular clouds (Federrath 2016).

To allow turbulence to fully develop and reach a steady state, we drive turbulence in the simulations without self-gravity for two turbulent crossing times, |$2\, t_\mathrm{turb}=L/(\mathcal {M}c_\mathrm{s}) = 2\, \mathrm{Myr}$| (Price & Federrath 2010; Federrath et al. 2010a). Once a steady state of turbulence is reached, we turn on self-gravity, which we define as time |$t = 0$|⁠. We allow the simulations to progress until a star formation efficiency (SFE) of 5 per cent is reached, i.e. when 5 per cent of the total cloud mass has been converted into stars, during which the turbulence driving is also sustained (Mac Low 1999). Throughout this phase of evolution, the overdense regions (analogous to dense cores) in the shocked gas become unstable and form stellar clusters.

With the numerical setup and initial conditions discussed above, we carry out three sets of simulations with different initial cloud virial parameters. For a homogeneous spherical cloud of radius |$R_{\mathrm{cl}}$|⁠, total mass M, and 3D turbulent velocity dispersion |$\sigma _v$|⁠, the cloud virial parameter is given by

where |$E_{\mathrm{turb}}=M\sigma _v^2/2$| and |$E_{\mathrm{grav}}=-3GM^2/(5R_{\mathrm{cl}})$|⁠, with |$R_{\mathrm{cl}} = L/2$| and |$M = \rho _0 L^3$|⁠. It is important to highlight that the cloud structure becomes highly inhomogeneous and the actual value of |$\alpha _{\mathrm{vir}}$| at the time when a steady state of turbulence is reached can be significantly different from the simple spherical approximation given by equation (1) of a homogeneous, gravitationally isolated cloud (Federrath & Klessen 2012). However, we refer to |$\alpha _{\mathrm{vir}}$| from equation (1) to simply denote the distinction between the three simulation models.

To investigate the impact of the virial parameter on the cloud evolution, we conduct three sets of simulations with |$\alpha _{\mathrm{vir}}= 0.0625, 0.125$|⁠, and 0.5 by varying the initial mean gas density |$\rho _0$| in the simulations (keeping the Mach number |$\mathcal {M}$| constant). We adjust the magnetic field strength |$B_0$| in each of the three simulation sets to keep the Alfvén Mach number |$\mathcal {M_\mathrm{A}}$| fixed. Maintaining constant values of |$\mathcal {M_\mathrm{A}}$| across the three models is crucial to ensure that the observed variations are independent of the magnetic fields (Mac Low et al. 1998; Krumholz & Federrath 2019; Beattie et al. 2021). We emphasize that rather than the absolute quantities like |$\sigma _v$| and |$B_0$|⁠, the relative (dimensionless) quantities like |$\mathcal {M}$| and |$\mathcal {M_\mathrm{A}}$| should be given greater consideration since the fundamental physical processes at play are better reflected in such relative (dimensionless) quantities. The mass-to-flux ratio |$\mu _{\mathrm{B}}$| is different in the three models because of the adjustments in |$\rho _0$| and |$B_0$|⁠, but the variations are not significant enough to affect the results. Table 1 lists initial values associated with important quantities in each model.

It is to be noted that changing |$\alpha _{\mathrm{vir}}$| without affecting any of the other dimensionless cloud parameters is not possible as many of them are mathematically and physically related. We choose to adjust |$\alpha _{\mathrm{vir}}$| by changing the mean gas density, which also changes the mean thermal Jeans mass (and the number of Jeans masses) but keeps the other crucial parameters like |$\mathcal {M}$| and |$\mathcal {M_{\mathrm{A}}}$| fixed (see Table 1). Alternatively, one could adjust |$\alpha _{\mathrm{vir}}$| by changing the velocity dispersion or the cloud size. However, changing the velocity dispersion would affect the Mach number, which would have distinct consequences on the cloud evolution (see Section 4). Changing the cloud size would vary the number of Jeans masses too. Moreover, since |$\alpha _{\mathrm{vir}}\propto L^{-2}$|⁠, the cloud size has to be increased to obtain lower virial parameters, which would make performing the corresponding simulations challenging due to the increase in the computational volume and therefore the computational expense.

3 RESULTS

For |$\alpha _{\mathrm{vir}}= 0.0625$|⁠, we carried out one simulation, which produced 517 sink particles (young stellar objects). For |$\alpha _{\mathrm{vir}}= 0.125$|⁠, we performed two simulations with different random seeds for the turbulence driving (Ornstein–Uhlenbeck process; see Section 2, but otherwise identical cloud properties), which generated a total of 520 sink particles. In the case of |$\alpha _{\mathrm{vir}}= 0.5$|⁠, we performed 10 simulations with different turbulence driving seeds, producing a total of 449 sink particles. In each of the models, the choice of the number of simulations performed is based on the aim to achieve a similar number of sink particles or similar IMF statistics among the three models. With the target of obtaining a statistically representative sample of the IMF, 10 simulations (turbulence seeds) were needed to produce a significant number of sinks or star particles (⁠|$> 400$|⁠), while for |$\alpha _{\mathrm{vir}}=0.125$| and 0.0625, only 2 and 1 simulation(s), respectively, were required to produce a similar overall number of sink particles. For |$\alpha _{\mathrm{vir}}= 0.125$| and 0.5, different quantities studied here are averaged (or compiled together) over the set of simulations in each of the models (2 and 10 simulations, respectively), while for |$\alpha _{\mathrm{vir}}= 0.0625$| it is from the single simulation performed.

Fig. 1 shows the mass-weighted gas density projection for the three simulation models at |$\mathrm{SFE}=5~{{\ \rm per\ cent}}$|⁠. For |$\alpha _{\mathrm{vir}}= 0.5$| (right panel), the high-density structures are localized and there are mainly three star-forming regions, while in the lower virial parameter cases of |$\alpha _{\mathrm{vir}}= 0.0625$| (left panel) and |$\alpha _{\mathrm{vir}}= 0.125$| (middle panel), the high-density regions are spread across the simulation box with a large network of star-forming regions. Since gravity is more dominant in the lower virial parameter models, it is able to amplify the anisotropies produced by turbulence more strongly, creating additional sheet-like and filamentary structures (Vázquez-Semadeni et al. 2007, 2019; Bonnell, Clark & Bate 2008; Federrath 2015; Matsumoto, Dobashi & Shimoikura 2015; Hennebelle & Grudić 2024). Thus, a higher fraction of the cloud can form stars. 3D visualizations of the respective simulations are shown in Section  A.

3.1 Gas density probability distribution function

Previous theoretical and numerical works find that the form of the density probability distribution function (PDF) of supersonic turbulence is approximately lognormal (Vazquez-Semadeni 1994; Padoan, Nordlund & Jones 1997; Kritsuk et al. 2007; Federrath et al. 2008; Federrath 2013; Hopkins 2013; Federrath & Banerjee 2015). In the presence of gravity, the PDF develops a power-law tail in the high-density regime (Klessen 2000; Ballesteros-Paredes et al. 2011; Kritsuk, Norman & Wagner 2011; Collins et al. 2012; Federrath & Klessen 2013; Girichidis et al. 2014; Matsumoto et al. 2015; Myers 2015; Schneider et al. 2015; Mocz et al. 2017; Padoan et al. 2017; Burkhart 2018). With time and increasing SFE, the slope of the power-law part (in logarithmic scale) becomes shallower (Federrath & Klessen 2013), and the density at which the PDF transitions from lognormal to a power-law decreases (Burkhart 2018).

Fig. 2 shows the PDFs of the volume-weighted logarithmic density contrast |$s = \mathrm{ln}(\rho /\rho _0)$| in our simulations with |$\alpha _{\mathrm{vir}} = 0.0625$| (green, dash–dotted), 0.125 (blue, solid), and 0.5 (red, dashed) at the moment the first sink particle forms (left panel) and at SFE = 5  per cent (right panel). The PDFs obtained from the simulations are plotted with circular markers and the curves denote the corresponding fits. Following the method in Khullar et al. (2021), we fit a lognormal + power-law function to the density PDFs, where the free parameters are the width of the lognormal part (⁠|$\sigma _\mathrm{s}$|⁠) and the slope of the power-law part (⁠|$\alpha _\mathrm{s}$|⁠). The form of the fit function is given by

For defining the parameters like the peak of the lognormal part (⁠|$s_0$|⁠), transition density (⁠|$s_\mathrm{t}$|⁠), and the normalization constants (⁠|$N_\mathrm{d}$| and |$k_\mathrm{d}$|⁠), we use the corresponding expressions derived in Khullar et al. (2021) based on constraints like continuity and differentiability at the transition density, mass conservation and normalization of the PDF to unity (see section 2.2 in Khullar et al. 2021 for more details). The derived parameter values of the fitted function are shown in the legend of Fig. 2.

At the time the first star particle forms (SFE = 0 per cent; left panel), the density PDF for the |$\alpha _{\mathrm{vir}} = 0.5$| model has a slightly higher fraction of dense gas (high s) than those with |$\alpha _{\mathrm{vir}} = 0.0625$| and |$\alpha _{\mathrm{vir}} = 0.125$|⁠. The reason for the relatively lower fraction in the simulations with |$\alpha _{\mathrm{vir}} = 0.0625$| and |$\alpha _{\mathrm{vir}} = 0.125$| is that stars start to form much earlier in these simulations (at |$t\sim 0.05-0.1\, \mathrm{Myr}$|⁠) compared to those with |$\alpha _{\mathrm{vir}} = 0.5$|⁠, where the first stars form at |$t\sim 0.2\, \mathrm{Myr}$|⁠. Therefore, the clouds corresponding to |$\alpha _{\mathrm{vir}} = 0.0625$| and |$\alpha _{\mathrm{vir}} = 0.125$| are comparatively less dynamically evolved when compared at the moment the first star forms. Conversely, at SFE = 5 per cent (right panel), the simulations with |$\alpha _{\mathrm{vir}} = 0.0625$| and |$\alpha _{\mathrm{vir}} = 0.125$| have a higher fraction of high-s regions than that for |$\alpha _{\mathrm{vir}} = 0.5$|⁠. This is also apparent when comparing the slope of the power-law part |$\alpha _\mathrm{s}$| of the density PDF between the simulation models (compare legends in Fig. 2). In the left panel, |$\alpha _\mathrm{s}$| in the model with |$\alpha _{\mathrm{vir}} = 0.5$| is shallower than that of the models with |$\alpha _{\mathrm{vir}} = 0.0625$| and |$\alpha _{\mathrm{vir}} = 0.125$|⁠, while in the right panel, it is steeper.

Fig. 3 shows the dependence of |$\alpha _\mathrm{s}$| on SFE. It is evident that |$\alpha _\mathrm{s}$| remains almost constant with time (SFE) in the case of |$\alpha _{\mathrm{vir}}= 0.5$|⁠, which agrees with Khullar et al. (2021) who find the same for |$\alpha _{\mathrm{vir}}=0.5-2$| (also, see Appel et al. 2022, 2023). This is because, in the case of high-virial parameters, stars primarily form in the high-density regions initially seeded by turbulence, i.e. the mass reservoir is pre-determined by the initial turbulence properties. On the other hand, |$\alpha _\mathrm{s}$| decreases (slopes become shallower) with SFE for |$\alpha _{\mathrm{vir}}= 0.0625$| and 0.125, which is because low-virial parameters imply that gravity is more efficient in enhancing the dense structures in the cloud. Once gravity has had sufficient time to amplify the anisotropies and develop filamentary channels, simulations with lower virial parameters develop overdensities, and consequently form stars at a higher rate, which is discussed next.

3.2 Star formation rate

The star formation rate is given by |$\mathrm{SFR} = M_*/\tau _{\mathrm{c}}$|⁠, where |$M_*$| is the cloud mass converted into stellar mass and |$\tau _{\mathrm{c}}$| is the average time-scale of the conversion process. SFR is generally referenced in its dimensionless form, i.e. in mass fraction per freefall time, given by (Krumholz & McKee 2005)

where |$M_{\mathrm{cl}}$| and |$t_{\mathrm{ff,\, \rho _0}}$| are the total mass of the cloud and the freefall time at the mean density |$\rho _0$|⁠, respectively.

Fig. 4 shows SFE (top panel) and |$\mathrm{SFR_{ff}}$| (bottom panel) as a function of time in our simulations with |$\alpha _{\mathrm{vir}}= 0.0625$| (dash–dotted), 0.125 (solid), and 0.5 (dashed). It is clear that, when considered at the same time, |$\mathrm{SFR_{ff}}$| is higher for lower virial parameters, with the highest |$\mathrm{SFR_{ff}}$| in the simulation with |$\alpha _{\mathrm{vir}}= 0.0625$|⁠. The offset between the |$\mathrm{SFR_{ff}}$| curves for |$\alpha _\mathrm{vir} = 0.0625$| and |$\alpha _\mathrm{vir} = 0.125$| is simply because the onset of star formation is earlier in |$\alpha _\mathrm{vir} = 0.0625$| than for |$\alpha _\mathrm{vir} = 0.125$|⁠. The lower the value of |$\alpha _{\mathrm{vir}}$|⁠, the faster gravity acts in rearranging gas to form stars. However, once the star formation gets underway, the rate of change of |$\mathrm{SFR_{ff}}$| is similar for |$\alpha _{\mathrm{vir}}= 0.0625$| and 0.125, which becomes more evident when |$\mathrm{SFR_{ff}}$| is compared at the same SFE. Fig. 5 shows the evolution of |$\mathrm{SFR_{ff}}$| with SFE in our three simulation models. Here, the |$\mathrm{SFR_{ff}}$| values for |$\alpha _{\mathrm{vir}}= 0.0625$| and 0.125 are similar, while |$\mathrm{SFR_{ff}}$| for |$\alpha _{\mathrm{vir}}= 0.5$| is lower. When |$\alpha _{\mathrm{vir}}$| is sufficiently low, the process is more strongly dictated by the gravitational interactions than other mechanisms and the SFR accelerates with time. In other words, when the virial parameter is sufficiently low, SFR is almost independent of the value of |$\alpha _{\mathrm{vir}}$|⁠. Hence, the rate of change of SFR is similar for |$\alpha _\mathrm{vir} = 0.0625$| and |$\alpha _\mathrm{vir} = 0.125$|⁠. However, in the case of |$\alpha _\mathrm{vir} = 0.5$|⁠, the turbulent support still plays a substantial role throughout the process, and thus the acceleration in SFR is not as significant. The |$\mathrm{SFR_{ff}}=0.30$|⁠, 0.29, and 0.14, averaged over the SFE range |$1-5~{{\ \rm per\ cent}}$|⁠, for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125$|⁠, and 0.5, respectively.

3.3 Initial mass function

Fig. 6 presents the mass distribution of sink particles (star + disc systems) obtained in the simulations with |$\alpha _{\mathrm{vir}}= 0.0625$| (histogram with dash–dotted edges), 0.0125 (histogram with solid edges), and 0.5 (histogram with dashed edges). The curves represent a modified version of the Chabrier (2005) IMF fitted to the simulation data using the Markov Chain Monte-Carlo (MCMC) sampler emcee of Foreman-Mackey et al. (2013) (see also Nam, Federrath & Krumholz 2021). The input to the sampling algorithm is the list of sink particle masses rather than the binned data, and hence the method has the advantage that the fit is not affected by the binning choice. The modified version is similar to the standard Chabrier (2005) IMF form except that, to take into consideration the finite mass in our computational domain, we include an exponential term that serves as a smooth cutoff in the high-mass regime of the Chabrier (2005) IMF,

where all masses are in units of |$\mathrm{M_\odot }$|⁠. Here, the appropriate values for the five free parameters |$\theta _{\mathrm{fit}} = (\sigma , \mathrm{log}\, M_0, \mathrm{log}\, M_{\mathrm{T}}, \Gamma , \mathrm{log}\, M_{\mathrm{cut}})$| are derived using MCMC sampling based on the simulated data. |$\sigma , M_0, M_{\mathrm{T}}$|⁠, and |$\Gamma$| are the standard deviation of the lognormal part, the peak mass, transition mass between the lognormal and power-law forms, and the slope of the power-law part, respectively. |$k_1$| and |$k_2$| are normalization factors, defined to ensure continuity at |$M_{\mathrm{T}}$|⁠. As a consequence of the exponential term in the power-law part, the IMF will be truncated at high masses. |$M_{\mathrm{cut}}$| characterizes the mass at which the exponential term begins to dominate and p describes how sharply the cutoff occurs around |$M_{\mathrm{cut}}$|⁠. We use |$p = 4$| to achieve a sufficiently sharp drop around |$M_{\mathrm{cut}}$|⁠. We refer the reader to section 3.2 in Mathew, Federrath & Seta (2023) for a more detailed description of the MCMC fitting technique used here.

The three stellar mass distributions are not substantially different, but they exhibit some distinctions. The derived values for the fit parameters corresponding to the simulated IMFs are listed in Table 2. The distributions for the lower virial parameters have a higher fraction of low-mass stars, which is expected since the large-scale turbulent support is weaker (or equivalently, gravity is comparatively more efficient), and therefore there is more fragmentation. The peak of the distributions for |$\alpha _{\mathrm{vir}}= 0.0625$| and |$\alpha _{\mathrm{vir}}= 0.125$| are comparable, however, the peak for |$\alpha _{\mathrm{vir}}= 0.5$| is at a relatively higher mass (see Table 2). Changing the virial parameter in the simulation from |$\alpha _{\mathrm{vir}}= 0.5$| to |$\alpha _{\mathrm{vir}}= 0.0625$|⁠, i.e. by a factor of 8, reduces the peak mass by a factor of |$\sim 2.4$|⁠. This shows that the IMF has a relatively weak, but systematic dependence on the cloud virial parameter. Lower |$\alpha _{\mathrm{vir}}$| suggests that the overdensities are more bound and stable against disruption by shocks (Federrath & Klessen 2012; Bertelli Motta et al. 2016), and hence more low-mass overdensities can collapse to form stars.

The slopes |$\Gamma$| of the fit at the high-mass end for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125,$| and 0.5 are |$0.90_{-0.30}^{+0.60}, 1.60_{-0.40}^{+0.70},$| and |$1.70_{-0.40}^{+0.90}$|⁠, respectively (note that the Salpeter slope would be |$\Gamma =1.35$| in this definition of |$\Gamma$|⁠). The slope for |$\alpha _{\mathrm{vir}}= 0.0625$| is shallower than that of the other |$\alpha _{\mathrm{vir}}$| values. However, since the high-mass range is narrow in our simulations, the error bars on the |$\Gamma$| estimates are large and therefore it is difficult to make conclusive remarks on the dependence of the high-mass regime of the IMF on the cloud virial parameter just based on |$\Gamma$|⁠.

The median masses of the IMFs at SFE = 5 per cent for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125,$| and 0.5 are |$\sim 0.26,\, 0.32\, $|⁠, and |$0.53\, \mathrm{M_\odot }$|⁠, respectively. The average masses at SFE = 5 per cent are |$\sim 0.59,\, 0.60\,$|⁠, and |$0.86\, \mathrm{M_\odot }$|⁠, respectively. We also evaluate |$M_{\mathrm{50}}$|⁠, which is, as outlined in Krumholz, Klein & McKee (2012), the median mass in a cumulative mass distribution, i.e. |$M_{\mathrm{50}}$| is the value below which 50 per cent of the total stellar mass is found. We find |$M_{\mathrm{50}} \sim 1.32,\, 1.17\,$|⁠, and |$1.56 \, \mathrm{M_\odot }$| at SFE = 5 per cent for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125,$| and 0.5, respectively. Fig. 7 shows the evolution of the median mass (top panel), the average mass (middle panel), and |$M_{50}$| (bottom panel) with SFE in the three |$\alpha _{\mathrm{vir}}$| models. The median mass remains relatively constant from |$\mathrm{SFE\sim 1}$| for the three models, while the average mass increases moderately with SFE. The increase in the average mass reflects the development of the high-mass end of the IMF with time. The median, mean, and |$M_{50}$| masses averaged over the SFE range |$1-5~{{\ \rm per\ cent}}$| for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125,$| and 0.5 are listed in Table 1.

3.4 Comparison with observational IMF

Fig. 8 depicts the comparison between our simulation IMFs and the observational IMF models. We compare our data with the system IMFs (unresolved close binaries) rather than the individual-star IMF since the accretion discs in our simulations are not fully resolved. All the observational models broadly agree with each other, though there are slight variations in the low-mass regime. While Parravano et al. (2011) propose a higher fraction of very low-mass stars or brown dwarfs than the Chabrier (2005) IMF, Da Rio et al. (2012) and Damian et al. (2021) suggest a lower fraction than Chabrier (2005). Such a difference also exists between the IMFs for our three simulation models, where the simulation corresponding to |$\alpha _{\mathrm{vir}}= 0.0625$| has the highest fraction of very low-mass objects and the simulations with |$\alpha _{\mathrm{vir}}= 0.5$| have the lowest fraction.

For comparisons between observations and simulations, it is also beneficial to measure the ratio between the number of stars in different mass ranges, which is more robust than deriving quantities such as the peak mass or the power-law slope. Observations find that the ratio of the number of substellar objects (⁠|$M \le 0.08\, \mathrm{M_\odot }$|⁠) to that with stellar masses (⁠|$0.15 < M \le 1.0\, \mathrm{M_\odot }$|⁠) is |$\sim 0.2$|⁠, i.e. one brown dwarf (substellar) is formed per every five late-type (subsolar) stars on average (Andersen et al. 2006, 2008; Thies & Kroupa 2007; Parravano et al. 2011; Kroupa et al. 2013). The corresponding ratios for the simulation models with |$\alpha _{\mathrm{vir}}= 0.0625, 0.125,$| and 0.5 are |$0.41, 0.25$|⁠, and 0.19, respectively. The ratio for |$\alpha _{\mathrm{vir}}= 0.5$| agrees with the typical value from observations, while that for |$\alpha _{\mathrm{vir}}= 0.0625$| is higher by a factor of |$\sim 2$|⁠.

3.5 Multiplicity

The multiplicity statistics is closely linked to the IMF and is a critical probe to test star formation theories. We identify multiple systems in our simulations based on the procedure employed in Bate (2009a), which we briefly summarize. In the list of stars (single objects) obtained from the simulations, we search for the closest gravitationally bound pair. The corresponding two stars are removed from the list and replaced by a binary object with the mass, centre-of-mass position, and centre-of-mass velocity of the removed bound pair. In the new list, we again find the closest bound pair. In the scenario where the pair consists of a single and a binary object, the pair is substituted by a triple object. The above algorithm is repeated until no new bound pairs can be found in the list or the only possible option is a quintuple. We rule out pairings that result in a quintuple or systems of higher order since such high-order multiple systems tend to be dynamically unstable and are likely to decay to lower-order systems as the cloud evolves. This procedure is the same as the one used in Mathew & Federrath (2021) and Mathew et al. (2023).

With the list of singles, binaries, triples, and quadruples derived from the above technique, we calculate the multiplicity fraction in different mass ranges. The multiplicity fraction (⁠|$\mathrm{ \mathrm{mf}}$|⁠) in a mass range is the ratio of the number of multiple systems to the total number of systems whose primary star falls within the specified mass range, i.e.

where |$S, B, T$|⁠, and Q represent the number of singles, binaries, triples, and quadruples, respectively.

We show |$\mathrm{mf}$| as a function of the primary mass for the three simulation models in Fig. 9. We see that |$\mathrm{mf}$| increases with primary mass, which is consistent with the established understanding (Bate 2012; Krumholz et al. 2012; Duchêne & Kraus 2013; Cunningham et al. 2018; Sharda, Federrath & Krumholz 2020; Offner et al. 2023). It is evident that, at subsolar masses, |$\mathrm{mf}$| for |$\alpha _{\mathrm{vir}}= 0.0625$| and |$\alpha _{\mathrm{vir}}= 0.125$| in each of the mass intervals are higher than that for |$\alpha _{\mathrm{vir}}= 0.5$|⁠. The multiplicity fraction over the subsolar range (primary mass |$< 1\, \mathrm{M_\odot }$|⁠) for |$\alpha _{\mathrm{vir}}= 0.0625, 0.125$|⁠, and 0.5 are |$\sim 0.27, 0.22,$| and 0.14, respectively. The higher |$\mathrm{mf}$| in lower |$\alpha _{\mathrm{vir}}$| cases might be because the local region in which the binaries form will be more bound and gas-rich in the case of low virial parameters. Therefore, the binaries are better shielded from dynamical interactions that result in binary decay (Rozner & Perets 2024). At supersolar masses, the derived |$\mathrm{mf}$| values are similar and |$\gtrsim 0.6$| for the three cases. High-mass binaries generally form in the most gas-rich regions within the cloud, which is why they can accrete gas rapidly. Hence, the high-mass binaries are relatively shielded from interactions even in the highest virial parameter case of |$\alpha _{\mathrm{vir}}= 0.5$|⁠. We note that, to confirm our interpretation of the higher |$\mathrm{mf}$| in lower |$\alpha _{\mathrm{vir}}$| cases, a comprehensive quantitative study is required, measuring the local virial parameter in cluster regions and analysing the evolution of binary orbits, which is beyond the scope of this work. Finally, we find that the fraction of singles across the full mass range is |$\sim 0.65$|⁠, irrespective of |$\alpha _{\mathrm{vir}}$|⁠. Previous numerical and observational works also yield similar values for the fraction of single stars (Lada 2006; Mathew & Federrath 2021; Rohde et al. 2021).

4 DISCUSSION

4.1 Comparison with previous numerical works

In our study, the mean thermal Jeans mass |$M_{\mathrm{J}}$| and the number of Jeans masses for the three simulations models are different since the initial mean density varies between them. Bate & Bonnell (2005) find that the IMF is linearly dependent on the mean thermal Jeans mass in their simulations. However, Bate (2009b) establish that the inclusion of radiative feedback removes the dependence on the thermal Jeans mass, finding that the IMFs are indistinguishable for simulations with mean Jeans masses of |$1\, \mathrm{M_\odot }$| and |$1/3\, \mathrm{M_\odot }$| (see also Krumholz et al. 2011). In our study, the mean Jeans mass in the simulations with |$\alpha _{\mathrm{vir}}= 0.0625$| and |$\alpha _{\mathrm{vir}}= 0.125$| also varies by a factor of |$\sim 3$|⁠. However, contrary to our case, |$\alpha _{\mathrm{vir}}$| is similar for the two simulation sets with different |$M_{\mathrm{J}}$| in Bate (2009b) as a result of using a different combination of cloud size and |$\mathcal {M}$| in each of their two sets. Therefore, a direct comparison between the two studies cannot be made.

Lee & Hennebelle (2018) performed a similar numerical experiment to Bate (2009b) by changing the gas density (and hence the thermal Jeans mass) while keeping |$\alpha _{\mathrm{vir}}$| fixed, by adjusting the cloud size and |$\mathcal {M}$|⁠. They find that the IMF is relatively flat in the |$0.1-1\, \mathrm{M_\odot }$| range at low densities, while at higher densities, the distribution starts to become more peaked, producing a power-law distribution with |$\Gamma \sim 3/4$| at masses greater than a few |$0.1\, \mathrm{M_\odot }$|⁠. However, stellar radiative heating was absent in their simulations, and therefore the influence of varying densities or thermal Jeans masses might be overestimated (Bate 2009b; Krumholz et al. 2011; Guszejnov et al. 2016; Federrath et al. 2017; Hennebelle et al. 2020; Mathew & Federrath 2020).

Ballesteros-Paredes et al. (2015) and Bertelli Motta et al. (2016) investigated the effect of the turbulent Mach number and found that the form of the IMF remains the same irrespective of |$\mathcal {M}$| in their simulations. This contradicts the analytical theories of the CMF or the IMF, which predict that the characteristic mass decreases with an increase in the Mach number (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Hopkins 2012). According to these theories, a higher Mach number results in stronger shocks, leading to greater density contrasts and higher fragmentation. The observed trends in Ballesteros-Paredes et al. (2015) and Bertelli Motta et al. (2016) arise because raising the Mach number through the increase of the velocity dispersion also increases the cloud virial parameter, as |$\alpha _{\mathrm{vir}}\propto \sigma _{\mathrm{v}}^2$|⁠. Higher virial parameters imply that overdensities are less susceptible to collapse, with only the sufficiently massive ones eventually collapsing (Bertelli Motta et al. 2016). On the other hand, Lee & Hennebelle (2018) observe slight variations in the IMF on changing the Mach number in their simulations by a factor of 5. Here again, the Mach number is varied by adjusting the velocity dispersion, which alters the virial parameter.

Guszejnov et al. (2022) find that a higher virial parameter results in a lower star formation rate, which agrees with our results. They find that variations in the virial parameter of the cloud do not significantly affect the form of the IMF, while in our study we see a weak dependence of the IMF on |$\alpha _{\mathrm{vir}}$|⁠. However, Guszejnov et al. (2022) also vary the virial parameter by adjusting the velocity dispersion, which in turn affects |$\mathcal {M}$|⁠. Since a higher Mach number is expected to enhance the formation of low-mass overdensities and a higher virial parameter is expected to suppress them, the net impact on the IMF is expected to be minimal as seen in the simulations of Ballesteros-Paredes et al. (2015), Bertelli Motta et al. (2016), and Lee & Hennebelle (2018). At very high-Mach number (or very high |$\sigma _\mathrm{v})$|⁠, the effect of the virial parameter (increased kinetic support) can dominate over the effect of the Mach number (density-enhancement via stronger shocks) and quench star formation altogether since |$\alpha _{\mathrm{vir}}\propto \sigma _{\mathrm{v}}^2$|⁠, while |$\mathcal {M} \propto \sigma _{\mathrm{v}}$|⁠. However, the details of the dependence are non-linear and depend on the region of parameter space considered (Federrath & Klessen 2012).

The numerical experiments carried out here are distinct from the numerical approaches discussed above. While the above works concentrate on the influence of the mean thermal Jeans mass (keeping the virial parameter constant) or the Mach number (keeping the mean thermal Jeans mass constant), i.e. in a sense the impact of the mean density or the velocity dispersion, we focus on the effect of the virial parameter, but keep the Mach number fixed. The mean density and hence the mean thermal Jeans mass varies between our simulation models, which could have some effect on the initial fragmentation in the cloud. However, since the turbulence in our simulations is supersonic (⁠|$\mathcal {M} = 5$|⁠), the kinetic energy dominates over the thermal energy and hence will govern most of the fragmentation in the early stages. Further, once the stars start to form, the protostellar heating increases the local temperature and Jeans mass in the star-forming cores on subparsec scales (Bate 2009b; Krumholz et al. 2011; Guszejnov et al. 2016; Federrath et al. 2017; Hennebelle et al. 2020; Mathew & Federrath 2020). Thus, the IMF variations that we observe between the three simulation models may be determined more by the variation in |$\alpha _{\mathrm{vir}}$| than by variations in the mean thermal Jeans mass.

Haugbølle, Padoan & Nordlund (2018) carried out a numerical experiment similar to our analysis, where they changed the virial parameter by varying the mean density, but keeping |$\mathcal {M}$| fixed. They find that lowering |$\alpha _{\mathrm{vir}}$| shifts the peak of the IMF to lower masses, which aligns with our findings. However, Haugbølle et al. (2018) focus on a regime with |$\alpha _{\mathrm{vir}}\sim 0.2-2$|⁠, while we concentrate on a lower |$\alpha _{\mathrm{vir}}$| regime (⁠|$0.06-0.5$|⁠). Haugbølle et al. (2018) study the |$\alpha _{\mathrm{vir}}$| range where either the gravitational and turbulent energies are comparable (⁠|$\alpha _{\mathrm{vir}}\sim 0.5$|⁠) or turbulence dominates (⁠|$\alpha _{\mathrm{vir}}\gtrsim 1$|⁠), while we study the |$\alpha _{\mathrm{vir}}$| range where either the gravitational and turbulent energies are comparable (⁠|$\alpha _{\mathrm{vir}}\sim 0.5$|⁠) or gravity dominates (⁠|$\alpha _{\mathrm{vir}}\lesssim 0.1$|⁠). Thus, the present work is complementary to the works of Bate & Bonnell (2005), Bate (2009b), Ballesteros-Paredes et al. (2015), Bertelli Motta et al. (2016), Lee & Hennebelle (2018), Haugbølle et al. (2018), and Guszejnov et al. (2022). Ultimately, star formation may depend on a combination of all of the dimensionless parameters (e.g. |$\mathcal {M}$|⁠, |$\alpha _{\mathrm{vir}}$|⁠, b, |$\mu _{\mathrm{B}}$|⁠, and |$\mathcal {M_{\mathrm{A}}}$|⁠), and future work is needed to quantify the effects of each of these parameters on the SFR and IMF.

4.2 Outlier IMFs

The high-gas density and large velocity dispersion measured in massive, elliptical early-type galaxies (ETGs) have been suggested as candidates for explaining their bottom-heavy nature (Chabrier, Hennebelle & Charlot 2014; Dib 2022). However, evaluation of the virial parameter in these regions is crucial to obtain conclusive results since understanding the respective |$\alpha _{\mathrm{vir}}$| regime is important to comprehend the net effect of turbulence. Large velocity dispersion or very high |$\alpha _{\mathrm{vir}}$| can quench star formation because the large-scale support by turbulence dominates over its ability to enhance fragmentation on small scales (Federrath & Klessen 2012; Chabrier et al. 2014; Brucy et al. 2024). Another possible candidate that can partly contribute to the bottom-heavy nature is the mode of turbulence driving in such regions. Since ETGs are considered to be formed in a starburst event associated with galaxy mergers, the turbulence driving is likely dominated by compressive modes (Federrath et al. 2008; Renaud et al. 2009), which can enhance the formation of low-mass stars (Chabrier et al. 2014; Mathew et al. 2023). Further, our results suggest that the high-virial parameters (⁠|$\alpha _{\mathrm{vir}}> 1$|⁠) measured in the central molecular zone (CMZ) (Myers, Hatchfield & Battersby 2022) contribute in part to the low-SFR and top-heavy IMFs (Figer et al. 1999; Kim et al. 2006; Lu et al. 2013; Hosek et al. 2019) observed in the embedded clouds, with the prevalence of solenoidal driving (Federrath et al. 2016; Rani et al. 2022) being another factor (see also Klessen, Spaans & Jappsen 2007; Haugbølle et al. 2018; Mathew et al. 2023). We note that, given the non-linear IMF dependence found in our study, to understand the extent of the contribution of the virial parameter to the top-heavy nature of the IMF observed in the CMZ, regimes where |$\alpha _{\mathrm{vir}}> 1$| have to be considered, which are not directly explored in this study.

4.3 Origin of the IMF peak

The weak dependence of the IMF on the cloud properties like the virial parameter and mode of turbulence driving (Padoan & Nordlund 2002; Hennebelle & Chabrier 2008; Schmidt et al. 2010; Hopkins 2012; Mathew et al. 2023) may account for why the region around the IMF peak is broad or resembles a plateau (see also Kroupa et al. 2024). When the net effect of the combination of cloud properties favours extensive fragmentation, the peak shifts to lower masses, and when the combination tends to suppress fragmentation, the peak shifts to higher masses. None the less, it still needs to be addressed why the broad peak or plateau exists in the subsolar or more specifically in the M-dwarf range. Some preferred candidates are the feedback mechanisms that emerge from the stars themselves, such as radiative heating and protostellar outflows, which have been found to be effective in self-regulation (Li & Nakamura 2006; Bate 2009b; Wang et al. 2010; Krumholz et al. 2011; Federrath et al. 2014, 2017; Frank et al. 2014; Guszejnov et al. 2016; Offner & Chaban 2017; Hennebelle et al. 2020; Mathew & Federrath 2020, 2021; Guszejnov et al. 2021; Lebreuilly et al. 2024), or generally, the thermodynamics in the immediate vicinity of stars (Li et al. 2003; Jappsen et al. 2005; Larson 2005; Elmegreen et al. 2008; Hennebelle et al. 2019; Colman & Teyssier 2020; Grudić & Hopkins 2023).

5 CONCLUSIONS

We conduct a series of star cluster formation simulations to investigate the impact of the cloud virial parameter |$\alpha _{\mathrm{vir}}$| on the formation process. The simulations incorporate gravity, magnetic fields, turbulence, protostellar radiative heating, and mechanical feedback in the form of jets/outflows. We have three sets of simulations corresponding to |$\alpha _{\mathrm{vir}}= 0.0625, 0.125$|⁠, and |$\alpha _{\mathrm{vir}}= 0.5$|⁠. We observe that the number of stars formed increases with decreasing |$\alpha _{\mathrm{vir}}$|⁠. In the case of high virial parameters, the formation of large-scale structures that host star clusters is primarily a consequence of turbulent shocks. By lowering the virial parameter, the relative influence of gravity in the large-scale structure building process increases. Gravity is able to enhance the anisotropies in the cloud and promote the formation of additional large-scale filaments, which eventually fragment to form stars (Vázquez-Semadeni et al. 2007, 2019; Bonnell et al. 2008; Mac Low 2013; Hennebelle & Grudić 2024).

We find that reducing the virial parameter from |$\alpha _{\mathrm{vir}}= 0.5$| to |$\alpha _{\mathrm{vir}}=0.125$| increases the star formation rate per freefall time |$\mathrm{SFR_{ff}}$| by a factor of |$\sim 2$|⁠. Reducing the virial parmeter further to |$\alpha _{\mathrm{vir}}= 0.0625$| does not result in an additional increase in |$\mathrm{SFR_{ff}}$|⁠, showing that |$\alpha _{\mathrm{vir}}$| has a non-linear effect on |$\mathrm{SFR_{ff}}$| (see Fig. 5). Similar |$\mathrm{SFR_{ff}}$| in simulations with |$\alpha _{\mathrm{vir}}= 0.0625$| and |$\alpha _{\mathrm{vir}}= 0.125$| imply that, at |$\alpha _{\mathrm{vir}}\sim 0.1$|⁠, gravity has transformed most of the anisotropies or overdensities in the turbulent cloud into stars, i.e. the maximum possible efficiency of star formation at the given Mach number of |$\mathcal {M} = 5$| is reached. Thus, lowering |$\alpha _{\mathrm{vir}}$| further would not affect |$\mathrm{SFR_{ff}}$| significantly. Nevertheless, the time at which star formation starts varies between the three models. The first stars form the latest in the simulations with |$\alpha _{\mathrm{vir}}= 0.5$| and the earliest in the simulation with |$\alpha _{\mathrm{vir}}= 0.0625$| (see Fig. 4).

Previous numerical studies of how the initial turbulence influences the IMF are primarily centred on the aspect of the velocity dispersion, with most of the studies finding that the IMF is nearly insensitive to the velocity dispersion or Mach number (see Section 4.1). However, analysing the effect of turbulence on the basis of the velocity dispersion alone is inadequate since changes in |$\sigma _{\mathrm{v}}$| change both |$\alpha _{\mathrm{vir}}$| and |$\mathcal {M}$| at the same time, each of which has a unique influence on the cloud dynamics, structure, SFR, and IMF. Since in our simulations, we keep the velocity dispersion and therefore |$\mathcal {M}$| fixed, our analysis is principally focused on the effect of changes in the virial parameter. Thus, our work complements previous studies on the influence of turbulence.

We find that the IMF has a weak dependence on the cloud virial parameter. On varying |$\alpha _{\mathrm{vir}}$| in the simulations from |$\alpha _{\mathrm{vir}}\sim 0.5$| to |$\alpha _{\mathrm{vir}}\lesssim 0.1$|⁠, the peak mass of the IMF and the median mass decreases by |$\sim 2$|⁠. We also find that the multiplicity fraction, |$\mathrm{mf}$|⁠, is higher for the lower cloud virial parameter cases, particularly in the subsolar range. |$\mathrm{mf}$| in the subsolar primary mass range for simulations with |$\alpha _{\mathrm{vir}}= 0.0625$| and 0.125 are higher than that for |$\alpha _{\mathrm{vir}}= 0.5$| by a factor |$\sim 2$| (see Fig. 9). Such a difference is important given that this is the mass range where the peak of the IMF is expected. Nevertheless, all three models reproduce the trend of increasing multiplicity fraction with mass and the |$\mathrm{mf}$| estimates in individual mass intervals broadly agree with the values derived in observational surveys.

We highlight that the non-linear dependence of the SFR and IMF on |$\alpha _{\mathrm{vir}}$| needs to be considered (see Figs 5 and 7). If |$\alpha _{\mathrm{vir}}$| is investigated solely within a limited range of either very low or high values, its impact on the SFR and IMF will appear to be negligible, simply because either gravity or turbulence dominates over the other. To reach conclusive results, the effect of the virial parameter should be examined across a broad range of |$\alpha _{\mathrm{vir}}$| values.

ACKNOWLEDGEMENTS

We thank the anonymous referee for their quick and constructive review of the work. We would like to thank Blakesley Burkhart, Patrick Hennebelle, No|$\acute{\text{e}}$| Brucy, and Mark Krumholz for discussions regarding SFR theory. SSM would like to acknowledge the useful discussions during ‘The Physics of Star Formation’ winter school held in 2024 February at Les Houches, France. CF acknowledges funding provided by the Australian Research Council (Discovery Project DP230102280), and the Australia-Germany Joint Research Cooperation Scheme (UA-DAAD). We further acknowledge 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 Australian National Computational Infrastructure (grant ek9), and the Pawsey Supercomputing Centre (project pawsey0810) in the framework of the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme. The simulation software flash was in part developed by the DOE-supported Flash Center for Computational Science at the University of Chicago and the University of Rochester. This work makes use of the yt-project (Turk et al. 2011) and colourmaps in the cmasher package (van der Velden 2020).

DATA AVAILABILITY

The data used in this article is available upon reasonable request to the authors.

REFERENCES

APPENDIX A: 3D VISUALISATION OF THE SIMULATIONS

Figs A1, A2, and A3 present 3D visualizations of the simulations shown in Fig. 1. Movies for these visualizations can be seen at https://sajaymathew.github.io/visualisations.html.