Constant rate coagulation

Background
The growth of planetesimals can be modelled as a function of collisions between smaller particles that lead to the distribution of mass over a range of object sizes. This can be modelled for a particle mass \( k \) and a constant collision rate using the coagulation equation:
\[ \frac{dn_k}{dt} = 0.5 \sum_{i+j = k} An_i n_j + \sum_{i=1}^\infty An_i n_k \]
Method
In order to model this growth, we can integrate this equation using a left Riemann sum in Python for sufficiently small time steps \( dt \).  To do this we loop over an array where each index represents a mass.  Taking the \( k^{th} \) element we add and subtract from it in accordance with the coagulation equation.

def coagulation(n, nalt, A, dt, nsteps):
    narr = zeros((nsteps, len(n)))
    for t in range(nsteps):
        for k in range(1,len(n)):
            added = 0.
            subtracted = 0.
            for i in range(1,len(n)):
                added = added + 0.5*A*n[i]*n[k-i]*dt
                subtracted = subtracted + n[k]*A*n[i]*dt
            dn_k = (added - subtracted)
            nalt[k] = nalt[k] + dn_k
        narr[t, ] = array(nalt)
        n = array(nalt)
    print n
    return narr

Results
The result of constant-rate coagulation is for the number of particles with smallest mass ( \( k=1 \)) to decrease with time while particles of larger mass increase in numbers until they all converge to a constant value.

Fig. 1—: The time evolution of several illustrative mass values with an initial mass \( n_1 =1 \).

This can also be seen by taking the ratio \( n_k/n_1 \) and graphing it against time:

Fig. 2—: \( n_k/n_1 \) as a function of time.  Note that each line eventually tends to 1.

We can also use this equation to see how the mass distribution evolves over time.  We predict that given a sufficient amount of time, the entire mass range will have a constant number of particles as in Figure 1.  Graphing our results a different way, we find that this is most likely true.

Fig. 3—: Number of particles at a particular mass for several illustrative times.
(time ranges from 0 to 100 in steps of size 0.1)

Fig. 4—: Total mass enclosed in particles size \( k \) or larger.

Discussion
Considering the coagulation equation, we note that there must be an infinite number of mass values in order to ideally show the evolution of particle sizes.  With a finite length of the array n, if there are fewer than \( 2^{\# steps} \) masses, the system eventually begins to lose mass as it is subtracted from the largest bin, but not added to any larger mass particles.  We can track this loss of mass using an initial  number of 1 in the \( k= 1 \) bin.  This problem can be potentially minimised for a greater number of time steps by using logarithmic mass bins.

Fig. 5—: Change in mass relative  to initial mass.  By the end  of this simulation we lose ~10%