Issue 
A&A
Volume 578, June 2015



Article Number  A18  
Number of page(s)  23  
Section  Interstellar and circumstellar matter  
DOI  https://doi.org/10.1051/00046361/201424798  
Published online  25 May 2015 
Substellar fragmentation in selfgravitating fluids with a major phase transition^{⋆}
Geneva Observatory, University of Geneva, 1290 Sauverny, Switzerland
email: andreas.fueglistaler@unige.ch
Received: 12 August 2014
Accepted: 22 March 2015
Context. The observation of various ices in cold molecular clouds, the existence of ubiquitous substellar, cold H_{2} globules in planetary nebulae and supernova remnants, or the mere existence of comets suggest that the physics of very cold interstellar gas might be much richer than usually envisioned. At the extreme of low temperatures (≲10 K), H_{2} itself is subject to a phase transition crossing the entire cosmic gas density scale.
Aims. This wellknown, laboratorybased fact motivates us to study the ideal case of a cold neutral gaseous medium in interstellar conditions for which the bulk of the mass, instead of trace elements, is subject to a gasliquid or gassolid phase transition.
Methods. On the one hand, the equilibrium of general nonideal fluids is studied using the virial theorem and linear stability analysis. On the other hand, the nonlinear dynamics is studied using computer simulations to characterize the expected formation of solid bodies analogous to comets. The simulations are run with a stateoftheart molecular dynamics code (LAMMPS) using the LennardJones intermolecular potential. The longrange gravitational forces can be taken into account together with shortrange molecular forces with finite limited computational resources, using supermolecules, provided the right scaling is followed.
Results. The concept of supermolecule, where the phase transition conditions are preserved by the proper choice of the particle parameters, is tested with computer simulations, allowing us to correctly satisfy the Jeans instability criterion for onephase fluids. The simulations show that fluids presenting a phase transition are gravitationally unstable as well, independent of the strength of the gravitational potential, producing two distinct kinds of substellar bodies, those dominated by gravity (planetoids) and those dominated by molecular attractive force (comets).
Conclusions. Observations, formal analysis, and computer simulations suggest the possibility of the formation of substellar H_{2} clumps in cold molecular clouds due to the combination of phase transition and gravity. Fluids presenting a phase transition are gravitationally unstable, independent of the strength of the gravitational potential. Arbitrarily small H_{2} clumps may form even at relatively high temperatures up to 400–600 K, according to virial analysis. The combination of phase transition and gravity may be relevant for a wider range of astrophysical situations, such as protoplanetary disks.
Key words: ISM: clouds / ISM: kinematics and dynamics / ISM: molecules / methods: analytical / methods: numerical / instabilities
Figures 33−44 are available in electronic form at http://www.aanda.org
© ESO, 2015
1. Introduction
Typically around 50% or more of the gas in spiral galaxies consists of H_{2}, inferred indirectly by CO or dust emission. Since the discovery of dark molecular hydrogen (Grenier et al. 2005; Langer et al. 2010; Planck Collaboration XIX 2011; Paradis et al. 2012), the estimated quantity of H_{2} in the Milky Way has essentially been doubled, effectively revealing the nature of some of the dark baryons. It is commonly admitted that the COrelated molecular hydrogen is present in relatively dense regions of the interstellar medium, molecular clouds with number density >10^{10} m^{3} and temperatures of 7–30 K (Draine 2011).
Even though H_{2} is by far the most abundant molecule (~90%), molecular clouds are mainly detected by CO emissions because of all the difficulties in detecting cold H_{2} (Bolatto et al. 2013). For example, H_{2} only starts emitting at temperatures >512 K. See Combes & Pfenniger (1997) for a review of several possible methods for detecting cold H_{2}. Because of these detection difficulties, the real quantity of H_{2} (as well as He, which shares similar properties of discreteness with H_{2}) in molecular clouds is still rather unknown, especially when the gas temperature is below ≲8 K down to the cosmic background temperature of 2.76 K.
The condensation properties of H_{2} relevant for molecular cloud conditions are well known from laboratory data (Air Liquide 1976). The phase diagram (Fig. 1) shows the domain of pressure conventionally attributed to molecular clouds. One has to keep in mind, however, that since molecular clouds are highly structured (Pfenniger & Combes 1994), and commonly observed in a state of supersonic turbulence (Elmegreen & Scalo 2004), large fluctuations in density and temperature must occur.
The presence of ice in the interstellar medium consisting of heavier molecules, such as H_{2}O, CO, CO_{2} and NH_{3} covering dust grains, is nowadays well documented (Allamandola et al. 1999). Figure 2 shows the location of the critical and triple points of abundant molecules existing in the ISM. H_{2} ice has been detected in the absorption band at 2.417 μm (Sandford & Allamandola 1993; Buch & Devlin 1994; Dissly et al. 1994), but the interpretation of this detection is that H_{2} is mixed within H_{2}Orich grains in conditions that are too warm to allow the bulk of H_{2} to condense (Kristensen et al. 2011).
Fig. 1 H_{2} phase diagram (bold line) and He (dotted line) in cold and low pressure conditions. 

Open with DEXTER 
Highresolution pictures of nearby planetary nebulae or remnants of supernova have shown the presence of substellar fragments (Walsh & Meaburn 1993). These very cold globules, or knots, each of the size of a few tens of AU, are at least as cold in the inner parts (~10 K) as molecular clouds, but much smaller. An important feature is that the apparent column density in these knots increases inwards as long as the resolution allows, or until the knot is optically thin (Burkert & O’dell 1998). If these trends extend to the centre, one can expect there would be much higher density at colder conditions. It would be ideal to eventually reach a regime where H_{2} could condense in liquid or solid form, especially because at high column density the medium blocks UV radiation and cosmic ray heating.
At the level of molecular clouds, even though their average properties are well separated from the H_{2} phase transition, these are only static properties that ignore the highly dynamical nature of supersonic turbulence observed in the interstellar medium, where fluctuations must be large. Lower temperatures can be reached with fast adiabatic decompression alone. An example of fast decompression is displayed in the Boomerang Nebula, which reaches a temperature of only 1 K because of a fast expanding wind (Sahai & Nyman 1997). Thus, one should expect, to be coherent with supersonic turbulence, that regions of expansion and compression must be common in the ISM.
The fragmentation of selfgravitating gas leading to collapse has been studied for over a century. The molecular cloud is normally represented as a selfgravitating, onephase fluid (i.e. pure gas), which is governed by the balance between gravity and gas pressure. The gravity of this fluid is directly proportional to its density, and the pressure to the adiabatic density derivative of the pressure, (∂P/∂ρ)_{s}. It has been shown by Jeans (1902) that if a perturbation with a long enough wavelength is introduced into the system, this perturbation will grow exponentially. With growing perturbation, the pressure of the fluid will not be able to withstand its gravity anymore, which will ultimately lead to gravitational collapse.
The approximation of molecular clouds as a onephase fluid is valid in many cases, but when considering very cold, highdensity regions, it is an oversimplification. In these condition, H_{2} will start creating ice, which has to be taken into account (Walker 2013). The dynamics of this kind of fluid presenting a phase transition is different from a onephase fluid, the most important difference being a very low value of (∂P/∂ρ)_{s} (Johnston 2014). But, as (∂P/∂ρ)_{s} is crucial for the stability of a fluid, a cold highdensity fluid presenting a phase transition is therefore expected to be very unstable.
Based on the observation of substellar globules in planetary nebulas and the dynamics of fluids presenting a phase transition, we may expect the presence of small, substellar H_{2} ice fragments in molecular clouds due to the fragmentation of cold highdensity regions. It is of great astronomical interest to study the nature of this substellar fragmentation, as the resulting bound objects may be too small to start nuclear fusion and, thanks to their very low temperature, are very difficult to detect. Some of the baryonic dark matter may actually consist of these fragments (Pfenniger & Combes 1994; Pfenniger et al. 1994).
Fig. 2 Critical points (bullet) and triple points (Y) of common molecules in the ISM. As for H_{2} in Fig. 1, the sublimation curves cross interstellar pressure conditions very steeply in pressure a few K below their triple point. For example, CO should essentially be frozen below ~ 20 K, so unable to emit rotation lines. 

Open with DEXTER 
In this article, we study the physics of a selfgravitating van der Waals fluid presenting a phase transition analytically and with simulations. In the analytic part (Sect. 2), the physics of a van der Waals fluid and the related LennardJones (thereafter LJ) potential are recalled. We calculate the virial theorem taking both the gravitational and the LJ potential into account. The virial analysis helps to characterize different types of fluids. The stability of a selfgravitating van der Waals fluid presenting a phase transition is then analyzed.
In Sect. 3, the molecular dynamics simulator LAMMPS (Plimpton 1995) is introduced. By proper scaling of physical constants the Coulomb force solver is used to calculate the gravitational force, and the short rang molecular force is calculated with the LJ force. We introduce the concept of supermolecules to enable us to perform simulations with a total mass high enough for the fluid to be selfgravitating.
The simulations performed are discussed in Sect. 4. First, the correctness of the supermolecule approach is tested. Second, onephase fluids (i.e. pure gas) are used to test the Jeans criterion. Third, simulations close to a phase transitions are performed, studying the properties of nongravitating fluids and fluids with a gravitational potential above and below the Jeans criterion.
2. Physics of a fluid presenting a phase transition
In this section, we describe the physics of a selfgravitating fluid presenting a phase transition.
2.1. van der Waals equation
The van der Waals equation is a classical equation of state (EOS) of a pure fluid presenting a first order phase transition. This EOS describes the macroscopic behaviour of a fluid in which at microscopic level the molecules are strongly repulsive at short distances and beyond that weakly attractive over a limited range.
The van der Waals equation (van der Waals 1910; Johnston 2014) is a modification of the idealgas law, taking the finite size of molecules and intermolecular interactions into account. It links pressure, density, and temperature as follows: (1)with a [Pa m^{6}] and b [m^{3}] being constants characteristic of the fluid. It can also be expressed in a reduced, dimensionless form as, (2)where P_{r} = P/P_{c}, n_{r} = n/n_{c}, and T_{r} = T/T_{c}. The parameters P_{c}, T_{c} and n_{c} are the values of the thermodynamic critical point (Kondepudi & Prigogine 1998; Carey 1999).
Figure 3 shows the phase diagram for a fluid with T = 0.9T_{c}. One can distinguish the gaseous phase (dotted line) and the solid phase (solid line); in between (dashed line), the fluid presents a phase transition. There are states where (∂P/∂ρ)_{s}< 0, which is thermodynamically unstable. The parts of the dashed curve where (∂P/∂ρ)_{s}> 0 are metastable because a lower entropy state is reached when the fluid splits into condensed and gaseous components. The fraction of condensed over gaseous phase grows from 0 to 1 from left to right.
Fig. 3 van der Waals phase diagram for a fluid with T_{r} = 0.9. Dotted line: gas phase; dashed line: phase transition; solid line: solid phase. 

Open with DEXTER 
Fig. 4 van der Waals phase diagram for a fluid with a temperature (from bottom to top) of T_{r} = 0.2, T_{r} = 0.4, T_{r} = 0.6, T_{r} = 0.8, T_{r} = 1.0, and T_{r} = 1.2. Solid line: van der Waals EOS with Maxwell construct, dotted line: original van der Waals EOS. 

Open with DEXTER 
Fig. 5 H_{2} and van der Waals phase diagrams. Solid line: H_{2} laboratory data; dotted line: van der Waals vapour curve derived with the Maxwell construct. 

Open with DEXTER 
When the two phases coexist, the van der Waals EOS is replaced by a constant pressure marked by a horizontal line. This constant pressure level is determined by the Maxwell equal area construct (ClerkMaxwell 1875), demanding a total zero P·v work for an adiabatic cycle between the fully gaseous to the fully condensed state.
Figure 4 shows the phase diagrams for fluids with a temperature from T = 0.2T_{c} to T = 1.2T_{c}. The dotted line shows the original van der Waals EOS whereas the solid line displays the modified law using the Maxwell construct. There is no Maxwell construct for T ≥ T_{c} as (∂P/∂ρ)_{s} is always ≥0. Lekner (1982) provides a parametric solution to the van der Waals and Maxwell construct, using Δs as variable.
The phase equilibrium pressure is almost identical with the coexistence curve of laboratory data for H_{2} over a wide range of pressure (Fig. 5).
Fig. 6 LennardJones potential. 

Open with DEXTER 
2.1.1. LennardJones forces
The van der Waals model of phase transition fluid is convenient for a continuum fluid description, but fails to represent all the phenomena around the phase transition, which is the reason for correcting it with the above mentioned Maxwell construct. Even with the Maxwell construct, however, a local thermal equilibrium is still supposed to hold. In astrophysical contexts, often even these assumptions cannot be granted. For instance, in supersonic turbulent situations, ubiquitous in the interstellar medium, thermal equilibrium cannot be satisfied since mechanical energy propagates faster than thermal energy through pressure. Thus any method using quantities, like temperature or pressure, implicitly assumes that a local thermal equilibrium is established, which makes their use in the supersonic turbulent regime uncertain.
A much less demanding model of fluid is provided by molecular dynamics, where the simplest molecule interactions close to the van der Waals model in equilibrium situations is provided by the LennardJones (LJ) intermolecular potential (Jones 1924). No local or global thermal equilibirum is required. The LJ potential consists of an attractive longrange term and a repulsive shortrange term (Fig. 6), (3)with r_{σ} = r/σ. Its minimum value, located at r_{σ} = 2^{1 / 6}, is −ϵ/m.
2.2. Virial theorem for a LennardJones gravitating fluid
2.2.1. LagrangeJacobi identity
The virial theorem (Clausius 1870) describes the statistical equilibrium of a system of interacting particles or fluid systems, relating the kinetic and potential energies. In the case of a selfgravitating LJ fluid, the LJ potential and gravity combine as a total potential Φ = Φ_{G} + Φ_{LJ}.
The virial theorem for this new potential can be derived using the LagrangeJacobi identity path, by taking the second time derivative of the moment of the polar inertia , we find, (4)where the first righthand term is twice the kinetic energy E_{kin}, and the second one is the virial term. For a system near a statistical equilibrium, both sides of this equation should oscillate around 0, so the respective time averages should vanish.
The LJ potential is a sum of two homogeneous functions^{1}Φ_{LJ} = Φ_{LJ,a} + Φ_{LJ,r} of degree −6 and −12 respectively, while gravity is of degree −1. So we can use Euler’s theorem of homogeneous functions to express the virial term as a sum of potential energies multiplied by minus the homogeneous degree. The LagrangeJacobi identity becomes (5)
2.2.2. Homogeneous sphere
For homogeneous, finite mass spheres at a given temperature, the individual terms of the LagrangeJacobi identity read,
where the constants c_{r} and c_{a} in the LJ terms have been calculated by straight summation over the 1.5 × 10^{12} nearest molecules, as given in Table 1 for both the simple cubic (SC) and hexagonal close packed (HCP) or facecentred cubic (FCC) lattices.
Lattice constants.
The virial equation, obtained by setting d^{2}I/ dt^{2} = 0 in the LagrangeJacobi identity, contains terms proportional to M except the gravity term, which is proportional to M^{5 / 3}. Dividing the equation by M we can separate M^{2 / 3}, yielding, (10)which is indeed positive if the sum of the terms inside the brackets are also positive. However the term proportional to c_{a} is negative, and when large introduces the possibility that no positive solution for M^{2 / 3} exists. Solving the bracket terms equal to zero for , we have then a quadratic equation for this term, which gives the solutions for the densities ρ_{0} for which M vanishes, (11)The + and − solutions are real nonnegative when the term inside the square root above is nonnegative. The maximum temperature at which M can vanish is thus, (12)Since the critical temperature for a phase transition in the absence of gravity is about ϵ/k_{B}, we see that T_{max} can be substantially larger than the critical temperature. With the values given in Table 1 for the constants c_{a} and c_{r}, the maximum temperature below which gravity combined with molecular forces enhances fragmentation are T_{max} = 414.3 K and 626.8 K for the SC and HCP/FCC lattices, respectively, which are 11.38 and 17.22 larger than ϵ/k_{B}.
The corresponding density ρ_{f} at which fragmentation can occur at arbitrarily small mass is (13)which is of the order of the individual molecule density m/σ^{3}.
Fig. 7 Minimum mass of isothermal equilibrium curves for H_{2} homogeneous spheres. 

Open with DEXTER 
Fig. 8 Virial equilibrium of gravitating LJ homogeneous isothermal spheres with H_{2} molecules arranged on a HCP/FCC lattice. The sphere mass as a function of density is plotted for temperatures ranging from τ ≡ T/T_{max} = 10^{2} to 10^{2}. 

Open with DEXTER 
The LagrangeJacobi identity therefore allows us to predict a density ρ_{f} and a maximum temperature T_{max} below which a homogeneous sphere made of LJ molecules can find a gravitational equilibrium at an arbitrarily small mass: in a way, this provides a simple model, an explanation, for the reason of forming substellar ice clumps out of a cold selfgravitating gas. For H_{2} molecules arranged as SC or HCP/FCC lattice, we find ρ_{f} = 73.8 kg m^{3} and 69.2 kg m^{3}, respectively. These densities are of the order or slightly below the solid or liquid H_{2} density (≈80 kg m^{3}).
At temperatures T>T_{max} the isothermal equilibrium curves have a minimum mass when dM/dρ  _{T} = 0, which is expressed as (14)where τ = T/T_{max} and (15)Figure 7 shows the minimum mass as a function of the temperature for HCP/FCC H_{2} homogeneous spheres. It is rising very steeply at T ≳ 627 K, but rising much slower after ~1000 K. Interestingly the mass range includes all the masses below terrestrial planet masses for temperatures below H_{2} dissociation.
2.2.3. Condensed bodies
Different condensed and uncondensed bodies can be identified using the LagrangeJacobi identity, as summarized in Table 2.
Condensed and uncondensed bodies in a LJ fluid using the LagrangeJacobi identity.
Figure 8 shows some M(ρ,T) equilibrium curves of gravitating homogeneous sphere whose molecules are arranged on a HCP/FCC lattice, which is able to represent the most compact sphere lattice at high density. At low density, on the left of the diagram, the equilibrium is principally fixed by the attractive part of molecules and their temperature, and the exact lattice structure is irrelevant. This is the domain of uncondensed molecular gas.
On the right, there is an accumulation curve at approximately ρ ≈ 100 kg m^{3} representing purely condensed H_{2} bodies mainly balancing gravity with the repulsive part of the molecular interaction. The area between the dashed line, connecting the minimum of the isotherms with T>T_{max}, and the accumulation curve is the domain of rocky planetoids.
At temperatures below T_{max} = 626.8 K, the equilibrium curves plunge to arbitrarily small masses along two regimes: the left vertical asymptotes represent the limit gaseous bodies, and the right asymptotes the condensed (solid or liquid) bodies called comets. The area of the comets lies between the dotted line, connecting the elbows of the isotherms with T<T_{max}, and the dashed line.
At high temperature and low density, to the left of both the dotted and dashed line and above the isotherm lines, lie the gaseous planetoids. These bodies’ equilibrium is fixed principally by gravity and temperature.
Of course real astrophysical bodies are not homogeneous, in practice in the small mass regime one can expect bodies made of a mixture of condensed state in the core surrounded by a less dense gaseous atmosphere, which could be calculated by solving the hydrostatic equilibrium, see Pfenniger (2004), where some models of isothermal bodies made of H_{2} containing a solid core and a gaseous atmosphere have been discussed.
2.3. Gravitational instability
The linearized wave equation for the isentropic density perturbation ρ_{per} of a selfgravitating fluid of density ρ is the following (Jeans 1902; Weinberg 1972): (16)Its solution is a superposition of modes of the form exp(i(k·x−ωt)), where ω is the frequency and k = 2π/λ the wavenumber, and where λ is the wavelength. This leads to the instability condition (Jeans 1902), (17)When fulfilled the modes are realexponential, therefore unstable. The classical Jeans criterion reads, (18)Therefore the effective gravitational instability is directly dependent on the EOS of the medium.
2.3.1. Ideal gas
In most textbooks about Jeans instability only the ideal gas is discussed. In the case of an ideal monoatomic gas, the pressure derivative term is (19)with the adiabatic index γ = 5 / 3. This value is always positive, which leads to the classical, ideal monoatomic gas Jeans criterion for instability, (20)
Fig. 9 van der Waals EOS, including the Maxwell construct. 

Open with DEXTER 
2.3.2. Fluid with phase transition
When considering a fluid presenting a phase transition the term (∂P/∂ρ)_{s} may differ from that of an ideal gas. As seen before, in the case of the van der Waals EOS, the Maxwell construct must be used, which may considerably change the combined stability of the fluid (Sect. 2.1). In the presence of a phase transition, the EOS gradient is strongly modified, in particular (∂P/∂ρ)_{T} = 0. Using (21)with c_{P} and c_{v} both finite values, we find (∂P/∂ρ)_{s} = 0. Therefore, a selfgravitating fluid in a phase transition is also automatically gravitationally unstable.
Figure 9 shows the 3D representation of the EOS including the Maxwell construct. Clearly the curved triangular region, the phase transition region (on the left) has a very different gradient than the almost idealgas region at low density and high temperature, or the condensed phase region at high density. It is remarkable that a temperature drop by one order of magnitude from the critical temperature leads to a drop in the critical pressure by 14 orders of magnitude. For example for H_{2} at 3 K the critical pressure is 6.4 × 10^{9} Pa.
3. Molecular dynamics
For all simulations, the LargeScale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is used (Plimpton 1995). The LAMMPS software is a state of the art and widely used single and multiprocessor code in chemistry, material sciences, and related fields. Its abilities to quickly compute short and longrange forces and the possibility to use a multitimescale integrator make it a suitable tool to perform our simulations.
3.1. Potential solver
The LAMMPS code has a wide range of force fields. For the simulations in this article we use only the LJ and selfgravitational potentials.
Since the straight calculation cost for exact pairwise interactions is O(N^{2}), approximate but still accurate methods are implemented that make the calculations possible at a much lower cost. A cutoff radius r_{c} is used for the shortrange forces. As the attractive part of the LJ potential drops with r^{6}, it is calculated only for neighbour particles within r_{c}. The neighbours for each particle are found using a Verlet neighbour list. This list is created with a radius of r_{n} = r_{c} + r_{s}, r_{s} being an extra skin distance to avoid the recalculation of the Verlet list at every timestep. This cutoff method is O(N) and therefore linearly scales with the number of particles.
The gravitational potential drops with r^{1} so the long range interactions cannot be ignored. It is calculated using the Particle^{3}Mesh (P^{3}M) method (Hockney & Eastwood 1981). The gravitational potential is split in shortrange and longrange parts in Fourier space (Ewald 1921): (22)where r_{s} is the splitting distance. The potential Φ_{SR} is calculated at the same time as the LJ potential, using the same Verlet list. The potential Φ_{LR} is calculated in Fourier space using a fast Fourier transform (FFT).
As LAMMPS is designed for the simulation of chemical substances not far from terrestrial conditions, and selfgravity is too weak to be of any influence, no specific selfgravity module is provided. For that reason, a tweak is used by using the Coulomb’s potential module for an electric field, (23)where C is the interaction constant, ϵ the dielectric constant, and q_{i,j} the charges of the molecules. As C is fixed in LAMMPS, to correctly calculate the gravitational potential we set the constants as follows:
3.2. Time integration
The time integration is done using the symplectic leapfrog scheme. The drift and kick operators, are applied according to the sequence at each elementary timestep. For shortrange force, the timestep is constant throughout the simulation and set as 10^{2}–10^{3}, which is the typical interaction timescale during nearestneighbour molecular interactions. This is reasonable since the repulsive interaction and the finite kinetic energy prevent strongly varying accelerations between particles.
Small changes of individual particle positions can significantly change the shortrange forces, which is why a very small timestep is required. But these small changes have almost no influence on longrange forces. For this reason, there is no need to calculate the longrange forces at every timestep.
In regard to shortrange calculations, longrange calculations are an order of magnitude more costly per step but always need the same amount of calculation time: creation of the densitymap with a fifth order interpolation procedure, a FFT solution of the potential, and the same interpolation procedure for finding the accelerations.
With the reversible reference system propagator algorithm (rRESPA), a multiple timescale integrator is available for LAMMPS that enables us to reduce the number of longrange force calculation (Tuckerman et al. 1992). It enables time integration in up to four hierarchical levels, but only two are needed for the simulations we present. The rRESPA timeintegrationscheme looks as follows: (28)where n_{SR} is the number of shortrange iterations and set as 10–100, K_{LR} is the kick operator using the longrange (FFT) accelerations, and K_{SR} the kick operator using the shortrange accelerations.
3.3. Supermolecules
The maximum number of particles that can be simulated on today’s supercomputers is of the order of 10^{9}–10^{10} particles. Even using such a huge number of molecules would only result in a fluid with a total mass of a few femtograms. It is obvious that gravity has no effect on this kind of fluid. For that reason, the concept of supermolecules is introduced. This concept is well established in galactic dynamics simulation, where superstars weigh typically 10^{4}−10^{6}M_{⊙}, and in cosmology where superWIMPS may weigh as much as ~10^{67} GeVc^{2} particles. Twobody relaxation or diffusion due to the low number of superparticles is negligible, provided the simulation time is not too large, depending on the specific problem. Practice has shown that this kind of an approximation is valid in galactic dynamics for instance if the simulation length is of order 100 dynamical times (Binney & Tremaine 2008).
The basic principle of the concept is that each supermolecule consists of η molecules, and its mass is therefore, (29)To have the same properties of a supermolecule fluid as for a normal fluid, every term of the virial Eq. (5) has to be transformed such that the ratios of the terms are invariant. The kinetic energy term reads (30)Since the total mass is constant, N_{M}m_{M} = N_{SM}m_{SM}, the velocity dispersion of molecules and supermolecules is also the same, (31)From the two LJ terms, Eqs. (7) and (8) remain the same if Solving these two equations and using Eq. (29), we derive The gravitational energy Eq. (9) does not change since it does not depend on supermolecule properties.
It is important to ensure that the gravitational force between two supermolecules remains small compared to the corresponding LJ force within the shortrange force cutoff radius. This is assured by setting F_{G}(r_{c}) ≪ F_{LJ}(r_{c}) with r_{c} being the cutoff radius. This leads to the following constraint: (36)where r_{c, SM} = (r_{c}/σ_{SM}). Using a cutoff radius of r_{c} = 4 σ and hydrogen supermolecules, we find a maximum supermolecule mass equal to 5.7 × 10^{6}M_{⊕}, meaning at least 1.7 × 10^{5} particles are required for simulating an Earthmass, virialized, H_{2}condensed body.
Fig. 10 Mean kinetic energy as a function of the maximum displacement in an ice clump. 

Open with DEXTER 
3.4. Ice clump detection
The LJ potential has its minimum at r_{m} = 2^{1 / 6}σ. If an ice clump is at absolute zero temperature, all bound molecules would have this distance from at least one other molecule. But, as the molecules in a clump are vibrating, the maximum binding distance between two molecules is generally larger than r_{m}.
The simplest clump consists of two molecules. Their mean kinetic energy as a function of the maximum displacement can be calculated numerically (Fig. 10).
The mean kinetic energy rises up to the maximum value of r_{σ,max} = 1.3625, after which it drops again. It can be assumed that all stable clumps have a binding distance below this value. The distance constraint for two molecules to be bound is therefore: (37)Using r_{σ,max} as the threshold distance, LAMMPS provides a list of clumps at fixed time intervals.
Fig. 11 LJ potential energy as a function of the number of supermolecules at t = 1τ. Solid line: SM15; dashed line: SM30; dotted line: SM60. See Table 3. 

Open with DEXTER 
Parameters of the supermolecule test simulations.
4. Simulations
4.1. Units
A fluid is defined by the number of supermolecules N_{SM}, the initial velocity distribution (temperature), the number density, and the strength of the gravitational potential. All simulation properties are moleculeindependent, the initial temperature and density are measured as a factor of the critical values T_{cr} and n_{cr}, and the gravitational potential is measured as the factor γ_{J} = G/G_{J} of the idealgas Jeans gravity G_{J}, defined as (38)with the box side L = (N_{SM}/n)^{1 / 3}. It is interesting to note that G_{J} is proportional to temperature and inversely proportional to density. For all simulations, the time unit is defined as the average time of a particle crossing the box (39)with .
The distance at which the gravitational and LJ forces are equal is denoted as (40)
4.2. Supermolecule concept tests
4.2.1. Potential energy
A fluid in a cubic box with periodic boundary conditions and initial constant density of n = 0.006 n_{cr} is simulated with LAMMPS. Three different initial Maxwellian velocity distributions with temperature T = 1.5, 3.0 and 6.0 T_{cr}, and a number of particles from N_{SM} = 30^{3}–200^{3} are used. Table 3 shows the parameters of the different simulations.
Figure 11 shows the LJ potential energy of those three fluids with respect to the number of supermolecules. One can see fluctuations in the simulations with low N_{SM}, but the result becomes reasonable when N_{SM}> 10^{4} and a satisfactory convergence is obtained if N_{SM}> 10^{5}. Therefore all the subsequent simulations are performed with N_{SM}> 10^{5}.
Fig. 12 Fraction of bound molecules as a function of the number of supermolecules at t = 5τ. Solid line: SM04^{a}, dashed line: SM06^{a}. ^{(a)} See Table 3. 

Open with DEXTER 
4.2.2. Cluster percentage
High density, low temperature fluids are simulated during 5τ at various number of supermolecules. These fluids will form comets and allow us to compare the final percentage of bound molecules. Table 3 shows the parameters of the different simulations.
Figure 12 shows the fraction of bound molecules with respect to the number of supermolecules. All simulations reach the same percentage, independent of the number of supermolecules.
4.2.3. Precipitation in an external field
An external acceleration is applied in the z direction to a fluid within a box with reflecting boundary conditions. The fluid is very cold and rather dense (see Table 3) and is therefore forming ice grains, or comets. Because of acceleration and Archimedes’s principle, the comets fall to the bottom and stay there, similar to snowfall on Earth.
The purpose of this simulation is to test the timescaling of the precipitation as function of the number of particles. Figure 13 shows the z −coordinate of the centre of mass of comets as a function of time. The curves are very similar, with a slight overestimation of the collapse time for the smallest simulation. Therefore the precipitation phenomenon timescale is not strongly dependent on the mass resolution.
Fig. 13 zCoordinate of the centre of mass of comets as a function of time of the SF simulation. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). See Table 3. 

Open with DEXTER 
4.3. Fluid with perturbation
To illustrate the behaviour of a homogeneous fluid perturbed by a plane sinusoidal wave, a velocity perturbation in the x direction is introduced with wavelength λ equal to the box side L. Starting with a Maxwellian distribution v_{i} of velocities for the homogeneous unperturbed case, the xcomponent of the velocities is perturbed in such a way as to conserve energy. The perturbed xvelocity component for each particle i reads, (41)with ω = 2π/L. The correction factor α ≤ 1 is used to keep the same kinetic energy in the x direction, with (42)β determines the strength of the perturbation, but is supposed to be small enough to be in the linear regime of perturbations. In the present simulations, β = 0.01.
Parameters of the onephase fluid simulations.
4.4. Onephase fluid
To study the reaction of a pure onephase idealgas fluid far from the phase transition but with the above planewave perturbation, the initial temperature is taken above the critical value (T_{0} = 1.25 T_{cr}) and the initial density well below the critical value (n_{0} = 10^{2}n_{cr}). Three different cases are studied: without gravity, weakly selfgravitating below the classical idealgas Jeans criterion with γ_{J} = 0.5, and sufficiently selfgravitating above the idealgas Jeans criterion with γ_{J} = 1.25. The simulations are performed with N_{SM} ranging from 50^{3} to 160^{3}. Table 4 shows the parameters of the different simulations.
Fig. 14 Temperature of the onephase fluid simulations OP0 (straight solid line) and OPGs as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). See Table 4. 

Open with DEXTER 
Fig. 15 Fraction of bound molecules in comets of the onephase fluids OP0 (straight solid line) and OPGs as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER 
4.4.1. Time evolution
As the simulations conserve total energy, the formation of comets, implying a decrease of potential energy, can be measured by the mean temperature (i.e., kinetic energy) of the system. This is equivalent to the release of latent heat when a gaseous fluid condenses.
Figures 14 and 15 display the temperature and the comet percentage evolutions of the onephase fluid simulations: OP0 without gravity and N_{SM} = 100^{3} supermolecules and the sufficiently selfgravitating fluid OPGs with N_{SM} ranging from 50^{3} to 160^{3}. The weakly selfgravitating fluid of simulation OPGw is not shown since it is identical to the no gravity case OP0.
For the fluid without gravity OP0 and the weakly selfgravitating fluid OPGw, no particular effect is observed. The percentage of comets rises to ≲2%, which can be attributed to the initial fluctuations in the distribution. The small perturbation introduced in the xdirection does not change the nature of the fluid; its temperature and comet percentage remains the same.
The sufficiently selfgravitating OPGs fluids, on the other hand, do change. Their temperatures and cluster percentages are rising steeply showing a substantial latent heat release. All simulations reach the same asymptotic temperature ~6T_{cr}, but different initial conditions (i.e., random seeds) yield slightly different behaviours during the collapse. This can also be observed in Fig. 16 where four identical parameter simulations but different random seeds are run with N_{SM} = 100^{3}, leading to a time difference of ~0.5τ for reaching the asymptotic upper value.
Fig. 16 Temperature of the onephase fluid simulations OPGs as a function of time. N_{SM} = 100^{3} with four different random seeds. 

Open with DEXTER 
Fig. 17 Density at x = (0.5 ± 0.05)L of the onephase fluids OP0 (straight solid line) and OPGs as a function of time. Bold line: analytic solution for ideal gas. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER 
Figure 17 shows the density increase around the centre of the perturbation in the range x = (0.5 ± 0.05)L and compares them to the analytic solution for an ideal gas. The simulations follow the idealgas solution with a slight dispersion up to t ≈ 2.5τ, where a density rise is not possible anymore because of the repulsive LJ force; the density declines for a short while because of the collapse rebound down to a minimum at t ≈ 2.7τ. This density local minimum corresponds with the break in the temperature increase visible in Fig. 14. Subsequently, the temperature growth is much slower.
Figure 18 shows the evolution of density in the simulations, and the predicted growth rate line for an ideal gas subject to Jeans’ instability (Sect. 2.3). All curve slopes correspond initially to the predicted Jeans’ growth rate. The initial steep growth and the small time delay are caused by the initial noise in the particle distribution. Thus, a certain time is needed for the perturbation to overcome the noise.
Fig. 18 Evolution of perturbation density. Straight dotted line: analytic solution for ideal gas. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER 
Fig. 19 3Dview of simulation OPGs with N_{SM} = 100^{3} at t = 4τ. 

Open with DEXTER 
Fig. 20 Distribution of y and zcoordinates of planetoid centre of mass on a 10 × 10 grid; all OPGs simulations at t = 2.5τ. 

Open with DEXTER 
Fig. 21 Temperature distribution of unbound molecules and comets as a function of comet mass of OPGs with N_{SM} = 100^{3} at t = 4τ. 

Open with DEXTER 
4.4.2. Comet mass distribution
Figure 33 shows snapshots and Fig. 34 the comet mass distribution at different times. During the temperature rise, the comets are distributed according to a power law, (43)and no planetoid is formed. The powerlaw index is very steep at the beginning with ξ_{c} ≈ −10, but quickly decreases as it reaches a value of ξ_{c} ≈ −2.5 at the end.
After the temperature increase break at t ≈ 2.7τ, the powerlaw distribution remains for small comets, but one bigger gaseous planetoid with over 1% of the total mass forms. It is shown in Fig. 19. The spherical body seen at t ≤ 3τ in the snapshots is uncondensed gas, which is why it does not figure in the comet mass distribution (see Sect. 3.4 for a discussion how condensed matter is identified).
A second power law: (44)on the right side of the diagram describes the comets and planetoids. Each body occurs typically only once with our limited N_{SM}, and ξ_{p} = 1. Equation (43) describes the mass distribution of small bodies where gravity is negligible whereas Eq. (44) describes large bodies where gravity is weak, but cannot be neglected.
The first highdensity plane is created thanks to the plane perturbation parallel to the yzplane at x = 0.5. One can observe in Fig. 33 at t = 2τ that two filaments form, along the y and zaxes, connecting the planetoid with itself thanks to the periodic boundary conditions. They are aligned with the mesh, but not primarily due to grid effects, but because these lines are the shortest distance between the planetoid and its periodic images. As can be seen in Fig. 39, diagonal filaments are also possible using the second shortest distance.
The position of the planetoid is, thanks to the planewave perturbation, situated in the x = 0.5L plane, the y and z positions are random and depend on the initial condition and vary with every simulation as shown in Fig. 20.
Figure 21 shows the temperature distribution in the comets as function of mass. As the number of large comets with the same number of supermolecules is small or even one, no error bars can be seen for the larger comets and the gaseous planetoid.
The broad Maxwellian velocity distribution is seen for the unbound molecules (M_{comet} = 10^{6}). The small comets (2 or more particles) clearly have a lower temperature and dispersion than the average temperature of the system. This is due to the fact that in order for two supermolecules to cluster together, a third one is needed, and this takes away some of their kinetic energy. The larger a comet is, the closer its temperature is to the system average temperature. On longer formation time the planetoids temperature tends towards the system average temperature, as expected.
Fig. 22 Mean temperature over four onephase fluid simulations OPGs with different random seeds as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER 
4.4.3. Scaling
Figure 22 shows the mean temperature value over four simulations with different random seeds. One can see that the two smallest simulations with N_{SM} = 50^{3} and 64^{3} are starting to collapse earlier than the other simulations, which are very similar. This can be explained by the fact that these simulations have very low x_{GLJ} values (Eq. (40)), 3.65 and 4.03 r_{σ}. The first value is in fact lower than the cutoff radius of 4 r_{σ} and fails to satisfy Eq. (36), the other one just slightly above it. In these simulations, the gravitational forces are strong even in shortrange intermolecular interactions, the simulations are therefore overemphasizing the gravitational effect.
Figure 35 shows the comet mass distribution of all OPGs simulations. For small comets, the percentage for a same number of supermolecules is the same for all N_{SM}. For example, at t = 5τ, all simulations have a fraction of ~10^{1} for comets consisting of two supermolecules. These small clumps should be called multimers as mentioned below. But, with increasing N_{SM}, the mass of a comet consisting of two supermolecules diminishes since M_{2 SM} = 2 /N_{SM}.
On the other hand, for large comets and the planetoid, the mass is invariant to N_{SM}. For example, the planetoid at t = 5τ has the same mass of M_{planetoid} ≈ 5 × 10^{2}M_{tot} for all N_{SM}.
The fact that the mass distributions are shifted to the left for larger N_{SM} is misleading, as shown in Fig. 36. Only the two largest simulations are shown, the simulation with N_{SM} = 128^{3} is shown normally, but for the simulation with N_{SM} = 160^{3} is downscaled to 160^{3}/ 2 ≈ 128^{3}, always two comet sizes are added together. As can be seen, when comparing the two simulation using the same scaling, the two mass distributions are in fact very similar.
It is interesting to look at the turning point, defined as the point where the mass sum of the biggest comets is larger than the mass sum of smaller comets. This point is interesting as it is the first indicator of the creation of a aggregate of molecules that start to be influenced by gravity. Figure 23 shows the time and the mass of the turning point. One can see than the time of appearance is independent of N_{SM}. The comet mass, on the other hand, clearly follows a power law with: (45)with ξ_{N} ≈ −1. Since N_{SM} = M_{tot}/M_{SM}, at turning point the number of supermolecules is always of order of 10. In other words, in critical conditions where the attractive molecular force adds up to gravity, the tendency to make a larger condensed body mass fraction already starts at about ten molecules.
Fig. 23 Time and size of the largest comet at the appearance of the turning point. 

Open with DEXTER 
4.4.4. Extrapolation to physical scale
Having studied the scaling behaviour of the simulations, we can attempt to extrapolate some properties for a fluid consisting of real molecules instead of supermolecules. Considering H_{2} molecules, the total mass of the OPGs simulations is 1.4 M_{⊕}, the box length equal to 2.3 × 10^{3} AU and τ = 1.5 × 10^{2} years.
Thus, in realistic conditions there would be a total amount of 2.5 × 10^{51} H_{2} molecules, the number of molecules per supermolecule η varying from 2.0 × 10^{46} for N_{SM} = 50^{3} to 3.2 × 10^{44} for N_{SM} = 160^{3}. So, even for the largest simulations, there is a difference in number of particles of more than 40 orders of magnitude with real situations. For that reason, the extrapolated data has to be treated with caution. As long as no other physics enters in the interval of scales, this is however not extraordinary in regards to common practices in simulation works, such as simulating stars with SPH particles. While a single SPH particle can at best represent a mass fraction of 10^{6}−10^{10} of the total, an SPH particle is supposed to behave as the smallest mass element in local thermal equilibrium, say a few 1000 protons over the stellar mass: ~10^{55}, which lies therefore well over 10^{40} times the SPH simulation resolving capacity.
As can be observed in Fig. 15, the fraction of bound molecules does not change when increasing the number of particles and will remain unchanged for a fluid consisting of molecules. The same is true for the size of the resulting planetoid, as can be seen in Fig. 35. Therefore, we can assume that after 5τ, the system will consist of ~80% unbound molecules and a planetoid with a mass of ~0.07 M_{⊕}.
Fig. 24 Extrapolation of the comet mass distribution for a real H_{2} molecular fluid at t = 5τ. 

Open with DEXTER 
At t = 5τ, there is a fraction of 10^{1} of two molecules aggregates i.e. dimers with a mass of 6.6 × 10^{27} kg. These multimolecules clumps should be called multimers instead of comets. As mentioned above the turning point occurs at ~2.4τ consistently for approximately 10 molecules, so its mass is about 20 m_{H} = 3.3 × 10^{26} kg.
Using the above values and the two power laws of Eqs. (43) and (44), Fig. 24 shows a schematic diagram showing the mass distribution at t = 5τ. While the mass fraction of the turning point multimers is tiny, the powerlaw distribution rapidly increases the mass in larger grain and cometsized bodies until they become a planetoid mass.
Fig. 25 Position of the LJ simulations in a corresponding van der Waals phase diagram. Solid line: Maxwell construct for the van der Waals EOS; dotted line: van der Waals EOS. Star: fluid presenting a phase transition; bullet: onephase fluid. 

Open with DEXTER 
4.5. Phase transition
We study several physical conditions close to a phase transition. Table 5 summarizes their properties. Three different densities, 10^{1}, 10^{2} and 10^{3}n_{cr}, and three different temperatures, 0.1, 0.2 and 0.3 T_{cr}, are chosen. As the onephase fluid studied in Sect. 4.4 also has a density of 10^{2}n_{cr}, three additional temperatures, 0.5, 0.7 and 0.9 T_{cr}, were used for this density to make a link between the onephase fluid and the fluids studied in this section.
Looking at the phase transition diagram (Fig. 25), all fluids with T ≤ 0.3T_{c} are on the Maxwell line. Those with ρ ≤ 10^{2}ρ_{c} are on the (∂P/∂ρ)_{s}> 0 part of the van der Waals phase diagram, which corresponds to metastable states on the verge of a deep phase transition. The fluids with ρ = 10^{1}ρ_{c} are on the (∂P/∂ρ)_{s} ≤ 0 part of the van der Waals phase diagram, which corresponds to unstable states in a phase transition. The fluids with T ≥ 0.5T_{c} are still in the stable regime.
Three different cases are studied: without gravity, strongly selfgravitating above the Jeans criterion, and weakly selfgravitating below the Jeans criterion.
Fig. 26 Temperature of nongravitating fluids as a function of time. Dotted line: n = 10^{1}n_{cr}; solid line: n = 10^{2}n_{cr}; dashed line: n = 10^{3}n_{cr}. All simulations with N_{SM} = 50^{3}. 

Open with DEXTER 
Parameters of the phase transition simulations.
Fig. 27 Closeup 3Dview of a medium sized comet and smaller aggregates. 

Open with DEXTER 
Fig. 28 Temperature distribution of unbound molecules and comets as a function of comet size of PT22^{a} without gravity and N_{SM} = 50^{3} at t = 100τ. ^{(a)} See Table 5. 

Open with DEXTER 
4.5.1. Fluid without gravity
To study the evolution of fluids presenting a phase transition without gravity, the x_{GLJ} value does not have to be considered and the number of supermolecules is set to N_{SM} = 50^{3}. Figure 26 shows the temperature evolution of some of these fluids. In the beginning, the molecules merge into small comets, which leads to a decrease of the potential energy and therefore a rise of the temperature. Once the temperature is high enough, the kinetic and potential energy reach a equilibrium and the fluid remains stable.
The number of comets formed is dependent on the density and inversely dependent on the temperature. The higher the number of formed comets, the longer it takes for the fluid to reach a stable regime. For example, the very cold fluids PT12 and PT22 only reach an asymptotic value around t ≈ 100τ.
The comets remain at moderate size as can be seen in the snapshots (Fig. 37) and in their mass distribution (Fig. 38): the number of supermolecules per comets is <10^{2}N_{SM}. The mass distribution follows the power law of Eq. (43), but at ~10^{4}, the mass distribution rises again, presenting a big number of comets between 10^{4} to 10^{2} supermolecules. Figure 27 shows a closeup 3D view of a mediumsized comet and smaller aggregates.
Figure 28 shows the temperature distribution as a function of the comet mass. As in simulation OPGs, the temperature for small bodies lies below average, but reaches average value for more massive bodies.
Fig. 29 Temperature of sufficiently selfgravitating fluids as a function of time. Dotted line: OPGs^{a}; dashed line: PT13^{a}; solid line (from left to right): PT32^{a}, PT52^{a}, PT72^{a}, PT92^{a}; dashdotted line: PT11^{a}. All simulations with N_{SM} = 100^{3}. ^{(a)} See Table 5 

Open with DEXTER 
4.5.2. Selfgravitating fluid above Jeans criterion
All sufficiently selfgravitating simulations use the same gravitational potential, G = 1.25 G_{J}(T = 1.25 T_{cr},n = 0.01 n_{cr}). As the temperatures and densities differ, G_{J} and therefore γ_{J} is different for each simulation, but this parameter is always >1.
As γ_{J}> 1 for all simulated fluids, they are unstable if perturbed. The lower the initial temperature of a fluid, the bigger is γ_{J}, which translates to a faster exponential growth. This can be seen in Fig. 29 where the fluids with a density of 10^{2}n_{cr} are starting to collapse one after another, from the lowest to the highest temperature.
The evolution of the fluids is identical to the fluid OPGs: an exponential rise of the temperature until it reaches a temperature break, after which the temperature rises only slowly and a single gaseous planetoid is formed. Figures 39 and 40 show snapshots and comet mass distribution, again very similar to the fluid OPGs. First, there is a formation of small comets following the power law of Eq. (43), and then the comets are massive enough to attract each other and merge into one big spherical planetoid (t> 1τ), following the power law of Eq. (44). The powerlaw index ξ_{c} also evolves in a similar fashion as in simulation OPGs, starting out very steep and reaching a final value of ≳−2.
Fig. 30 Temperature of weakly selfgravitating fluids as a function of time. Dotted line: PT22 and PT32, with N_{SM} = 50^{3} and without gravity; dashed line: PT22, γ_{J} = 0.5 and N_{SM} = 80^{3} and 100^{3}; solid line PT32, γ_{J} = 0.5 and N_{SM} = 80^{3} and 100^{3}; dashdotted line PT52, γ_{J} = 0.5 and N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 31 Temperature distribution of unbound molecules and comets as a function of comet size of PT32 with γ_{J} = 0.5 and N_{SM} = 100^{3} at t = 200τ. 

Open with DEXTER 
4.5.3. Selfgravitating fluid below Jeans criterion
To compare weakly selfgravitating fluids in a phase transition (on the Maxwell line) with outofphasetransition fluids (below the Maxwell line), γ_{J} is set to 0.5 for all fluids. This means that the individual G value is different for each simulation.
As shown in Sect. 4.5.1, the perturbation of the weakly selfgravitating onephase fluid OPGw does not create any effect if γ_{J}< 1. Figure 30 shows the temperature as a function of time of the selfgravitating fluids PT52, PT32, and PT22, and compares them to the nongravitating fluids. The outofphasetransition fluid PT52 remains unchanged and does not differ from the nongravitating fluid, identical to OPGw.
The temperatures of the fluids in a phase transition PT32 and PT22 are exponentially growing and differ clearly from the nongravitating ones. In the beginning, the LJ forces dominate the gravitational forces and the fluids behave like nongravitating fluids. PT22 reacts much faster than PT32, thanks to its low initial temperature and fluctuations due to the Maxwellian velocity distribution,and forms many small and mediumsized comets (see Sect. 4.5.1). Once these comets are formed, their mass is high enough to attract each other with gravity, forming a single “rocky planetoid”.
This can be seen in the snapshots (Fig. 41) and comet mass distribution (Fig. 42). During the first 15τ, the fluid is identical to the lowtemperature, nongravitating fluid (compare to Fig. 38). Only then does gravity start to play a role, one can see the formation of a planetoid and how it swallows the small and mediumsized comets during its growth.
With its higher temperature, PT32 forms almost no small or mediumsized comets, which can be seen in the snapshots (Fig. 43) and comet mass distribution (Fig. 44) for the first ≈50τ. Only after 50τ does gravity show its effect with the appearance of a mediumsized comet, which is too massive to fit into the power law ξ_{1} (see Fig. 43, t = 50τ). From then on, this comet attracts supermolecules and thus grows in size until at t ≈ 125τ, a single big rocky planetoid is formed.
Figure 31 shows the temperature distribution as a function of the comet mass of PT32. The temperature of the comets containing few supermolecules lies below the average, but the planetoids temperature is the same as the average.
PT22, PT32 and PT52 all scale very well. The two fluids in a phase transition, PT22 and PT32, have an almost identical timescale for the exponential growth for N_{SM} = 80^{3} and 100^{3}. The slight difference is due to the initial random seed. As x_{GLJ} depends on γ_{J} and the temperature, both much lower than for OPGs, its value is clearly above the cutoff radius, which is 7.63 and 8.34 for PT22 and 7.03 and 7.69 for PT32.
PT52 shows no difference for any number of supermolecules, with N_{SM} ranging from 50^{3} to 100^{3}.
Fig. 32 Density of planetoid as a function of the radius. Dashed line: OPGs with N_{SM} = 100^{3} at t = 4τ, the gaseous planetoid consists of 41559 supermolecules (0.042 M_{tot}). Dashdotted line: PT32, with γ_{J} = 0.5 and N_{SM} = 100^{3} at t = 200τ; the “rocky planetoid” consists of 74208 supermolecules (0.074 M_{tot}). Solid lines: hydrogen, higher value: solid; lower value: liquid. Dotted line: ρ_{f} of LagrangeJacobi identity, higher value: SC; lower value: HCP/FCC. 

Open with DEXTER 
4.6. Planetoid densities
Figure 32 shows the densities of the planetoids as a function the radius of the simulations OPGs and PT32 with γ = 0.5γ_{J}, and compares them to hydrogen laboratory data and values found in the LagrangeJacobi identity (Sect. 2.2).
The density of the OPGs gaseous planetoid is below that of a conventional solid or liquid. As their temperature is above the critical value (T> 6 T_{cr}), supermolecules are not able to condense without the aid of gravity, and because of their high kinetic energy, are vibrating at much higher amplitudes. This leads to more space between the supermolecules and therefore a lower density compared to supermolecules below the critical temperature.
The density of the gaseous planetoid drops with increasing radius and hence does not qualify as a homogeneous sphere, but its average density is below the ρ_{f} value of the LagrangeJacobi identity (Eq. (13)), in accordance with Table 2 and Fig. 8 for gaseous planetoids.
The density of the PT32 rocky planetoid lies between the liquid and solid phase. There is no continuous density drop as for the OPGs simulation, the density remains stable up to almost the outer radius and the body can be approximated as a homogeneous sphere. Its density lies above ρ_{f} in accordance with Table 2 and Fig. 8 for rocky planetoids.
5. Conclusions
We used analytic methods (Sect. 2) and computer simulations (Sects. 3 and 4) to study substellar fragmentation of fluids presenting a phase transition. The motivating astrophysical context are molecular clouds where H_{2} forms the bulk of the gravitating mass and is not very far from condensation conditions.
5.1. Analytic results
The study of the virial theorem, using the gravitational and the LennardJones potential energies, has shown that there is a maximum temperature below which a fluid can fragment at arbitrary small masses. This temperature is an order of magnitude above the critical temperature. This shows, granted the right circumstances, that comets can form at a temperatures an order of magnitude larger than the critical temperature. In the case of H_{2}, this maximum temperature lies in the range of 400–600 K, depending on the solid crystalline structure.
A van der Waals fluid can be in three different states: gaseous, solid/liquid, or in a phase transition where the two phases coexist. The latter is defined as lying on the line of the Maxwell construct, where (∂P/∂ρ)_{s} = 0. A fluid is gravitationally unstable if an introduced perturbation has a wavelength above a certain value, which depends on (∂P/∂ρ)_{s}. Since this quantity vanishes for a fluid presenting a phase transition, such a fluid is also gravitationally unstable.
5.2. Simulation results
We performed simulations using the stateoftheart molecular dynamics simulator LAMMPS. We used supermolecules to simulate gravitational and molecular forces together with a computationally tractable number of particles.
5.2.1. Supermolecules
We tested the supermolecule concept thoroughly. We achieved good scaling for nongravitational effects if the number of supermolecules is large enough (≳10^{5}). The sticking point is to ensure that the molecular LJ forces are dominating the gravitational forces in closerange interactions. This is achieved by setting the number of supermolecules to be large enough so that the distance at which the two forces are equal is above the cutoff radius 4 σ_{SM}. For an H_{2} fluid this sets a maximum supermolecule mass equal to 5.7 × 10^{6}M_{⊕}, limiting the total mass that can be simulated by the available computing power.
In principle, the supermolecule concept should not perfectly reproduce the time dependence of diffusive properties, as bigger particles introduce faster relaxation. But our experiments using different resolutions spanning several orders of magnitude did not reveal important modifications in the timings of major collapse and asymptotic evolution state. This enables us to study the fragmentation of large bodies, which we call comets and planetoids.
5.2.2. Onephase fluid
We applied a plane sinusoidal perturbation into a onephase fluid with a temperature above the critical value. The results reproduce the idealgas Jeans instability: no collapse is seen for conditions below the Jeans criterion, while an exponential growth of the perturbation is observed for conditions above it.
As a result, the temperature rises, and small and mediumsized comets form, some of which later merge into one big planetoid. An interesting observation is the mass distribution of these comets, which follows a power law for the small and mediumsized comets, while the planetoid and the largest comets follow a different power law.
5.2.3. Phase transition fluid
We simulated fluids with temperatures below the critical value and close to the effective phase transition for three different cases: without gravity, sufficiently, and weakly selfgravitating (above and below the idealgas Jeans criterion).
Because of the Maxwellian velocity distribution fluctuations, the nongravitating fluids form small and mediumsized comets until the potential and kinetic energy reach an equilibrium; thereafter they remain statistically stable. In the absence of any longrange force, no planetoid forms.
The selfgravitating fluids with a gravitational potential above the idealgas Jeans criterion do not react differently to a plane sinusoidal perturbation from a onephase fluid: the perturbation grows exponentially and its temperature rises, and small and mediumsized comets form, which ultimately merge into one planetoid.
The analysis predicts that fluids presenting a phase transition are unstable even if the gravitational potential alone is below the idealgas Jeans threshold. The performed simulations of weakly selfgravitating fluids in the phase transition regime did reproduce the prediction. Because of the weak nature of the gravitational potential, the timescale to see this reaction is two orders of magnitude longer than for sufficiently selfgravitating fluids, but, in the end, a big central planetoid forms anyway. While the outofphase transition fluids do not amplify perturbations, those in a phase transition indeed amplify them.
This study is general enough to be relevant for many astrophysical situations. In the case of H_{2} as main component, it concerns situations where temperature may drop well below 10 K, such as in cold molecular gas in the outer galactic disks, in expanding winds of planetary nebulae or supernovae, where cold substellar mass condensations are known to form, or in protoplanetary disks. Other molecules besides H_{2}, such as CO, can condense even sooner, but as their mass fraction is low they should not significantly perturb the gravitational balance, while He should not condense at all and should keep a minimal amount of gaseous phase^{2}. The precipitation of formed comets and planetoids in a gravitational field, separating the gaseous and condensed phases, is also a fascinating aspect that is as yet impossible to capture with traditional hydrodynamical codes, but possible with the molecular dynamics approach. We will pursue further simulation work, including more specific astrophysical applications.
In Safa & Pfenniger (2008) it was calculated that the mixture of He and H_{2} at cosmic abundance does not change the conclusion regarding the phase transition of H_{2} alone, as the species do not remain mixed when H_{2} alone condenses.
Acknowledgments
This work is supported by the STARFORM Sinergia Project funded by the Swiss National Science Foundation. We thank the LAMMPS team for providing a powerful open source tool to the scientific community. We thank the referee for a thorough reading of the manuscript and constructive comments, which substantially improved the paper.
References
 Air Liquide 1976, Gas Encyclopedia (Editor Elsevier) [Google Scholar]
 Allamandola, L. J., Bernstein, M. P., Sandford, S. A., & Walker, R. L. 1999,Space Sci. Rev., 90, 219 [Google Scholar]
 Binney, J., & Tremaine, S. 2008, Galactic Dynamics, 2nd edn. (Princeton: University Press) [Google Scholar]
 Bolatto, A. D., Wolfire, M., & Leroy, A. K. 2013, ARA&A, 51, 207 [NASA ADS] [CrossRef] [Google Scholar]
 Buch, V., & Devlin, J. P. 1994, ApJ, 431, L135 [NASA ADS] [CrossRef] [Google Scholar]
 Burkert, A., & O’dell, C. R. 1998, ApJ, 503, 792 [NASA ADS] [CrossRef] [Google Scholar]
 Carey, V. P. 1999, Statistical Thermodynamics and Microscale Thermophysics (Cambridge University Press) [Google Scholar]
 Clausius, R. 1870, Ann. Phys., 217, 124 [CrossRef] [Google Scholar]
 ClerkMaxwell, J. 1875, Nature, 11, 357 [NASA ADS] [CrossRef] [Google Scholar]
 Combes, F., & Pfenniger, D. 1997, A&A, 327, 453 [NASA ADS] [Google Scholar]
 Dissly, R. W., Allen, M., & Anicich, V. G. 1994, ApJ, 435, 685 [NASA ADS] [CrossRef] [Google Scholar]
 Draine, B. T. 2011, Physics of the Interstellar and Intergalactic Medium (Princeton: University Press) [Google Scholar]
 Elmegreen, B. G., & Scalo, J. 2004, ARA&A, 42, 211 [NASA ADS] [CrossRef] [Google Scholar]
 Ewald, P. P. 1921, Ann. Phys., 369, 253 [CrossRef] [Google Scholar]
 Grenier, I. A., Casandjian, J.M., & Terrier, R. 2005, Science, 307, 1292 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Hockney, R. W., & Eastwood, J. W. 1981, Computer Simulation Using Particles (New York: McGrawHill) [Google Scholar]
 Jeans, J. H. 1902, Roy. Soc. London Philos. Trans. Ser. A, 199, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Johnston, D. C. 2014, in Advances in Thermodynamics of the van der Waals Fluid, Chap. 4 (Morgan & Claypool Publishers) [Google Scholar]
 Jones, J. E. 1924, Roy. Soc. London Proc. Ser. A, 106, 463 [NASA ADS] [CrossRef] [Google Scholar]
 Kondepudi, D. K., & Prigogine, I. 1998, Modern thermodynamics: from heat engines to dissipative structures (John Wiley) [Google Scholar]
 Kristensen, L. E., Amiaud, L., Fillion, J.H., Dulieu, F., & Lemaire, J.L. 2011, A&A, 527, A44 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Langer, W. D., Velusamy, T., Pineda, J. L., et al. 2010, A&A, 521, L17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lekner, J. 1982, AIP, 50, 161 [Google Scholar]
 Paradis, D., Dobashi, K., Shimoikura, T., et al. 2012, A&A, 543, A103 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pfenniger, D. 2004, in Dark Matter in Galaxies, eds. S. Ryder, D. Pisano, M. Walker, & K. Freeman, IAU Symp., 220, 241 [Google Scholar]
 Pfenniger, D., & Combes, F. 1994, A&A, 285, 94 [NASA ADS] [Google Scholar]
 Pfenniger, D., Combes, F., & Martinet, L. 1994, A&A, 285, 79 [NASA ADS] [Google Scholar]
 Planck Collaboration XIX. 2011, A&A, 536, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Plimpton, S. 1995, J. Comput. Phys., 117, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Safa, Y., & Pfenniger, D. 2008, Eur. Phys. J. B, 66, 337 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sahai, R., & Nyman, L.A. 1997, ApJ, 487, L155 [NASA ADS] [CrossRef] [Google Scholar]
 Sandford, S. A., & Allamandola, L. J. 1993, ApJ, 409, L65 [NASA ADS] [CrossRef] [Google Scholar]
 Tuckerman, M., Berne, B. J., & Martyna, G. J. 1992, J. Chem. Phys., 97, 1990 [NASA ADS] [CrossRef] [Google Scholar]
 van der Waals, J. D. 1910, Koninklijke Nederlandse Akademie van Wetenschappen, Proc. Ser. B Phys. Sci., 13, 1253 [Google Scholar]
 Walker, M. A. 2013, MNRAS, 434, 2814 [NASA ADS] [CrossRef] [Google Scholar]
 Walsh, J. R., & Meaburn, J. 1993, The Messenger, 73, 35 [NASA ADS] [Google Scholar]
 Weinberg, S. 1972, Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity [Google Scholar]
Online material
Fig. 33 Snapshot time sequence of the simulation OPGs with N_{SM} = 100^{3}, showing the smooth formation of a planetoid. The slice size is 0.025 × 1 × 1 L^{3} at x = (0.5 ± 0.0125)L. 

Open with DEXTER 
Fig. 34 Snapshot comet massdistribution sequence of the simulation OPGs with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 35 Snapshot comet massdistribution sequence of all OPGs simulations with N_{SM} equal to the triangle pointing down: 50^{3}, triangle pointing up: 64^{3}, square: 80^{3}, diamond: 100^{3}, circle: 128^{3}, star: 160^{3}. 

Open with DEXTER 
Fig. 36 Snapshot comet massdistribution summation sequence of OPGs simulations. Reference with N_{SM} = 128^{3}: circle; sum of 4 comet sizes per star: 100^{3}. 

Open with DEXTER 
Fig. 37 0.2 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.1)L of simulation PT22 with N_{SM} = 50^{3}, and without gravity. 

Open with DEXTER 
Fig. 38 Comet massdistribution snapshots of simulation PT22 with N_{SM} = 50^{3}, and without gravity. 

Open with DEXTER 
Fig. 39 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of sufficiently selfgravitating simulation PT13 with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 40 Comet massdistribution snapshots of sufficiently selfgravitating simulation PT13 with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 41 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of weakly selfgravitating simulation PT22 with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 42 Comet mass distribution of weakly selfgravitating simulation PT22 with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 43 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of weakly selfgravitating simulation PT32 with N_{SM} = 100^{3}. 

Open with DEXTER 
Fig. 44 Comet mass distribution of weakly selfgravitating simulation PT32 with N_{SM} = 100^{3}. 

Open with DEXTER 
All Tables
Condensed and uncondensed bodies in a LJ fluid using the LagrangeJacobi identity.
All Figures
Fig. 1 H_{2} phase diagram (bold line) and He (dotted line) in cold and low pressure conditions. 

Open with DEXTER  
In the text 
Fig. 2 Critical points (bullet) and triple points (Y) of common molecules in the ISM. As for H_{2} in Fig. 1, the sublimation curves cross interstellar pressure conditions very steeply in pressure a few K below their triple point. For example, CO should essentially be frozen below ~ 20 K, so unable to emit rotation lines. 

Open with DEXTER  
In the text 
Fig. 3 van der Waals phase diagram for a fluid with T_{r} = 0.9. Dotted line: gas phase; dashed line: phase transition; solid line: solid phase. 

Open with DEXTER  
In the text 
Fig. 4 van der Waals phase diagram for a fluid with a temperature (from bottom to top) of T_{r} = 0.2, T_{r} = 0.4, T_{r} = 0.6, T_{r} = 0.8, T_{r} = 1.0, and T_{r} = 1.2. Solid line: van der Waals EOS with Maxwell construct, dotted line: original van der Waals EOS. 

Open with DEXTER  
In the text 
Fig. 5 H_{2} and van der Waals phase diagrams. Solid line: H_{2} laboratory data; dotted line: van der Waals vapour curve derived with the Maxwell construct. 

Open with DEXTER  
In the text 
Fig. 6 LennardJones potential. 

Open with DEXTER  
In the text 
Fig. 7 Minimum mass of isothermal equilibrium curves for H_{2} homogeneous spheres. 

Open with DEXTER  
In the text 
Fig. 8 Virial equilibrium of gravitating LJ homogeneous isothermal spheres with H_{2} molecules arranged on a HCP/FCC lattice. The sphere mass as a function of density is plotted for temperatures ranging from τ ≡ T/T_{max} = 10^{2} to 10^{2}. 

Open with DEXTER  
In the text 
Fig. 9 van der Waals EOS, including the Maxwell construct. 

Open with DEXTER  
In the text 
Fig. 10 Mean kinetic energy as a function of the maximum displacement in an ice clump. 

Open with DEXTER  
In the text 
Fig. 11 LJ potential energy as a function of the number of supermolecules at t = 1τ. Solid line: SM15; dashed line: SM30; dotted line: SM60. See Table 3. 

Open with DEXTER  
In the text 
Fig. 12 Fraction of bound molecules as a function of the number of supermolecules at t = 5τ. Solid line: SM04^{a}, dashed line: SM06^{a}. ^{(a)} See Table 3. 

Open with DEXTER  
In the text 
Fig. 13 zCoordinate of the centre of mass of comets as a function of time of the SF simulation. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). See Table 3. 

Open with DEXTER  
In the text 
Fig. 14 Temperature of the onephase fluid simulations OP0 (straight solid line) and OPGs as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). See Table 4. 

Open with DEXTER  
In the text 
Fig. 15 Fraction of bound molecules in comets of the onephase fluids OP0 (straight solid line) and OPGs as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER  
In the text 
Fig. 16 Temperature of the onephase fluid simulations OPGs as a function of time. N_{SM} = 100^{3} with four different random seeds. 

Open with DEXTER  
In the text 
Fig. 17 Density at x = (0.5 ± 0.05)L of the onephase fluids OP0 (straight solid line) and OPGs as a function of time. Bold line: analytic solution for ideal gas. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER  
In the text 
Fig. 18 Evolution of perturbation density. Straight dotted line: analytic solution for ideal gas. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER  
In the text 
Fig. 19 3Dview of simulation OPGs with N_{SM} = 100^{3} at t = 4τ. 

Open with DEXTER  
In the text 
Fig. 20 Distribution of y and zcoordinates of planetoid centre of mass on a 10 × 10 grid; all OPGs simulations at t = 2.5τ. 

Open with DEXTER  
In the text 
Fig. 21 Temperature distribution of unbound molecules and comets as a function of comet mass of OPGs with N_{SM} = 100^{3} at t = 4τ. 

Open with DEXTER  
In the text 
Fig. 22 Mean temperature over four onephase fluid simulations OPGs with different random seeds as a function of time. N_{SM} = 50^{3} and 64^{3} (dotted line), 80^{3} and 100^{3} (dashed line), 128^{3} and 160^{3} (solid line). 

Open with DEXTER  
In the text 
Fig. 23 Time and size of the largest comet at the appearance of the turning point. 

Open with DEXTER  
In the text 
Fig. 24 Extrapolation of the comet mass distribution for a real H_{2} molecular fluid at t = 5τ. 

Open with DEXTER  
In the text 
Fig. 25 Position of the LJ simulations in a corresponding van der Waals phase diagram. Solid line: Maxwell construct for the van der Waals EOS; dotted line: van der Waals EOS. Star: fluid presenting a phase transition; bullet: onephase fluid. 

Open with DEXTER  
In the text 
Fig. 26 Temperature of nongravitating fluids as a function of time. Dotted line: n = 10^{1}n_{cr}; solid line: n = 10^{2}n_{cr}; dashed line: n = 10^{3}n_{cr}. All simulations with N_{SM} = 50^{3}. 

Open with DEXTER  
In the text 
Fig. 27 Closeup 3Dview of a medium sized comet and smaller aggregates. 

Open with DEXTER  
In the text 
Fig. 28 Temperature distribution of unbound molecules and comets as a function of comet size of PT22^{a} without gravity and N_{SM} = 50^{3} at t = 100τ. ^{(a)} See Table 5. 

Open with DEXTER  
In the text 
Fig. 29 Temperature of sufficiently selfgravitating fluids as a function of time. Dotted line: OPGs^{a}; dashed line: PT13^{a}; solid line (from left to right): PT32^{a}, PT52^{a}, PT72^{a}, PT92^{a}; dashdotted line: PT11^{a}. All simulations with N_{SM} = 100^{3}. ^{(a)} See Table 5 

Open with DEXTER  
In the text 
Fig. 30 Temperature of weakly selfgravitating fluids as a function of time. Dotted line: PT22 and PT32, with N_{SM} = 50^{3} and without gravity; dashed line: PT22, γ_{J} = 0.5 and N_{SM} = 80^{3} and 100^{3}; solid line PT32, γ_{J} = 0.5 and N_{SM} = 80^{3} and 100^{3}; dashdotted line PT52, γ_{J} = 0.5 and N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 31 Temperature distribution of unbound molecules and comets as a function of comet size of PT32 with γ_{J} = 0.5 and N_{SM} = 100^{3} at t = 200τ. 

Open with DEXTER  
In the text 
Fig. 32 Density of planetoid as a function of the radius. Dashed line: OPGs with N_{SM} = 100^{3} at t = 4τ, the gaseous planetoid consists of 41559 supermolecules (0.042 M_{tot}). Dashdotted line: PT32, with γ_{J} = 0.5 and N_{SM} = 100^{3} at t = 200τ; the “rocky planetoid” consists of 74208 supermolecules (0.074 M_{tot}). Solid lines: hydrogen, higher value: solid; lower value: liquid. Dotted line: ρ_{f} of LagrangeJacobi identity, higher value: SC; lower value: HCP/FCC. 

Open with DEXTER  
In the text 
Fig. 33 Snapshot time sequence of the simulation OPGs with N_{SM} = 100^{3}, showing the smooth formation of a planetoid. The slice size is 0.025 × 1 × 1 L^{3} at x = (0.5 ± 0.0125)L. 

Open with DEXTER  
In the text 
Fig. 34 Snapshot comet massdistribution sequence of the simulation OPGs with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 35 Snapshot comet massdistribution sequence of all OPGs simulations with N_{SM} equal to the triangle pointing down: 50^{3}, triangle pointing up: 64^{3}, square: 80^{3}, diamond: 100^{3}, circle: 128^{3}, star: 160^{3}. 

Open with DEXTER  
In the text 
Fig. 36 Snapshot comet massdistribution summation sequence of OPGs simulations. Reference with N_{SM} = 128^{3}: circle; sum of 4 comet sizes per star: 100^{3}. 

Open with DEXTER  
In the text 
Fig. 37 0.2 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.1)L of simulation PT22 with N_{SM} = 50^{3}, and without gravity. 

Open with DEXTER  
In the text 
Fig. 38 Comet massdistribution snapshots of simulation PT22 with N_{SM} = 50^{3}, and without gravity. 

Open with DEXTER  
In the text 
Fig. 39 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of sufficiently selfgravitating simulation PT13 with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 40 Comet massdistribution snapshots of sufficiently selfgravitating simulation PT13 with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 41 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of weakly selfgravitating simulation PT22 with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 42 Comet mass distribution of weakly selfgravitating simulation PT22 with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 43 0.025 × 1 × 1 L^{3} snapshots at x = (0.5 ± 0.0125)L of weakly selfgravitating simulation PT32 with N_{SM} = 100^{3}. 

Open with DEXTER  
In the text 
Fig. 44 Comet mass distribution of weakly selfgravitating simulation PT32 with N_{SM} = 100^{3}. 

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.