Free Access
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/0004-6361/201424798
Published online 25 May 2015

© ESO, 2015

1. Introduction

Typically around 50% or more of the gas in spiral galaxies consists of H2, 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 H2 in the Milky Way has essentially been doubled, effectively revealing the nature of some of the dark baryons. It is commonly admitted that the CO-related molecular hydrogen is present in relatively dense regions of the interstellar medium, molecular clouds with number density >1010 m-3 and temperatures of 730 K (Draine 2011).

Even though H2 is by far the most abundant molecule (~90%), molecular clouds are mainly detected by CO emissions because of all the difficulties in detecting cold H2 (Bolatto et al. 2013). For example, H2 only starts emitting at temperatures >512 K. See Combes & Pfenniger (1997) for a review of several possible methods for detecting cold H2. Because of these detection difficulties, the real quantity of H2 (as well as He, which shares similar properties of discreteness with H2) 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 H2 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 H2O, CO, CO2 and NH3 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. H2 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 H2 is mixed within H2O-rich grains in conditions that are too warm to allow the bulk of H2 to condense (Kristensen et al. 2011).

thumbnail Fig. 1

H2 phase diagram (bold line) and He (dotted line) in cold and low pressure conditions.

Open with DEXTER

High-resolution 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 H2 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 H2 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 self-gravitating gas leading to collapse has been studied for over a century. The molecular cloud is normally represented as a self-gravitating, one-phase 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 one-phase fluid is valid in many cases, but when considering very cold, high-density regions, it is an over-simplification. In these condition, H2 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 one-phase 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 high-density 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 H2 ice fragments in molecular clouds due to the fragmentation of cold high-density 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).

thumbnail Fig. 2

Critical points (bullet) and triple points (Y) of common molecules in the ISM. As for H2 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 self-gravitating 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 Lennard-Jones (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 self-gravitating 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 super-molecules to enable us to perform simulations with a total mass high enough for the fluid to be self-gravitating.

The simulations performed are discussed in Sect. 4. First, the correctness of the super-molecule approach is tested. Second, one-phase 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 non-gravitating 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 self-gravitating 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 ideal-gas law, taking the finite size of molecules and intermolecular interactions into account. It links pressure, density, and temperature as follows: (1)with a [Pa m6] and b [m3] being constants characteristic of the fluid. It can also be expressed in a reduced, dimensionless form as, (2)where Pr = P/Pc, nr = n/nc, and Tr = T/Tc. The parameters Pc, Tc and nc 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.9Tc. 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.

thumbnail Fig. 3

van der Waals phase diagram for a fluid with Tr = 0.9. Dotted line: gas phase; dashed line: phase transition; solid line: solid phase.

Open with DEXTER

thumbnail Fig. 4

van der Waals phase diagram for a fluid with a temperature (from bottom to top) of Tr = 0.2, Tr = 0.4, Tr = 0.6, Tr = 0.8, Tr = 1.0, and Tr = 1.2. Solid line: van der Waals EOS with Maxwell construct, dotted line: original van der Waals EOS.

Open with DEXTER

thumbnail Fig. 5

H2 and van der Waals phase diagrams. Solid line: H2 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 (Clerk-Maxwell 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.2Tc to T = 1.2Tc. 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 TTc 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 H2 over a wide range of pressure (Fig. 5).

thumbnail Fig. 6

Lennard-Jones potential.

Open with DEXTER

2.1.1. Lennard-Jones 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 Lennard-Jones (LJ) intermolecular potential (Jones 1924). No local or global thermal equilibirum is required. The LJ potential consists of an attractive long-range term and a repulsive short-range term (Fig. 6), (3)with rσ = r/σ. Its minimum value, located at rσ = 21 / 6, is ϵ/m.

2.2. Virial theorem for a Lennard-Jones gravitating fluid

2.2.1. Lagrange-Jacobi 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 self-gravitating 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 Lagrange-Jacobi identity path, by taking the second time derivative of the moment of the polar inertia , we find, (4)where the first right-hand term is twice the kinetic energy Ekin, 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 functions1Φ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 Lagrange-Jacobi identity becomes (5)

2.2.2. Homogeneous sphere

For homogeneous, finite mass spheres at a given temperature, the individual terms of the Lagrange-Jacobi identity read,

where the constants cr and ca in the LJ terms have been calculated by straight summation over the 1.5 × 1012 nearest molecules, as given in Table 1 for both the simple cubic (SC) and hexagonal close packed (HCP) or face-centred cubic (FCC) lattices.

Table 1

Lattice constants.

The virial equation, obtained by setting d2I/ dt2 = 0 in the Lagrange-Jacobi identity, contains terms proportional to M except the gravity term, which is proportional to M5 / 3. Dividing the equation by M we can separate M2 / 3, yielding, (10)which is indeed positive if the sum of the terms inside the brackets are also positive. However the term proportional to ca is negative, and when large introduces the possibility that no positive solution for M2 / 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 non-negative when the term inside the square root above is non-negative. 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 ϵ/kB, we see that Tmax can be substantially larger than the critical temperature. With the values given in Table 1 for the constants ca and cr, the maximum temperature below which gravity combined with molecular forces enhances fragmentation are Tmax = 414.3 K and 626.8 K for the SC and HCP/FCC lattices, respectively, which are 11.38 and 17.22 larger than ϵ/kB.

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.

thumbnail Fig. 7

Minimum mass of isothermal equilibrium curves for H2 homogeneous spheres.

Open with DEXTER

thumbnail Fig. 8

Virial equilibrium of gravitating LJ homogeneous isothermal spheres with H2 molecules arranged on a HCP/FCC lattice. The sphere mass as a function of density is plotted for temperatures ranging from τT/Tmax = 10-2 to 102.

Open with DEXTER

The Lagrange-Jacobi identity therefore allows us to predict a density ρf and a maximum temperature Tmax 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 self-gravitating gas. For H2 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 H2 density (80 kg m-3).

At temperatures T>Tmax the isothermal equilibrium curves have a minimum mass when dM/dρ | T = 0, which is expressed as (14)where τ = T/Tmax and (15)Figure 7 shows the minimum mass as a function of the temperature for HCP/FCC H2 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 H2 dissociation.

2.2.3. Condensed bodies

Different condensed and uncondensed bodies can be identified using the Lagrange-Jacobi identity, as summarized in Table 2.

Table 2

Condensed and uncondensed bodies in a LJ fluid using the Lagrange-Jacobi 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 H2 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>Tmax, and the accumulation curve is the domain of rocky planetoids.

At temperatures below Tmax = 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<Tmax, 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 H2 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 self-gravitating 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 real-exponential, 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)

thumbnail 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 cP and cv both finite values, we find (∂P/∂ρ)s = 0. Therefore, a self-gravitating 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 ideal-gas 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 H2 at 3 K the critical pressure is 6.4 × 10-9 Pa.

3. Molecular dynamics

For all simulations, the Large-Scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) is used (Plimpton 1995). The LAMMPS software is a state of the art and widely used single and multi-processor code in chemistry, material sciences, and related fields. Its abilities to quickly compute short and long-range forces and the possibility to use a multi-timescale 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 self-gravitational potentials.

Since the straight calculation cost for exact pairwise interactions is O(N2), approximate but still accurate methods are implemented that make the calculations possible at a much lower cost. A cut-off radius rc is used for the short-range forces. As the attractive part of the LJ potential drops with r-6, it is calculated only for neighbour particles within rc. The neighbours for each particle are found using a Verlet neighbour list. This list is created with a radius of rn = rc + rs, rs being an extra skin distance to avoid the recalculation of the Verlet list at every time-step. This cut-off 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 Particle3-Mesh (P3M) method (Hockney & Eastwood 1981). The gravitational potential is split in short-range and long-range parts in Fourier space (Ewald 1921): (22)where rs 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 self-gravity is too weak to be of any influence, no specific self-gravity 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 qi,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 time-step. For short-range force, the time-step is constant throughout the simulation and set as 10-210-3, which is the typical interaction timescale during nearest-neighbour 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 short-range forces, which is why a very small time-step is required. But these small changes have almost no influence on long-range forces. For this reason, there is no need to calculate the long-range forces at every time-step.

In regard to short-range calculations, long-range calculations are an order of magnitude more costly per step but always need the same amount of calculation time: creation of the density-map 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 long-range 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 time-integration-scheme looks as follows: (28)where nSR is the number of short-range iterations and set as 10100, KLR is the kick operator using the long-range (FFT) accelerations, and KSR the kick operator using the short-range accelerations.

3.3. Super-molecules

The maximum number of particles that can be simulated on today’s supercomputers is of the order of 1091010 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 super-molecules is introduced. This concept is well established in galactic dynamics simulation, where super-stars weigh typically 104−106M, and in cosmology where super-WIMPS may weigh as much as ~1067 GeVc-2 particles. Two-body relaxation or diffusion due to the low number of super-particles 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 super-molecule consists of η molecules, and its mass is therefore, (29)To have the same properties of a super-molecule 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, NMmM = NSMmSM, the velocity dispersion of molecules and super-molecules 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 super-molecule properties.

It is important to ensure that the gravitational force between two super-molecules remains small compared to the corresponding LJ force within the short-range force cut-off radius. This is assured by setting FG(rc) ≪ FLJ(rc) with rc being the cut-off radius. This leads to the following constraint: (36)where rc, SM = (rc/σSM). Using a cut-off radius of rc = 4 σ and hydrogen super-molecules, we find a maximum super-molecule mass equal to 5.7 × 10-6M, meaning at least 1.7 × 105 particles are required for simulating an Earth-mass, virialized, H2-condensed body.

thumbnail 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 rm = 21 / 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 rm.

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.

thumbnail Fig. 11

LJ potential energy as a function of the number of super-molecules at t = 1τ. Solid line: SM15; dashed line: SM30; dotted line: SM60. See Table 3.

Open with DEXTER

Table 3

Parameters of the super-molecule test simulations.

4. Simulations

4.1. Units

A fluid is defined by the number of super-molecules NSM, the initial velocity distribution (temperature), the number density, and the strength of the gravitational potential. All simulation properties are molecule-independent, the initial temperature and density are measured as a factor of the critical values Tcr and ncr, and the gravitational potential is measured as the factor γJ = G/GJ of the ideal-gas Jeans gravity GJ, defined as (38)with the box side L = (NSM/n)1 / 3. It is interesting to note that GJ 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. Super-molecule 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 ncr is simulated with LAMMPS. Three different initial Maxwellian velocity distributions with temperature T = 1.5, 3.0 and 6.0 Tcr, and a number of particles from NSM = 3032003 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 super-molecules. One can see fluctuations in the simulations with low NSM, but the result becomes reasonable when NSM> 104 and a satisfactory convergence is obtained if NSM> 105. Therefore all the subsequent simulations are performed with NSM> 105.

thumbnail Fig. 12

Fraction of bound molecules as a function of the number of super-molecules at t = 5τ. Solid line: SM04a, dashed line: SM06a. (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 super-molecules. 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 super-molecules. All simulations reach the same percentage, independent of the number of super-molecules.

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 zcoordinate of the centre of mass of comets as a function of time. The curves are very similar, with a slight over-estimation of the collapse time for the smallest simulation. Therefore the precipitation phenomenon timescale is not strongly dependent on the mass resolution.

thumbnail Fig. 13

z-Coordinate of the centre of mass of comets as a function of time of the SF simulation. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (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 vi of velocities for the homogeneous unperturbed case, the x-component of the velocities is perturbed in such a way as to conserve energy. The perturbed x-velocity 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.

Table 4

Parameters of the one-phase fluid simulations.

4.4. One-phase fluid

To study the reaction of a pure one-phase ideal-gas fluid far from the phase transition but with the above plane-wave perturbation, the initial temperature is taken above the critical value (T0 = 1.25 Tcr) and the initial density well below the critical value (n0 = 10-2ncr). Three different cases are studied: without gravity, weakly self-gravitating below the classical ideal-gas Jeans criterion with γJ = 0.5, and sufficiently self-gravitating above the ideal-gas Jeans criterion with γJ = 1.25. The simulations are performed with NSM ranging from 503 to 1603. Table 4 shows the parameters of the different simulations.

thumbnail Fig. 14

Temperature of the one-phase fluid simulations OP0 (straight solid line) and OPGs as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line). See Table 4.

Open with DEXTER

thumbnail Fig. 15

Fraction of bound molecules in comets of the one-phase fluids OP0 (straight solid line) and OPGs as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (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 one-phase fluid simulations: OP0 without gravity and NSM = 1003 super-molecules and the sufficiently self-gravitating fluid OPGs with NSM ranging from 503 to 1603. The weakly self-gravitating 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 self-gravitating 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 x-direction does not change the nature of the fluid; its temperature and comet percentage remains the same.

The sufficiently self-gravitating 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 ~6Tcr, 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 NSM = 1003, leading to a time difference of ~0.5τ for reaching the asymptotic upper value.

thumbnail Fig. 16

Temperature of the one-phase fluid simulations OPGs as a function of time. NSM = 1003 with four different random seeds.

Open with DEXTER

thumbnail Fig. 17

Density at x = (0.5 ± 0.05)L of the one-phase fluids OP0 (straight solid line) and OPGs as a function of time. Bold line: analytic solution for ideal gas. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (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 ideal-gas 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.

thumbnail Fig. 18

Evolution of perturbation density. Straight dotted line: analytic solution for ideal gas. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line).

Open with DEXTER

thumbnail Fig. 19

3D-view of simulation OPGs with NSM = 1003 at t = 4τ.

Open with DEXTER

thumbnail Fig. 20

Distribution of y- and z-coordinates of planetoid centre of mass on a 10 × 10 grid; all OPGs simulations at t = 2.5τ.

Open with DEXTER

thumbnail Fig. 21

Temperature distribution of unbound molecules and comets as a function of comet mass of OPGs with NSM = 1003 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 power-law 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 power-law 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 NSM, 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 high-density plane is created thanks to the plane perturbation parallel to the yz-plane at x = 0.5. One can observe in Fig. 33 at t = 2τ that two filaments form, along the y- and z-axes, 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 plane-wave 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 super-molecules 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 (Mcomet = 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 super-molecules 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.

thumbnail Fig. 22

Mean temperature over four one-phase fluid simulations OPGs with different random seeds as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (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 NSM = 503 and 643 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 xGLJ values (Eq. (40)), 3.65 and 4.03 rσ. The first value is in fact lower than the cut-off 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 short-range 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 super-molecules is the same for all NSM. For example, at t = 5τ, all simulations have a fraction of ~10-1 for comets consisting of two super-molecules. These small clumps should be called multimers as mentioned below. But, with increasing NSM, the mass of a comet consisting of two super-molecules diminishes since M2 SM = 2 /NSM.

On the other hand, for large comets and the planetoid, the mass is invariant to NSM. For example, the planetoid at t = 5τ has the same mass of Mplanetoid ≈ 5 × 10-2Mtot for all NSM.

The fact that the mass distributions are shifted to the left for larger NSM is misleading, as shown in Fig. 36. Only the two largest simulations are shown, the simulation with NSM = 1283 is shown normally, but for the simulation with NSM = 1603 is downscaled to 1603/ 2 ≈ 1283, 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 NSM. The comet mass, on the other hand, clearly follows a power law with: (45)with ξN ≈ −1. Since NSM = Mtot/MSM, at turning point the number of super-molecules 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.

thumbnail 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 super-molecules. Considering H2 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 × 1051 H2 molecules, the number of molecules per super-molecule η varying from 2.0 × 1046 for NSM = 503 to 3.2 × 1044 for NSM = 1603. 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 1040 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.

thumbnail Fig. 24

Extrapolation of the comet mass distribution for a real H2 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 multi-molecules 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 mH = 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 power-law distribution rapidly increases the mass in larger grain- and comet-sized bodies until they become a planetoid mass.

thumbnail 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: one-phase 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-3ncr, and three different temperatures, 0.1, 0.2 and 0.3 Tcr, are chosen. As the one-phase fluid studied in Sect. 4.4 also has a density of 10-2ncr, three additional temperatures, 0.5, 0.7 and 0.9 Tcr, were used for this density to make a link between the one-phase fluid and the fluids studied in this section.

Looking at the phase transition diagram (Fig. 25), all fluids with T ≤ 0.3Tc 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.5Tc are still in the stable regime.

Three different cases are studied: without gravity, strongly self-gravitating above the Jeans criterion, and weakly self-gravitating below the Jeans criterion.

thumbnail Fig. 26

Temperature of non-gravitating fluids as a function of time. Dotted line: n = 10-1ncr; solid line: n = 10-2ncr; dashed line: n = 10-3ncr. All simulations with NSM = 503.

Open with DEXTER

Table 5

Parameters of the phase transition simulations.

thumbnail Fig. 27

Close-up 3D-view of a medium sized comet and smaller aggregates.

Open with DEXTER

thumbnail Fig. 28

Temperature distribution of unbound molecules and comets as a function of comet size of PT2-2a without gravity and NSM = 503 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 xGLJ value does not have to be considered and the number of super-molecules is set to NSM = 503. 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 PT1-2 and PT2-2 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 super-molecules per comets is <10-2NSM. 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 super-molecules. Figure 27 shows a close-up 3D view of a medium-sized 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.

thumbnail Fig. 29

Temperature of sufficiently self-gravitating fluids as a function of time. Dotted line: OPGsa; dashed line: PT1-3a; solid line (from left to right): PT3-2a, PT5-2a, PT7-2a, PT9-2a; dash-dotted line: PT1-1a. All simulations with NSM = 1003. (a) See Table 5

Open with DEXTER

4.5.2. Self-gravitating fluid above Jeans criterion

All sufficiently self-gravitating simulations use the same gravitational potential, G = 1.25 GJ(T = 1.25 Tcr,n = 0.01 ncr). As the temperatures and densities differ, GJ 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-2ncr 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 power-law index ξc also evolves in a similar fashion as in simulation OPGs, starting out very steep and reaching a final value of −2.

thumbnail Fig. 30

Temperature of weakly self-gravitating fluids as a function of time. Dotted line: PT2-2 and PT3-2, with NSM = 503 and without gravity; dashed line: PT2-2, γJ = 0.5 and NSM = 803 and 1003; solid line PT3-2, γJ = 0.5 and NSM = 803 and 1003; dash-dotted line PT5-2, γJ = 0.5 and NSM = 1003.

Open with DEXTER

thumbnail Fig. 31

Temperature distribution of unbound molecules and comets as a function of comet size of PT3-2 with γJ = 0.5 and NSM = 1003 at t = 200τ.

Open with DEXTER

4.5.3. Self-gravitating fluid below Jeans criterion

To compare weakly self-gravitating fluids in a phase transition (on the Maxwell line) with out-of-phase-transition 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 self-gravitating one-phase fluid OPGw does not create any effect if γJ< 1. Figure 30 shows the temperature as a function of time of the self-gravitating fluids PT5-2, PT3-2, and PT2-2, and compares them to the non-gravitating fluids. The out-of-phase-transition fluid PT5-2 remains unchanged and does not differ from the non-gravitating fluid, identical to OPGw.

The temperatures of the fluids in a phase transition PT3-2 and PT2-2 are exponentially growing and differ clearly from the non-gravitating ones. In the beginning, the LJ forces dominate the gravitational forces and the fluids behave like non-gravitating fluids. PT2-2 reacts much faster than PT3-2, thanks to its low initial temperature and fluctuations due to the Maxwellian velocity distribution,and forms many small- and medium-sized 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 low-temperature, non-gravitating 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 medium-sized comets during its growth.

With its higher temperature, PT3-2 forms almost no small- or medium-sized 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 medium-sized comet, which is too massive to fit into the power law ξ1 (see Fig. 43, t = 50τ). From then on, this comet attracts super-molecules 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 PT3-2. The temperature of the comets containing few super-molecules lies below the average, but the planetoids temperature is the same as the average.

PT2-2, PT3-2 and PT5-2 all scale very well. The two fluids in a phase transition, PT2-2 and PT3-2, have an almost identical timescale for the exponential growth for NSM = 803 and 1003. The slight difference is due to the initial random seed. As xGLJ depends on γJ and the temperature, both much lower than for OPGs, its value is clearly above the cut-off radius, which is 7.63 and 8.34 for PT2-2 and 7.03 and 7.69 for PT3-2.

PT5-2 shows no difference for any number of super-molecules, with NSM ranging from 503 to 1003.

thumbnail Fig. 32

Density of planetoid as a function of the radius. Dashed line: OPGs with NSM = 1003 at t = 4τ, the gaseous planetoid consists of 41559 super-molecules (0.042 Mtot). Dash-dotted line: PT3-2, with γJ = 0.5 and NSM = 1003 at t = 200τ; the “rocky planetoid” consists of 74208 super-molecules (0.074 Mtot). Solid lines: hydrogen, higher value: solid; lower value: liquid. Dotted line: ρf of Lagrange-Jacobi 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 PT3-2 with γ = 0.5γJ, and compares them to hydrogen laboratory data and values found in the Lagrange-Jacobi 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 Tcr), super-molecules 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 super-molecules and therefore a lower density compared to super-molecules 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 Lagrange-Jacobi identity (Eq. (13)), in accordance with Table 2 and Fig. 8 for gaseous planetoids.

The density of the PT3-2 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 H2 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 Lennard-Jones 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 H2, this maximum temperature lies in the range of 400600 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 state-of-the-art molecular dynamics simulator LAMMPS. We used super-molecules to simulate gravitational and molecular forces together with a computationally tractable number of particles.

5.2.1. Super-molecules

We tested the super-molecule concept thoroughly. We achieved good scaling for non-gravitational effects if the number of super-molecules is large enough (105). The sticking point is to ensure that the molecular LJ forces are dominating the gravitational forces in close-range interactions. This is achieved by setting the number of super-molecules to be large enough so that the distance at which the two forces are equal is above the cut-off radius 4 σSM. For an H2 fluid this sets a maximum super-molecule mass equal to 5.7 × 10-6M, limiting the total mass that can be simulated by the available computing power.

In principle, the super-molecule 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. One-phase fluid

We applied a plane sinusoidal perturbation into a one-phase fluid with a temperature above the critical value. The results reproduce the ideal-gas 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 medium-sized 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 medium-sized 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 self-gravitating (above and below the ideal-gas Jeans criterion).

Because of the Maxwellian velocity distribution fluctuations, the non-gravitating fluids form small- and medium-sized comets until the potential and kinetic energy reach an equilibrium; thereafter they remain statistically stable. In the absence of any long-range force, no planetoid forms.

The self-gravitating fluids with a gravitational potential above the ideal-gas Jeans criterion do not react differently to a plane sinusoidal perturbation from a one-phase fluid: the perturbation grows exponentially and its temperature rises, and small- and medium-sized 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 ideal-gas Jeans threshold. The performed simulations of weakly self-gravitating 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 self-gravitating fluids, but, in the end, a big central planetoid forms anyway. While the out-of-phase 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 H2 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 H2, 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 phase2. 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.


1

By definition, a homogeneous function f of degree k satisfies f(λx) = λkf(x).

2

In Safa & Pfenniger (2008) it was calculated that the mixture of He and H2 at cosmic abundance does not change the conclusion regarding the phase transition of H2 alone, as the species do not remain mixed when H2 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

Online material

thumbnail Fig. 33

Snapshot time sequence of the simulation OPGs with NSM = 1003, showing the smooth formation of a planetoid. The slice size is 0.025 × 1 × 1 L3 at x = (0.5 ± 0.0125)L.

Open with DEXTER

thumbnail Fig. 34

Snapshot comet mass-distribution sequence of the simulation OPGs with NSM = 1003.

Open with DEXTER

thumbnail Fig. 35

Snapshot comet mass-distribution sequence of all OPGs simulations with NSM equal to the triangle pointing down: 503, triangle pointing up: 643, square: 803, diamond: 1003, circle: 1283, star: 1603.

Open with DEXTER

thumbnail Fig. 36

Snapshot comet mass-distribution summation sequence of OPGs simulations. Reference with NSM = 1283: circle; sum of 4 comet sizes per star: 1003.

Open with DEXTER

thumbnail Fig. 37

0.2 × 1 × 1 L3 snapshots at x = (0.5 ± 0.1)L of simulation PT2-2 with NSM = 503, and without gravity.

Open with DEXTER

thumbnail Fig. 38

Comet mass-distribution snapshots of simulation PT2-2 with NSM = 503, and without gravity.

Open with DEXTER

thumbnail Fig. 39

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of sufficiently self-gravitating simulation PT1-3 with NSM = 1003.

Open with DEXTER

thumbnail Fig. 40

Comet mass-distribution snapshots of sufficiently self-gravitating simulation PT1-3 with NSM = 1003.

Open with DEXTER

thumbnail Fig. 41

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of weakly self-gravitating simulation PT2-2 with NSM = 1003.

Open with DEXTER

thumbnail Fig. 42

Comet mass distribution of weakly self-gravitating simulation PT2-2 with NSM = 1003.

Open with DEXTER

thumbnail Fig. 43

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of weakly self-gravitating simulation PT3-2 with NSM = 1003.

Open with DEXTER

thumbnail Fig. 44

Comet mass distribution of weakly self-gravitating simulation PT3-2 with NSM = 1003.

Open with DEXTER

All Tables

Table 1

Lattice constants.

Table 2

Condensed and uncondensed bodies in a LJ fluid using the Lagrange-Jacobi identity.

Table 3

Parameters of the super-molecule test simulations.

Table 4

Parameters of the one-phase fluid simulations.

Table 5

Parameters of the phase transition simulations.

All Figures

thumbnail Fig. 1

H2 phase diagram (bold line) and He (dotted line) in cold and low pressure conditions.

Open with DEXTER
In the text
thumbnail Fig. 2

Critical points (bullet) and triple points (Y) of common molecules in the ISM. As for H2 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
thumbnail Fig. 3

van der Waals phase diagram for a fluid with Tr = 0.9. Dotted line: gas phase; dashed line: phase transition; solid line: solid phase.

Open with DEXTER
In the text
thumbnail Fig. 4

van der Waals phase diagram for a fluid with a temperature (from bottom to top) of Tr = 0.2, Tr = 0.4, Tr = 0.6, Tr = 0.8, Tr = 1.0, and Tr = 1.2. Solid line: van der Waals EOS with Maxwell construct, dotted line: original van der Waals EOS.

Open with DEXTER
In the text
thumbnail Fig. 5

H2 and van der Waals phase diagrams. Solid line: H2 laboratory data; dotted line: van der Waals vapour curve derived with the Maxwell construct.

Open with DEXTER
In the text
thumbnail Fig. 6

Lennard-Jones potential.

Open with DEXTER
In the text
thumbnail Fig. 7

Minimum mass of isothermal equilibrium curves for H2 homogeneous spheres.

Open with DEXTER
In the text
thumbnail Fig. 8

Virial equilibrium of gravitating LJ homogeneous isothermal spheres with H2 molecules arranged on a HCP/FCC lattice. The sphere mass as a function of density is plotted for temperatures ranging from τT/Tmax = 10-2 to 102.

Open with DEXTER
In the text
thumbnail Fig. 9

van der Waals EOS, including the Maxwell construct.

Open with DEXTER
In the text
thumbnail Fig. 10

Mean kinetic energy as a function of the maximum displacement in an ice clump.

Open with DEXTER
In the text
thumbnail Fig. 11

LJ potential energy as a function of the number of super-molecules at t = 1τ. Solid line: SM15; dashed line: SM30; dotted line: SM60. See Table 3.

Open with DEXTER
In the text
thumbnail Fig. 12

Fraction of bound molecules as a function of the number of super-molecules at t = 5τ. Solid line: SM04a, dashed line: SM06a. (a) See Table 3.

Open with DEXTER
In the text
thumbnail Fig. 13

z-Coordinate of the centre of mass of comets as a function of time of the SF simulation. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line). See Table 3.

Open with DEXTER
In the text
thumbnail Fig. 14

Temperature of the one-phase fluid simulations OP0 (straight solid line) and OPGs as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line). See Table 4.

Open with DEXTER
In the text
thumbnail Fig. 15

Fraction of bound molecules in comets of the one-phase fluids OP0 (straight solid line) and OPGs as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line).

Open with DEXTER
In the text
thumbnail Fig. 16

Temperature of the one-phase fluid simulations OPGs as a function of time. NSM = 1003 with four different random seeds.

Open with DEXTER
In the text
thumbnail Fig. 17

Density at x = (0.5 ± 0.05)L of the one-phase fluids OP0 (straight solid line) and OPGs as a function of time. Bold line: analytic solution for ideal gas. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line).

Open with DEXTER
In the text
thumbnail Fig. 18

Evolution of perturbation density. Straight dotted line: analytic solution for ideal gas. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line).

Open with DEXTER
In the text
thumbnail Fig. 19

3D-view of simulation OPGs with NSM = 1003 at t = 4τ.

Open with DEXTER
In the text
thumbnail Fig. 20

Distribution of y- and z-coordinates of planetoid centre of mass on a 10 × 10 grid; all OPGs simulations at t = 2.5τ.

Open with DEXTER
In the text
thumbnail Fig. 21

Temperature distribution of unbound molecules and comets as a function of comet mass of OPGs with NSM = 1003 at t = 4τ.

Open with DEXTER
In the text
thumbnail Fig. 22

Mean temperature over four one-phase fluid simulations OPGs with different random seeds as a function of time. NSM = 503 and 643 (dotted line), 803 and 1003 (dashed line), 1283 and 1603 (solid line).

Open with DEXTER
In the text
thumbnail Fig. 23

Time and size of the largest comet at the appearance of the turning point.

Open with DEXTER
In the text
thumbnail Fig. 24

Extrapolation of the comet mass distribution for a real H2 molecular fluid at t = 5τ.

Open with DEXTER
In the text
thumbnail 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: one-phase fluid.

Open with DEXTER
In the text
thumbnail Fig. 26

Temperature of non-gravitating fluids as a function of time. Dotted line: n = 10-1ncr; solid line: n = 10-2ncr; dashed line: n = 10-3ncr. All simulations with NSM = 503.

Open with DEXTER
In the text
thumbnail Fig. 27

Close-up 3D-view of a medium sized comet and smaller aggregates.

Open with DEXTER
In the text
thumbnail Fig. 28

Temperature distribution of unbound molecules and comets as a function of comet size of PT2-2a without gravity and NSM = 503 at t = 100τ. (a) See Table 5.

Open with DEXTER
In the text
thumbnail Fig. 29

Temperature of sufficiently self-gravitating fluids as a function of time. Dotted line: OPGsa; dashed line: PT1-3a; solid line (from left to right): PT3-2a, PT5-2a, PT7-2a, PT9-2a; dash-dotted line: PT1-1a. All simulations with NSM = 1003. (a) See Table 5

Open with DEXTER
In the text
thumbnail Fig. 30

Temperature of weakly self-gravitating fluids as a function of time. Dotted line: PT2-2 and PT3-2, with NSM = 503 and without gravity; dashed line: PT2-2, γJ = 0.5 and NSM = 803 and 1003; solid line PT3-2, γJ = 0.5 and NSM = 803 and 1003; dash-dotted line PT5-2, γJ = 0.5 and NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 31

Temperature distribution of unbound molecules and comets as a function of comet size of PT3-2 with γJ = 0.5 and NSM = 1003 at t = 200τ.

Open with DEXTER
In the text
thumbnail Fig. 32

Density of planetoid as a function of the radius. Dashed line: OPGs with NSM = 1003 at t = 4τ, the gaseous planetoid consists of 41559 super-molecules (0.042 Mtot). Dash-dotted line: PT3-2, with γJ = 0.5 and NSM = 1003 at t = 200τ; the “rocky planetoid” consists of 74208 super-molecules (0.074 Mtot). Solid lines: hydrogen, higher value: solid; lower value: liquid. Dotted line: ρf of Lagrange-Jacobi identity, higher value: SC; lower value: HCP/FCC.

Open with DEXTER
In the text
thumbnail Fig. 33

Snapshot time sequence of the simulation OPGs with NSM = 1003, showing the smooth formation of a planetoid. The slice size is 0.025 × 1 × 1 L3 at x = (0.5 ± 0.0125)L.

Open with DEXTER
In the text
thumbnail Fig. 34

Snapshot comet mass-distribution sequence of the simulation OPGs with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 35

Snapshot comet mass-distribution sequence of all OPGs simulations with NSM equal to the triangle pointing down: 503, triangle pointing up: 643, square: 803, diamond: 1003, circle: 1283, star: 1603.

Open with DEXTER
In the text
thumbnail Fig. 36

Snapshot comet mass-distribution summation sequence of OPGs simulations. Reference with NSM = 1283: circle; sum of 4 comet sizes per star: 1003.

Open with DEXTER
In the text
thumbnail Fig. 37

0.2 × 1 × 1 L3 snapshots at x = (0.5 ± 0.1)L of simulation PT2-2 with NSM = 503, and without gravity.

Open with DEXTER
In the text
thumbnail Fig. 38

Comet mass-distribution snapshots of simulation PT2-2 with NSM = 503, and without gravity.

Open with DEXTER
In the text
thumbnail Fig. 39

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of sufficiently self-gravitating simulation PT1-3 with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 40

Comet mass-distribution snapshots of sufficiently self-gravitating simulation PT1-3 with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 41

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of weakly self-gravitating simulation PT2-2 with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 42

Comet mass distribution of weakly self-gravitating simulation PT2-2 with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 43

0.025 × 1 × 1 L3 snapshots at x = (0.5 ± 0.0125)L of weakly self-gravitating simulation PT3-2 with NSM = 1003.

Open with DEXTER
In the text
thumbnail Fig. 44

Comet mass distribution of weakly self-gravitating simulation PT3-2 with NSM = 1003.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text 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 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.