Planetesimal formation

dustpylib.planetesimals.formation contains several prescriptions for planetesimal formation. They are briefly described here. For more detailed discussions, please have a look at their respective publications, which are referenced here.

Drążkowska et al. (2016)

In Drążkowska et al. (2016) only pebbles contribute to planetesimal formation. Pebbles are defined as particles with sizes above a certain cristical Stokes number \(\mathrm{St}_\mathrm{crit}\). In the default case \(\mathrm{St}_\mathrm{crit}=0.01\).

\(\Sigma_\mathrm{peb} = \Sigma_\mathrm{dust} \left( \mathrm{St} \geq \mathrm{St}_\mathrm{crit} \right)\)

If the midplane pebbles-to-gas ratio is above a critical threshold, planetesimal formation is triggered. In the default case, this critical ratio is \(1\).

If planetesimal formation is triggered, a fraction \(\zeta\) of the pebbles is converted into planetesimals per orbit, where \(\zeta=0.01\) in the default case.

Implementing this into DustPy would work as follows. First we define the hyperparameters of the method. These are the critical pebbles-to-gas ratio, the critical Stokes number and \(\zeta\).

p2g_crit = 1.
St_crit = 0.01
zeta = 0.01

Then we create a simulation object and initialize it.

from dustpy import Simulation

s = Simulation()

s.initialize()

In the next step we add an updater to the external sources of the dust surface density, that removes dust if planetesimal formation is triggered with the the method above. Here we can import the drazkowska2016() function, which returns the source term of dust due to planetesimal formation. If will be negative if planetesimal formation is triggered.

from dustpylib.planetesimals.formation import drazkowska2016

def S_ext(s):
    return drazkowska2016(
        s.grid.OmegaK,
        s.dust.rho,
        s.gas.rho,
        s.dust.Sigma,
        s.dust.St,
        p2g_crit=p2g_crit,
        St_crit=St_crit,
        zeta=zeta
    )

s.dust.S.ext.updater = S_ext

Then we are going to add a group for the planetesimals and a field to store their surface density.

import numpy as np

s.addgroup("planetesimals", description="Planetesimal quantities")
s.planetesimals.addfield("Sigma", np.zeros_like(s.gas.Sigma), description="Surface density of planetesimals [g/cm²]")

We are going to let DustPy integrate this field over time to get the evolution of the planetesimal surface density. Therefore we need to define a derivative of the planetesimal surface density. The derivative is simply the negative sum over the external source terms defined above.

def dSigma_planetesimals(s, t, Sigma_planetesimals):
    return -s.dust.S.ext.sum(-1)

This function is added to the differentiator of the field.

s.planetesimals.Sigma.differentiator = dSigma_planetesimals

In the next step we have to create an integration instruction of this field. We are going to integrate the planetesimal surface density with a simple explicit 1st-order Euler scheme.

from simframe import Instruction
from simframe import schemes

instruction = Instruction(schemes.expl_1_euler, s.planetesimals.Sigma, description="Planetesimals: explicit 1st-order Euler")

This instruction is added to the existing integration instructions.

s.integrator.instructions.append(instruction)

s.integrator.instructions

[Instruction (Dust: implicit 1st-order direct solver),
 Instruction (Gas: implicit 1st-order direct solver),
 Instruction (Planetesimals: explicit 1st-order Euler)]

We can now update the simulation object.

s.update()

The simulation is now ready to go and can be started with s.run().

Note, that this setup will not create any planetesimals, since the conditions for planetesimal formation will never be fulfilled without any mechanism to concentrate dust above the given threshold.

Schoonenberg et al. (2018)

The prescription of Schoonenberg et al. (2018) is very similar to Drążkowska et al. (2016). But instead of only considering pebbles, all particles contribute to planetesimal formation.

If the midplanet dust-to-gas ratio is above a critical threshold, planetesimal formation is triggered. The critical threshold is \(1\) in the default case. As soon as planetesimal formation is triggered, a fraction of \(\zeta\) of the dust particles is converted to planetesimal per settling time scale

\(t_\mathrm{sett} = \frac{1}{\mathrm{St}\Omega_\mathrm{K}}\).

Since the settling time scale is inversly proportional to the Stokes number of the particles, the contribution of small particles to planetesimal formation is supressed. In the default case the planetesimal formation efficiency is \(\zeta=0.1\)

Implementation into DustPy is identical to the example above. First, we define the hyperparameters of the method.

d2g_crit = 1.
zeta = 0.1

Then we create a simulation object and initialize it.

from dustpy import Simulation

s = Simulation()

s.initialize()

In the next step we add an updater to the external sources of the dust surface density, that removes dust if planetesimal formation is triggered with the the method above. Here we can import the schoonenberg2018() function, which returns the source term of dust due to planetesimal formation. If will be negative if planetesimal formation is triggered.

from dustpylib.planetesimals.formation import schoonenberg2018

def S_ext(s):
    return schoonenberg2018(
        s.grid.OmegaK,
        s.dust.rho,
        s.gas.rho,
        s.dust.Sigma,
        s.dust.St,
        d2g_crit=d2g_crit,
        zeta=zeta
    )

s.dust.S.ext.updater = S_ext

Then we are going to add a group for the planetesimals and a field to store their surface density.

import numpy as np

s.addgroup("planetesimals", description="Planetesimal quantities")
s.planetesimals.addfield("Sigma", np.zeros_like(s.gas.Sigma), description="Surface density of planetesimals [g/cm²]")

We are going to let DustPy integrate this field over time to get the evolution of the planetesimal surface density. Therefore we need to define a derivative of the planetesimal surface density. The derivative is simply the negative sum over the external source terms defined above.

def dSigma_planetesimals(s, t, Sigma_planetesimals):
    return -s.dust.S.ext.sum(-1)

This function is added to the differentiator of the field.

s.planetesimals.Sigma.differentiator = dSigma_planetesimals

In the next step we have to create an integration instruction of this field. We are going to integrate the planetesimal surface density with a simple explicit 1st-order Euler scheme.

from simframe import Instruction
from simframe import schemes

instruction = Instruction(schemes.expl_1_euler, s.planetesimals.Sigma, description="Planetesimals: explicit 1st-order Euler")

This instruction is added to the existing integration instructions.

s.integrator.instructions.append(instruction)

s.integrator.instructions

[Instruction (Dust: implicit 1st-order direct solver),
 Instruction (Gas: implicit 1st-order direct solver),
 Instruction (Planetesimals: explicit 1st-order Euler)]

We can now update the simulation object.

s.update()

The simulation is now ready to go and can be started with s.run().

Note, that as above this setup will not create any planetesimals, since the conditions for planetesimal formation will never be fulfilled without any mechanism to concentrate dust above the given threshold.

Miller et al. (2021)

The presciption of Miller et al. (2021) is very similar to Schoonenberg et al. (2018). But instead of a hard threshold in midplane dust-to-gas ratio it employs a smoooth transition. The probability that planetesimal formation is triggered is given by

\(\mathcal{P} = \frac{1}{2} \left[ 1 + \tanh \left( \frac{\log\left( \varepsilon \right) - \log\left( \varepsilon_\mathrm{crit} \right)}{n} \right) \right]\),

where \(\varepsilon\) and \(\varepsilon_\mathrm{crit}\) are the midplane dust-to-gas ratio and its critical value at which this transition occurs (the default is \(\varepsilon_\mathrm{crit}=1\)) and \(n\) is a smoothness parameters (the default is \(n=0.03\)).

Setting up this prescription in DustPy works identical to the examples above. First, we define the hyperparameters of the method.

d2g_crit = 1.
n = 0.03
zeta = 0.1

Then we create a simulation object and initialize it.

from dustpy import Simulation

s = Simulation()

s.initialize()

In the next step we add an updater to the external sources of the dust surface density, that removes dust if planetesimal formation is triggered with the the method above. Here we can import the miller2021() function, which returns the source term of dust due to planetesimal formation. If will be negative if planetesimal formation is triggered.

from dustpylib.planetesimals.formation import miller2021

def S_ext(s):
    return miller2021(
        s.grid.OmegaK,
        s.dust.rho,
        s.gas.rho,
        s.dust.Sigma,
        s.dust.St,
        d2g_crit=d2g_crit,
        n=n,
        zeta=zeta
    )

s.dust.S.ext.updater = S_ext

Then we are going to add a group for the planetesimals and a field to store their surface density.

import numpy as np

s.addgroup("planetesimals", description="Planetesimal quantities")
s.planetesimals.addfield("Sigma", np.zeros_like(s.gas.Sigma), description="Surface density of planetesimals [g/cm²]")

We are going to let DustPy integrate this field over time to get the evolution of the planetesimal surface density. Therefore we need to define a derivative of the planetesimal surface density. The derivative is simply the negative sum over the external source terms defined above.

def dSigma_planetesimals(s, t, Sigma_planetesimals):
    return -s.dust.S.ext.sum(-1)

This function is added to the differentiator of the field.

s.planetesimals.Sigma.differentiator = dSigma_planetesimals

In the next step we have to create an integration instruction of this field. We are going to integrate the planetesimal surface density with a simple explicit 1st-order Euler scheme.

from simframe import Instruction
from simframe import schemes

instruction = Instruction(schemes.expl_1_euler, s.planetesimals.Sigma, description="Planetesimals: explicit 1st-order Euler")

This instruction is added to the existing integration instructions.

s.integrator.instructions.append(instruction)

s.integrator.instructions

[Instruction (Dust: implicit 1st-order direct solver),
 Instruction (Gas: implicit 1st-order direct solver),
 Instruction (Planetesimals: explicit 1st-order Euler)]

We can now update the simulation object.

s.update()

The simulation is now ready to go and can be started with s.run().

Note, that as above this setup will not create any planetesimals, since the conditions for planetesimal formation will never be fulfilled without any mechanism to concentrate dust above the given threshold.