Users guide, Part II: Parametrized coagulation simulations

1. Coagulation of naked grains

V. Ossenkopf, MPG AG “Dust in Star-Forming Regions” Jena

Abstract

SIMNAKED is a program for the simulation of the coagulation of dust grains without ice mantles. It solves the coagulation equation by an Eulerian integration in conditions typical for the interior of a cold molecular cloud core. This means that the physical effects of thermal and turbulent motions, relative grain velocities due to accidential grain asymmetries, grain charges and a finite critical velocity for grain sticking are implemented. The dust particles are characterized by two parameters - the compact volume and an internal density parameter. The evolution of the dust distribution over these two parameters is followed at any fixed point in the cloud. No dust transport is incorporated. The resulting dust distributions can be stored after arbitrary time intervals. Any resulting distribution can be used as starting distribution of further computations. The theoretical basis for the simulation program is discussed in Ossenkopf (1993, 1994).

1 Basic algorithms

1.1 Analytical formulae and approximations used

The physical background, the astrophysical conditions and the principal details of the program are published by Ossenkopf (1993,1994).

The goal of the program is to integrate the coagulation equation:

\frac{∂f(V, r, t)}{∂t} = −f(V, r, t) ∫ dV′ ∫ dr′ f(V′, r′, t)Q(V, r, V′, r′)
+1/2 ∫ dV′ ∫ dr′ ∫ dV′′ ∫ dr′′f(V′, r′, t)f(V′′, r′′, t)
Q(V′, r′, V′′, r′′) δ(V − Γ_{V} (V′, V′′))δ(r − Γ_{r}(r′, r′′))  (1)

where f denotes the number density distribution. The particles are characterized by two parameters: the volume V and an internal density parameter r. Q is the coagulation probability and Γ are the parameters of the newly formed particles.

The set of new parameters for a coagulation product is computed under the assumption that a nondestructive ballistic coagulation happens. The new compact volume of a coagulation product is simply given as

Γ_{V} = V_{1} + V_{2}  (2)

The new density parameter r is computed from

Γ_{r} = [\frac{r_{1}V_{1} + r_{PCA}(V_{1}n_{sub,2}/V_{2})V_{2}}{Γ_{V}} ] (\frac{V_{1}}{Γ_{V}} )^{ex}  (3)

where the indices 1 and 2 denote the larger and the smaller collision partner. The quantity r_{PCA}(n_{sub}) is the density parameter for PCA particles consisting of n_{sub} subgrains of equal size

r_{PCA}(n_{sub}) = [a_{0} exp(−b/n_{sub}^{γ})]^{−1}  (4)
a_{0} = 15.17, b = 2.855, γ = 0.0963 for n_{sub} ≤ 1500;
a_{0} = 4.70, b = 9.016, γ = 0.50328 for n_{sub} > 1500;

and inversely the number of subgrains of a particle is computed from

n_{sub}(r) = (\frac{ln(a_{0} × r)}{b})^{−1/γ}  (5)

with the first parameter set of Eq. (4) for r ≥ 0.27 and the second set for r < 0.27.

The exponent of the motion in the V − r−diagram ex is computed from

ex = (\frac{d ln r(n_{sub})}{d ln V})_{CCA}  (6)

+ ((\frac{d ln r(n_{sub})}{d ln V})_{PCA}− (\frac{d ln r(n_{sub})}{d ln V})_{CCA}) × (\frac{a_{tooth,1}/a_{tooth,2} − 1}{a_{tooth}(n_{sub}, r_{PCA}(n_{sub}))/a_{tooth}(1, 1) − 1})

The coefficients are given as

(\frac{d ln r(n_{sub})}{d ln V})_{CCA} = 1/3 − 0.05 − \frac{0.30072}{ln^{2}(n_{sub}) + 0.30072 ln(n_{sub})}  (7)
(\frac{d ln r(n_{sub})}{d ln V})

PCA

=

\frac{bγ}{n_{sub}^{γ}}

(8)

where b and γ are to be taken from Eq. (4).

The toothing radii a_{tooth} of the particles are computed from

a_{tooth}(V, r) = 1.276 \frac{V^{1/3}}{r^{0.63}} (1 − 0.216r) .  (9)

where the volume V has to be taken in volume units V_{0} of the smallest considered grains (a = 5 nm) and the toothing radius is obtained in units of the radius of these grains.

The coefficients for the exponent of the motion of the particles in the V − r−diagram (Eq. (7,8)) are not computed during the integration of the coagulation equation but a table for n_{sub} between 1 and 50000 is created at the start of the program which is only read by the routine computing the new parameters of a coagulation product.

A separate treatment is performed for collisions between compact grains. Here, the analytic approximations for the behaviour of PCA and CCA particles are not yet valid but a direct geometric approach is possible. We obtain

Γ_{r} = \frac{Γ_{V}^{2/3}}{A_{1} + 2/3 A_{2}\sqrt{\frac{2a_{tooth,2}}{a_{tooth,1} + a_{tooth,2}}}}  (10)

The letter A stands for the average projected area of a particle which is exactly

A = \frac{V^{2/3}}{r}  (11)

(A and V in surface and volume units, A_{0} and V_{0} respectively.) The indices 1 and 2 again denote the larger and smaller particle.

The second important quantity in the coagulation equation is the coagulation probability Q. It is given as

Q(r, V, r′, V′) = ∫ dv_{rel} p_{v_{rel}}(V, r, V′, r′) A_{coll}(r, V, r′, V′, v_{rel}) v_{rel} S(V, r, V′, r′, v_{rel}) .  (12)

where the collision cross section A_{coll} and the sticking efficiency S may explicitely depend on the relative velocity v_{rel}. p_{v_{rel}} is the normalized distribution of relative velocities between two grains. All relative velocities are counted positive.

For the sticking efficiency S, we assumed a Heaviside function in the form

S(V, r, V′, r′) = \{ \begin{matrix}1 & for v_{rel}(V, r, V′, r′) ≤ v_{stick}(V, r, V′, r′) \\ 0 & for v_{rel}(V, r, V′, r′) > v_{stick}(V, r, V′, r′) \end{matrix} .  (13)

For the critical sticking velocity Chokshi et al. (1993) derived

v_{stick} = 4.39\frac{γ^{5/6}a_{r}^{2/3}}{E^{1/3}m_{r}^{1/2}}  (14)

where γ is the surface energy per unit area and E the Young’s modulus of the material. The quantities a_{r} and m_{r} denote the reduced radius and mass of the two grains. The program uses a simplification by replacing the reduced quantities by the lower one of the quantities of the two particles.

v_{stick} = 0.147 \frac{min(a_{sub,1}, a_{sub,2})^{2/3}}{min(V_{1}ρ_{1}, V_{2}ρ_{2})^{1/2}} × \frac{cm^{1/3}g^{1/2}}{s}  (15)

where the E and γ values of ice are assumed. The radius of the subgrains in a cluster a_{sub} is determined by its total compact volume V and the number of subgrains n_{sub}

a_{sub} = (\frac{3V}{4πn_{sub}} )^{1/3} .  (16)

The average number of subgrains is computed from Eq. (5) assuming a PCA structure. This can lead to artificially small subgrains, so that we set a_{sub} to the lowest radius of the original particles whenever Eq. (16) gives lower values.

The integral in Eq. (12) is treated in two ways. In a high velocity regime, the distribution of relative velocities p(v_{rel}) covers a range where the collision cross section A_{coll} does not depend on v_{rel}, but which contains the sticking velocity. Then, Q is obtained from

Q(r, V, r′, V′) = A_{coll}(r, V, r′, V′) ∫_{0}^{v_{stick}} dv_{rel} p_{v_{rel}}(V, r, V′, r′) v_{rel}  (17)

In a low velocity regime S is always 1 but A_{coll} depends on v_{rel}. We compute Q from

Q(r, V, r′, V′) = ⟨v_{rel}(r, V, r′, V′)⟩⟨A_{coll}(r, V, r′, V′)⟩_{v}  (18)

with

⟨v_{rel}(r, V, r′, V′)⟩ = ∫ dv_{rel} p_{v_{rel}}(V, r, V′, r′) v_{rel}  (19)
⟨A_{coll}(r, V, r′, V′)⟩_{v} = \frac{∫ dv_{rel} p_{v_{rel}}(V, r, V′, r′) v_{rel} A_{coll}(r, V, r′, V′, v_{rel})}{⟨v_{rel}(r, V, r′, V′)⟩}  (20)

The high velocity regime is applied when the average relative velocity ⟨v_{rel}(r, V, r′, V′)⟩ is above 0.4 times the sticking velocity v_{stick}(V, r, V′, r′). Otherwise, the low velocity equation is applied.

The criterion for the high velocity regime is only fulfilled in case of large turbulent velocities. Their velocity distribution is almost Maxwellian:

p(v_{rel}) = \frac{4}{α^{3}π^{1/2}} v_{rel}^{2} exp (−v_{rel}^{2}/α^{2})  (21)

with

α = \sqrt{π}/2 × ⟨v_{rel}(r, V, r′, V′)⟩ .  (22)

Consequently, the integral on the right hand side of Eq. (17) can be computed

⟨v_{rel}⟩_{eff} = ⟨v_{rel}⟩ × (1 − [1 + \frac{v_{stick}^{2}}{α^{2}} ] exp (−\frac{v_{stick}^{2}}{α^{2}} ))  (23)

For relative velocities ⟨v_{rel}⟩ > 5 v_{stick}, the fraction of sticking grains is vanishing so that ⟨v_{rel}⟩_{eff} is set to 0.

Three different processes contribute to the relative grain velocities. The average relative velocity ⟨v_{rel}⟩ between two grains is determined from the relative velocities produced by these processes by

⟨v_{rel}⟩ = \sqrt{∑_{i=1}^{3}⟨v_{rel,i}⟩^{2}} .  (24)

The first considered velocity is the thermal one

⟨v_{therm}⟩ = \sqrt{\frac{8kT (V_{1}ρ_{1} + V_{2}ρ_{2})}{π V_{1}ρ_{1}V_{2}ρ_{2}}} .  (25)

where we took T = 20K, so that kT = 2.76132e-15 g cm^{2}/s^{2}.

The second one is given by turbulent motions. The relative velocity between two grains is approximately given by the surface to mass ratio of the larger one:

⟨v_{turb}⟩ = max (\sqrt{\frac{V ρ}{A}}) \sqrt{\frac{f_{v}^{3}}{f_{l}}} (\frac{kT G}{n_{H} m_{H}^{2}} )^{1/4}  (26)

where G is the gravitational constant, n_{H} the hydrogen density and m_{H} the hydrogen mass, i.e. the constant quantities give (kT G/m_{H}^{2})^{1/4} = 2.84875e6 cm^{5/4}/(s g^{1/2}). The filling factors f_{v} and f_{l} are assumed to be 1.

There is a cut-off of turbulent relative velocities for particles with average velocities below

⟨v_{crit}⟩ = \frac{f_{v}}{Re^{1/4}} \sqrt{\frac{2kT}{μ m_{H}}} .  (27)

Here, we need the Reynolds number Re = L_{J}n_{H}σf_{v}/μ, the relative molecule mass μ, the Jeans length L_{J} and the H_{2} cross section σ. Altogether, we have ⟨v_{crit}⟩ = 1900.59 cm^{5/8}s^{−1} /n_{H}^{1/8}. The effective turbulent relative velocities of grains where Eq. (26) gives lower velocities are 0.

The third type of velocities comes from accidential grain asymmetries. The relative velocity between two grains is given by

⟨v_{asymm}⟩ = \sqrt{(v_{0,1} ⟨|γ|⟩_{1})^{2} + (v_{0,2} ⟨|γ|⟩_{2})^{2}}  (28)

with v_{0} × ⟨|γ|⟩ = (1 − \frac{π}{4} ) × \sqrt{\frac{2kT_{g}}{πμ}} (1 − \sqrt{\frac{T_{d}}{T_{g}}} ) δα × \sqrt{\frac{2}{πn}} .  (29)

Here, T_{d} and T_{g} are the dust and gas temperature, δα is the variation of the accomodation coefficient, μ is the average molecule mass, and n is the number of exceptional sites on a grain. The prefactor (1 − π/4) comes from the precession of the grains in a magnetic field.

The total number of exceptional surface sites and their variation of the accomodation coefficient relative to the average value α are given by

n = n_{H_{2}} + n_{sub} − 1  (30)
δα = \frac{0.03n_{H_{2}} + 0.15A_{sub}/A_{0}(n_{sub} − 1)}{n_{H_{2}} + A_{sub}/A_{0}(n_{sub} − 1)}  (31)
with n_{H_{2}} = 10\frac{A}{A_{0}} ,  (32)

where n_{sub} is the number of subgrains in a cluster, A_{sub} = πa_{sub}^{2} is the mean projected area of the subgrains and A_{0} is this area for compact 5-nm grains.

The temperature ratio between gas and dust temperature is obtained from the equation

2An_{g}v_{g}αkT_{g}(1 − T_{d}/T_{g}) = 4AσT_{d}^{4}⟨Q(T_{d}, a)⟩_{Pl} ,  (33)

where the quantity σ denotes the Stefan-Boltzmann constant and ⟨Q(T_{d}, a)⟩_{Pl} is the Planck-averaged emission or absorption efficiency of the dust particles. Since no exact solution is possible, we have approximated it by a simple function:

1 − \sqrt{T_{d}/T_{g}} = \frac{1 − X^{1/12}/(1 + X)}{1 + 2 X}  (34)

where X = \frac{n_{g}v_{g}αkT_{g}}{2σ Q_{0}(3V/4π)^{1/3}} × \frac{1}{T_{g}^{6}}

The temperature dependent absorption efficiency is assumed to vary by Q_{0} = 2 10^{−5} μm^{−1} K^{−2} and for the average accomodation coefficient α a value 0.3 was chosen. Altogether, this gives X = 0.032461 cm^{4}× n_{H}/V^{1/3}.

Whether asymmetry velocities may really occur depends on their amount relative to the thermal velocity. Only for particles with asymmetry energies above 3/2 kT, the process is efficient. This is due to the equivalence of rotational and translational asymmetry velocities and the fact that a suprathermal rotation is required for the exitation of translational asymmetry velocities. is necessary to stabilize a dust particle For the whole distribution of velocities, the average asymmetry velocitiy is reduced by the factor (1 − f_{ex}) with

f_{ex} = Erf (\sqrt{\frac{3kT_{g}}{πρV}} \frac{1}{v_{0} × ⟨|γ|⟩}) = Erf (\sqrt{\frac{3}{2}} \frac{v_{therm}}{v_{0} × ⟨|γ|⟩}) .  (35)

The error function is replaced by a linear function Erf(x) ≈ 2/\sqrt{π} x. For f_{ex} > 0.9, the contribution of asymmetry velocities is completely neglected.

The influence of electric charges is expressed by the integral (20) since charges on the dust grains do not change their average relative velocities but their collision cross sections ⟨A_{coll}⟩_{v} depending on the relative velocity v_{rel}. The integral has to be computed for different distributions of relative velocities. This is expressed in four different routines computing the average collision cross section between two particles depending on the relative velocity regime and the charge distribution on the grains.

The first case is the computation of the collision cross section with the neglect of grain charges, so that Eq. (20) gives exactly the geometric collision cross section and the velocity distribution plays no role. The collision cross section of two neutral particles is computed from the projected area of the larger one and the ratio of the toothing radii

A_{coll} = A_{1} (1 − [\frac{a_{tooth,2}}{a_{tooth,1}} ]^{2}) + (4π − 8.3 (1 − [\frac{a_{tooth,1}}{a_{tooth,2}} ]^{1.22})) a_{tooth,2}^{2}  (36)

where the subscripts 1 and 2 denote the particles with the larger and the smaller toothing radii, respectively. This value of the geometric corss section is always used in the high velocity regime where the coagulation rate is computed from Eq. (17). Moreover, the other routines for the computation of the collision cross section also use this value for neutral grains since they compute the changes relative to this value.

The second case assumes a thermal distribution of relative velocities, i.e. a Maxwellian velocity distribution with an average energy kT. This case is applied when the asymmetry velocities are negligible (see above) and the square of the ratio between the relative thermal and turbulent velocities is above 19.

The third case consideres a general Maxwellian velocity distribution. This is characterized by an effective temperature T_{eff} = T ×(1+⟨v_{turb}⟩^{2}/⟨v_{therm}⟩^{2}). It is applied when only the asymmetry velocities are negligible. The splitting between the second and the third case is done due to the possibility to precompute more relevant quantities in case of a fixed effective temperature before the coagulation integration is carried out.

The fourth case is the most general but most time consuming case. It represents the superposition of a Maxwellian and a Gaussian velocity distribution.

In all cases which take the grain charging into account, the total average collision cross section has to be computed as

⟨A_{coll}⟩ = ∑_{Z_{1}} ∑_{Z_{2}} f_{Z_{1}}f_{Z_{2}}⟨A_{coll}(Z_{1}, Z_{2})⟩  (37)

where the sum runs over all possible charging states Z_{1} and Z_{2} of the collision partners and the quantities f_{Z_{1}} and f_{Z_{2}} denote the fraction of particles in each state. ⟨A_{coll}(Z_{1}, Z_{2})⟩ is the average collision cross section according to Eq. (20) for definite particle charges Z_{1} and Z_{2}.

For all cases the distribution of electric charges on the grains, i.e. the fraction of grains in each charging state, is computed in the same way. At first two quantities characterizing the charges in the gas phase are computed depending on the gas density. We have approximated the numerical values for the ratio between the electron and the ion density in the gas pahse obtained by Nishi et al. (1991) by

μ := \frac{n_{e}S_{e}}{n_{i}S_{i}} (\frac{m_{i}}{m_{e}} )^{1/2} = \{ \begin{matrix}126.8 & if n_{H} < 10^{6.5} cm^{−3}; \\ 126.8 × (n_{H}/10^{6.5} cm^{−3})^{−4/5} & if 10^{6.5} cm^{−3} < n_{H} < 10^{9} cm^{−3}; \\ 1.268 & if n_{H} > 10^{9} cm^{−3}. \end{matrix}  (38)

The quantity ψ which is defined as the solution of 1 − ψ = μ exp(ψ) is approximated by

ψ = 0.04555 − 1.6229 lg(μ)  (39)

For every value of ψ, a table of charge distributions is computed, which contains the fraction of particles in the charging states +1, 0, -1, -2, -3, -4 depending on the average particle charge. According to Draine & Sutin (1988) Gaussian distributions with the FWHM

σ = \sqrt{(1 − ψ)⟨−Z⟩/(2 − ψ)(−ψ)} − 0.12⟨−Z⟩^{−3/2}  (40)

are assumed. The small correction term outside of the sqare root is introduced to obtain a smooth transition to the behaviour at average particle charges ⟨−Z⟩ below 1.

The average particle charge is computed from the electron and ion currents onto the grains. They are determined by the ratio of the collision cross sections of a dust grain for infinitely small particles with different charges. We obtain

⟨−Z⟩ = f_{−1} − ψ × τ − 0.2 ,  (41)

with

f_{−1} = [1 + \frac{J(−1)}{J(0)μ}]^{−1}

f_{0} = [1 + \frac{J(0)μ}{J(−1)}]^{−1}
f_{+1} = f_{−1}/μ^{2}  (42)

and

J(0) = (1 + a_{add}/a_{tooth})^{2}  (43)
J(−1) = 1 + 1/τ  (44)

where

a_{add} = a_{sub} [\sqrt{1 + \sqrt{\frac{π}{2τ_{sub}}}} − 1]  (45)

and

τ = \frac{4πϵ_{0} a_{tooth} kT}{e^{2}}  (46)
τ_{sub} = \frac{4πϵ_{0} a_{sub} kT}{e^{2}}  (47)

If ⟨−Z⟩ is lower than 1, the fraction of particles in each charging state is not taken from the table but given by the quantities f_{−1}, f_{0} and f_{+1}. The fraction of particles in all other charging states is negligible in this case. The quantities a_{tooth} and a_{sub} denote the toothing radius and the average subgrain radius of a given (coagulated) particle. The quantities a_{add}, τ and τ_{sub} will also be used to determine the change of the grain collision cross sections due to particle charges.

In the computation of the average collision cross section all three routines which consider the particle charges follow the same structure. At first the cross section for collisions with one neutral and one charged partner is computed. Here, the radius difference a_{add} is computed by which the radius of the smaller particle is virtually enlarged due to the induced dipole interaction. From the projected area of the larger particle and the newly computed radius of the smaller one the collision cross section is computed using Eq. (36).

This part consists of two steps. It is checked whether the low or the high x_{lim} formulae have to be taken by comparing the minimal distance from the induced dipole a_{min} to the minimal distance between two grains given by their radii a_{tooth}. For the three velocity regimes, a_{min} is given by

a_{min} = \{ \begin{matrix}a_{sub} [\sqrt{1 + \sqrt{\frac{1}{2τ_{sub}}}} − 1] \\ a_{sub} [\sqrt{1 + \sqrt{\frac{1}{2τ_{sub}} \frac{T}{T_{eff}}}} − 1] \\ a_{sub} [\sqrt{1 + \sqrt{\frac{1}{2 × 0.7145τ_{sub}}}} − 1] \end{matrix}  (48)

This has to be computed twice since both partners may be neutral. The quantities a_{tooth} and τ have to be taken for the charged particle and a_{sub} and τ_{sub} for the neutral particle.

If the influence of the induced dipole is dominant (a_{min} > a_{tooth}, low x_{lim}-regime), the additional radius by which the surface of the smaller particle has to be shifted outwards is given by

a_{add} = a_{x} − a_{tooth}  (49)

with

a_{x} = \{ \begin{matrix}a_{sub} [\sqrt{1 + \sqrt{\frac{π}{2τ_{sub}}}} − 1] \\ a_{sub} [\sqrt{1 + \sqrt{\frac{π}{2τ_{sub}} \frac{T}{T_{eff}}}} − 1] \end{matrix}  (50)

for the two Maxwellian velocity distributions. In case of the superposition of Maxwellian and Gaussian distributions, we have

a_{x} = \sqrt{(a_{tooth} + a_{sub})^{2} + \frac{a_{tooth}^{2} + 2a_{tooth}a_{sub}}{g} \{ \begin{matrix}3.5449 + \sqrt{g} − g & for g < 1.56046 \\ 2.4228 + ln(g) & for g > 1.56046 \end{matrix}} − a_{sub} (51)

with

r = \frac{v_{therm}^{2} + v_{turb}^{2}}{4v_{asymm}^{2}}  (52)
g = \frac{2τ_{sub}(a_{tooth}^{2} + 2a_{tooth}a_{sub})^{2}}{ra_{sub}^{4}}  (53)

In the high x_{lim}-regime (a_{min} > a_{tooth}) the additional radius a_{add} can be written as

a_{x} = −a_{sub} + \sqrt{(a_{tooth} + a_{sub})^{2} + \frac{a_{tooth}^{2} + 2a_{tooth}a_{sub}}{g}}  (54)

when one takes g from Eq.(53) and substitutes r by 1, by T/T_{eff}, and by r(2.3533 − ln(r)) with r from Eq. (52) for the three velocity distributions.

Then, this radius a_{add} is added to the toothing radius of the smaller collision partner and the resulting collision cross section A_{coll}(Z_{1} × Z_{2} = 0) is computed from Eq. (36).

In the collision of partners none of which is neutral, the geometric collision cross section is reduced or enhanced by a factor J(Z_{1}, Z_{2}) = A_{coll}(Z_{1}, Z_{2})/A_{coll}(neutral) given as

⟨J(Z_{1} × Z_{2} < 0)⟩ = 1 − \frac{1}{4πϵ_{0}} \frac{Z_{1}Z_{2}e^{2}}{(a_{tooth,1} + a_{tooth,2})⟨E⟩}
⟨J(Z_{1} × Z_{2} > 0)⟩ = exp (−\frac{1}{4πϵ_{0}} \frac{Z_{1}Z_{2}e^{2}}{(a_{tooth,1} + a_{tooth,2})⟨E⟩}) .  (55)

according to the charge situation for the two Maxwellian velocity distributions. In case of the superposition of Gaussian and Maxwellian velocities we have

⟨J(Z_{1} × Z_{2} < 0)⟩ = exp(r) + 5.2778r [0.26424 − 0.63212f_{Z}] + rf_{Z}Ei(−r)  (56)
⟨J(Z_{1} × Z_{2} > 0)⟩ = exp(rf_{Z}) + rf_{Z}Ei(−rf_{Z})  (57)
with  f_{Z} = \frac{1}{4πϵ_{0}} \frac{Z_{1}Z_{2}e^{2}}{(a_{tooth,1} + a_{tooth,2})kT} .  (58)

and r from Eq. (52).

This is treated in the next part of the routines to compute the integral (20). At first, the contributions from the terms with both neutral particles, one neutral particle and a charged one (calculated above), and particles with different signs of the charges to the average collision cross section are added up. The contribution for particles with different signs of the charge may be computed straight forward, since only one positive charge may occur. The additional cross section due to the attractive force is proportional to the particle charge. Hence, we obtain the contributions

f_{0}^{(1)}f_{0}^{(2)}A_{coll}(neutral)  (59)

f_{0}^{(1)} ∑_{Z≠0} |Z|f_{Z}^{(2)} (A_{coll}(Z_{1} = 0) − A_{coll}(neutral)) + ∑_{Z≠0} f_{Z}^{(2)}A_{coll}(neutral)

+f_{0}^{(2)} ∑_{Z≠0} |Z|f_{Z}^{(1)} (A_{coll}(Z_{2} = 0) − A_{coll}(neutral)) + ∑_{Z≠0} f_{Z}^{(1)}A_{coll}(neutral)  (60)

f_{+1}^{(1)} ∑_{Z<0} f_{Z}^{(2)} (1 + |Z|/(τ̃_{1} + τ̃_{2})) A_{coll}(neutral)

+f_{+1}^{(2)} ∑_{Z<0} f_{Z}^{(1)} (1 + |Z|/(τ̃_{1} + τ̃_{2})) A_{coll}(neutral)  (61)

to the sum (37). Herein, the quantity τ̃ is exactly the τ given by Eq. (46) for the thermal velocities, it is τ × T_{eff}/T for the general Maxwell distribution of velocities, and it is 0.23406τ/r with r from Eq.

(52) for the superposition of Gaussian and Maxwellian velocities.

Finally, the contribution of collisions where both particles are negatively charged are added. Starting from the factor ⟨J(Z_{1} × Z_{2} = 1)⟩ = exp (−1/(τ̃_{1} + τ̃_{2})), a series of contributions with growing Z_{1} × Z_{2}- values is iteratively created from ⟨J(Z_{1}×Z_{2} = i+1)⟩ = ⟨J(Z_{1}×Z_{2} = i)⟩×⟨J(Z_{1}×Z_{2} = 1)⟩ and summed up to the total collision cross section for the whole charge distribution on both collision partners.

For the superposition of Gaussian and Maxwellian velocities two different approximations are used to avoid a direct computation of the error function Ei during the integration. For r/(τ_{1} + τ_{2}) < 1.7 the quantity τ̃ is τ/(2.6r) and the prefactor for the geometric collision cross section is corrected to A_{coll}(neutral) × (1 + r^{2}/(τ_{1} + τ_{2})^{2}). For r/(τ_{1} + τ_{2}) ≥ 1.7 we use τ̃ = τ/(1.13328r) and A_{coll}(neutral) × 0.31837. Due to their low number and the reduced collision cross section the contribution from collisions between two positively charged grains is neglected.

1.2 Numerical details

To solve the coagulation equation numerically a binning method was chosen. According to the twodimensional description of the dust aggregates, the numerical integration was carried out on a grid with 200 points in V -direction and 40 points in r-direction. The points were chosen to be equidistant in the logarithm of the particle volume V and equidistant in the density parameter r.

The contribution of a newly formed particle with coordinates V and r somewhere between the grid points is distributed over the next four neighbouring grid points.

f(V_{i}, r_{j}) = (1 − \frac{V − V_{i}}{V_{i+1} − V_{i}} ) (1 − \frac{r − r_{j}}{r_{j+1} − r_{j}} )  (62)
f(V_{i+1}, r_{j}) = \frac{V − V_{i}}{V_{i+1} − V_{i}} (1 − \frac{r − r_{j}}{r_{j+1} − r_{j}} )  (63)
f(V_{i}, r_{j+1}) = (1 − \frac{V − V_{i}}{V_{i+1} − V_{i}} ) \frac{r − r_{j}}{r_{j+1} − r_{j}}  (64)
f(V_{i+1}, r_{j+1} = \frac{V − V_{i}}{V_{i+1} − V_{i}} \frac{r − r_{j}}{r_{j+1} − r_{j}}  (65)

This allows the treatment of collisions between particles with very different volumes but slightly smooths out the flanks of the distribution.

The coordinate grid is dynamically adjusted to the size of the current distribution. It is tested whether the product of the collision of two particles from that point of the distribution with the highest volume and the lowest density parameter would still fall into the range covered by the coordinate grid. Whenever this fails, the grid is either shifted or when this is not possible it is streched by a factor 2. A grid streching by non-integer numbers has turned out to produce artefacts. Since all density parameters r have to be above 0, the r scale of the initial distribution has to be chosen carefully so that the streching may cover the whole parameter range considered but does not contain 0 or negative values. Whenever an extension of the grid to r-values below 0 would occur, coagulation products with r values below the smallest possible one in the current grid are deposited at this lowest value.

The integration stepwidth is constrained by the requirement that the relative change of the number density of the smallest particles per step is less than a user defined accuracy.

All quantities which are not dimensionless by their own are normalized. The unit number density is 1 per H-atom in the gas, the unit volume is 4π/3(5nm)^{3} = 5.23599e-19 cm^{3}, the unit surface is 7.85398e-13 cm^{2}. The lengths are either treated in units of \sqrt{π} × 5nm = 8.86227e-7 cm or in cm.

2 User interface

2.1 Terminal input

Since all physical processes are defined in the program depending on only one parameter, the gas density, the user has to input this physical quantity and some technical quantities which determine the integration.

At first, the file of the starting distribution has to be given. Since all resulting files are numbered, an offset number can be given when starting from a result file of a former simulation so that the numbering is continued. Then, the user has to define the time interval of the simulation after which a resulting distribution shall be stored. When the integration time has reached this number the distribution is written to the disk and the simulation continues. In this way, a series of files with consecutive numbers is produced. As fourth parameter the user has to define the gas density in the form of hydrogen atoms per cubic centimetre. The next parameter determines the relative accuracy of the integration. It limits the integration step size by the condition that the number of dust grains with the smallest volume may only change by this relative amount during one integration step. I recommend values of 0.002… 0.005. At last the first letters of the output file names have to be defined. The program automatically appends a two digits number for the numbering of consecutive distributions and the extension “.dis”.

2.2 Data files

The program reads and writes the distribution of dust grains f(V, r) as data files. All distribution files have the same structure so that output files may be reused as input files. The files are ASCII files with all numbers in E11.5 format. They consist of a three lines header followed by the field of number densities for all parameter points.

The first header line contains the coagulation time. In other versions of the program it may additionally contain the gas density or the ice molecule density in the gas phase. The second and the third line give the first and the last value of the parameters lg(V) and r in the distribution, i.e. they define the coordinate grid. The body of the file contains the number densities at all 200×40 parameter points. Herein, the density parameter varies in an inner loop and the logarithmic volume in an outer loop. The data are organized in lines with 10 numbers.

3 System requirements

The program is written in standard FORTRAN 77 so that it can be compiled on all machines providing an f77 compiler. It has been tested on a wide variety of machines under Unix and VMS.

The data files of the resulting distributions contain 96863 Byte each. A complete run produces 99 output files, so that together with the input file at least 10 MB hard disk space are necessary to run a complete simulation.

The main storage has to be able to keep 29 fields with 200×40 real*4 numbers and some additional variables what is in total about 1.1 MB. Together with the program code, this gives a minimum storage requirement of 1.6–1.8 MB depending on the compiler efficiency.

4 Limitations and bugs

The applicability of the program is restricted to the physical conditions it is made for. There is no sense in computing a coagulation with the formulae above for gas densities below 10^{4} and above 10^{12} cm^{−3} or for astrophysical environments which considerably differ from cold molecular cloud cores.

The initial range of r values should allow an extension close to zero by streching the grid by factors 2, 4, 8 or so on. An appropriate choice is e.g. r = 1 … 0.501.

There is no error handling when the user selects too high values for the integration stepwidth so that negative points in the number density appear. In this case, the program aborts with an error message.

5 Description of the different parts of the program

5.1 Main program and program control(SIMNAKED, EINGABEN)

These routines represent the outer shell of the program. They call the routines for data input and output, read the control data from the terminal, call the initialization routines, control the integration loop and the storage conditions.

SIMNAKED is the main program. It calls at first the routines for file and terminal input (FREAD, EINGABEN) and then the initialization routines (PRODDA, INITFELDS). Then, the integration loop is started. It contains the test of the correct grid size (NEWGRID), the computation of the appropriate step width (EFFT), and the real one-step integration (INTEGR). Both a time counter and a counter for the number of integration steps run. After a selected time interval or after 150 integration steps, the distribution is stored and the integration continues. The program terminates when 99 distribution files have been stored or when the total integration time exceeds 3.2 10^{12} s.

EINGABEN reads the control data from the terminal. The user has to define the time interval of the simulation after which a resulting distribution is to be stored. Then, the gas density has to be given in the form of hydrogen atoms per cubic centimetre. The next parameter determines the relative accuracy of the integration. It is the fraction by which the most abundant particles may change their number density within one integration step. At last the first letters of the output file names have to be defined. The program automatically appends a two digits number for the numbering of consecutive distributions and the extension “.dis”.

5.2 File operations (FREAD, FSAVE)

The file input/output routines manage the transformation of the distributions between the ASCII file format and the internal representation. They create or resolve the file header and they contain the direct file operations.

FREAD asks for the file name of the initial dust distribution and it reads the coordinate grid data and the distribution from the file. A renormalization of the distribution is done so that the most abundant dust particles are treated with a number density 1. The normalization constant is given back to the main program so that the integration routine can use it. FSAVE is the inverse procedure. It removes the normalization and it creates a file name from the first letters for the output files, the file number and the extension “.dis”. Then, the header and the distribution is written to the output file.

5.3 Initialization of auxiliary quantities for the integration (INITFELDS, INITDIS, PRODDA, DA1, DA2, DAR)

Here, quantities or functions are computed which will be used in the integration kernel but which do not change during the integration so that they may be precomputed. Since the real integration routine has to compute in the order of 200 × 40 × 200 × 40 interaction terms, most subterms which only depend on the parameters of one collision partner are computed before the integration routine is called (only 200 × 40 computations), they are stored and directly used in the integration routine.

This is mainly done in the routine INITFELDS. It computes for every parameter point the quantities: r = V^{2/3}/A, n_{sub} = V/⟨V_{sub}⟩, V, v_{therm}^{2}, v_{turb}, a_{sub}, \sqrt{V ρ}, 0.147 × a_{sub}^{2/3}, a_{tooth}, (a_{tooth})^{1.22}, πa_{tooth}^{2}, τ, π^{2}a_{sub}^{4}/(2τ_{sub}), \sqrt{π/(2τ_{sub})}, \sqrt{1/(2τ_{sub})}, a_{x}(thermal), a_{min}(thermal), a_{min}(Maxwell+Gauss), f_{0}, ⟨−Z⟩, ∑_{Z<0} f_{Z}, ∑_{Z<0} |Z|f_{Z}, f_{+1}, A, 2/3A\sqrt{2\sqrt{π}a_{tooth}}, (v_{0}⟨γ⟩)^{2}, f_{ex}. Additionally, it calls the routine INITDIS which computes for a given ψ a table of charge distributions depending on the average particle charge ⟨−Z⟩. The table contains f_{+1}, f_{0}, f_{−1}, f_{−2}, f_{−3} and f_{−4}. The tenfold average charge is directly used as table index. The table is used to compute ∑_{Z<0} f_{Z} and ∑_{Z<0} |Z|f_{Z} within INITFELDS and the interactions between negatively charged particles within APROB, APROM, and APROG.

The initialization routine PRODDA produces four tables which contain the r-value for the PCA curve, the exponents of the PCA and CCA curve, and the quantity

((\frac{d ln r(n_{sub})}{d ln V})_{PCA}− (\frac{d ln r(n_{sub})}{d ln V})_{CCA}) × (\frac{1}{a_{tooth}(n_{sub}, r_{PCA}(n_{sub}))/a_{tooth}(1, 1) − 1})

depending on the number of subgrains which is used as index. All tables contain entries for 1 to 50000 subgrains. The functions DA1, DA2, and DAR make the tables 3, 4, and 1 available in the form of function calls. They convert a real argument to the next table index and return the value from that table point.

5.4 Integration routine (INTEGR, EFFT)

These routines compute the interaction between all particles and assign the newly formed particles to the points of the distribution. They represent the integration kernel of the program.

INTEGR is the real integration routine which computes the interactions between all particles. In a combined loop running over all parameter points of two collision partners the coagulation probability Q and the parameters of a coagulation product Γ_{V} and Γ_{r} are computed by calling the routines QCOLL and NEU. Then, the position of the coagulation product within the coordinate grid is determined by INDEXM and INDEXR and its contribution is distriubuted over the next four parameter points with the weighting factors determined by AUFM2 and AUFR2. So a second distribution of coagulation products is created. The particles removed by the coagulation are not immediately substracted from the distribution but stored in a third distribution of removed particles. This gives the same treatment for all particle combinations independent from their appearance in the loop. After the loop all three distributions are summed up to obtain the new distribution of dust particles. It is tested whether negative number densities occur. If this happens, the integration stepwidth was too large, i.e. the chosen accuracy is insufficient, and the program terminates.

EFFT is the function to compute the integration stepwidth depending on the user defined accuracy. Analogously to INTEGR, the interaction between particles is computed, but here, only the interaction between the smallest particles of the distribution and all other particles is computed. Then, the time stepwidth is computed which provides a relative change of the number density of the smallest particles by the chosen accuracy. The stepwidth is rounded to a two digits mantissa+decimal exponent. This stepwidth is then used by INTEGR.

5.5 Interaction routines (NEU, QCOLL, ANEUT, APROB, APROM, APROG)

These routines contain the functions in the integration kernel. NEU computes the Γ-values and QCOLL gives the coagulation rate using the collision cross sections for the different velocity regimes implemented as ANEUT, APROB, APROM, and APROG.

In the subroutine NEU, the lg(V) and the r values of the product of the coagulation of two particles with parameters given by two points of the coordinate grid are calculated. At first the volumes are added up. Then, it is tested whether the larger particle is compact and if so the formula for compact particles (Eq. (10)) is applied. Otherwise, the total number of subgrains in both particles is computed, the r_{PCA}(V_{1}n_{sub,2}/V_{2})-value and the Exponent ex are determined, and finally the new r parameter is calculated from Eq. (3). All operations are separated for the case that particle 1 is larger than particle 2 and vice versa.

The function QCOLL determines the coagulation rate depending on the parameters of two particles according to two points of the coordinate grid. At first the function ANEUT is called to determine the geometric collision cross section since this value is used in all velociity regimes. Then, the different relative velocities are compared and the function for the average collision cross section is chosen that is corresponding to the velocity regime.

In the first case, the turbulent velocities are above 0.4 times the stickiung velocity. Here, the geometric collision cross section A_{coll}(neutral) is used and the effective relative velocity is reduced according to Eq. (23). For average turbulent velocities above 5×v_{stick} the coagulation rate is zero. In the second case the asymmetry velocities are negligible due to f_{ex} > 0.9 (see Eq. 35). If T/T_{eff} > 0.95 the average collision cross section for a thermal velocity distribution as computed by APROB is used. Otherwise, the collision cross section for a general Maxwellian velocity distribution with T_{eff} as computed by APROM is used. In the fourth case the thermal and turbulent velocities are negligible with respect to the asymmetry velocities. Here, the average collision cross section for the superposition of a Maxwellian and a Gaussian velocity distribution (APROG) is used and simply multiplied by the relative velocity. In the fifth and most general case non of the velocities is negligible. Here, we have to combine the collision cross section for the superposition (APROG) for that part of the dust particles which obtains suprathermal relative velocities with the collision cross section for the Maxwellian velocity distribution (either APRAB or APROM) for the fraction f_{ex} of particles which are excluded from the asymmetry velocities.

The function ANEUT computes the geometric collision cross section of two particles using Eq. (36). The functions APROB, APROM, APROG compute the average collision cross section for thermal, Maxwellian and a superposition of Maxwellian and Gaussian velocities. They use the value of ANEUT. All three routines have the same structure. At first the additional radius for the collision between neutral and charged particles is computed depending on the ratio of the particle radii (Eq. (50)). This value is added to the toothing radius of the smaller particle and Eq. (36) is used to compute the effective collision cross section for particles with this corrected radius. This is done twice since both particles may be neutral or charged. Then, the contributions from the terms with both neutral particles, one neutral particle and a charged one, and particles with different signs of the charges to the average collision cross section are added up (Eq. (59–61)). Finally, the contributions of collisions where both particles are negatively charged are added. A series of contributions with growing Z_{1} × Z_{2}-values is iteratively created and summed up to the total collision cross section for the whole charge distribution on both collision partners.

5.6 Transformation between grid coordinates and physical quantities (WM,WR, INEXM,INDEXR, AUFM2, AUFR2)

These routines compute the real parameters lg(V) and r from the grid coordinates and vice versa the grid coordinates from the parameters. In the second case a parameter characterizing the distance of a point from the next grid point is computed which will be used to determine the weighting factors for the distribution of one point over the next four grid points.

The functions WM and WR compute the logarithmic volume and the density parameter for a given grid index. The subroutines INDEXM and INDEXR compute the next lower grid index for a given parameter point in lg(V) or r, respectively. Additionally, they give the distance of the point from that grid point relative to the grid spacing. AUFM2 and AUFR2 are subroutines that compute the factors for the distribution of a point over the next two neighbouring grid points in lg(V) and r depending on the distance value computed by INDEXM and INDEXR, respectively. These factors are chosen in such a way that the total volume and the total density parameter remain conserved in this splitting process.

5.7 Grid treatment (EFFGRID, NEWGRID, SHIFTGRID, STRECHGRID)

These routines adjust the coordinate grid to the dust distribution. They compute the coordinate range covered by the distribution and cut the distribution at negligible number densities.

EFFGRID computes the maximum number density of the distribution and defines the distribution cut-off at 10^{−16} of this value. Only particles with higher number densities are treated by INTEGR. The routine tests in which parameter range the non-negligible distribution lies. By this way the effective coordinate grid is restricted so that all operations may run on a subgrid of the total 200×40 parameter grid.

NEWGRID computes the parameters of a coagulation product in the collision of two particles with the highest volume and the lowest density parameter of the current subgrid. If the coordinates of this particle would fall out of the coordinate grid a change of the grid will be performed. At first it is tested whether the new coordinates may be covered by a simple shift of the grid. This is done by SHIFTGRID. If this is not possible, the grid is streched by a factor 2 by STRECHGRID. Both operations are done separately for the logarithimic volume scale and the density parameter scale. SHIFTGRID and STRECHGRID consist of two parts which are alternatively called - one for the lg(V)- and one for the r-scale. At first, the new parameter range is computed and then the data points are assigned to their new coordinate points.

5.8 Used data structures and the data exchange between different parts

All computation operations necessary for the kernel integration which are limited to one collision partner are tried to be carried out before the integration so that auxiliary quantities are created by INITFELDS which are used in the integration kernel routine QCOLL. They are stored in COMMON blocks. These are

HELPNEU V – V in V_{0} VVDIS – n_{sub} = V/⟨V_{sub}⟩

RR – r = V^{2/3}/A

HELPVEL VBROWN – v_{therm}^{2} in (cm/s)^{2} VTURB – v_{turb} in cm/s

WURZ – \sqrt{V ρ} in g^{1/2}
RSUB2 – 0.147 × a_{sub}^{2/3} in cm^{2/3}

HELPAREA AP – A in cm^{2} RTOOTH – \sqrt{π}a_{tooth} in cm RV2 – πa_{tooth}^{2} in cm^{2} RVD – (\sqrt{π}a_{tooth})^{1.22} in cm^{1.22}

HELPCHRG ERV – τ F0 – f_{0} FMINUS – ∑_{Z<0} f_{Z} FPLUS – f_{+1} CHRG – ⟨−Z⟩ for ⟨−Z⟩ > 1, 1 otherwise FCH – ∑_{Z<0} |Z|f_{Z}

DIPOLE RSUB – \sqrt{π}a_{sub} in cm EHSUB – π^{2}a_{sub}^{4}/(2τ_{sub}) in cm^{4} RVADD – \sqrt{π}a_{x}(thermal) in cm RVMIN – \sqrt{π}a_{min}(thermal) in cm RVMINA – \sqrt{π}a_{min}(Maxwell+Gauss) in cm

COMPNEU AN – A in A_{0}

AR – 2/3A\sqrt{2\sqrt{π}a_{tooth}} in cm^{5/2}

ASYMM VBS – (v_{0}⟨γ⟩)^{2} in (cm/s)^{2}

FEX – f_{ex}

ASYMCH HADD – \sqrt{π/(2τ_{sub})}

HMIN – \sqrt{1/(2τ_{sub})}

CHRGES DIS – table of charge distributions with f_{+1}, f_{0}, f_{−1}, f_{−2}, f_{−3}, f_{−4} for every index i =

⟨−Z⟩ × 10

Additional COMMON blocks are used for the main data characterizing the dust distribution. They are accessed by most routines directly. We have FELD containing the two dimensional field for the dust distribution f(V, r) SIZE with the first and the last coordinate value in lg(V), the first and the last coordinate value

in r, and three auxiliary quantities used by the transformation routines: 199/(lg(V_{max})− (lg(V_{min})), 39/(r_{max} − r_{min}), and (V_{max}/V_{min})^{1/199}.

GRENZE with the cut-off value for the distribution and the four grid coordinates of the current

subgrid as computed by EFFGRID

All other data are exchanged as parameters in subroutine calls.

6 References

Chokshi A., Tielens A.G.G.M., Hollenbach D. 1993, ApJ 407, 806 Draine B.T., Sutin B. 1988, ApJ 320, 803 Nishi R., Nakano T., Umebayashi T., ApJ 368, 181 Ossenkopf V. 1993, A&A 280, 617 Ossenkopf V. 1994, Ph.D thesis, Friedrich-Schiller-Universität Jena