Issue 
A&A
Volume 574, February 2015



Article Number  A130  
Number of page(s)  11  
Section  Astrophysical processes  
DOI  https://doi.org/10.1051/00046361/201424705  
Published online  06 February 2015 
Shell instability of a collapsing dense core
^{1}
Laboratoire AIM, ParisSaclay, CEA/IRFU/SAp – CNRS – Université Paris
Diderot,
91191
GifsurYvette Cedex,
France
email:
eva.ntormousi@cea.fr
^{2}
LERMA, UMR CNRS 8112, École Normale Supérieure, 75231
Paris Cedex,
France
Received: 29 July 2014
Accepted: 7 November 2014
Aims. Understanding the formation of binary and multiple stellar systems largely comes down to studying the circumstances under which a condensing core fragments (or not) during the first stages of the collapse. However, both the probability of fragmentation and the number of fragments seem to be determined to a large degree by the initial conditions. In this work we explore this dependence by studying the fate of the linear perturbations of a homogeneous gas sphere, both analytically and numerically.
Methods. In particular, we investigate the stability of the wellknown homologous solution that describes the collapse of a uniform spherical cloud. One problem that arises in such treatments is the mathematical singularity in the perturbation equations, which corresponds to the location of the sonic point of the flow. This difficulty is surpassed here by explicitly introducing a weak shock next to the sonic point as a natural way of connecting the subsonic to the supersonic regimes. In parallel, we perform adaptive mesh refinement (AMR) numerical simulations of the linear stages of the collapse and compare the growth rates obtained by each method.
Results. With this combination of analytical and numerical tools, we explore the behavior of both axisymmetric and nonaxisymmetric perturbations. The numerical experiments provide the linear growth rates as a function of the core’s initial virial parameter and as a function of the azimuthal wave number of the perturbation. The overlapping regime of the numerical experiments and the analytical predictions is the situation of a cold and large cloud, and in this regime the analytically calculated growth rates agree very well with the ones obtained from the simulations.
Conclusions. The use of a weak shock as part of the perturbation allows us to find physically acceptable solutions to the equations for a continuous range of growth rates. The numerical simulations agree very well with the analytical prediction for the most unstable cores, while they impose a limit of a virial parameter of 0.1 for core fragmentation in the absence of rotation.
Key words: stars: formation / ISM: clouds
© ESO, 2015
1. Introduction
By now it is well understood that stars form out of the gravitational condensation of cold and dense molecular cores, which are very structured in density and in velocity (Larson 1981; Falgarone & Phillips 1990). In parallel, it has also been argued that most stars live in binaries or multiple systems (e.g., Duquennoy & Mayor 1991; Goodwin et al. 2007). It is therefore natural to assume that owing to these internal velocities and density enhancements, a core typically fragments at a certain stage of its collapse and eventually produces two or more stars.
From an observational perspective, the matter can be elucidated by registering the degree of multiplicity at different stages of the core contraction. Most data for early stellar multiplicity are available for Class I/II objects, in other words, for protostars with almost no trace of their initial envelope left (Ghez et al. 1993; Patience et al. 2002; Duchêne et al. 2004, 2007). But, to clarify the dependence on the initial physical properties and internal structure of a core, one must allude to even earlier phases in the collapse (Class 0 objects), when the protostars are still largely embedded in their parental core (Andre et al. 1993, 2000). Although challenging because of the obscuration from the envelope, there have been numerous attempts to constrain stellar multiplicity that early in the collapse stage of a core (Chen et al. 2008; Girart et al. 2009; Duchêne et al. 2004, 2007; Jørgensen et al. 2007). The results show an increase in the number of fragments from the Class 0 to the Class I stage, which is evidence that suggests that fragmentation is favored either in the Class 0 stage or very shortly afterward (Maury et al. 2010).
There have been also numerous theoretical studies, both of the collapse process and of the fragmentation during collapse. Ebert (1955) and Bonnor (1956) independently derived an analytic criterion for a static isothermal sphere to be stable against gravitational collapse. In their numerical simulations of collapsing, initially uniform isothermal spheres, Larson (1969) and Penston (1969) found a onedimensional, selfsimilar solution (called LP solution or flow in the following). If x is the selfsimilar variable of the profile, this solution describes a cloud with a flat density profile that turns into a x^{2} dependence toward infinity and a velocity profile proportional to x. It was later discovered by Shu (1977) that this solution is a member of a whole family of selfsimilar collapse solutions, many of which contain critical points. In that paper, he arrived at the conclusion that an initially centrally condensed cloud will collapse from the inside out, establishing an expansion wave. This behavior was integrated into a more complete picture of the collapse and the behavior of selfsimilar collapse solutions by Whitworth & Summers (1985).
The stability of the collapse solutions is understood to a much lesser extent. For the case of a static, uniform, presureless spherical cloud, Hunter (1964) showed that there is an unstable shell mode that grows like (t_{0} − t)^{1}, where t is the time and t_{0} the time at which collapse occurs.
Hanawa & Matsumoto (1999, HM99 in the following) performed a stability analysis of the LP flow using a shooting method. In this type of analysis, one starts integrating the equations from one boundary and varies some parameters until the solution at the other boundary matches the desired conditions there. They came to the conclusion that the solution is stable overall, with the exception of a slowly growing l = 2 mode (where l the azimuthal wave number of a spherical harmonic perturbation) with a growth rate of 0.354. This very weak instability is consistent with the analysis by Ori & Piran (1988), who derived a stability criterion for selfsimilar isothermal collapse flows that is based on the gradient of the radial velocity close to the critical point. According to that criterion, a homogeneous isothermal sphere is an unstable configuration that could naturally converge to a (generally much stabler) LarsonPenston type flow.
While quite interesting, this result poses a few questions. One is that the growth rate is rather slow, so it is unclear whether perturbations can really grow sufficiently (see Sect. 5). Another issue is that only the l = 2 mode was found to be unstable, since HM99 report the nonexistence of eigenvectors for higher values of l. Finally, no total eigenvector has been calculated for the spherical mode. Altogether, this suggests that the stability analysis of a collapsing cloud is not yet complete.
On the other hand, there is abundant literature on numerical simulations of the fragmentation of rotating and/or turbulent cores. Boss & Bodenheimer (1979) studied the fate of an l = 2 perturbation on a collapsing core, and Burkert & Bodenheimer (1993) repeated the experiment, finding that a filament connecting the fragments should develop. (This experiment has been repeated many times since and is used extensively as a code testing tool.) A criterion for fragmentation in the presence of rotation was provided by the semianalytical work of Tohline (1981), quantifying thermal support with the virial parameter α and the rotational versus gravitational energy in the cloud with the corresponding parameter β. The quantity αβ has since been used repeatedly to delimit the conditions for fragmentation (Miyama et al. 1984; Hachisu & Eriguchi 1984; Tsuribe & Inutsuka 1999; Tohline 2002, for instance), and it was found that αβ< 0.1−0.2 typically leads to fragmentation.
At the same time, increasingly more complex direct simulations of the collapse and fragmentation processes have been performed, including the internal and external thermal pressure of the core, rotation (Myhill & Kaula 1992; Cha & Whitworth 2003; Matsumoto & Hanawa 2003; Hennebelle et al. 2004; Machida et al. 2008; Commerçon et al. 2008), magnetic fields (Banerjee & Pudritz 2006; Price & Bate 2007; Hennebelle & Teyssier 2008; Boss & Keiser 2013), and turbulence (Klessen et al. 1998; Klessen & Burkert 2001; Offner et al. 2008; Bate 2009; Joos et al. 2013), which seem to affect the number of fragments and their separation. Girichidis et al. (2011) conducted a parameter study in which they varied the initial density profile of the core and the level of turbulence. They provide evidence of the strong influence of the mean radial density profile to the result of the collapse.
Given the abundance and complexity of the existing models for fragmentation, it is somewhat surprising that there is so little to be said for the stability of a solution that is as simple as the homologous, uniform, isothermal collapse. It is clear that this solution has limited applicability to real cores, which are sometimes far from isothermal spheres and do host turbulent motions and magnetic fields, which are all potentially important and complex effects that can substantially alter the behavior of the solution. However, a study of its linear stability offers important insight into the principal mechanisms causing fragmentation, since it does not suffer from the complexity of nonlinear effects and can contribute to the overall understanding of core collapse.
In the present paper we rectify by performing a linear stability analysis of this flow. Our study shows that the homologous solution is indeed unstable and we obtain the corresponding growth rates. This analysis is presented in Sects. 2 and 3. In the fourth section we present numerical simulations of a uniform dense core subject to a spherical perturbation and we look for the range of parameters for which it is unstable. We measure the growth rate of the shell mode (a spherically symmetric perturbation), as well as higherorder spherical harmonic perturbations and find them to be in good agreement with the result of our linear analysis. The fifth section discusses the implications of our results and proposes further interpretation. The sixth section concludes the paper.
2. Stability analysis: the spherical case
2.1. Equations and selfsimilarity
We investigate the stability of a collapsing isothermal sphere against linear perturbations. Given the symmetry of the problem, we start with the equations of hydrodynamics in spherical coordinates in one dimension: (1)where ρ the gas mass density, u_{r} the radial velocity, and ψ the gravitational potential. It is well known (Larson 1969; Penston 1969; Shu 1977) that this system admits selfsimilar solutions of the form (2)where the density ℛ, the radial velocity , and the gravitational potential, Φ are functions of the selfsimilar variable X, C_{s} is the sound speed, and G the gravitational acceleration. Among the various solutions known in the literature the simplest is the homologous solution (3)which describes a cloud with a uniform density R_{0} and a velocity field proportional to X with a slope α_{0}. Clearly, like other selfsimilar profiles, this solution is only valid within certain spatial and temporal limits, as both the mass and the velocity diverge for large values of X. In addition, in a real collapse situation there is a rarefaction wave propagating outwards (Tsuribe & Inutsuka 1999; Truelove et al. 1997), which breaks the selfsimilarity. In order then for this solution to remain valid the gas has to be cold enough for the collapse to happen before the rarefaction wave can reach the center. Then Eqs. (3) can be used to describe the inner parts of the cloud.
The simplicity of this profile and the fact that it describes the center of a cold core quite accurately make it a popular initial condition for numerical simulations (Larson 1969; Price & Bate 2007; Commerçon et al. 2008; Boss & Keiser 2013, for example). It is therefore both interesting and relevant to perform a stability analysis of this solution and to confront the analytical results with numerical simulations.
We first look for perturbations of the form (4)where R_{0}, U_{0}, and φ_{0} denote the equilibrium state and δR(X), δU(X), and δφ(X) the perturbed quantities. The quantity σ in the exponent is assumed to be positive and represents the growth rate of the perturbation. Inserting these expressions into Eqs. (1), we get with primes denoting the spatial derivatives. Combining Eqs. (5) and (7), one can show that (8)where K is a constant chosen to ensure the continuity of δφ. From Eqs. (3), (5), (6), we get (9)which together with Eqs. (6) and (8) allows us to solve the system. Once the boundary conditions have been specified, a numerical integration is carried out by a standard RungeKutta scheme.
2.2. Inner boundary and the critical point
The boundary conditions for the perturbation at X ≃ 0 are given by (10)where cst is a constant. So by simply specifying a constant value for R_{1} we can integrate the system of Eqs. (6) and (8) for different values of σ.
One complication to this otherwise straightforward approach is that the system possesses a critical point at X_{crit} = 3, where the gas becomes supersonic with respect to the selfsimilar profile. The presence of such a boundary poses a restriction to our numerical solutions, since only the ones which cross it can be considered physically acceptable.
It is easy to show that the system (5)–(7) has two exact solutions that resemble these boundary conditions and are given by While mathematically acceptable, these solutions offer little physical insight: they merely describe a variation of the mean density. They are, however, mathematically acceptable, since they remain small at any X with respect to the perturbed solutions.
Apart from these trivial forms, the integrations of the system of Eqs. (6) and (8) we performed for different values of σ gave no other solutions able to cross the critical point.
2.3. Shock conditions
From the above discussion it becomes clear that a more general form of the perturbation should exist which allows for physically meaningful solutions to the corresponding perturbation equations and which does not suffer from the issue of critical point crossing.
To begin with, one can imagine a core whose size is much larger than the sonic radius. In this case the perturbation can be restricted to the supersonic parts, where X is much greater than the critical value. For a cloud this large we are essentially in the same regime as the cold cloud described by the homologous collapse solution (Eqs. (3)), where gravitational collapse happens before thermal pressure effects have had time to break the selfsimilarity. We expect then that the behavior of the perturbations will be similar to the general t^{1} growth found by Hunter (1964) for the case of a cold cloud.
In order to integrate the solution for the outer parts of the cloud, we need a physically meaningful boundary condition at the vicinity of the critical point. We thus seek solutions that include a shock. In the vicinity of the critical point, the velocity of the fluid is by definition transonic with respect to the selfsimilar profile, so an arbitrary weak shock can naturally happen in this area. Note that, although the spatial derivatives of the density and velocity fields are infinite at the shock, all fields remain small if the shock is weak. Therefore, such a discontinuity indeed generates a linear perturbation.
Let the shock be located at X_{shock} = X_{crit} + δX where δX> 0. This is the location where the unperturbed solution, given by δR = 0, δU = 0 is connected to a state given by the RankineHugoniot (RH) conditions, which in the frame of reference of the shock are (14)where V_{shock} is the velocity of the shock and subscripts 1 and 2 mark the pre and post shock quantities, respectively.
The velocity of the flow with respect to the selfsimilar profile is given by α_{0}X + X = X/ 3. Thus at X_{shock}, the Mach number can be expressed as ℳ = X_{shock}/ 3 = 1 + δX/ 3. In the frame of the shock, the gas from the subsonic region is entering supersonically into the shock and therefore constitutes the preshock medium whose density is amplified.
Since we consider that the subsonic region is unperturbed, the expressions for ρ_{1} and u_{r,1} that enter the jump conditions in Eqs. (14) are simply given by Eqs. (2) and (3).
By replacing we get (15)At the limit of a weak shock, we can expand these relations and get (16)As the perturbation develops it moves toward larger δX, so the density increases and the shock becomes stronger.
Since the perturbed density remains zero up to the position of the shock, the constant K in Eq. (8) must be chosen in such a way that δφ′(X_{shock}) = 0.
2.4. Limit at large radii
In order to identify the physically relevant solutions, the asymptotic behavior at large X must also been known. From Eqs. (5)–(7), it is easy to infer that there must be an asymptotic behavior of the form where the exponent v in Eqs. (17) and (18) should be chosen such that the asymptotic form satisfies Eqs. (5)–(7). By inserting expressions (17) and (18) into Eqs. (5)–(7) and dropping the thermal pressure term, which is negligible at large X, we obtain two linear equations for R_{∞} and U_{∞}, which depend on σ and v. By demanding that the determinant of this system be zero, we get the two solutions for v: (19)We thus recover the exact homologous solution mentioned above, since when σ = 1 or −2/3, we get v = −1, which means U ∝ X and R = cst. On closer inspection we also see that the first branch (v = 3σ−4) diverges most rapidly. Since mathematically it is required that v> −1 for the solutions to behave well towards infinity, it follows that σ> 1. When 1 <σ< 4/3, the velocity diverges at infinity but still remains everywhere negligible with respect to the perturbed solution.
2.5. Numerical integration
The shock is now essentially the lowX boundary for the numerical integration, so we only need to specify its position δX from the critical point in order to define the boundary conditions. This value must be small in order to ensure linearity.
As shown by Eqs. (16), the exact values of δU and δR are not important (they depend on t_{0}, which can be freely chosen) but their ratio is fixed.
Starting from a location close to zero (specifically X_{s} = 10^{2}), and using a step dX = 10^{3}, we integrate the perturbation equations up to X = 100. We have checked that the solutions are well resolved and converged.
Fig. 1
Density (top panel) and velocity (bottom panel) fields of the solutions for various values of σ. For values of X smaller than the position of the shock (X<X_{sh} ≃ 3), the density and the velocity are unperturbed, i.e. δR = δV = 0. The solid lines show the solutions for the system of Eqs. (5)–(7) For comparison, solutions to the same equations but without the thermal pressure term are shown with the dotted lines. 

Open with DEXTER 
Figure 1 shows the density and the velocity fields for various values of the growth rate, σ. The solid lines show the solutions of Eqs. (5)–(7) and the dotted lines show the solutions for the same system of equations when the thermal pressure is suppressed. The latter case corresponds to a cold cloud and is given for comparison.
The velocity and the density fields vary rapidly in the vicinity of the shock, X ≃ 3. The velocity, which is negative, increases and the density decreases. At X ≃ 4−5, the behavior of the density changes: it starts decreasing much less rapidly and it quickly reaches the asymptotic behavior expressed by Eqs. (17) and (18).
The behavior of the velocity is slightly more complex. For σ< 4/3, it decreases for X> 4−5 as expected from the asymptotic form (Eqs. (17) and (18)). For 4/3 <σ< 1.6, the velocity eventually decreases and tends towards zero after a small increase. Again, this is in good agreement with the expected asymptotic profile. In contrast, for larger values of σ, the behavior is identical, only the velocity becomes positive just after the shock at X ≃ 3−5, which means that the perturbation is expanding.
Based on this form of the perturbations, we believe that the regime 1 <σ< 4/3 corresponds to the most physically relevant solutions. For smaller σ, the velocity diverges at large X and is therefore not really a perturbation with respect to the unperturbed solution. When σ> 4/3, the velocity quickly tends to zero, meaning that it is very localized around the shock. Finally, when σ> ≃ 1.7, the velocity becomes positive. This situation probably corresponds to somehow artificial perturbations that are unlikely to represent a physically relevant case.
This conclusion remains almost identical for the cold case since for σ< 4/3 the velocity diverges toward large X while for σ> ≃ 1.8, the density becomes negative at X> ≃ 5.
To summarize, we find that linear perturbations that include a weak shock in the vicinity of the sonic point present physically acceptable behavior if the growth rate is larger than unity. As the growth rate becomes larger than 4/3, the form of the eigenvectors suggests that the corresponding perturbation is unlikely to occur: the velocity field rapidly converges toward zero, implying that the perturbation is isolated in the sonic region.
3. Stability analysis: the nonaxisymmetric modes
The instability of a shelllike structure is an important result, but if we care about the formation of multiple fragments we must also consider nonspherically shaped perturbations. In this section we present the stability analysis for nonspherically symmetric modes. Since the method is practically the same as the spherical case, here we will only highlight the differences.
3.1. Equations
The fluid equations for a selfgravitating gas in spherical coordinates and three dimensions are: (20)where u_{r}, u_{θ}, and u_{φ} denote the velocities along each of the spherical coordinates r, θ and φ. To study the stability of the homologous solution with respect to nonspherical modes, we look for perturbations of the form (HM99): (21)where V and W are the angular velocities and δV(X), δW(X) are the corresponding perturbations in the selfsimilar frame and are the usual spherical harmonics. As in the study of HM99, the equations for δV and δW are identical so the perturbations are given the same amplitude, δV(X). Replacing these expressions into Eqs. (20), we get Thus, the system of Eqs. (22)–(25) consists of three firstorder ordinary differential equations and one secondorder ordinary differential equation. Like in the spherical case, we solve it using a standard RungeKutta method.
3.2. Inner boundary and shock condition in the nonaxisymmetric case
The boundary conditions of Eqs. (22)–(25) at X ≃ 0 are given by (26)where (27)for the system to have a solution. Equations (22)–(25) admit a family of solutions that resemble these boundary conditions and are given by Eqs. (26) and (28)The situation here resembles what we found for the spherical problem: while these solutions are mathematically acceptable, they diverge at large X so they do not represent physical solutions. Nonetheless, they will play an important role in finding a solution, as we will show later.
It is easy to see that this system also admits a critical point located at X_{crit} = 3. Again, the only solutions that satisfy the inner boundary conditions and cross the critical point are the trivial ones stated by Eqs. (28), which for l> 2 diverge at infinity.
Following the same line of thought, we look for solutions that present a shock at the vicinity of the critical point, X_{shock} = X_{crit} + δX where δX> 0. The unperturbed solution given by δR = 0, δU = 0 is connected to a postshock state given by the RankineHugoniot conditions. The same applies to the other variables δV, δφ and δφ′. They are all assumed to be zero in the inner part before the shock. These variables however are continuous so they are also equal to zero immediately after the shock.
The RankineHugoniot conditions are identical to the spherical case. The same calculations then lead to (29)This implies that the surface of the shock itself is spherical only at the zeroth order. At the first order, the shock surface is described by a spherical harmonic, .
3.3. Limit at large radii for nonaxisymmetric modes
From Eqs. (22)–(25), it is easy to infer that the asymptotic behavior towards infinity is (30)By inserting these expressions into Eqs. (22)–(25) we obtain a fourthorder polynomial whose four roots are: (31)The two first roots are the same ones we obtained before. The third branch, v = −l + 1, leads to a strong divergence with l, with the exception of the value l = 2 for which it leads to the same asymptotic behavior as the solution stated by Eqs. (3). This asymptotic form is associated to the solution stated by Eqs. (28) and it is, generally speaking, unphysical.
However, if this branch is linearly combined with the solution (28) in such a way that their asymptotic behaviors compensate, the new solution can be made to present a physically acceptable asymptotic behavior described by the root v = 3σ−4, like in the spherical case. It is again required that v> −1, so we must have σ> 1. Also, the condition that the velocity go to 0 at large X implies again that σ> 4/3.
To summarize, in order to get physically meaningful perturbations of the exact solution stated by Eqs. (3), we need to combine the solution stated by Eqs. (28) with the solution obtained by applying the RankineHugoniot conditions at the vicinity of the critical point in such a way that the linear combination does not diverge too strongly.
If S_{l} is the desired solution and φ_{1} is the value of the linear combination in Eqs. (28), then before the shock the solution is given by Eqs. (26) and (28). At the shock position the solution must satisfy the jump relations stated by Eqs. (29). The difference with the spherical case is that here u_{r,1} is not zero, but given by Eqs. (28). This leads to (32)while (33)Thus S_{l} is a linear combination of the two above solutions. This is made clearer by writing (34)where . We see that the first terms of the righthand side of Eqs. (34) are identical to Eqs. (29), while the second terms are compatible with Eqs. (28).
It is interesting that, unlike in the spherically symmetric problem, here the perturbations cannot vanish in the subsonic region. Since gravity is a nonlocal force, these perturbations are induced inside the subsonic region by the nonaxisymmetric density distribution of the supersonic region.
Strictly speaking, since the spherical harmonics also take negative values, there are locations where expressions (29) and (34) place the shock in the subsonic regions of the core. In order for the shock to always be outside the critical point or, in other words, in order for δX to always be positive, these expressions should contain a spherically symmetric term. This can be achieved by combining the nonaxisymmetric perturbations with the spherically symmetric perturbation described by Eq. (13) for σ = 1, so essentially adding an extra term, constant in (θ, φ), in Eqs. (29)–(34) that describe the RankineHugoniot conditions. We have omitted such a term for the sake of simplicity, but the subsequent analysis would remain identical. Indeed, we have verified that the solutions depend only weakly on the inner boundary conditions.
3.4. Results for the nonaxisymmetric modes
As in the spherical case, here as well we specify a position for the shock δX small enough to ensure linearity. We use a step dX = 10^{3} and integrate up to X = 10^{4}. At the high X limit we divide the potential by X^{l + 1}, which yields the parameter φ_{1} (see Eq. (33)). To obtain S_{l} then we subtract the unphysical solution (28) from the result of the numerical integration and recover the expected asymptotic behavior for values of X smaller than 10^{4}. An artifact of this subtraction is that, when X ≃ 10^{4} the potential goes abruptly to zero. The integration to larger values of X compared to the spherical problem enables us to discard the diverging term while maintaining a physical behavior for a large enough range of X values.
Fig. 2
From top to bottom, density R, radial velocity U, azimuthal velocity V and gravitational potential φ for the l = 2 mode and for a series of growth rates. 

Open with DEXTER 
Fig. 3
Same as Fig. 2 for the l = 3 mode. 

Open with DEXTER 
Fig. 4
Same as Fig. 2 for the l = 4 mode. 

Open with DEXTER 
Figures 2–4 show the density R, radial and azimuthal velocities U and V and the gravitational potential φ for the azimuthal wave numbers l = 2, 3 and 4 respectively and for a series of growth rates σ. These profiles show the solutions described by Eqs. (28) in the inner subsonic regions, then a shock in the neighborhood of the critical point X_{c} ≃ 3 and finally a behavior at large X consistent with the asymptotic analysis above and similar to the spherically symmetric case.
Although not identical, the three modes behave similarly to the spherically symmetric example. The values of the various fields tend to be smaller for larger l because the gravitational potential has stiffer spatial variations in the azimuthal directions, which lowers the resulting force inside the subsonic regions. Like in the simple onedimensional case, here as well higher values of σ lead to sharply peaked profiles around the shock position. This appears to suggest that the physically relevant values of σ are between 1 and ≃1.5, but at the same time, there does not seem to be an obvious way of deciding whether a particular growth rate is expected or preferred. For this reason and also for the sake of verifying the reliability of the present approach, we perform numerical simulations, which are presented in the following section.
4. Stability of collapsing cores in numerical simulations
The above analysis is complemented with a series of numerical experiments, consisting of collapsing cores with different initial ratios of gravitational to thermal energy (virial parameter α), and different initial perturbations. These simulations are a very efficient way to witness the behavior of a linearly perturbed core in the full complexity that threedimensionality entails, meaning the inclusion of the azimuthal and vertical velocity components. It also allows us to probe the limit for fragmentation as the thermal support of the core increases.
Fig. 5
Contours of the log n_{H} in units of m_{H} cm^{3} along the yz plane. From top to bottom: l = 0, l = 1, l = 2 and l = 3 modes for virial parameter α = 0.006. The units along each axis are parsecs and t_{ff} stands for freefall time. 

Open with DEXTER 
4.1. Code and initial conditions
The simulations are performed with the publicly available Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002), which solves the equations of hydrodynamics on a Cartesian grid with a secondorder Godunov scheme and comes with wide range of options for the simulated physical processes. In this case we take advantage of its AMR capabilities by requiring that the Jeans length be always resolved with at least 20 grid cells. An isothermal equation of state is used throughout.
The initial condition of each simulation is a core of uniform density ρ_{0}, equal to 4 × 10^{18} g cm^{3}, and uniform temperature T, equal to 10 K. To this profile we add a linear perturbation, given by (35)where k_{r} = 0.25 is the radial wavenumber of the perturbation, r is the distance from the core center, L_{c} is the core outer radius and ϵ_{1}, ϵ_{2} the amplitudes of each component of the perturbation, which, according to the case, have values from 0.05 to 0.1. The core resides in a uniform medium of pressure and density 100 times smaller than the core’s interior.
For all these simulations we employ periodic boundary conditions; when selfgravity is involved, such boundaries can potentially affect the evolution of the system, so the core is placed at the center of the computational volume, at a distance of two core radii from each box side. This was shown to be enough to avoid boundary effects by convergence tests we performed for the largest core.
4.2. Stability as a function of cloud parameters
In this study of the core stability, we vary only two parameters: the azimuthal wave number l of the spherical harmonics and the virial parameter of the core, which gives the thermal over gravitational energy ratio in the core and is defined as (36)where c_{s} is the sound speed in the core interior. Here the temperature of the core is held constant and α is varied by changing the radius R of the core. We have not varied the radial wavenumber of the perturbation, which means the l = 0 case is a single shell around the center of the core.
Some examples of the core behavior are shown in Figs. 5 and 6, where we show density contour plots of a cut along the yz plane in the middle of the simulation box for different times of each simulation. For low virial parameters (largest cores, α< 0.15) the density peaks are more and more enhanced with respect to the background, which, as the collapse continues, can lead to the formation of as many objects as the peaks of the initial perturbation. In the intermediate virial parameter regime (0.15 <α< 0.6), the perturbations are erased by the expanding rarefaction wave and the core collapses into one single fragment. In the large virial parameter regime (0.6 <α, smallest cores) the core is no longer held by gravity and reexpands under the influence of its thermal pressure.
Fig. 6
Contours of the log n_{H} in units of m_{H} cm^{3} along the yz plane. From top to bottom: l = 2 mode for virial parameters α = 0.1, α = 0.2 and α = 0.8. Like in Fig. 5, units along each axis are parsecs and t_{ff} stands for freefall time. 

Open with DEXTER 
4.3. Growth rates
The growth rate of the initial perturbation is measured as the rate of change of δρ/ρ_{0}, where δρ = ρ_{max} − ρ_{0}. Here ρ_{max} is defined as the maximum density in the core and ρ_{0} as the central density of the core at each instant. The time is normalized as t′ = 1 − t/t_{ff}, where t the actual time in the simulation and t_{ff} the initial freefall time of the core, t_{ff} = 3π/32Gρ_{0}, where ρ_{0} the initial density of the core. In practice, we plot the logarithm of these quantities and calculate a leastsquare fit to the linear regime of the plot. The slope of this fit is the negative of the growth rate we are after. To make it clearer, if the denote the growth rate with σ, like in the analytical treatment of the previous sections, then ^{(}δR/R_{0}^{)} = ^{(}δR/R_{0}^{)}_{0}t^{′σ} (see Eq. (4) for the definition of the quantities, with the exception of t′, which here denotes the simulation time). By taking the logarithm of this relation one can clearly see that σ becomes the slope of a linear relation between log ^{(}δR/R_{0}^{)} and log (t′). But since in our notation time is negative, we must also take the negative of the slope as the growth rate of the instability.
Some growth rates thus calculated are shown plotted in Figs. 7–9.
Fig. 7
Growth rates for α = 0.006 mode for azimuthal wave numbers l. 

Open with DEXTER 
Following a simple line of thought, we expect the highest growth rates to appear for the smallest l, since the azimuthal flow of material is then focused on a smaller number of fragments, which can therefore grow faster. At the limit of very high l on the other hand, one expects to recover σ values similar to the l = 0 case, since a very large number of fragments in the azimuthal direction is roughly equivalent to a continuous shell structure.
Figure 7 shows the growth rates for the case of the largest core, with a virial parameter of α = 0.006. For the lowl regime (approximately 0 <l< 10), the growth rates inferred from the simulations are above the value of 1, in excellent agreement with the analytical prediction of the previous sections. Also, the highest growth rates are observed for l = 1 and l = 2, in accordance with the simple prediction just above. But for very high values of l, instead of tending towards the value of σ for l = 0 the growth rates in the simulations fall below unity. This disagreement can probably be attributed to a lack of angular resolution.
Figures 8 and 9 show growth rates for the same azimuthal wavenumber (l = 0 and l = 2, respectively), but for different virial parameters α. Both figures show a similar pattern: a decrease in the growth rate as the virial parameter increases, which can be understood as an increasing thermal stability of the core. As thermal pressure becomes important, the rarefaction wave is increasingly more efficient in weakening the perturbation, to the point that, for large enough α, it does not grow anymore.
At this point we should mention that it was impossible to locate the weak shock in the simulations, due to the very limited resolution in the central parts of the core (the region with X< 3 was only resolved with three to five cells). This makes a direct comparison of the eigenvectors very challenging. Nonetheless, the excellent agreement of the growth rates obtained with two such different methods is very encouraging.
5. Discussion and implications
Fig. 8
Growth rates of the l = 0 mode for various virial parameters α. 

Open with DEXTER 
Fig. 9
Growth rates of the l = 2 mode for various virial parameters α. 

Open with DEXTER 
The most important result of the present work is that a homologously collapsing cloud is prone to a shell instability due to its selfgravity, even for nonaxisymmetric perturbations with large azimuthal wave numbers. This instability develops relatively fast when the cloud has little thermal support, showing growth rates typically larger than 1 and up to 1.2–1.3 for l = 2. These growth rates decrease when the thermal support increases and eventually reach zero, meaning that the cloud is not expected to undergo fragmentation, at least not through the shelllike mode studied here.
It is important to contrast this result with the analysis of HM99 for the LP solution. In their study, only the l = 2 mode is unstable and it grows like ≃(t_{0} − t)^{0.354}. Since the density ρ is proportional to (t_{0} − t)^{2}, this implies that the perturbation grows like ≃ρ^{0.175}. Therefore, a perturbation starting with an amplitude ϵ ≃ 0.1 should become nonlinear (ϵ ≃ 1), only when the density within the cloud has grown by 10^{1/0.175} ≃ 5 × 10^{5}. This is such a significant density increase that this instability may never happen. Take for example a dense core with a peak density of 10^{5} cm^{3}. According to the HM99 growth rate, the density perturbation should become nonlinear only when the density is about 5 × 10^{10} cm^{3}. At these densities the cloud is not expected to be isothermal anymore since dust is already opaque to its radiation.
In contrast, if the core is cold enough, the shell mode may grow as ≃(t_{0} − t)^{1} and therefore as ≃ρ^{0.5} This would imply that a 10% amplitude perturbation in a core of mean particle density of 10^{5} cm^{3} becomes nonlinear for a density of about 10^{7}cm^{3}, which is entirely reasonable.
Since the LP flow also describes a very fast collapse, the question arises why do the two types of solutions present such different fragmentation properties? We believe the answer relies on the tidal forces, which differ much in the two cases (see also Jog (2013) for an analysis of the effect of tidal forces on the Jeans stability criterion).
To illustrate this, let us consider a core whose density profile is ρ = Ar^{− p} and a perturbation δρ of characteristic size δr. A fluid particle located at the edge of the perturbation, i.e. at a distance δr from the perturbation center, is subject to the perturbed gravitational force, (37)where δM = 4π/ 3δr^{3}ρ, as well as to the tidal forces, that is the gradient of the gravitational force produced by the mean density. The gravitational acceleration within the core is a_{g} = GM(r) /r^{2} with M = (4π/ 3)(A/ (3−p))r^{3−p} and (38)Thus the tidal force is (39)The fluid particle at the edge of the perturbation feels a total gravitational force equal to (40)Thus the tidal force modifies the effective gravity of the perturbation. We can identify the following regimes for the values of the exponent p in the equilibrium profile:

If p = 0, then δF_{tot} = 4/3dF_{sg}, which means that the tidal force enhances the gravitational instability because it compresses the perturbation.

If p = 1, then δF_{tot} = dF_{sg} and the tidal force has no effect.

If p = 1.5, then δF_{tot} = 2/3dF_{sg}, in other words the tidal force works against total gravity because it tends to shear the fluid elements apart.

Finally, if p = 2, then δF_{tot} = 0 and the tidal force cancels the gravitational force of the perturbation. In this case, only nonlinear or nonaxisymetric perturbations can develop.
This short analysis suggests that the homologously collapsing, uniform density solution is more prone to the gravitational shell instability that the LarsonPenston solution, which presents an r^{2} density profile at large radii. This result is also in excellent agreement with both the analytical and the numerical findings presented in this work.
Our results concerning the role of the equilibrium profile in the fate of the perturbations are also compatible with the numerical simulations performed by Girichidis et al. (2011). These authors performed a series of simulations of massive turbulent cores, in which they vary the initial density profiles of the clouds. In particular, they found that when the density is initially uniform or given by a BonnorEbert density profile, the cloud fragments in many objects. On the contrary, if the density profile is initially proportional to r^{2}, only one or very few fragments form. This is entirely consistent with the above interpretation of the tidal forces.
In terms of understanding the origin of stellar multiplicity, where do the above results leave us? In short, our study complements a large volume of previous literature on the subject of core fragmentation, which so far has been mostly dealing with rotating or turbulent environments. This, of course, is justified by the very structured kinematics typically observed in prestellar cores (Goodwin et al. 2007, for a partial review), which very often do show irrefutable rotation signatures (Goodman et al. 1993).
And although there is little doubt that rotation does lead to fragmentation, here we suggest a “thermal” type of fragmentation that does not require any initial rotation. In a large and cold enough cloud, if its density profile is flat, perturbations will grow to become nonlinear, eventually leading to fragmentation. The fragments, unlike in cases with rotation, which produce tens to hundreds of AU separations, (Hennebelle et al. 2004; Commerçon et al. 2008, for example) will be located at distances of more than 1000 AU during the first core collapse, with a possibility of migrating inward at later stages and, with the onset of rotation, become binaries, like observed, for example in the simulations of Hennebelle & Teyssier (2008) for magnetized rotating environments.
Although limited in their application due to the idealized character of the setup, the results of this study are useful, both for predicting the behavior of initially flat, quiescent cores at a stage of their collapse almost inaccessible observationally, but also, very importantly, fot deciphering the behavior of more complex, dynamical models.
6. Conclusions
We have calculated the form of linear perturbations for the selfsimilar flow that describes the homologous collapse of an isothermal sphere. The calculation was done by introducing a weak shock at the sonic point of the flow and integrating the corresponding perturbation equations from that location to infinity with a RungeKutta scheme. The problem was treated in one dimension for the case of a spherical shell and in three dimensions for nonaxisymmetric perturbations. Results of this analysis follow.

Using a shock at the critical point of the selfsimilar profile allows the calculation of physically meaningful perturbations.

In the spherical case, physically acceptable solutions form a continuum with growth rates σ> 1.

The nonspherical modes exhibit growth rates similar to the shell mode. This implies that it is the shell mode that drives the instability, with the higher order modes growing on top of it and providing the multiplicity of the fragments.
In parallel, we performed direct numerical simulations of a collapsing uniform isothermal sphere with a linear perturbation, varying the degree of thermal support (quantified by the ratio of thermal to gravitational energy, or virial parameter α) and the shape of the initial perturbation. For the case of smallest thermal support we tried (α = 0.006), the growth rates estimated from the simulations are in perfect accordance to the analytical expectations.
By increasing the degree of thermal support, we find that the density peaks grow slower and slower with respect to the background, while for larger values of the virial parameter the density contrast decreases with time. This suggests that eventually, nonrotating cores with 0.15 <α< 0.6 would collapse to a single object, while for larger values of α they would reexpand.
Our results suggest the existence of a “thermal” type of fragmentation for cold, large clouds that does not require the presence of initial rotation and that can produce very widely separated fragments at the early stages of the collapse.
Acknowledgments
We are grateful to Philippe André, Anaelle Maury, Tomoyuki Hanawa and Shuichiro Inutsuka for their comments and suggestions. We also thank the anonymous referee for useful comments that helped improve this manuscript. This work is largely a result of our participation in the International Summer Institute for Modeling in Astrophysics (ISIMA) during the summer of 2011. We especially thank Pascale Garaud for organizing the meeting and for useful discussions. This research has received funding from the European Research Council under the European Community’s Seventh Framework Programme (ERC Grant Agreement “ORISTARS”, No. 291294 and FP7/20072013 Grant Agreement No. 306483).
References
 Andre, P., WardThompson, D., & Barsony, M. 1993, ApJ, 406, 122 [NASA ADS] [CrossRef] [Google Scholar]
 Andre, P., WardThompson, D., & Barsony, M. 2000, in Protostars and Planets IV (Tucson: University of Arizona Press), 59 [Google Scholar]
 Banerjee, R., & Pudritz, R. E. 2006, ApJ, 641, 949 [NASA ADS] [CrossRef] [Google Scholar]
 Bate, M. R. 2009, MNRAS, 397, 232 [NASA ADS] [CrossRef] [Google Scholar]
 Bonnor, W. B. 1956, MNRAS, 116, 351 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P., & Bodenheimer, P. 1979, ApJ, 234, 289 [NASA ADS] [CrossRef] [Google Scholar]
 Boss, A. P., & Keiser, S. A. 2013, ApJ, 764, 136 [NASA ADS] [CrossRef] [Google Scholar]
 Burkert, A., & Bodenheimer, P. 1993, MNRAS, 264, 798 [NASA ADS] [CrossRef] [Google Scholar]
 Cha, S.H., & Whitworth, A. P. 2003, MNRAS, 340, 91 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Launhardt, R., Bourke, T. L., Henning, T., & Barnes, P. J. 2008, ApJ, 683, 862 [NASA ADS] [CrossRef] [Google Scholar]
 Commerçon, B., Hennebelle, P., Audit, E., Chabrier, G., & Teyssier, R. 2008, A&A, 482, 371 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duchêne, G., Bouvier, J., Bontemps, S., André, P., & Motte, F. 2004, A&A, 427, 651 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duchêne, G., Bontemps, S., Bouvier, J., et al. 2007, A&A, 476, 229 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Ebert, R. 1955, ZAp, 37, 217 [NASA ADS] [Google Scholar]
 Falgarone, E., & Phillips, T. G. 1990, ApJ, 359, 344 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Ghez, A. M., Neugebauer, G., & Matthews, K. 1993, AJ, 106, 2005 [NASA ADS] [CrossRef] [Google Scholar]
 Girart, J. M., Rao, R., & Estalella, R. 2009, ApJ, 694, 56 [NASA ADS] [CrossRef] [Google Scholar]
 Girichidis, P., Federrath, C., Banerjee, R., & Klessen, R. S. 2011, MNRAS, 413, 2741 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, A. A., Benson, P. J., Fuller, G. A., & Myers, P. C. 1993, ApJ, 406, 528 [NASA ADS] [CrossRef] [Google Scholar]
 Goodwin, S. P., Kroupa, P., Goodman, A., & Burkert, A. 2007, in Protostars and Planets V (Tucson: University of Arizona Press), 133 [Google Scholar]
 Hachisu, I., & Eriguchi, Y. 1984, A&A, 140, 259 [NASA ADS] [Google Scholar]
 Hanawa, T., & Matsumoto, T. 1999, ApJ, 521, 703 [NASA ADS] [CrossRef] [Google Scholar]
 Hennebelle, P., & Teyssier, R. 2008, A&A, 477, 25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hennebelle, P., Whitworth, A. P., Cha, S.H., & Goodwin, S. P. 2004, MNRAS, 348, 687 [NASA ADS] [CrossRef] [Google Scholar]
 Hunter, C. 1964, ApJ, 139, 570 [NASA ADS] [CrossRef] [Google Scholar]
 Jog, C. J. 2013, MNRAS, 434, L56 [NASA ADS] [CrossRef] [Google Scholar]
 Joos, M., Hennebelle, P., Ciardi, A., & Fromang, S. 2013, A&A, 554, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jørgensen, J. K., Bourke, T. L., Myers, P. C., et al. 2007, ApJ, 659, 479 [NASA ADS] [CrossRef] [Google Scholar]
 Klessen, R. S., & Burkert, A. 2001, ApJ, 549, 386 [NASA ADS] [CrossRef] [Google Scholar]
 Klessen, R. S., Burkert, A., & Bate, M. R. 1998, ApJ, 501, L205 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1969, MNRAS, 145, 271 [NASA ADS] [CrossRef] [Google Scholar]
 Larson, R. B. 1981, MNRAS, 194, 809 [NASA ADS] [CrossRef] [Google Scholar]
 Machida, M. N., Omukai, K., Matsumoto, T., & Inutsuka, S.I. 2008, ApJ, 677, 813 [NASA ADS] [CrossRef] [Google Scholar]
 Matsumoto, T., & Hanawa, T. 2003, ApJ, 595, 913 [NASA ADS] [CrossRef] [Google Scholar]
 Maury, A. J., André, P., Hennebelle, P., et al. 2010, A&A, 512, A40 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Miyama, S. M., Hayashi, C., & Narita, S. 1984, ApJ, 279, 621 [NASA ADS] [CrossRef] [Google Scholar]
 Myhill, E. A., & Kaula, W. M. 1992, ApJ, 386, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Offner, S. S. R., Klein, R. I., & McKee, C. F. 2008, ApJ, 686, 1174 [NASA ADS] [CrossRef] [Google Scholar]
 Ori, A., & Piran, T. 1988, MNRAS, 234, 821 [NASA ADS] [CrossRef] [Google Scholar]
 Patience, J., White, R. J., Ghez, A. M., et al. 2002, ApJ, 581, 654 [NASA ADS] [CrossRef] [Google Scholar]
 Penston, M. V. 1969, MNRAS, 144, 425 [NASA ADS] [CrossRef] [Google Scholar]
 Price, D. J., & Bate, M. R. 2007, MNRAS, 377, 77 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, F. H. 1977, ApJ, 214, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Teyssier, R. 2002, A&A, 385, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tohline, J. E. 1981, ApJ, 248, 717 [NASA ADS] [CrossRef] [Google Scholar]
 Tohline, J. E. 2002, ARA&A, 40, 349 [NASA ADS] [CrossRef] [Google Scholar]
 Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179 [NASA ADS] [CrossRef] [Google Scholar]
 Tsuribe, T., & Inutsuka, S.I. 1999, ApJ, 526, 307 [NASA ADS] [CrossRef] [Google Scholar]
 Whitworth, A., & Summers, D. 1985, MNRAS, 214, 1 [NASA ADS] [CrossRef] [Google Scholar]
All Figures
Fig. 1
Density (top panel) and velocity (bottom panel) fields of the solutions for various values of σ. For values of X smaller than the position of the shock (X<X_{sh} ≃ 3), the density and the velocity are unperturbed, i.e. δR = δV = 0. The solid lines show the solutions for the system of Eqs. (5)–(7) For comparison, solutions to the same equations but without the thermal pressure term are shown with the dotted lines. 

Open with DEXTER  
In the text 
Fig. 2
From top to bottom, density R, radial velocity U, azimuthal velocity V and gravitational potential φ for the l = 2 mode and for a series of growth rates. 

Open with DEXTER  
In the text 
Fig. 3
Same as Fig. 2 for the l = 3 mode. 

Open with DEXTER  
In the text 
Fig. 4
Same as Fig. 2 for the l = 4 mode. 

Open with DEXTER  
In the text 
Fig. 5
Contours of the log n_{H} in units of m_{H} cm^{3} along the yz plane. From top to bottom: l = 0, l = 1, l = 2 and l = 3 modes for virial parameter α = 0.006. The units along each axis are parsecs and t_{ff} stands for freefall time. 

Open with DEXTER  
In the text 
Fig. 6
Contours of the log n_{H} in units of m_{H} cm^{3} along the yz plane. From top to bottom: l = 2 mode for virial parameters α = 0.1, α = 0.2 and α = 0.8. Like in Fig. 5, units along each axis are parsecs and t_{ff} stands for freefall time. 

Open with DEXTER  
In the text 
Fig. 7
Growth rates for α = 0.006 mode for azimuthal wave numbers l. 

Open with DEXTER  
In the text 
Fig. 8
Growth rates of the l = 0 mode for various virial parameters α. 

Open with DEXTER  
In the text 
Fig. 9
Growth rates of the l = 2 mode for various virial parameters α. 

Open with DEXTER  
In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.