© The Authors 2024

1 Introduction

Molecular clouds are the privileged environment for the formation of young stars. These stars are rarely born in isolation (Lada & Lada 2003) as they emerge in clusters or pair together in multiple systems that can extend from a few AU to several hundred (Duchêne & Kraus 2013) or thousand AU (Joncour et al. 2017). In addition,the newly born stars are characterised by an initial mass function (IMF) that represents the probability density function (PDF) of their mass distribution. Among star formation processes operating in molecular clouds, some models favour local stochastic processes to regulate (i) the mass of the stars, (ii) their spatial distribution, and (iii) the multiplicity of stellar groups through dynamical interactions between protostars themselves, for example coalescence (Dib et al. 2017); dynamical ejection or capture (Reipurth et al. 2014); stellar outflows (Adams & Fatuzzo 1996); or competitive accretion (Bonnell et al. 2001). These local stochastic processes at small scale would erase residual imprints of large-scale physical processes structuring both the cloud and the pre-stellar cores, regarding the IMF and multiplicity of stellar systems.

Alternatively, it has been proposed that the stellar IMF may be directly inherited from the mass distribution of the parental dense gas cores at a larger scale (Motte et al. 1998; Sadavoy et al. 2010; Offner et al. 2014). The increasing number of observations at different spatial resolutions and the improved performance of recent instruments have revealed that these dense structures of gas are actually made up of several buried substructures (Offner et al. 2012; Pokhrel et al. 2018; Thomasson et al. 2022). The existence of such substructures is often interpreted as the result of core fragmentation, although their physical properties may depend on the instrumental resolution (Louvet et al. 2021), which complicates their full characterisation.

Consequently, the fragmentation of molecular clouds all the way down to dense cores that may in turn subfragment to constitute the birth sites of young stellar systems is one of the driving processes shaping both the IMF and young stellar groups. In that case, the spatial distribution of young stars would simply mimic the structure of the cloud. To support this idea, multiplicity analysis of ultra-wide stellar systems in Taurus (Joncour et al. 2017) suggests a fragmentation cascade scenario for their formation. Further studies unravel the presence of dense stellar Nested Elementary STructures (NESTs) in star-forming regions (Joncour et al. 2018; González et al. 2021) whose origins may be attributed to the cloud fragmentation into gas cores clusters that in turn may fragment further into young stars. Hence, these stellar groups are suspected to be the pristine vestiges of star formation, and the emergence of these groups may be the result of physical mechanisms that regulate both the hierarchical fragmentation and the molecular cloud’s structure at larger scale. Hereafter we refer to a local concentration of stars, gravitationally bound or not, within a cloud as a stellar structure or stellar system.

Structures of different morphologies such as filaments, gas clumps, or pre-stellar cores may emerge at different spatial scales inside molecular clouds depending on the relevant physical processes shaping these scales (e.g. turbulence, magnetic field, rotation, or even stellar feedback). In that context, molecular clouds are characterised by their multi-scale nature. To investigate the origin of stellar systems through the fragmentation of the ISM, every process that influences the structure of the cloud at their relevant scale needs to be accounted for. In this work we focus on the fragmentation of gas clumps considering turbulence and the thermodynamics of gravitational collapse.

Quasi-static gravitational instabilities (Jeans 1902) constitute the basis of our current understanding concerning the local collapse of self-gravitating gas clumps into young stars. Throughout their evolution, molecular clouds condense into self-gravitating structures from which stars eventually form. A self-gravitating clump may collapse monolithically (Larson 1969; Shu 1977) under isothermal pressure free condition (Tohline 1980), but may also fragment in order to form multiple dense structures (Hoyle 1953; Guszejnov et al. 2018b; Vázquez-Semadeni et al. 2019). We explore in this work the possibility that fragmentation of molecular clouds shapes the stellar properties directly from cloud physical properties.

Provided that a self-gravitating clump neither loses mass during its evolution nor accretes mass from its neighbouring environment, the number of Jeans masses contained in the initial structure increases as it contracts quasi-statically (Hoyle 1953; Hunter 1962). New gravitationally unstable regions then appear within the initial structure, causing its subfragmentation and so on (Guszejnov et al. 2018a; Vázquez-Semadeni et al. 2019). The collapse of each fragment may eventually end as the gas structure becomes opaque to its own radiation and stabilises. This stage is known as the first Larson core (Hoyle 1953; Larson 1969), and the following adiabatic phase prevents any other fragmentation events. However, to achieve total stabilisation, the structure should exert an effective pressure force large enough to support its own mass. Thus, the final mass and size properties of these structures should depend on the growth rate of the instability during collapse as mass can be dynamically added or removed, but also on the degree of instability of the initial structure (Tohline 1980).

Hierarchical fragmentation of clumps within other clumps has been analytically described by Hoyle (1953) using geometric sequences tracing the evolution of the mass of the fragments after series of fragmentation events. Since then, stochastic models of fragmentation were developed for solely to derive the stellar IMF (Larson 1973; Elmegreen & Mathieu 1983) or to study the transfer of the angular momentum throughout many fragmentation events (Zinnecker 1984). Even though it has not been done, these discrete and stochastic approaches are adapted to investigate the multiplicity of the resulting stellar groups as it is possible to simply count the number of stars at the end of the process. However, these models were not designed to account for the physical mechanisms or the thermodynamics that structure the parental cloud from which young stars emerge.

On the other hand, more recent analytical approaches based on the interplay between supersonic turbulence and gravity in the interstellar medium account for the physical processes shaping the continuous multi-scale structure of molecular clouds (Hennebelle & Chabrier 2008; Hopkins 2012; Guszejnov & Hopkins 2015). In these gravo-turbulent models, self-gravitating structures are formed locally within the turbulent cloud through small density enhancements that in turn collapse or eventually subfragment. In this framework, the fragmentation and collapse of a structure is modulated by the balance between (i) gravity that triggers local collapse within a structure, (ii) turbulence that supports the collapse at large scale, but also induce local density enhancements that may collapse, and (iii) thermal pressure. Although these gravo-turbulent models provide a continuous description of a cloud, they do not consider the discrete nature of fragmentation so the multiplicity of stellar groups remains unknown.

Our aim is then to connect the continuous and diffuse structure of the cloud that regulates the fragmentation processes with the result of this fragmentation (i.e. a collection of countable and discrete distribution of stars). Based on a combined stochastic and gravo-turbulent approach, we introduce a framework to model hierarchical fragmentation and monitor the multi-scale structure of a cloud shaped by hierarchical fragmentation. Implementing a set of physical processes at different scales, this model allows us to predict the characteristic mass of last fragmenting pre-stellar structures as well as the characteristic clustering of the pristine stellar groups produced through hierarchical fragmentation.

The paper is organised as follows. In Sect. 2 we establish a general stochastic model that describes the successive hierarchical subfragmentation of local regions within a cloud. In Sect. 3 we implement our model within the framework of the gravo-turbulent fragmentation, and assess the fragmentation properties. Then, in Sect. 4 we constrain our model to existing observational data and evaluate its current limits before concluding this work in Sect. 5.

2 Description of the model

The following stochastic model aims to describe the multi-scale fragmentation of any extended region of a cloud and to provide an analytical framework to investigate the outcome of cloud fragmentation in terms of fragments mass and multiplicity. In this model, a dense structure represents a multi-scale object whose substructures depend on the local physical processes that regulates its fragmentation until star formation.

2.1 A stochastic and geometrical multi-scale model

To introduce our model of fragmentation, we define an extended structure as a subregion within the cloud that undergoes local gravitational collapse and may fragment into several parts with smaller spatial extent, called children. Then these children become the parents of the next generation and subfragment further. This operation is repeated for increasingly smaller scales (see Fig. 1). The main challenge of such a model is to account for the discrete nature of fragmentation, which results in the formation of point-like and countable stars, but also on the continuous and multi-scale nature of the interstellar medium from which these stars form. Therefore, the model used to describe this cascade of fragmentation has to consider the duality between both the discrete product of fragmentation and the continuous origin of its media. Hence, we discretise the spatial dimension into several scales Rl, associated with the fragmentation level l. The level l + 1 defines an additional fragmentation step at lower spatial scales (i.e. R_l > R_l+1), such that hierarchical fragmentation cascades from large fragments to smaller fragments. The design of the following model ressembles the fragmentation model of Larson (1973) that stochastically described the mass hierarchy resulting from fragmentation. In this work we consider the spatial hierarchy resulting from fragmentation. This different approach allows us to model the multi-scale structure of the interstellar medium and implement relevant physical processes at specific spatial scales (see Sect. 3.1).

Due to non-isothermal processes on small scales (e.g. the first Larson core, Larson 1969) the fragmentation cascade results in the formation of children that are eventually unable to subfrag-ment but instead collapse to form individual stars. We define by Rstop the scale beyond which no further fragmentation occurs. At this point, one fragment produces exactly one single young star. We choose to represent this process with a multi-layered network of connected nodes (Fig. 1) in which each level l is populated by fragments whose typical size corresponds to the spatial scale R_l (Fig. 1) while the genealogy between fragments is represented by directed edges from parents to children.

As we consider the geometrical extension of fragments, the maximum number of children of size R_l+1 that can be inserted inside a parent of size R_l is limited by the filling volume of the parent. In fact, a parent of scale R_l cannot contain too many children of size R_l+1 enclosed in its borders, unless these children overlap or merge together. To allow the existence of distinct children between two successive levels, a fragmentation event can only be triggered at a sufficiently high scaling ratio r, where $r=\frac{R_l}{R_{l+1}}>1$. Assuming that the children do not intersect, the maximum number of children n_l,max that can exist within a parent in three dimensions, without prior about the geometrical shape, is the ratio of the volume of the parent at scale R_l to the volume of its children at scale R_l+1: $$n_{l,\max }~r^3.$$(1)

Between two successive levels l and l + 1 , the fragmentation outcome (i.e. the children masses and multiplicity) is regulated both by the number of produced children and by the fraction of the parental mass that each child inherits (Fig. 2). We define the discrete probability density function p_l(n_l) as the probability that at a level l, a parent of mass M_l fragments into n_l children such that $${\displaystyle \sum _{n_l=1}^{n_l\text{ max }}{p}_l}\left(n_l\right)=1,$$(2)

where n_l,max is the maximum number of children possibly generated at the level l + 1 given the geometrical constraints of Eq. (1). The expected number of fragments $\overline{n_l}$ produced at the level l + 1 by one parent at the level l is given by $$\overline{n_l}={\displaystyle \sum _{n_l=1}^{n_{l,\max }}{n}_l}\times p_l\left(n_l\right).$$(3)

The expected number of fragments $\overline{n_l}$ can be computed from any underlying probability distribution p_l(n_l). As we aim to evaluate the general characteristics of fragmentation throughout the spatial scales, $\overline{n_l}$ is the only relevant parameter in the following.

Next, a parent at level l splits its mass M_l between its nl children at level l + 1. We can define the mass efficiency ϵ_l as the fraction of mass inherited by the children from their parent, $${ϵ}_l=\frac{{\displaystyle \sum _{i=1}^{n_l}{M}_{l+1,i}}}{M_l},$$(4)

where the index 1 ⩽ i ⩽ n_l points to the i-th child produced. The mass reservoir ϵ_lM_l is then partitioned between the n_l children with respect to a partition function ψ_l,i describing the fraction of mass the i-th child inherits from the reservoir, such that $${\displaystyle \sum _{i=1}^{n_l}{\psi }_{l,i}}=1.$$(5)

Consequently, the i-th child inherits a mass M_l+1,i = ψ_l,i∊_lM_l and we can verify that Eq. (4) holds (see Fig. 2). In this work we don’t impose additional constraints on the partition function ψ_l,i. Next, we focus on the expected outcome of fragmentation by computing the average number of children produced at any scale R_l (Sect. 2.2) and their average mass (Sect. 2.3). In particular we introduce two parameters that describe the continuous spatial variations in the number of fragments and their mass.

Fig. 1 Schematics of the hierarchical fragmentation process of a cloud with its network representation. Each spatial scale R is associated with a discrete level of fragmentation l. The fragmentation cascades from large to small scales, down to a level l = L corresponding to a scale Rstop beyond which we assume the objects no longer fragment, but collapse into a single star of scale R*.

Fig. 2 Microscopic description of our hierarchical fragmentation model. Left: fragmentation procedure of a single parent contained in a multi-scale structure of size R0. Any parent at a level l, of scale Rl and mass Ml, can be selected in order to evaluate the number nl of fragments it produces and then determine the fraction of its mass ϵlψl,i·distributed to the i-th child. The size of the arrow represents the amount of mass inherited by a child. Right: parent’s probability pl(nl) to produce nl children before transferring a proportion ϵlψl,i of its mass to each of its children.

2.2 Average number of fragments produced

2.2.1 The fragmentation rate ϕ(R)

Considering a cloud of size R₀, at any level l a parent has the probability p_l (n_l) to fragment into n_l children. Hence, all the parents n_par on this level l produce on average a total of 〈N_tot(R_l+1)〉 children, corresponding to $$\langle N_{\text{tot}}\left(R_{l+1}\right)\rangle ={\displaystyle \sum _{n_{\text{par}}}{\overline{n}}_l}.$$(6)

The 〈⋅〉 operator corresponds here to the average amongst all the fragmentation outcomes that may occur within a structure of initial size R₀. Thus, for any scale R there is on average a total amount of 〈N_tot(R)〉 children laying within the structure of size R₀ and 〈N_tot(R)〉 represents the whole collection of children of size R produced by every single parent localised within their hosting structure. We define the fragmentation rate ϕ(R) as the variation of the fragments number per unit of logarithmic size: $$\varphi (R)=-\frac{\text{d\hspace{0.17em}ln}\langle N_{\text{tot}}(R)\rangle }{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}.$$(7)

The minus sign is set such that ϕ(R) > 0 characterises an increase of the number of fragments at lower spatial scale (i.e. an actual fragmentation). With this definition, the fragmentation rate statistic represents the average variation of fragments number inside a multi-scale structure at each size R, rather than the local variation within a specific parent of size R < R₀. In addition, we can show that all the parents of scale R_l ≡ R experiencing a fragmentation rate ϕ(R) produce on average $\langle {\overline{n}}_l\rangle$ fragments at the next scale R_l+1 according to $$\langle {\overline{n}}_l\rangle =\exp \left[-{\displaystyle {\int }_{R_l}^{R_{l+1}}\varphi }\left(R^{\prime }\right)\text{d\hspace{0.17em}ln\hspace{0.17em}}{R}^{\prime }\right].$$(8)

The previous equation connects the statistic associated with the fragmentation rate ϕ(R) with the average multiplicity statistic of individual parents $\langle {\overline{n}}_l\rangle$. The fragmentation rate ϕ(R) is the first parameter of our model and aims to quantify the amount of fragments produced during the fragmentation cascade.

2.2.2 Stellar density

As fragmentation eventually ends at a stoping scale Rstop, we can compute the mean stellar density 〈n,(R₀)〉 as the average number of newly born stars ${\langle N_{*}\rangle }_{R_0}$ contained within a region of any size R₀ per unit of volume, $$\langle n_{*}\left(R_0\right)\rangle =\frac{{\langle N_{*}\rangle }_{R_0}}{V\left(R_0\right)},$$(9)

where V(R₀) is the volume of the region.

Assuming an object of scale R_stop within a region of size R₀ produces a single star, then $\langle N_{\text{tot }}\left(R_{\text{stop }}\right)\rangle ={\langle N_{*}\rangle }_{R_0}$. Then, we can express the average number of newly born stars ${\langle N_{*}\rangle }_{R_0}$ as a function of the fragmentation rate ϕ(R) by integrating Eq. (7), $${\langle N_{*}\rangle }_{R_0}=N_{\text{tot }}\left(R_0\right)\exp \left[{\displaystyle {\int }_{R_{\text{stop }}}^{R_0}\varphi }\left(R^{\prime }\right)\text{d\hspace{0.17em}ln\hspace{0.17em}}{R}^{\prime }\right],$$(10)

where N_tot(R₀) = 1 because we compute the stellar multiplicity of a single structure of size R₀. Here the size R₀ of the region from which we average the number of young stars is to be considered as a variable while R_stop is constant. The average number of newly born stars ${\langle N_{*}\rangle }_{R_0}$ is usefull to predict the multiplicity of stellar groups given their spatial extension (e.g. Sect. 4.1.1) The mean stellar density can also be used to compute the mean separation d between stars in a system as $\overline{d}~{\langle n_{*}\rangle }^{-1/3}$ regardless of their actual spatial distribution.

2.3 Average mass of the fragments

2.3.1 The mass transfer rate ξ(R)

To obtain a second parameter based on the mass of the fragments, we introduce the fragment average formation mass efficiency 〈ℰ(R)〉 at any scale R as $$\langle ℰ(R)\rangle =\frac{\langle M_{\text{tot}}(R)\rangle }{M_0},$$(11)

where M₀ = M(R₀) is the mass of the fragment at the initial scale R₀ and 〈M_tot(R)〉 is the average total mass of the fragments of scale R resulting from every possible fragmentation outcome. Similarly as Eq. (7) we define the mass transfer rate ξ(R) as the log-derivative of 〈ℰ(R)〉 with respect to the size R: $$\xi (R)=-\frac{\text{d\hspace{0.17em}ln}\langle ℰ(R)\rangle }{\text{d\hspace{0.17em}ln}R}=-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}\langle M_{\text{tot}}(R)\rangle }{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}.$$(12)

The mass transfer rate ξ(R) quantifies the variation of the total mass that ends up in the whole collection of fragments of size R. If ξ(R) < 0 then the total mass contained in the fragments of size R is lower than the mass of the fragments of size R + dR. This means that the parents lose a fraction of their mass reservoir during the fragmentation process. The mass dispersal of a fragment into its environment can be done through outflows or radiative feedback (Kim et al. 2018) or any random turbulent motion that could unbind fluid particles from the fragment and partially destroy it (e.g. Heitsch et al. 2006; Camacho et al. 2016; Lu et al. 2020). On the contrary, if ξ > 0 the mass reservoir of the parents is smaller than the children reservoir. Thus, some mass is injected into the children, for example through accretion. Finally, if the mass of the parents is similar to the mass used to form the whole collection of children, we have ξ = 0.

The mass transfer rate ξ(R) only takes into account the net outcome of the interplay between accretion and mass dispersal. Hence, the mass transfer rate does not discriminate the individual contributions of these two effects. Given the definition of the mass transfer rate ξ(R) (Eq. (12)), we can express the average formation efficiency 〈ℰ(R)〉 as $$\langle ℰ(R)\rangle ={ℰ}_0\exp \left[{\displaystyle {\int }_R^{R_0}\xi }\left(R^{\prime }\right)\text{d\hspace{0.17em}ln\hspace{0.17em}}{R}^{\prime }\right],$$(13)

with ℰ₀ = 1 as the total mass of fragments of size R₀ is equivalent to the mass M₀ of the initial structure. At any discrete scale R_l+1, we can write $$\langle ℰ\left(R_{l+1}\right)\rangle =\langle ℰ\left(R_l\right)\rangle \exp \left[{\displaystyle {\int }_{R_{l+1}}^{R_l}\xi }\left(R^{\prime }\right)\text{d\hspace{0.17em}ln\hspace{0.17em}}{R}^{\prime }\right].$$(14)

Substituing 〈ℰ(R_l)〉 and 〈ℰ(R_l+1)〉 by their integral expression Eq. (13), we can express the average efficiency 〈ϵ_l〉 between two discrete levels of fragmentation by $$\langle {ϵ}_l\rangle =\exp \left[{\displaystyle {\int }_{R_{l+1}}^{R_l}\xi }\left(R^{\prime }\right)\text{d\hspace{0.17em}ln\hspace{0.17em}}{R}^{\prime }\right].$$(15)

The mass transfer rate ξ(R) represents the second parameter of our model based on the total mass of the fragments produced in a given region. Its sign is determined by the dominant process between accretion and mass-loss processes (see Table 1).

Next, we consider the mass of a single fragment at scale R. Among all possible fragmentation outcomes, a randomly selected fragment has on average a mass 〈M(R)〉. This average mass is defined considering every object that might be produced within every fragmentation outcome. Hence, the expected mass 〈M(R)〉 one fragment have at scale R considering every fragmentation outcome can be written as $$\langle M(R)\rangle =\frac{\langle M_{\text{tot}}(R)\rangle }{\langle N_{\text{tot}}(R)\rangle },$$(16)

that can be rewritten as $$\langle M(R)\rangle =M_0\frac{\langle ℰ(R)\rangle }{\langle N_{\text{tot }}(R)\rangle }.$$(17)

By taking the log-derivative of the Eq. (17) and injecting the definition of both the fragmentation rate ϕ(R) (Eq. (7)) and the mass transfer rate ξ(R) (Eq. (12)) it can be shown that $$\frac{\text{d }\ln \langle M(R)\rangle }{\text{d }\ln R}=\varphi (R)-\xi (R).$$(18)

This relation implies that the mass efficiency of individual fragment takes into account two components. The first is given by the fragmentation rate ϕ(R) which is linked to the amount of children that have to share the same parental mass reservoir. The second is given by the mass transfer rate ξ(R) that holds information concerning the amount of mass lost or obtained during the fragmentation process.

2.3.2 Volumetric mass density

Given the size R of a fragment and its average mass 〈M(R)〉, one can compute the average mass density 〈ρ(R)〉 of spherical fragments as $$\langle \rho (R)\rangle =\frac{\langle M(R)\rangle }{k{R}^3},$$(19)

where k = 4π/3 is the geometrical constant if the scale R is defined as the radius of the structure while k = π/6 if it represents its diameter. By taking the log-derivation with respect to the size and using Eq. (18) we can show that $$\frac{\text{d }\ln \langle \rho \rangle }{\text{d }\ln R}=\varphi (R)-\xi (R)-3.$$(20)

The dependence of the average mass density with the fragmentation rate ϕ(R) and the mass transfer rate ξ(R) implies that the fragmentation outcome is determined by the physical description of the fragments. In particular, a density fluctuation, induced for example by turbulence, can grow if the condition ϕ(R) − ξ(R) − 3 < 0 is satisfied. In that case, this density enhancement generates a gravitational instability and the structure can locally collapses. If the internal heating and cooling processes of the collapsing structure is modelled with an equation of state (EOS), then the Eq. (20) may provide a mean to study the consequences of the EOS on fragmentation. We investigate in more details in Sect. 3.2.2 the effect of a polytropic EOS on the fragmentation rate.

2.4 Three modes of fragmentation

A structure may fragment hierarchically or collapse monolithi-cally without fragmentation but it may also disperse under the influence of internal support. We evaluate the required conditions a parent needs to (i) fragment into several pieces, (ii) collapse to a single child or (iii) disperse into its environment. According to the definition of the fragmentation rate ϕ(R) (Eq. (7)), a parent associated with ϕ = 0 only forms on average a single child (see Sect. 2.2.1). The outcome of this situation is similar to that associated with a monolithic collapse (i.e. the formation of a single object) (Shu 1977). However, there may be other pairs (ϕ, ξ) that allow single fragment formation.

To qualitatively investigate the conditions describing the different modes under which fragments can be formed, we introduce the free-fall time t_ff that represents the characteristic time with which a gas structure collapses under its own gravity. We consider that a child emerging in a parent whose free-fall time is t_ff,l exists if its own free-fall time t_ff,l+1 is smaller than that of its parent. If the parent collapses faster than its children, the net product of fragmentation is equivalent to a monolithic collapse with virtual children that could merge during the collapse of the parent. Hence, a child exists if $$\frac{t_{\text{ff},l+1}}{t_{\text{ff},l}}<1.$$(21)

For spherical fragments the free-fall time depends on the mass density ρ such that $$t_{\text{ff}}=\sqrt{\frac{3\pi }{32G\rho }}\propto {\rho }^{-1/2}.$$(22)

Then, a child exists within its parent as an actual fragment if the parent density is lower than the density of the possible child: $$\frac{{\rho }_l}{{\rho }_{l+1}}<1.$$(23)

As a consequence, the density of the fragments needs to increase as fragmentation proceeds, which leads to the following condition for the existence of the children according to Eq. (20): $$\xi (R)+3>\varphi (R).$$(24)

By definition hierarchical fragmentation also needs ϕ(R) > 0 such that multiple fragments are formed (Eq. (7)) and we get the condition for hierarchical fragmentation with children of increasing density: $$\xi +3>\varphi >0.$$(25)

Then, we investigate the hierarchical fragmentation process (ϕ > 0) under which children are formed with decreasing density (i.e. when their parent free-fall time is higher). First, we consider $$\varphi ⩾\xi +3⩾0.$$(26)

In this case, the children collapse less rapidly than their parent as they are more diffuse. The children material disperse and the parental structure collapses, so the net product of fragmentation is the formation of a single object. The children formed through fragmentation are called virtual as they do not exist anymore after one parental free-fall time. Since the children are virtual, the outcome of this condition can be modelled using an effective fragmentation rate ϕ_eff = 0. This scenario implies that the parental structure is not necessarily substructured across the scales for all ϕ > 0 but may be observed as collapsing monolithically with dispersed children (red structure in Fig. 3).

Finally, we consider the condition $$\varphi ⩾0>\xi +3,$$(27)

under which hierarchical fragmentation is also prevented and the mass transfer rates satisfies ξ < −3. In this condition, the children also collapse less rapidly than their parent as they are more diffuse. Nonetheless, the mass loss by the collapsing parent is not compensated by its volume reduction so its density also decreases as it contracts (Eq. (20) with ϕ_eff = 0). Therefore, any structure with ξ < −3 eventually disperses into its environment with its virtual children as their density decreases. The net product of fragmentation at smaller scale is then equivalent to a fragmentation disruption (ϕ < 0) as fragments progressively disappear throughout the spatial scales (black structure in Fig. 3).

In a nutshell, we distinguish three fragmentation modes (Fig. 3) in the case of ϕ > 0 with this qualitative study based on the free-fall time:

• a dispersal mode if ξ < −3 as the children and the parent densities decrease with their respective collapse;

• an effective monolithic collapse if ϕ ⩾ ξ + 3 ⩾ 0 as the children disperse into their denser collapsing parent;

• a hierarchical fragmentation mode if ξ(R) + 3 > ϕ(R) > 0 as the children are denser than their parental structure.

In addition, we also recall that fragmentation is disrupted for ϕ < 0 while the monolithic collapse can be represented by ϕ = 0.

It is important to note that the modes of fragmentation we extracted are based on the qualitative study of the free-fall time that does not account for thermal nor turbulent support. As turbulent support is more effective for larger structures (see Sect. 3.1.1), their free-fall time should be higher than the free-fall time of their children. In that case, the hierarchical condition should be less restrictive (i.e. some of the ϕ ⩾ ξ + 3 ⩾ 0 solutions could be hierarchical). Similarily for thermally supported fragments, if the children have a higher internal pressure than their parents (for example because of opacity heating), we may expect the hierarchical condition to be more restrictive and some of the 0 < ϕ < ξ + 3 solutions may not be hierarchical.

Fig. 3 Different fragmentation modes in a fragmentation rate ϕ vs mass injection rate ξ diagram. The hierarchical, monolithic, and mass dispersal modes are shown respectively as blue, red, and black areas and each respective sketch is above the figure. When ϕ − ξ − 3 < 0, the children collapse before their parent (hierarchical). At lower ξ the children do not have time to collapse before their parent does and a single fragment is effectively formed (monolithic). For ξ < −3 this single fragment progressively disperses as it gets more and more diffuse during the collapse (dispersal). In the hatched region ϕ < 0, fragmentation is disrupted.

3 Hierarchical model applied to a gravo-turbulent framework

We introduced in Sect. 2 the general structure of our fragmentation model. Hereafter we implement our model into a physical background using a gravo-turbulent framework (Sect. 3.1) in order to find analytical solution of the fragmentation rate ϕ(R) and determine the conditions that allow a cloud to hierarchically fragment (Sect. 3.2). Then, we investigate the implications of this gravo-turbulent fragmentation on how fragmentation behave with respect to the spatial scales (Sect. 3.3).

3.1 Theoretical framework

3.1.1 Turbulent energy cascade

Turbulent kinetic energy cascades along spatial scales (Kolmogorov 1941). The power spectrum of turbulent energy is scale-free such that the resulting fragmentation process is inherently hierarchical until the turbulent cascade is terminated. In that perspective, we may assume that the turbulent velocity of molecular clouds V_rms could be modelled with any function of the spatial scale so that we can introduce an exponent η(R) characterising the local variations of V_rms(R) as $$\eta (R)=\frac{d\ln V_{\text{rms}}}{d\ln R}.$$(28)

If such a function exists, the turbulent velocity V_rms(R) is scale-free if η(R) is constant, which implies $$V_{\text{rms}}=V_0{\left(\frac{R}{R_0}\right)}^{\eta }.$$(29)

This relation is similar to the Larson law (Larson 1981) for non-thermal motions with V₀ = 1 km s⁻¹, R₀ = 1 pc and η = 0.38, while η = 1/3 corresponds to a Kolmogorov cascade (Kolmogorov 1941). Here the scale R corresponds to the length of the cloud. This scale-free relation has been interpreted as the consequence of a cloud remaining in virial equilibrium but may also arise from a Kolmogorov cascade of turbulent energy (Kolmogorov 1941). Although Larson’s law is empirically valid for sizes R > 0.1 pc (Hennebelle & Falgarone 2012), we assume in the following that Eq. (29) remains valid for sizes R < 0.1 pc. This assumption is reasonable at large scales but may be questionable at smaller scales. This assumption neglects the additional non-thermal support on the fragmentation outcome at smaller scale, for example the magnetic support. In fact, the Eq. (29) may not be verified for non-virialised collapsing structures (Traficante et al. 2018). As pressure support could be underestimated, we expect to overestimate the number of fragments produced with our model (e.g. Sect. 4.2). As long as the function V_rms(R) exists, it could be modelled with any η(R) and could be adapted for higher non-thermal support by changing the value of the power index η. The choice of using the scale-free Larson law should be considered as purely illustrative as we could assess the influence of higher velocities on the fragmentation properties.

3.1.2 ρ-PDF induced by turbulence

In turbulent flow, the density p of a fluid particle contained in a cloud of mean density $\overline{\rho }$ fluctuates stochastically. In a gravo-turbulent fragmentation, the random field of density fluctuations in an isothermal medium has a log-normal statistic (Passot & Vázquez-Semadeni 1998; Padoan & Nordlund 2002). Thus, the density contrast $\delta =\ln (\rho /\overline{\rho })$ follows a normal distribution P(δ): $$P(\delta )=\frac{1}{\sqrt{2\pi }\sigma }\exp \left[-\frac{{\left(\delta +\frac{{\sigma }^2}{2}\right)}^2}{2{\sigma }^2}\right],$$(30)

with a variance σ² = ln (1 + b²ℳ²) where b² ≈ 0.25 and ℳ is the turbulent Mach number of the cloud.

The choice of describing the field of density fluctuations with a normal distribution is a valid solution for an isothermal supersonic medium. Under the assumption that turbulent velocity can be modelled by the Larson law (see Sect. 3.1.1), the former condition is always satisfied for spatial scales above R ~ 0.1 pc at cloud temperature T ~ 10 K. Below R ~ 0.1 pc, adiabatic conditions arise when the fragments are dense enough such that they become opaque to their own radiation. At these small scales the turbulence is subsonic according to the Larson law. In such a subsonic regime the amplitude of density fluctuations is reduced due to the low level of turbulence. According to Passot & Vázquez-Semadeni (1998), we can consistently employ the log-normal ρ-PDF given Eq. (30) under this small density fluctuations approximation.

A collapsing structure develops over time a powerlaw radial density profile ρ(r) ∝ r⁻² where r is the distance to the centre of the structure (Shu 1977). This radial profile may bias the ρ-PDF since the densest regions concentrate within the centre of a structure while the edges are more diffuse. As a consequence, the ρ-PDF develops a power law tail at high densities (Kritsuk et al. 2011). Therefore, this ρ-PDF would no longer be log-normal and one should consider the density profile of a collapsing structure. However, our aim is to evaluate the condition for an instability to develop in the future. If a structure is characterised by this radial profile, the central instability has already grown and no additional fluctuation has appeared in the meantime to disrupt the emergence of a ρ(r) = r⁻² profile. If other fluctuations have appeared in the meantime, subfragments may actually exist. For a more detailed description of this mechanism, we would have to consider the dynamics of the formation of fragments which goes beyond the scope of our quasi-static model.

3.1.3 Collapse threshold

A random density fluctuation may lead to a local collapse within the cloud if the resulting density is sufficient to trigger a gravitational instability. The critical value δ_c beyond which such an instability can grow has been formally expressed by Hennebelle & Chabrier (2008) in the framework of Jeans’ instability considering the turbulent support of a non-magnetised isothermal star-forming collapsing structure $${\delta }_c=\text{\hspace{0.17em}ln\hspace{0.17em}}\left[\frac{{\pi }^{5/3}}{6^{2/3}}\frac{3C_s^2+V_{\text{rms}}^2}{3G{R}^2\overline{\rho }}\right],$$(31)

where G = 6.67 × 10⁻¹¹ m³kg⁻¹s⁻² is the gravitational constant and C_s = (k_BT/µm_H)^1/2 is the sound velocity in which k_B ≈ 1.38 × 10⁻²³ JK⁻¹ the Boltzmann constant, T is the cloud mean temperature, µ = 2.33 the mean molecular weight of the cloud and m_H = 1.67 × 10⁻²⁷ kg the hydrogen mass. Then, a volume V(R) contained in a larger structure collapses if a fluctuation δ ⩾ δ_c arises. The probability 𝒫(δ ⩾ δ_c) for this volume to collapse is then characterised by the complementary cumulative distribution of the PDF P(δ). Considering the ρ-PDF given by Eq. (30), we obtain after integration $$𝒫\left(\delta ⩾{\delta }_c\right)=\frac{1}{2}\left(1-\mathrm{erf}\left[\frac{{\delta }_c+\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}\right]\right),$$(32)

where erf(x) is the error function.

In the following, fragmentation is considered to be the result of the local conditions within the cloud. Thus, we assume that the variance σ(R)² = ln(1 + 0.25ℳ(R)²) is a function of the spatial scale R through the local mach number of the fragment we consider. This local Mach number is defined as the ratio between the turbulent velocity V_rms(R) and the sound velocity C_s. Hence, the complementary cumulative function of Eq. (32) represents the probability that any unstable fluctuation of size R emerges within a structure of density p with a variance σ² determined by its local mean Mach number.

3.1.4 Equation of state

During the collapse of a star-forming structure, the fragments density eventually reaches an opacity limit that induces internal heating and increases the fragment thermal pressure. This heating can be modelled with an EOS regulated by an adiabatic index γ that describes the heating rate of a fragment (Masunaga et al. 1998). Here the index γ represents the heat capacities ratio of a specific chemical specie. Adiabaticity is reached when the fragment density ρ reaches the density ρ_L ~ 10⁻¹³ g/cm⁻³ of the first Larson core (Larson 1969). For molecular hydrogen, the adiabatic index γ = 5/3 for temperatures T < 100 K. The relationship between the temperature T of a structure and its density ρ is described by the following EOS (Lee & Hennebelle 2018): $$T(\rho )=T_0\left[1+{\left(\frac{\rho }{{\rho }_L}\right)}^{\gamma -1}\right],$$(33)

with T₀ ≈ 10 K for a dense clump. We also define the local polytropic EOS at scale R $$\frac{\text{d }\ln T}{ \text{d }\ln \rho }=p-1,$$(34)

where p is the polytropic index describing locally the heating rate of a fragment of density ρ(R) at scale R. Thus, for densities ρ > ρ_L, we have asymptotically T ∝ ρ^γ−1 with a local variation modelled with p = γ. On the contrary, for densities ρ < ρ_L temperature is constant and we have p = 1 (i.e. an isothermal EOS). The sound speed C_s is a function of temperature C_s(T) such that C_s depends on the density of the fragments. Thus, the collapsing criteria δ_c can be expressed as a function of size R and density ρ only.

3.2 An analytical solution for the fragmentation rate ψ(R)

In the following, we assume an isotropic fragmentation, meaning that the probability of finding fragments in specific parts of a parent does not depend on their location within their parent. Since magnetic field usually breaks spherical symmetries, we consider non-magnetised collapsing structures. We also consider the fragmentation as an ergodic process such that the evolution over time of both ϕ(R) and ξ(R) yields to the same variations as the average variation taken over a whole cloud. We also consider that fragmentation only depends on the local properties of the parents such that they do not interact with each other. Finally, we assume that no fragments can be formed in the inter-fragments media at scale R_l+1 without a direct parental structure at scale R_l so fragmentation is entirely hierarchical.

3.2.1 Procedure to compute the fragmentation rate

Our aim is to compute the fragmentation rate ϕ(R) using the ρ-PDF we introduced in Sect. 3.1.2. The maximum number of fragments a structure can contain corresponds roughly to its filling factor V₀ /V(R), where V(R) is the volume of a child with characteristic size R. Considering a fluctuation δ(R), each one of these children has a probability 𝒫(δ ⩾ δ_c) of collapsing. The average number of fragments of size R is then $$\langle N(R)\rangle =N\left(R_0\right){\left(\frac{R}{R_0}\right)}^{-3}\mathcal{P}\left(\delta ⩾{\delta }_c\right),$$(35)

where N(R₀) = 1 as we start from a single cloud. To obtain the fragmentation rate ϕ(R), we use its definition given by Eq. (7) and take the log-derivative of the equation above with respect to the logarithm of the size, so we obtain $$\varphi (R)=3-\frac{\text{d }\ln \mathcal{P}\left(\delta ⩾{\delta }_c\right)}{\text{d }\ln R}.$$(36)

This expression can be computed analytically for any ρ-PDF provided it is differentiable. Since δc depends on the average density $\overline{\rho }$ (Eq. (31)) of the parents, it is necessary to solve simultaneously Eq. (20) which describes the variation of the average density of the fragments. Provided we set the initial conditions on the size R₀ of a structure and its initial mean density ${\overline{\rho }}_0$ to define a starting collapsing threshold δ_c(R₀, ρ₀), we can compute the fragmentation rate by solving the following coupled system: $$\{\begin{array}{l}\varphi (R,\overline{\rho })=3-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}\left(\delta ⩾{\delta }_c(R,\overline{\rho })\right)}{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}\hfill \\ \frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}\overline{\rho }}{\text{d\hspace{0.17em}ln\hspace{0.17em}}R}=\varphi (R,\overline{\rho })-\xi (R)-3\hfill \end{array}.$$(37)

The analytical derivation of the fragmentation rate $\varphi (R,\overline{\rho })$ under all of our considerations (see Sect. 3.1) is given in Appendix A. Expanding Eq. (36) we find that $$\varphi (R,\overline{\rho })=\frac{3}{1-A_{\text{th}}}-(\xi (R)+3)\frac{A_{\text{th}}}{1-A_{\text{th}}}+\frac{2A_{\text{turb}}}{1-A_{\text{th}}},$$(38)

where $A_{\text{th}}(R,\overline{\rho }),A_{\text{turb }}(R,\overline{\rho })$ are associated with the thermodynamics of the collapse (i.e. the EOS) and the turbulent energy cascade (i.e. in this case the scale-free Larson law) respectively, with $$\{\begin{array}{l}{A}_{\text{th}}=A_{\text{p}}\left[1-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}T}{ \text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }\left({\widehat{C}}_s^2+\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }\right)\right]\hfill \\ A_{\text{turb }}=A_{\text{p}}\left[1-\eta \left({\widehat{V}}_{\text{rms }}^2-\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }\right)\right],\hfill \\ A_{\text{p}}=-\frac{\exp \left[-u^2\right]}{\mathcal{P}\left(\delta ⩾{\delta }_c\right)\sigma \sqrt{2\pi }}\hfill \end{array}$$(39)

where $A_{\text{p}}(R,\overline{\rho })$ ponderates both $A_{\text{th}}(R,\overline{\rho })$ and $A_{\text{turb }}(R,\overline{\rho })$ with the fluctuation statistic. The reduced variables ${\widehat{C}}_s^2,{\widehat{V}}_{\text{rms}}^2,{\widehat{ℳ}}^2$ are given by $$\{\begin{array}{l}{\widehat{C}}_s^2=\frac{3C_s^2}{3C_s^2+V_{\text{rms}}^2}\hfill \\ {\widehat{V}}_{\text{rms}}^2=\frac{V_{\text{rms}}^2}{3C_s^2+V_{\text{rms}}^2}\hfill \\ {\widehat{ℳ}}^2=\frac{b^2{ℳ}^2}{1+b^2{ℳ}^2}\hfill \end{array},$$(40)

while the parameters u and υ are reduced centred variables: $$\{\begin{array}{l}u(R,\overline{\rho })=\frac{{\delta }_c(R,\overline{\rho })+\sigma {(R)}^2/2}{\sigma (R)\sqrt{2}}\hfill \\ v(R,\overline{\rho })=\frac{{\delta }_c(R,\overline{\rho })-\sigma {(R)}^2/2}{\sigma (R)\sqrt{2}}\hfill \end{array}.$$(41)

The fragmentation rate expressed in Eq. (38) consists of the addition of three terms. From left to right, the first one reflects both the isotropy of fragmentation and the spheroidal geometry of the fragments we consider. The second one characterises the dependancy of ϕ(R) with the mass transfer rate ξ(R) since the density of fragments inherently couples them. The third term is associated with the turbulent energy cascade that contributes only in the supersonic limit as this term cancels out in a subsonic flow. All these terms are weighted by both the density fluctuation statistic and the local EOS. In particular, as the fragmentation rate sensitively depends on the EOS, heating processes are a dominant and determining factors in the fragmentation outcomes as we show in Sect. 3.3.4.

In Eq. (38) the mass transfer rate ξ(R) can be represented by any function. Without any suitable model to predict ξ(R), we assume that ξ is a constant so that we can test in the following different fragmentation outcomes as a function of the mass transfer rates. We have for the mean efficiency $$\langle ℰ(R)\rangle ={ℰ}_0{\left(\frac{R}{R_0}\right)}^{-\xi },$$(42)

with ℰ₀ = 1 as the total mass of fragments of size R₀ is equivalent to the mass M0 of the initial structure.

3.2.2 Condition for a hierarchical fragmentation

In Sect. 2.4 we established qualitatively different fragmentation modes depending on the numerical values of both the fragmentation rate ϕ and mass transfer rate ξ (see Fig. 3). We can compute analytically using Eq. (38) the solutions that determine these modes under subsonic and supersonic conditions. A parent hierarchically fragments if the condition ξ + 3 > ϕ(R) > 0 is satisfied. In the subsonic limit in which C_s ≫ V_rms, we can show (see Appendix B) that $$\varphi (R,\overline{\rho })\approx \xi +3-\frac{2\left(1+\eta \frac{{\delta }_c}{2}\right)}{1-\frac{\text{d\hspace{0.17em}ln\hspace{0.17em}}T}{ \text{d\hspace{0.17em}ln\hspace{0.17em}}\rho }\left(1+\frac{{\delta }_c}{2}\right)}.$$(43)

To simplifiy Eq. (43) we assume that the mass reservoir is conserved between the scales so that ξ = 0. We also assume marginally unstable parents such that ${\rho }_c~\overline{\rho }$ and δ_c ~ 0. Then, using the local EOS (see Eq. (34)), Eq. (43) can be simplified to $$\varphi \approx \frac{4-3p}{2-p}.$$(44)

In these conditions, hierarchical fragmentation occurs when the condition 3 > ϕ(R) > 0 is satisfied, which means when the polytropic index p satisfies $$p<\frac{4}{3}.$$(45)

In order to collapse monolithically we can also show that the condition ϕ = 0 implies (under the same approximations) that p = 4/3. Beyond p = 4/3, the children produced disperse their mass. This result is similar to the traditional analysis of Jeans fragmentation that derives a critical Jeans mass $M_J\propto {\rho }^{\frac{3}{2}\left(p-\frac{4}{3}\right)}$. Hence, p = 4/3 is a critical point below which the Jeans instabilities grow as the fragments contract quasi-statically. The Eq. (43) generalises the Jeans analysis in the subsonic limit by taking into account the amount of mass added or removed from the fragments (ξ), the stability of the collapsing structure (δ_c) but also the decay of the turbulent cascade (η).

In the supersonic limit (C_s ≪ V_rms), ϕ → 3 for any ξ. Hence, any loss of mass ξ < 0 during the fragmentation may likely induce a monolithic collapse as the pair (ϕ, ξ) would satisfy the monolithic condition ϕ ⩾ ξ + 3 ⩾ 0 (see Fig. 3).

3.3 Variations of the fragmentation rate

3.3.1 Turbulent vs. thermal support

To assess the variations of the fragmentation rate ϕ(R) with the size R of the fragments we consider in the following a cloud of size R₀ = 10 pc and mass M₀ = 10⁴ M_⊙, so the initial mean density $\overline{{\rho }_0}~{10}^{-21}$ g/cm³ with a cloud radius R₀/2. The fragmentation rate is computed using an isothermal EOS with T = 10 K and for different mass transfer rates ξ = constants varying between 1 and −1 (Fig. 4).

We distinguish two regimes of fragmentation. The first one is a fragmentation driven by turbulence at large scales while the second one is driven by gravity at smaller scales. The transition between these two regimes occurs for spatial scale around R ~ 0.1 pc that is generally designated as the characteristic sonic length (Vázquez-Semadeni et al. 2003; Ballesteros-Paredes 2004; Palau et al. 2015) from which turbulence becomes negligible compared to thermal support.

In the framework of our model, this transition can be understood as the moment when the collapse threshold δ_c reaches a maximum (Fig. 5). Since δ_c < 0 at R₀, the cloud is initially unstable (see Fig. 5). The threshold δ_c eventually becomes positive while the turbulence is still supersonic. If R > 0.1 pc then δ_c increases for decreasing R and it is more difficult for a fluctuation to be gravitationnally unstable because of strong turbulence support. Once R < 0.1 pc, turbulent support has sufficiently decayed not to be dominant anymore compared to thermal support. Hence, the threshold δ_c decreases as R decreases since a random fluctuation is more and more likely to induce a local collapse. The maximum of the function δ_c(R) for which this transition occurs can be analytically computed by taking the derivative with respect to the variable ln(R). Then we can show that this derivative cancels out when $$\frac{2\eta {ℳ}^2}{3+{ℳ}^2}=\varphi -\xi -1.$$(46)

For the values of ξ we considered, this transition occurs when ϕ − ξ ≈ 1.1−1.4. Using the analytical solution given Eq. (46) for η = 0.38, we find that super-to-subsonic transition occurs when ℳ ≈ 0.7–1.8 such that the Mach number needs to be of the order of unity.

In the supersonic limit, the fragmentation rate asymptotically tends to ϕ ≈ 3 (see Sect. 3.2.2) while in the subsonic limit the fragmentation rate is defined by Eq. (43). In addition, the fragmentation is scale-free if ϕ is constant, which is the case in the supersonic limit, as we may expect from a scale-free turbulence cascade. This transition between supersonic and subsonic fragmentation necessarily breaks the turbulent scale-free fragmentation. In the subsonic case, the fragmentation is scale-free if the structure is marginally unstable (Eq. (44)) and if ξ is constant. In fact, the fragmentation rate varies slowly in the subsonic regime for R < 0.01 pc (Fig. 4) and can be considered as locally scale-free.

Nonetheless, the scale-free subsonic fragmentation cannot be reached for all mass transfer rate. If ξ < −0.5 fragmentation falls into a disrupted regime as ϕ(R) < 0 (Fig. 4). In that case, fragmentation stops at scale R_stop > 10⁻² pc. Using Eq. (42) we find that fragmentation is prevented if the last fragments are produced with amass efficiency 〈ℰ(R_stop)〉 < 3%. Hence, we should expect large-scale mass efficiencies to be higher than the 3% order of magnitude in order to produce star-forming structures.

Fig. 4 Fragmentation rate ϕ as a function of spatial scale R for isothermal EOS and different values of mass transfer rate ξ with cloud initial size R0 = 10 pc and initial density ρ0 = 10−21 g/cm3. The hatched area represents the ϕ ⩽ 0 space in which there is no hierarchical fragmentation. The grey vertical line separates turbulent (supersonic) fragmentation from thermal (subsonic) fragmentation.

Fig. 5 Contrast density threshold ${\delta }_c=\text{\hspace{0.17em}ln\hspace{0.17em}}\left({\rho }_c/\overline{\rho }\right)$ as a function of spatial scale R for isothermal EOS and different values of mass transfer rate ξ with cloud initial size R0 = 10 pc and initial density ρ0 = 10−21 g/cm3. The grey vertical line separates turbulent (supersonic) fragmentation from thermal (subsonic) fragmentation.

3.3.2 Self-regulation of the fragmentation

We may check the influence of the cloud’s initial conditions on the evolution of the fragmentation rate ϕ(R). Since the important physical property of a structure for the fragmentation is its mean density $\overline{\rho }$, we choose to vary the initial mass of the cloud. In addition, we assume an isothermal EOS, a cloud initial size R₀ = 10 pc and a mass transfer rate ξ = 0. The fragmentation rate ϕ(R) is computed for different cloud initial masses M₀ = 10³−10⁴−10⁵−10⁶ M_⊙, corresponding respectively to mean densities $\overline{{\rho }_0}\approx {10}^{-22}-{10}^{-21}-{10}^{-20}-{10}^{-19}$ g/cm⁻³ with a cloud radius of R_0/2. The relative variations of the fragmentation rate are all similar even though more diffuse clouds have smaller fragmentation rates (Fig. 6a).

In the subsonic regime associated with a thermal support, we notice that the fragment average masses and average densities are similar no matter the initial density (respectively Figs. 6b and c). Indeed, as the spatial scale decreases the mass and density slowly converge towards the same numerical values. The asymptotic trend is analytically understandable within the framework of our model. On small scales, there is a convergence towards marginally stable structures as δ_c → 0, which induces convergence of ϕ − ξ → 1 (Eq. (43)). On the other hand, the difference ϕ − ξ corresponds to the local variation of average mass (Eq. (18)). Hence, if ϕ − ξ is constant, then $$\langle M\rangle =M_0{\left(\frac{R}{R_0}\right)}^{\varphi -\xi }.$$(47)

With ϕ − ξ = 1 we have 〈M〉 ∝ R. This scaling relation is typical of the critical unstable Bonnor-Ebert isothermal gas spheres embedded in a pressurised medium (Bonnor 1956) as $$M_{BE}=2.4R_{BE}\frac{C_s^2}{G},$$(48)

where R_BE is the critical Bonnor–Ebert radius of a spheroidal isothermal gas clump. In addition, under the same conditions, the average density of the fragments can be written as $$\langle \rho \rangle ={\rho }_0{\left(\frac{R}{R_0}\right)}^{\varphi -\xi -3}.$$(49)

With ϕ − ξ = 1 we have 〈ρ〉 ∝ R⁻². This scaling relation is what to be expected as a result from collapsing self-gravitating structures (Shu 1977). Both of these relations suggest that fragmentation is a mechanism tending to naturally produce marginally stable structures that collapse quasi-statically in order to form stars, regardless of the cloud initial stability (i.e. its initial density). Even though these structures have not actually reached a stage of marginal stability with δ_c ≈ 0, these asymptotic behaviours are reasonable approximations on a small scale R < 0.01 pc (Figs. 6b and c).

Nevertheless, the similar convergence on both the mass and density numerical values suggests that fragmentation is also a self-regulated process. The more massive a parent is, the higher its fragmentation rate. As it fragments more, the average mass of the children is reduced accordingly to the value of ϕ − ξ, as they need to share their parent mass. A less massive parent tends to fragment less, meaning its fragmentation rate is lower than that of a more massive parent. Thus, the average mass of each children decreases less rapidly. The two parents eventually form children with the same mass on average. As a consequence, these children fragments identically as they have the same fragmentation rate.

Fig. 6 Variations of the fragmentation rate ϕ (a), mean fragment mass 〈M〉 ∝ Rϕ−ξ (b), and mean fragment density 〈ρ〉 ∝ Rϕ−ξ−3 (c) as a function of their spatial scale R for different initial cloud densities ρ0 with an isothermal EOS using ξ = 0. The legend of (a) is common to the three plots. The scaling relation M ∝ R in (b) indicates the asymptotic limits of a marginally unstable Bonnor-Ebert structure (Bonnor 1956). The scaling relation ρ ∝ R−2 in (c) represents the density profile of a self-gravitating isothermal sphere (Shu 1977).

3.3.3 The structure of the cloud

In the supersonic regime associated with turbulent fragmentation, the average mass 〈M〉 evolves like a power law as a function of the spatial scale R (Fig. 6b), reflecting the fractal structure of the cloud above 0.1 pc as a consequence of turbulence (Elmegreen 1997; Hennebelle & Falgarone 2012; Kritsuk et al. 2013). However, the fragmentation rate, which traces the local slope as 〈M〉 ∝ R^ϕ considering ξ = 0, varies between ϕ(R) = 2 − 3 for R > 0.1 pc (Fig. 6a). Hence, the apparent fractal-ity of the cloud seen from the mass-size relation at large scales is an artifact of the low dynamical spatial range at which this fractality is considered. Nonetheless, this variation range is consistent with mass-size relation fractal index derived for example by Larson (1981) and Myers (1983) (M ∝ R²), Traficante et al. (2018) (M ∝ R^2.38) and Elmegreen & Falgarone (1996) (between M ∝ R^2.3−3.7). The variability of these slopes may result either from a difference in the mass transfer rates (ξ) within the clouds and/or as a difference in the cloud initial density $\overline{\rho }$ (Fig. 6).

3.3.4 Characteristic scale as a stop for fragmentation

To model the increasing opacity of the fragments as their density increases, we use the polytropic EOS given by Eq. (33) with an adiabatic index γ = 5/3. We consider a cloud of size R₀ = 10 pc and mass M₀ = 10⁴ M_⊙, so the initial mean density $\overline{{\rho }_0}~{10}^{-21}$ g/cm³. The fragmentation rate is computed for different mass transfer rates ξ = constants varying between 1 and −1 (Fig. 7a).

For R > 10⁻² pc, the fragmentation rate follows isothermal behaviour (Fig. 4a) as the cloud is still under isothermal conditions with T = 10 K. However, unlike the case of pure isothermal EOS, ϕ(R) eventually becomes negative and goes into a disrupted regime as thermal pressure overcomes gravity. The spatial scales R_stop at which ϕ(R) = 0 so fragmentation ends are defined as the stopping scales. For each mass transfer rate we can attribute one different stopping scale that is typically at R_stop ~ 40–100 AU for ξ ⩾ −0.5 and R_stop ~ 1 kAU for ξ = −1. These characteristic scales of a few tens of AU represent the scales below which the fragmentation of a gaseous clump is prevented and set the initial conditions for the formation of the first Larson core (Lee & Hennebelle 2018). However, at these scales the presence of a proto-planetary disk is expected. Such a disk structure is fundamentally different from a spherical gas structure we have considered and have a different fragmentation process. The consequences of disk fragmentation are discussed more specifically in Sect. 4.3.

The density of the last fragmenting structures formed of size Rstop is represented in Fig. 7b while their associated temperature is represented in Fig. 7c. For ξ ⩾ 0.5 the last fragmenting structures are adiabatic, setting the initial conditions for the formation of the first Larson core in these structures. The last fragments associated with −0.5 < ξ < 0.5 have started to heat up without reaching the adiabatic regime, that is the initial conditions for the first Larson core. It is necessary for these structures to get denser in order to host stars. Hence, they need to accrete from their environment and increase their mass transfer rate, otherwise any pre-stellar core of size R < R_stop is expected to disperse. Then, if ξ < −0.5 the fragments are still isothermal and they are unlikely to form the first Larson core.

Fig. 7 Theoretical variations in the fragmentation rate φ (a), the mean fragment density 〈ρ〉 (b) as a function of the spatial scale R for different mass transfer rates ξ. These profiles were computed using an adiabatic EOS (c), a cloud initial size R0 = 10 pc with an initial density ρ0 = 10−21 g/cm3. Red triangle, black dot, blue triangle, and blue square indicate the position of the last fragmenting structures for mass transfer rates ξ = 1,0.5,0, −0.5, −1, respectively. (a) The black dotted line is associated with an isothermal EOS with a polytropic index p = 1. The transition between turbulent and thermal fragmentation is delimited by the grey line. (b) The black dashed line indicates the density of the first Larson core (Larson 1969) when p = 5/3. (c) Evolution of temperature T as a function of fragments density.

4 Discussion

4.1 Observational constraints

4.1.1 Number of YSOs in stellar groups

Using the fragmentation rate solutions we computed within a gravo-turbulent framework in the case of an adabatic EOS (Sect. 3.3.4), we can evaluate the average number of stars 〈N_*〉_R located within spherical stellar group of size R using Eq. (10). As the fragmentation rate ϕ(R) cancels out for sizes R_stop, we assume below this point a one-to-one correspondence between the prestellar cores and young stars such that at R_stop the number of stars ${\langle N_{*}\rangle }_{R_{\text{stop }}}=1$.

The average number of stars 〈N_*〉_R can be directly compared to the NESTs as defined by Joncour et al. (2018) as spatially extended stellar structures, that are local concentrations of stars within the cloud. We show (Fig. 8) that the spatial extension of NESTs of size R > 20 kAU is compatible with our model considering mass transfer rates ξ ~−0.75 for an initial cloud density ρ₀ = 10⁻²⁰ g/cm⁻³. In a star formation efficiency perspective, this means that a 0.1 pc clump produces the last fragments at scale R_stop ≈ 800 AU with a 〈ℰ(R_stop)〉 ~ 9% efficiency. For low clouds initial density ρ₀ > 10⁻²³ g/cm⁻³, NESTs of size R > 20 kAU remain compatible with our model for ξ < −0.25, for average efficiencies 〈ℰ(R_stop)〉 < 20% at stopping scales R_stop ≈ 40 AU.

In addition, the number of stars within NESTs saturates around ${\langle N_{*}\rangle }_R^{\text{NESTs }}~5$, between R = 10 kAU and R = 20 kAU. We can interprete this plateau with our model as the spatial properties of the smaller NESTs are compatible with higher mass transfer rates, suggesting that ξ may not be a constant across the spatial scales but instead increases as the size of the fragments decreases. The plateau at ${\langle N_{*}\rangle }_R^{\text{NESTs }}~5$ implies that ϕ = 0 at least within the R = 5–20 kAU range of scales. Hence, the collapse should be monolithic at these scales. However, given the multiplicity of these NESTs, we expect that fragmentation resumes on a spatial scales R < 5 kAU. The actual multiplicity of the fragments that could have generated these NESTs over this spatial range appears compatible with our model at scales R > 20 kAU considering constant mass transfer rate while the mass transfer rate should be increasing for R < 20 kAU suggesting a better formation efficiency or a more efficient accretion of mass at smaller scales, or both. Although, the NESTs multiplicity is well reproduced at all scales for diffuse cloud (ρ₀ = 10⁻²³ g/cm³) considering a constant mass transfer rate ξ = −0.25.

Fig. 8 Average number of stars 〈N*〉R located within a region of size R for different initial cloud density ρ0 as predicted by our model of fragmentation. The solid, dashed, dot-dashed, and dotted curve represent our model for different mass transfer rates ξ = 0, −0.25, −0.5, −0.75, respectively. The triangles mark the points below which gas clumps do not fragment anymore. The crosses represent the number of YSOs in NESTs (Joncour et al. 2018) as a function of their size, as observed in the Taurus cloud.

4.1.2 Observational constraints on the fragmentation rate

This model of fragmentation can be constrained using multi-scale analysis of molecular clouds. If one is able to count the number of fragments at different spatial scales, the fragmentation rate may be measured. Such a counting was performed using a combination of Spitzer and Herschel 70 - 160 - 250 - 350 and 500 µm images of NGC 2264 star-forming region (Thomasson et al. 2022). In this work, we define a fractality coefficient F as the average number of fragments a single parent produces between two scales separated by a ratio r = R_l/R_l+1 = 2. Using the definition of ϕ(R) given by Eq. (7) we can define a two dimensional fragmentation rate ϕ_2D describing the variation of the average number of fragments 〈N_2D〉 observed in the two dimensional projection of a molecular cloud. Because the data sample the spatial scales into discrete levels of fragmentation (e.g. Fig. 1), we cannot compute the ‘instantaneous’ fragmentation rate at scale R but instead the average fragmentation rate $\overline{{\varphi }_{2D}}$ integrated over the range of scales R that separate the largest fragment to the smallest one $$\overline{{\varphi }_{2D}}=-\frac{\text{Δ\hspace{0.17em}ln}\langle N_{2D}\rangle }{\text{Δ\hspace{0.17em}ln\hspace{0.17em}}R}.$$(50)

By definition of the fractality coefficient, there is a production of ∆ ln〈N_2D〉 = ln ℱ fragments for every reduction of scales ∆ ln R = −ln 2, so $$\overline{{\varphi }_{2D}}=\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}ℱ}{\text{\hspace{0.17em}ln\hspace{0.17em}}2}.$$(51)

Subsequently we can compute the error on the fragmentation rate associated with the measure of ℱ as $$\overline{{\varphi }_{2D}}=\frac{\text{\hspace{0.17em}ln\hspace{0.17em}}ℱ}{\text{\hspace{0.17em}ln\hspace{0.17em}}2}\pm \frac{{\sigma }_{ℱ}}{ℱ\text{\hspace{0.17em}ln\hspace{0.17em}}2}.$$(52)

With a measurement ℱ = 1.45 ± 0.12 (Thomasson et al. 2022), we obtain $\overline{{\varphi }_{2D}}=0.54\pm 0.12$. To compute the equivalent three dimensional fragmentation rate $\overline{{\varphi }_D}$ we need to account for projection effects that may coincidentally align two structures on the line of sight and reduce the number of fragments we can detect. To obtain this three dimensional fragmentation rate, ellipsoidal geometries are randomly sampled and fragmented accordingly to a fragmentation rate $\overline{{\varphi }_{3D}}$. The resulting ellipsoidal fragments are randomly placed inside their parents and the fragmentation process is repeated for at smaller scales. Between two levels of fragmentation, we used a scaling ratio r = 1.5, similar to the scaling ratios defined in Thomasson et al. (2022) (more details in Appendix C). We then project the obtained fragments into a 2D plane and filter the aligned ones. We finally count the amount of fragments remaining after the projection in order to derive the fragmentation rate $\overline{{\varphi }_{2D}}$. The error induced by the projection is given in Fig. 9 in which we estimate a variation $f=\frac{\overline{{\varphi }_{2D}}-\overline{{\varphi }_{3D}}}{\overline{{\varphi }_{3D}}}\approx -0.30\pm 0.10$. From the relation $$\overline{{\varphi }_{3D}}=\frac{\overline{{\varphi }_{2D}}}{1+f},$$(53)

we can compute the uncertainty ${\sigma }_{\overline{{\varphi }_{3D}}}$ with a quadratic error propagation $${\sigma }_{\overline{{\varphi }_{3D}}}=\sqrt{{\left(\frac{{\sigma }_{\overline{{\varphi }_{2D}}}}{1+f}\right)}^2+{\left(\frac{\overline{{\varphi }_{2D}}}{{(1+f)}^2}{\sigma }_f\right)}^2}.$$(54)

With $\overline{{\varphi }_{2D}}=0.54,{\sigma }_{\overline{{\varphi }_{2D}}}=0.12$ f = −0.3 and σ_f = 0.1, we obtain the three dimensional average fragmentation rate $$\overline{{\varphi }_{3D}}=0.77\pm \mathrm{0.20.}$$(55)

As this fragmentation rate has been computed within the 1.4–26 kAU range, we can compute the mean fragmentation rate we obtain in the gravo-turbulent model (Sect. 3.2) for ξ ∊ [−1; 1] with $${\overline{\varphi }}_{3D_{\text{model }}}=\frac{1}{\text{\hspace{0.17em}ln\hspace{0.17em}}{R}_b-\text{\hspace{0.17em}ln\hspace{0.17em}}{R}_a}{\displaystyle {\int }_{\text{\hspace{0.17em}ln\hspace{0.17em}}{R}_b}^{\text{\hspace{0.17em}ln\hspace{0.17em}}{R}_a}\varphi }(R)\text{d\hspace{0.17em}ln\hspace{0.17em}}R,$$(56)

where R_a = 1.4 kAU and R_b = 26 kAU. We numerically computed the Eq. (56) with the gravo-turbulent model with a R₀ = 10 pc cloud of initial mass M₀ ≈ 10⁴ M_⊙ to simulate NGC 2264 (Dahm 2008) corresponding to initial cloud density ρ₀ ≈ 10⁻²¹ g/cm³. We determine that the measured fragmentation rate is compatible with our theoretical value for mass injection rates ξ ∊ [−0.49; −0.15] (see Fig. 10), meaning that between 7 and 40% of the mass of a 0.1 pc gas structure is used to form its stars. In NGC 2264 hierarchical fragmentation is characterised by a loss of fragment mass with ξ < 0, for example via molecular jets as it has been observed in the fragmented hubs of the cloud (Maury et al. 2009; Cunningham et al. 2016). Since we define no proper physical model to describe ξ we are not able to quantify the mass ejected per unit time for better comparison. Nonetheless, the range of scales (R < 0.1 pc) and the numerical value of the fragmentation rate (ϕ ~ 0.77) suggests that the subfragmentation of NGC 2264 corresponds to an isothermal fragmentation (Fig. 7), dominated by gravity with a marginal contribution from turbulence.

Furthermore, the loss of fragment mass with ξ < 0 evaluated in NGC 2264 is not incompatible with the NESTs properties from Taurus (see Sect. 4.1.1) despite the fact that NGC 2264 is a dense star-forming region while Taurus is a more diffuse cloud. We may hypothesise that ξ < 0 at scales R > 10 kAU is essentially the rule within molecular clouds regardless of their densities. As we show in Sect. 3.3.1, these scales characterised by a negative mass transfer rate are also characterised by turbulent fragmentation. Hence, the mode of fragmentation (here turbulent) could set the stellar multiplicity rather than the cloud initial conditions.

Fig. 9 Relative error caused on the fragmentation rate by the projection of 3D ellipsoidal fragments onto a 2D plane. The 2D fragmentation rate $\overline{{\varphi }_{2D}}=0.54\pm 0.12$ (red area) is associated with a relative error ~−30% ± 10 with the 3D fragmentation rate $\overline{{\varphi }_{3D}}$. The dashed grey line indicates a 0% error.

Fig. 10 Theoretical average fragmentation rate $\overline{{\varphi }_{3D_{model}}}$ integrated between Ra = 1.4 kAU and Rb = 26 kAU as a function of the injected mass rate ξ. The red horizontal line and grey uncertainty band corresponds to the fragmentation rate $\overline{{\varphi }_{3D}}=0.77\pm 0.20$ measured in NGC 2264. The arrows indicate the ξ boundaries between which our model is compatible with the measurement for different cloud initial densities ρ0.

4.2 Influence of a magnetic field

Throughout our model we have considered a non-magnetised cloud. To qualitatively assess the local influence of a magnetic field, we consider that the presence of a magnetic field brings an additional support to counteract gravity, and therefore prevents fragmentation. To evaluate the variation induced by the magnetic field, or any mechanism that tends to stabilise a structure through additional support, we may add a constant velocity component V_add to Eq. (31), so $${\delta }_c=\text{\hspace{0.17em}ln\hspace{0.17em}}\left[\frac{{\pi }^{5/3}}{6^{2/3}}\frac{3C_s^2+V_{\text{rms}}^2+V_{\text{add }}^2}{3G{R}^2\overline{\rho }}\right].$$(57)

To qualitatively estimate the additional support arising from the magnetic field, we introduce the Alfven velocity $$v_A=\sqrt{\frac{B^2}{{\mu }_0\rho }},$$(58)

where B is the magnetic field strength, µ₀ = 1.26 × 10⁻⁶ T m/A is the vacuum magnetic permittivity and ρ ~ 4 × 10⁻¹⁴ g/cm⁻³ the fragment density at R ~ 40–100 AU. We computed the fragmentation rate ϕ(R) for V_add < 500 m/s (Fig. 11) which corresponds to an equivalent magnetic fields B < 30 mG in order of magnitude (see e.g. Biersteker et al. 2019; Vlemmings et al. 2019), assuming υ_A ~ V_add. This additional support shifts the R_stop from 40 AU towards larger scales, up to 100 AU as the fragmentation rate decreases. Nevertheless, the magnetic field in proto-planetary disks of tens of AU is typically about B ~ 1 mG (Brauer et al. 2017; Lesur et al. 2023). Considering this magnetic strength, the additional support is negligible as υ_A ~ 3 m/s. We suppose here that the support Vadd is constant over the whole R = 40 AU − 10 pc range even though at scales R > 0.1 pc the magnetic field is negligible (B ~ 10–30µG, Crutcher et al. 2010; Pattle et al. 2023). Therefore, the variation of R_stop caused by the magnetic-field should be taken as a maximum variation.

In addition to this support, accounting for the magnetic field in our model would break the isotropic assumption. As a consequence, the effective parental volume that can be filled by the children would be reduced. This is reflected in the volume ratio V(R)/V₀ = (R/R0)⁻³ in Eq. (35). Under the influence of a magnetic field, fragmentation tends to be inhibited in the perpendicular direction of the magnetic field lines, reshaping the geometries of the structures. Because of this anisotropy, the size of fragments may vary more quickly in one direction compared to another (e.g. an ellipsoid that fragments along its length and the fragments have the same width of their parent). As a result, the power index −3 should decrease. As fragmentation is geometrically constrained by the effective parental volume, ϕ should decrease. Hence, even without introducing an additional support, the fragmentation rate can be reduced by any anisotropy.

Fig. 11 Influence of additional support Vadd on both the fragmentation stopping scale Rstop (blue) and the density of the last fragmenting structures < ρ(Rstop) > at this stopping scale (red). The three vertical greydashed lines indicates the estimated equivalent magnetic field strengths from their Alfven velocities at ρ ~ 4 × 10−14 g/cm−3 using Eq. (58).

4.3 Disk fragmentation

Below the stopping scale R_stop ~ 40–100 AU, we expect the first Larson core to emerge and later the formation of a protoplanetary disk. Hence, the collapse criteria of a gaseous sphere δ_c we used to compute the fragmentation rate (Eq. (31)) should not be valid anymore. Instead, one should consider the instability condition of a rotating disk that may be magnetised. The instability of such a disk may generate additional fragmentation events that result in the formation of binaries. Hence, to obtain the actual multiplicity of the close stellar systems and the mass of their stars it would be necessary to account for these fragmentation events.

At first approximation, we may assume that the prestellar core turns into a single star with a one-to-one correspondence. We can define the efficiency between the core and the star ϵ_CS representing the fraction of the prestellar core mass that effectively end up in the star it hosts. However, if multiple stars are formed because of disk fragmentation, for example a binary, this efficiency needs to be corrected by the mass ratio q of the binary (see Sect. 2.1). This mass partition would reduce the individual mass of the stars compared to a similar core that would only produce a single star. Hence, disk fragmentation would tend to produce less massive stars.

Furthermore, the systematic formation of a binary would systematically multiply the number of stars produced within gaseous fragments by a factor of 2. Thus, the spatial and multiplicity properties of the NESTs would become compatible for mass transfer rates ξ < −0.75 for scales R > 0.1 pc suggesting that supersonic turbulence produces fragments with low mass efficiencies. As an order of magnitude, a R₀ = 10 pc cloud produces R = 0.1 pc fragments with an efficiency <β> = (0.1/10)^−ξ ~ 3% with ξ = −0.75. Under such circumstances the fragments do not reach their adiabatic phase as their density is not high enough, assuming ξ = −0.75. This suggests that the mass injection rate ξ increases for R < 0.1 pc in order for the fragments to reach their adiabatic phase and effectively form a star.

4.4 The mass of the stars

In Sect. 3.3.2 we have shown that fragmentation is self-regulated at small scales, suggesting that the final properties of stars are marginally determined by the initial conditions of the cloud. Provided that this self-regulation is sufficiently fast to occur, the young stars or even pre-stellar cores should possess a characteristic average mass independant of their cloud. Under adiabatic conditions, hierarchical fragmentation leads to fragments of size R_stop with characteristic mass M ≈ 4 × 10⁻³ M_⊙. These low mass structures are consistent with fragments built under adi-abatic conditions when they reach their opacity limit (Boyd & Whitworth 2005; Lee & Hennebelle 2018). Nonetheless, the mass of the first Larson core that supposedly forms inside these fragments is about 10⁻² M_⊙ (Larson 1969; Lee & Hennebelle 2018), while the characteristic stellar mass is 0.1 M_⊙, corresponding to the ‘peak’ of the IMF. This characteristic stellar mass can emerge within the framework of a single step fragmentation process (i.e. without hierarchy) in which filaments can produce fragments of mass 0.1–10 M_⊙ (Larson 1985). Sim-ilarily, the thermal Jeans mass of ~0.1 pc clumps is about 0.1–1 M_⊙ assuming no hierarchical substructures (Larson 2005; Alves et al. 2007). Hence, the introduction of a hierarchical process inherently build less massive fragments.

Using our model and the observational constraints we derived, we propose a scenario that may reconcile the mass discrepancies between the stars and the pre-stellar structures. For scales R > 0.1 pc, we extract a mass transfer rate ξ < 0 (Sect. 4.1). At these scales, turbulence is in a supersonic regime and the fragments are not completely bound because of turbulent dissipation. As a result, they lose a fraction of their mass during their formation so ξ < 0. On scales of R < 0.1 pc, turbulence enters a subsonic regime. The fragments at these scales are more bound than at larger scales so they conserve their mass better. This mechanism would result in a progressive increase of the mass transfer rate ξ(R) as R decreases and as the turbulent support vanishes, while these fragments finally accrete their environment to subsequently form stars.

Alternatively, hierarchical fragmentation may have transient phases in which fragments collapse monolithically in a range of spatial scales before subfragmenting hierarchically again. This solution would produce more massive fragments as fewer fragments share the mass reservoir. However, this solution is speculative and the physical mechanism that would slow down the hierarchical process remains unkown.

5 Conclusion

5.1 Summary

We introduced a probabilistic multi-scale fragmentation model to describe the multi-scale structure of molecular clouds down to star formation. With this model we aim to predictthe stellar multiplicity of the systems formed through hierarchical fragmentation as well as the mass of the fragments produced. This model connects the discrete nature of fragmentation, which results in the formation of countable stars, with the continuous, multi-scale structure of the molecular clouds from which young stars inherit their mass. The discrete nature of fragmentation process is represented as a network organised into different spatial scales populated by fragments interconnected according to their genealogical affiliation. This architecture provides a framework to perform a multi-scale analysis that can be compared to the continous model of fragmentation. Our model introduces two continuous parameters that determine the multiplicity and mass of fragments at any spatial scale R:

1. the fragmentation rate ϕ(R), which controls the amount of fragments produced at each scale and the multiplicity of young stars in stellar systems;

2. the mass transfer rate ξ(R), characterises the efficiency of fragment formation.

The fragmentation rate and mass transfer rate were used to compute the average number of fragments produced, their average mass, their average mass density, and their average stellar density. We also associated three fragmentation modes considering the numerical values of the fragmentation rate ϕ and the mass transfer rate ξ:

1. a dispersal mode if ξ < −3, or by definition if ϕ > 0;

2. an effective monolithic collapse mode if ϕ ⩾ ξ + 3 ⩾ 0, or if ϕ = 0;

3. a hierarchical fragmentation mode if 0 <ϕ<ξ + 3.

Then we implemented a gravo-turbulent fragmentation in our model to analytically compute the fragmentation rate ϕ, evaluate its variation along the spatial scales, and compare it with observational measurement. The mass transfer rate ξ was linked to the mass budget of fragment formation that represents the balance between accretion processes (Bonnell et al. 2001; Longmore et al. 2014), which add mass to the fragments, and outflow processes (Krumholz et al. 2012; Clark & Whitworth 2021) or turbulence (Camacho et al. 2016), which deplete the fragments of their mass. Although this balance was not explicitly modelled, it is necessary if one needs to quantitatively constrain the mass transfer rate ξ with respect to any mechanism that locally modifies the mass contained within the fragments. We highlighted three major results regarding the cloud structure, the fragmentation property, and the resulting star formation scenario.

Cloud structure. We find that molecular clouds do not systematically exhibit fractal structure as the fragmentation rate is not scale-free in the general case. We highlight two characteristic spatial scales affecting fragmentation. The first is the typical sonic scale R ≈ 0.1 pc, from which fragmentation shifts from a turbulent regime to a thermal regime. The second scale R_stop ~ 40 AU marks the point beyond which a structure becomes optically opaque to its own radiation field, preventing further fragmentation from occurring, while allowing the first Larson core to form. Under supersonic turbulence the cloud may exhibit fractal properties at large scales. Similarily in the subsonic regime, scale-free Jeans-like fragmentation can be retrieved for marginally unstable fragments with an isothermal EOS.

Fragmentation inherent property. This gravo-turbulent model of fragmentation also revealed that fragmentation is a self-regulated process. At similar spatial scales, comparatively denser structures tend to fragment more, producing more diffuse fragments that subsequently subfragment less and form denser fragments. As a consequence, fragmentation progressively loses the information about the initial conditions of the cloud. At small scales, fragments eventually reach marginal instability as mass evolves linearly with their spatial scale M ∝ R, characteristic of marginally unstable isothermal Bonnor-Ebert structures.

Star formation scenario. In our model, hierarchical fragmentation produces pre-stellar structures of typical mass M ~ 10⁻³ M_⊙, for which mass accretion is necessary to produce stars of characteristic mass M ~ 10⁻¹ M_⊙. On the other hand, the comparison of our model with the Taurus NESTs’ multiplicity and size (Joncour et al. 2018) indirectly suggests that fragments are produced more efficiently at low scales, which implies that the mass transfer rate ξ(R) < 0 increases as the spatial scale R decreases. These measurements are compatible with a sequential star formation scenario in the early stages of fragmentation at large scales, dense structures are not formed efficiently as random turbulent motion could unbind fluid particles from the fragments material (ξ < 0). Then, as thermal instability becomes the dominant driver for fragmentation, turbulence mass dispersal is progressively reduced (ξ increases). Finally, pre-stellar structures accrete their surrounding material (ξ > 0) at small scales without fragmenting so that stars can form.

5.2 Potential of this model

The main purpose of our model is to test hierarchical fragmentation as the mechanism setting the discrete properties of stars with respect to their diffuse environment at larger scale. Under mass-loss conditions for the fragments formation (i.e. ξ < 0), the fragmentation cascade can explain the origin of NESTs (Joncour et al. 2018) and the subfragmentation of gas clumps observed in NGC 2264 (Thomasson et al. 2022). However, further investigation is needed in order to fully map the transition between the dense structures located within molecular clouds and the young stars they form through fragmentation. In particular, we did not quantify how stellar IMF may rise from the fragmentation of cores composing a CMF (Larson 1973; Elmegreen & Mathieu 1983; Guszejnov & Hopkins 2015, Thomasson in prep.), nor did we assess how these stellar objects are dynamically bound in order to form stable multiple systems. In addition, the model we implement focuses on the isotropic fragmentation of roundish gas clumps, but our framework can be adapted and improved to investigate the fragmentation of alternative geometries such as filaments, galactic disks, or proto-planetary disks. To describe different objects, different turbulent or magnetic prescriptions should be implemented so the model can be reprocessed in order to derive different fragmentation rate functions.

In order to obtain better constraints on the fragmentation rate, multi-scale analysis should be performed on a wide range of spatial scales, for example within different molecular cloud environments. We can also challenge the model’s theoretical prediction, given a set of implemented physical processes, using hydrodynamical or magneto-hydrodynamical numerical simulations of molecular clouds down to core scale. These simulated data would provide useful theoretical constraints regarding the fragmentation rate. As these simulations describe the evolution of a cloud it could be possible to evaluate the temporal dependency of both the fragmentation rate and the mass transfer rate.

Although some observations have highlighted a correlation between the mass efficiency and the column density of dense clumps (Louvet et al. 2014), no information is provided concerning the multi-scale behaviour of the mass efficiency. As the fragmentation rate, numerical simulations could be employed in order to measure the mass transfer rate. Hence, we could build a semi-empirical law to constrain our fragmentation model and have a mean for comparison with indirect observations. In a more general perspective, measurements of the fragmentation rate and mass transfer rate could both be used as a proxy to directly compare, under the same theoretical basis, the multi-scale properties of the structures observed in real molecular clouds, as we did in this work with NGC 2264, with those that emerge from the numerical simulations.

Acknowledgements

The authors wish to acknowledge the anonymous referee whose useful comments helped us to clarify this work. B.T. carried out this project under equal supervision of I.J. and E.M. who both provided comments and suggestions on this manuscript. F.L., F.M., M.G. and T.N. provided helpful discussion and comments on the paper. We thank Patrick Hennebelle for in-depth discussions which greatly improved the quality of this work. BT has received financial support from the European Research Council synergy grant “ECOGAL” (Grant 855130), the “StarFormMapper” project funded by the European Union’s Horizon 2020 Research and Innovation Action programme (Grant 687528) and from the French Agence Nationale de la Recherche (ANR) through the project “COSMHIC” (ANR-20-CE31-0009).

Appendix A Calculation of the fragmentation rate

This Appendix aims to detail the formula for the fragmentation rate $\varphi (R,\overline{\rho })$ given in Sect. 3.2.1. We recall here the Eq.36 which determines the fragmentation rate as a function of the probability 𝒫(δ > δ_c) that a density fluctuation δ induces a local collapse with respect to a threshold fluctuation δ_c: $$\varphi (R)=3-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}\left(\delta ⩾{\delta }_c\right)}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R},$$(A.1)

with $$\mathcal{P}\left(\delta ⩾{\delta }_c\right)=\frac{1}{2}\left(1-\mathrm{erf}\left[\frac{{\delta }_c+\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}\right]\right).$$(A.2)

Our aim is to compute $\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$. Let $u=\frac{{\delta }_c+\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}$ such that $$\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=-\frac{1}{2\mathcal{P}}\frac{d\mathrm{erf}(u)}{du}\frac{du}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}.$$(A.3) $$\text{ As }\frac{d\mathrm{erf}(u)}{du}=\frac{2}{\sqrt{\pi }}\exp \left[-u^2\right]\text{ we obtain }$$ $$\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=-\frac{1}{\mathcal{P}}\frac{\exp \left[-u^2\right]}{\sqrt{\pi }}\frac{du}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}.$$(A.4)

Expanding the derivative $\frac{du}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ we finally obtain $$\frac{du}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{1}{\sigma \sqrt{2}}\left[\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}-\sqrt{2}v\frac{d\sigma }{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}\right],$$(A.5)

where we let $v=\frac{{\delta }_c-\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}$

Let $A_{\text{p}}\left({\delta }_c(R);\sigma (R)\right)=-\frac{\exp \left[-u^2\right]}{{\mathcal{P}}_{\sigma }\sqrt{2\pi }}$ so we can write $$\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=A_{\text{p}}\left[\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}-\sqrt{2}v\frac{d\sigma }{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}\right].$$(A.6)

A.1 Critical threshold

First we compute the derivative of δc appearing in the expres-sion of $\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ of Eq. A.6. We use the threshold determined by Hennebelle & Chabrier (2008) that consider both thermal and turbulent support: $${\delta }_c=\text{\hspace{0.17em}ln\hspace{0.17em}}\left[\frac{{\pi }^{5/3}}{6^{2/3}}\frac{3C_s^2+V_{\text{rms}}^2}{3G{R}^2\overline{\rho }}\right],$$(A.7)

where G = 6 67 × 10⁻¹¹ m³kg⁻¹s⁻² is the gravitational constant, C_s is the sound velocity, $\overline{\rho }$ the parent mean density and V_rms the turbulent velocity.

We can compute the derivative of $\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ appearing in Eq. A.6 that sets the fragmentation rate as $$\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{3\frac{d{C}_s^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}+\frac{d{V}_{\text{rms}}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}-\frac{3C_s^2+V_{\text{rms}}^2}{3G{R}^2\overline{\rho }}\frac{d\left(3G{R}^2\overline{\rho }\right)}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}}{3C_s^2+V_{\text{rms}}^2}.$$(A.8)

With $$\frac{1}{3G{R}^2\overline{\rho }}\frac{d\left(3G{R}^2\overline{\rho }\right)}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\left(3G{R}^2\overline{\rho }\right)}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\varphi -\xi -1,$$(A.9)

we obtain $$\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{3\frac{d{C}_s^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}+\frac{d{V}_{\text{rms}}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}}{3C_s^2+V_{\text{rms}}^2}+\xi -\varphi +1.$$(A.10)

A.1.1 Derivative of $C_s^2$

At temperature T the sound velocity can be expressed as $$C_s^2=\frac{k_{B}T}{\mu m_H}.$$(A.11)

where k_B = 1.38 × 10⁻²³ J/K, µ = 2.33 et m_H = 1.67 × 10⁻²⁷ kg respectively the Boltzmann constant, mean molecular weight and hydrogen mass. By computing $$\frac{1}{C_s^2}\frac{d{C}_s^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}{C}_s^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R},$$(A.12)

and considering the Eq 20 that express the log-derivation of density, we obtain $$\frac{d{C}_s^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=C_s^2\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }(\varphi -\xi -3).$$(A.13)

A.1.2 Derivative of $V_{\text{rms}}^2$

The quadratic velocity resulting from Larson’s laws (Larson 1981) is expressed as $$V_{\text{rms}}^2=V_0^2{\left(\frac{R}{1\text{\hspace{0.17em}pc}}\right)}^{2\eta }.$$(A.14)

where V₀ = 1000 m/s and η = 0.38 represents respectively the velocity dispersion of a region of size R = 1 pc and the power index characterising the turbulent cascade across the spatial scales. By computing $$\frac{1}{V_{\text{rms }}^2}\frac{d{V}_{\text{rms }}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}{V}_{\text{rms }}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R},$$(A.15)

we obtain $$\frac{d{V}_{\text{rms}}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=2\eta V_{\text{rms}}^2.$$(A.16)

A.1.3 Derivative of δ_c

Injecting both the derivatives of sound and turbulent velocities Eqs. A.13 and A.16 within the expression of $\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ Eq. A.10, we finally obtain $$\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=(\xi -\varphi +3)\left[1-{\widehat{C}}_s^2\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\right]-2\left[1-\eta {\widehat{V}}_{\text{rms}}^2\right],$$(A.17)

where we defined the reduced variables ${\widehat{C}}_s^2$ and ${\widehat{V}}_{\text{rms }}^2$: $$\{\begin{array}{l}{\widehat{C}}_s^2=\frac{3C_s^2}{3C_s^2+V_{\text{rms }}^2}\hfill \\ {\widehat{V}}_{\text{rms }}^2=\frac{V_{\text{rms }}^2}{3C_s^2+V_{\text{rms }}^2}=1-{\widehat{C}}_s^2\hfill \end{array}.$$(A.18)

A.2 Variance of density fluctuations

Next we compute the derivative of σ(R) appearing in the expression of $\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ of Eq. A.6 with $${\sigma }^2=\text{\hspace{0.17em}ln\hspace{0.17em}}\left[1+b^2{ℳ}^2\right],$$(A.19)

where b² ≈ 0.25. We compute $$\frac{d{\sigma }^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{d{ℳ}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}\frac{b^2}{1+b^2{ℳ}^2}.$$(A.20)

A.2.1 Derivative of ℳ²

The Mach number ℳ is defined as the ratio between the turbulent velocity to the sound velocity $$ℳ=\frac{V_{\text{rms }}}{C_s}.$$(A.21)

By computing $\frac{1}{{ℳ}^2}\frac{d{ℳ}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ and substituing both the derivatives of sound and turbulent velocities by Eqs. A.13 and A.16 obtain $$\frac{d{ℳ}^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}={ℳ}^2\left[2\eta -\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }(\varphi -\xi -3)\right].$$(A.22)

A.2.2 Derivative of σ²

Injecting the derivatives of the Mach number Eq. A.22 within the expression of $\frac{d{\sigma }^2}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}$ Eq. A.20, we finally obtain $$\frac{d\sigma }{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{{\widehat{ℳ}}^2}{2\sigma }\left[2\eta +\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }(\xi -\varphi +3)\right],$$(A.23)

where we defined the reduced variables ${\widehat{ℳ}}^2$: $${\widehat{ℳ}}^2=\frac{b^2{ℳ}^2}{1+b^2{ℳ}^2}.$$(A.24)

A.3 The fragmentation rate

We now have the following equations: $$\{\begin{array}{l}\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=(\xi -\varphi +3)\left[1-{\widehat{C}}_s^2\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\right]-2\left[1-\eta {\widehat{V}}_{\text{rms}}^2\right]\hfill \\ \frac{d\sigma }{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=\frac{{\widehat{ℳ}}^2}{2\sigma }\left[2\eta +\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }(\xi -\varphi +3)\right]\hfill \\ \frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}=A_{\text{p}}\left[\frac{d{\delta }_c}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}-\sqrt{2}v\frac{d\sigma }{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}\right]\hfill \\ \varphi (R)=3-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}\mathcal{P}}{d\text{\hspace{0.17em}ln\hspace{0.17em}}R}\hfill \end{array}.$$(A.25)

We showed that both the critcal threshold δ_c and the standard deviation σ depend on the fragmentation rate ϕ but also on the mass transfer rate ξ. To obtain the actual function that describes the fragmentation rate at a given R, one needs to isolate ϕ(R). After refactoring the set of equations (Eqs. A.25) and isolating the fragmentation rate we finally find $$\varphi (R)=\frac{3}{1-A_{\text{th}}}-\frac{A_{\text{th}}(\xi +3)}{1-A_{\text{th}}}+\frac{2A_{\text{turb}}}{1-A_{\text{th}}},$$(A.26)

where we let: $$\{\begin{array}{l}{A}_{\text{p}}=-\frac{\exp \left[-u^2\right]}{\mathcal{P}{\sigma }_{}\sqrt{2\pi }}\hfill \\ A_{\text{th}}=A_{\text{p}}\left[1-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\left({\widehat{C}}_s^2+\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }\right)\right]\hfill \\ A_{\text{turb }}=A_{\text{p}}\left[1-\eta \left({\widehat{V}}_{\text{rms }}^2-\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }\right)\right]\hfill \end{array},$$(A.27)

and: $$\{\begin{array}{l}u=\frac{{\delta }_c+\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}\hfill \\ v=\frac{{\delta }_c-\frac{{\sigma }^2}{2}}{\sigma \sqrt{2}}\hfill \end{array}.$$(A.28)

Appendix B Asymptotic limits

B.1 Subsonic limit

First we assume that V_rms ≪ C_s such that turbulent support is negligible. In this limit, the Mach number ℳ² goes to 0. By using the expansion ln(1 + x) ≈ x, we also compute the asymptotic limit of density fluctuations variance σ²: $${\sigma }^2\approx b^2{\left(\frac{V_{\text{rms}}}{C_s}\right)}^2,$$(B.1)

which goes to 0 in the subsonic limit. We can also compute in this limit the following asymptotic relations: $$\{\begin{array}{l}{\lim }_{\sigma \to 0}u=\frac{{\delta }_c}{\sigma \sqrt{2}}\hfill \\ {\lim }_{\sigma \to 0}v=\frac{{\delta }_c}{\sigma \sqrt{2}}\hfill \\ {\lim }_{ℳ\to 0}{\widehat{C}}_s^2=1\hfill \\ {\lim }_{ℳ\to 0}{\widehat{V}}_{\text{rms}}^2=0\hfill \\ {\lim }_{ℳ\to 0}{\widehat{ℳ}}^2\approx {\sigma }^2\approx b^2{\left(\frac{V_{\text{rms}}}{C_s}\right)}^2\hfill \end{array}.$$(B.2)

We may also compute the asymptotic limit of the quantity $\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }$ involved within the coefficients A_th and A_turb as it will become useful in the following: $$\underset{\sigma \to 0}{\lim }\frac{v{\widehat{ℳ}}^2}{\sqrt{2}\sigma }=\frac{{\delta }_c}{2}.$$(B.3)

Next we can compute the asymptotic limit of Ap as given in Eq. A.27. In the subsonic limit the variable u goes to +∞ if δ_c > 0 or goes to −∞ if δ_c < 0. Assume first that δ_c > 0, we computed the case δ_c < 0 in Appendix B.3. With δ_c > 0 $$\underset{u\to +\infty }{\lim }\mathcal{P}\left(\delta ⩾{\delta }_c\right)=\frac{1}{2}\left(1-\underset{u\to +\infty }{\lim }erf(u)\right).$$(B.4)

The expansion for u → +∞ of the error function is $$erf(u)\approx 1-\frac{\exp \left[-u^2\right]}{u\sqrt{\pi }}.$$(B.5)

In that case we can write 𝒫(δ ⩾ δ_c) as $$\underset{u\to +\infty }{\lim }\mathcal{P}\left(\delta ⩾{\delta }_c\right)=\frac{1}{2}\frac{\exp \left[-u^2\right]}{u\sqrt{\pi }},$$(B.6)

and we can compute the asymptotic limit of A_p $$\underset{\sigma \to 0}{\lim }{A}_{\text{p}}=-\frac{{\delta }_c}{\sigma },$$(B.7)

which goes to −∞ for δ_c > 0. Then we can compute the limits of A_th and A_turb given Eq. A.27: $$\{\begin{array}{l}{\lim }_{\sigma \to 0}{A}_{\text{th}}=A_{\text{p}}\left[1-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\left(1+\frac{{\delta }_c}{2}\right)\right]\hfill \\ {\lim }_{\sigma \to 0}{A}_{\text{turb }}=A_{\text{p}}\left[1+\eta \frac{{\delta }_c}{2}\right]\hfill \end{array}.$$(B.8)

Using these asymptotic limits, the fragmentation rate ϕ(R) given Eq. A.26 can finally be written as $$\varphi (R)\approx \xi +3-\frac{2\left(1+\eta \frac{{\delta }_c}{2}\right)}{1-\frac{d\text{\hspace{0.17em}ln\hspace{0.17em}}T}{d\text{\hspace{0.17em}ln\hspace{0.17em}}\rho }\left(1+\frac{{\delta }_c}{2}\right)}.$$(B.9)

B.2 Supersonic limit

Next we assume that V_rms ≫ C_s such that turbulent support dominates. In this limit, the Mach number ℳ² goes to +∞ and so is the density fluctuations variance σ². We can compute in the supersonic limit the following asymptotic relations: $$\{\begin{array}{l}{\lim }_{\sigma \to +\infty }u=\frac{\sigma }{2\sqrt{2}}\hfill \\ {\lim }_{\sigma \to +\infty }v=-\frac{\sigma }{2\sqrt{2}}\hfill \\ {\lim }_{ℳ\to +\infty }{\widehat{C}}_s^2=0\hfill \\ {\lim }_{ℳ\to +\infty }{\widehat{V}}_{\text{rms }}^2=1\hfill \\ {\lim }_{ℳ\to +\infty }{\widehat{ℳ}}^2=1\hfill \end{array}.$$(B.10)