text:
# Untitled Document

A SURVEY OF NUMERICAL SOLUTIONS TO THE

COAGULATION EQUATION

Man Hoi Lee

Department of Physics, University of California, Santa Barbara, CA 93106, USA

E-mail: mhlee@europa.physics.ucsb.edu

ABSTRACT

We present the results of a systematic survey of numerical solutions to the coagulation equation for a rate coefficient of the form $`A_{ij}(i^\mu j^\nu +i^\nu j^\mu )`$ and monodisperse initial conditions. The results confirm that there are three classes of rate coefficients with qualitatively different solutions. For $`\nu 1`$ and $`\lambda =\mu +\nu 1`$, the numerical solution evolves in an orderly fashion and tends toward a self-similar solution at large time $`t`$. The properties of the numerical solution in the scaling limit agree with the analytic predictions of van Dongen and Ernst. In particular, for the subset with $`\mu >0`$ and $`\lambda <1`$, we disagree with Krivitsky and find that the scaling function approaches the analytically predicted power-law behavior at small mass, but in a damped oscillatory fashion that was not known previously. For $`\nu 1`$ and $`\lambda >1`$, the numerical solution tends toward a self-similar solution as $`t`$ approaches a finite time $`t_0`$. The mass spectrum $`n_k`$ develops at $`t_0`$ a power-law tail $`n_kk^\tau `$ at large mass that violates mass conservation, and runaway growth/gelation is expected to start at $`t_{\mathrm{crit}}=t_0`$ in the limit the initial number of particles $`n_0\mathrm{}`$. The exponent $`\tau `$ is in general less than the analytic prediction $`(\lambda +3)/2`$, and $`t_0=K/[(\lambda 1)n_0A_{11}]`$ with $`K=1`$$`2`$ if $`\lambda \mathrm{}>1.1`$. For $`\nu >1`$, the behaviors of the numerical solution are similar to those found in a previous paper by us. They strongly suggest that there are no self-consistent solutions at any time and that runaway growth is instantaneous in the limit $`n_0\mathrm{}`$. They also indicate that the time $`t_{\mathrm{crit}}`$ for the onset of runaway growth decreases slowly toward zero with increasing $`n_0`$.

1. INTRODUCTION

Smoluchowski’s coagulation equation is the mean-rate equation that describes the evolution of the mass spectrum of a collection of particles due to successive mergers. It is widely used for modeling growth in many fields of science. Examples include planetesimal accumulation, mergers in dense clusters of stars, coalescence of interstellar dust grains, and galaxy mergers in astrophysics, aerosol coalescence in atmospheric physics, colloids, and polymerization and gelation (see, e.g., Drake 1972, Ernst 1986, Jullien and Botet 1987, Lee 1993, 2000, and references therein).

If the masses of the particles are integral multiples of a minimum mass $`m_1`$, the coagulation equation is written in discrete form as

$$\frac{dn_k}{dt}=\frac{1}{2}\underset{i+j=k}{}A_{ij}n_in_jn_k\underset{i=1}{\overset{\mathrm{}}{}}A_{ki}n_i$$
$`(1)`$
where $`n_k`$ is the number of particles of mass $`m_k=km_1`$ in a volume $`V`$ and $`A_{ij}`$ is the rate coefficient (or coagulation kernel) for mergers between particles of mass $`m_i`$ and $`m_j`$.1 We can interpret $`n_k`$ as the concentration (i.e., the number of particles per unit volume) if we replace $`A_{ij}`$ by $`A_{ij}^{}=VA_{ij}`$. In equation (1), it is assumed that the merging of two particles of mass $`m_i`$ and $`m_j`$ results in one particle of mass $`m_i+m_j`$. The coagulation equation can also be written in continuous form as

$$\begin{array}{cc}\hfill \frac{dn(m)}{dt}& =\frac{1}{2}_0^m𝑑m^{}A_{m^{},mm^{}}n(m^{})n(mm^{})\hfill \\ & n(m)_0^{\mathrm{}}𝑑m^{}A_{m,m^{}}n(m^{})\hfill \end{array}$$
$`(2)`$
where $`n(m)dm`$ is the number of particles of mass between $`m`$ and $`m+dm`$ and $`A_{m,m^{}}`$ is the rate coefficient for mergers between particles of mass $`m`$ and $`m^{}`$.

Examples of the rate coefficient $`A_{ij}`$ as a function of $`m_i`$ and $`m_j`$ (or equivalently $`i`$ and $`j`$) that arise in various problems can be found in the references cited above. Most rate coefficients used in the literature are homogeneous functions of degree $`\lambda `$, i.e., $`A_{ai,aj}=a^\lambda A_{ij}`$. The exponent $`\lambda `$ specifies the mass dependence of the probability of merger for two particles of comparable mass ($`ij`$). It is also useful to classify $`A_{ij}`$ according to the exponents $`\mu `$ and $`\nu `$ for the merger between a light particle and a heavy particle: $`A_{ij}i^\mu j^\nu `$ for $`ij`$ and $`\mu +\nu =\lambda `$ (see, e.g., Ernst 1986). For example, $`A_{ij}i+j`$ has $`\mu =0`$, $`\nu =1`$, and $`\lambda =1`$.

For a few simple rate coefficients and monodisperse initial conditions (i.e., $`n_0`$ particles of mass $`m_1`$ at $`t=0`$), there are exact analytic solutions to the discrete form of the coagulation equation (Trubnikov 1971, Ernst 1986, and references therein). The analytic solutions for $`A_{ij}`$ constant and $`i+j`$ show orderly evolution of a smooth mass spectrum at all times, and they agree with the results from Monte Carlo simulations of the merger process in the limit $`n_0\mathrm{}`$ (with $`n_0/V`$ fixed). These two cases are examples of orderly growth. The analytic solution for $`A_{ij}ij`$ develops a power-law tail with $`n_kk^{5/2}`$ at large $`k`$ as $`tt_0=1/(n_0A_{11})`$. This power-law tail violates mass conservation because it implies a nonzero mass flux to the infinite-mass bin. In this case, the results from Monte Carlo simulations in the limit $`n_0\mathrm{}`$ agree with the solution to the coagulation equation at $`tt_0`$, but they show a transition from a smooth mass spectrum to a smooth spectrum plus a massive runaway particle at $`t=t_0`$ (Spouge 1985, Wetherill 1990). The runaway particle (gel) acquires a mass much larger than that of the other particles (sol) in the system and becomes detached from the smooth mass spectrum of the rest of the particles at $`t>t_0`$. This phenomenon is known as runaway growth in the astrophysics literature, and the transition is considered to be the gelation transition in studies of polymerization and gelation.

For most rate coefficients, there are no exact analytic solutions to the coagulation equation. However, there are extensive analytic results on the asymptotic properties of the solutions for $`A_{ij}`$ with $`\nu 1`$ (see, e.g., review by Ernst 1986). It is important to note that some of these analytic results (such as the shape of the mass spectrum at small and large mass) are derived based on assumptions (such as self-similar evolution) that have not been verified. Nevertheless, the analytic results indicate that there are qualitatively two types of solutions to the coagulation equation if $`\nu 1`$:

(1) if $`\nu 1`$ and $`\lambda 1`$, the solution shows orderly growth at all times;

(2) if $`\nu 1`$ and $`\lambda >1`$, the solution develops in a finite time $`t_0`$ a power-law tail at large mass that violates mass conservation.

For $`A_{ij}`$ with essentially $`\nu <1`$, it has been proved that a solution to the coagulation equation exists for all times (including $`t>t_0`$ if $`\lambda >1`$) and that the coagulation equation is the limit of finite system (whether or not the runaway particle and the other particles are allowed to interact at $`t>t_0`$ if $`\lambda >1`$) (Leyvraz and Tschudi 1981, Spouge 1985, Bak and Heilmann 1994, Jeon 1998). (For $`\nu =1`$ and $`\lambda >1`$, we have the example $`A_{ij}ij`$, where the coagulation equation needs to be modified for $`t>t_0`$ if there is sol-gel interaction; Ziff et al 1983, Bak and Heilmann 1994.)

Several authors have investigated the properties of the solutions to the coagulation equation for $`A_{ij}`$ with $`\nu >1`$, using series expansion of the mass spectrum $`n_k(t)`$ about $`t=0`$ and moments of the mass spectrum (McLeod 1962, Hendriks et al 1983, Ernst et al 1984, van Dongen 1987a). The results suggest that

(3) if $`\nu >1`$, there are no self-consistent solutions that conserve mass at any time.

An alternative to the analytic approach is to solve the coagulation equation numerically. In Paper I (Lee 2000), a numerical code that can yield accurate solutions to the discrete form of the coagulation equation, equation (1), with a reasonable number of numerical mass bins was developed. The numerical code was used to study solutions to the coagulation equation for $`A_{ij}`$ that are limiting cases for gravitational interaction. We considered geometric or gravitational focusing dominated cross-section, mass-independent or equipartition velocity dispersion, and the power-law index of the mass-radius relation $`\beta =1/3`$ (for planetesimals) or $`2/3`$ (for stars). For the two cases with geometric cross-section and $`\beta =1/3`$, which have $`\nu 1`$ and $`\lambda 1`$, the mass spectrum evolves in an orderly fashion and tends towards a self-similar solution at large time. For the remaining cases, which have $`\nu >1`$, the numerical mass spectrum shows, after some evolution, an exponential drop in an intermediate mass range and a power-law tail of the form $`n_kk^\nu `$ (or $`n(m)m^\nu `$) at large mass. This mass spectrum is not self-consistent because the power-law tail implies a mass flux2 As we pointed out in Paper I, in these cases, a massive particle grows mainly by accumulating low-mass particles because of the much larger number of low-mass particles. So the growth rate of a massive particle is $`\dot{m}=𝑑m^{}A_{m,m^{}}n(m^{})m^{}m^\nu `$, since the integral is dominated by the range $`m^{}m`$. Hence $`n\dot{m}`$ at the high-mass end of the mass spectrum is non-zero and independent of $`m`$ if $`nm^\nu `$. However, we did not point out that the mass flux from particles of mass $`m^{}m`$ to particles of mass $`m^{}>m`$ is $`F_mmn(m)\dot{m}(m)+_m^{m_{\mathrm{max}}}𝑑m^{}n(m^{})\dot{m}(m^{})`$. With $`n\dot{m}`$ independent of $`m`$, $`F_mn\dot{m}m_{\mathrm{max}}`$, which is independent of $`m`$ but increases with $`m_{\mathrm{max}}`$. We have verified this by an explicit evaluation of the mass flux (equation ) for the numerical solutions.  and, if $`1<\nu 2`$, a cumulative mass that diverge with the maximum particle mass, $`m_{\mathrm{max}}`$, included in the computational grid. The time at which the power-law tail develops decreases toward zero as the numerical parameter $`n_{\mathrm{min}}`$ decreases (see Section 3 for the definition of $`n_{\mathrm{min}}`$). Thus the numerical results strongly suggest that there are no self-consistent solutions to the coagulation equation at any time if $`\nu >1`$. We also considered a case with $`\beta =0`$ as an example with $`\nu 1`$ and $`\lambda >1`$, and its mass spectrum develops a power-law tail that violates mass conservation in a finite time $`t_0`$. We discussed a simplified merger problem that illustrates the qualitative differences in the solutions to the coagulation equation for the three classes of $`A_{ij}`$. The results in Paper I (and the analytic results cited above) strongly suggest that there are two types of runaway growth. For $`A_{ij}`$ with $`\nu 1`$ and $`\lambda >1`$, runaway growth starts at a finite time $`t_{\mathrm{crit}}=t_0`$, the time at which the coagulation equation solution begins to violate mass conservation, in the limit $`n_0\mathrm{}`$. For $`A_{ij}`$ with $`\nu >1`$, runaway growth is instantaneous in the limit $`n_0\mathrm{}`$, and there are indications (since decreasing $`n_{\mathrm{min}}`$ is similar to increasing $`n_0`$) that the time $`t_{\mathrm{crit}}`$, in units of $`1/(n_0A_{11})`$, for the onset of runaway growth decreases slowly toward zero with increasing initial number of particles $`n_0`$. Recent Monte Carlo simulations have shown that the time for all particles to coalesce into a single particle decreases as a power of the logarithm of $`n_0`$ if $`\nu >1`$ (Malyshkin and Goodman 2001; see also Spouge 1985, Jeon 1999).

The study in Paper I was focused on the rate coefficient for gravitational interaction, and the range of $`\mu `$ and $`\nu `$ studied was limited. In particular, $`\mu `$ was limited to $`0`$ and $`\pm 1/2`$, and the region $`\mu >0`$ and $`\nu <1`$ was not studied at all. Other authors have obtained numerical solutions to the coagulation equation for rate coefficients that arise in specific problems. However, we are not aware of any study that has systematically surveyed the properties of numerical solutions as a function of $`\mu `$ and $`\nu `$ (and $`\lambda `$) and compared them to the analytic results on the asymptotic properties. For such a study, it is important that the numerical code used can follow the evolution of the mass spectrum accurately for a long time. After the computation for this paper was nearly complete, it came to our attention that Krivitsky (1995) has obtained numerical solutions to the continuous form of the coagulation equation, equation (2), for $`A_{m,m^{}}(mm^{})^{\lambda /2}`$ (which has $`\mu =\nu =\lambda /2`$) and $`A_{m,m^{}}(m+m^{})^\nu `$ (which has $`\mu =0`$). For $`A_{m,m^{}}(mm^{})^{\lambda /2}`$ with $`\lambda 1`$, Krivitsky found that the numerical solutions are self-similar at large time, but that unlike the analytic result, the asymptotic behavior of the scaling function at small mass is not a power law. As we shall see, the latter result is incorrect because Krivitsky did not evolve the numerical solutions for a sufficiently long time to see the true asymptotic behavior at small mass. The scaling function does in fact approach the analytically predicted power law at small mass, but in a damped oscillatory fashion that was not known previously. It is unlikely that Krivitsky could have solved the coagulation equation accurately for the necessary amount of time because the numerical code used by Krivitsky does not conserve mass. For $`A_{m,m^{}}(m+m^{})^\nu `$ with $`\nu >1`$, Krivitsky found that the mass spectrum develops a slowly decreasing tail at very small time. Only the numerical solution for $`\nu =2`$ was shown. Its evolution is qualitatively similar to that found in Paper I for $`\nu >1`$, but it is not clear that the tail is power-law in nature because the maximum particle mass included in Krivitsky’s computational grid is not large enough. It was also not demonstrated that the numerical solution is not sensitive to the other numerical parameters.

In this paper we present the results of a systematic survey of numerical solutions to the coagulation equation. The purpose of this survey is (1) to confirm that there are three classes of rate coefficients with qualitatively different solutions to the coagulation equation and that the boundaries of these three classes are as stated above; (2) to investigate, in the cases where self-consistent solutions exist, whether the solutions approach self-similar solutions as $`t\mathrm{}`$ or $`tt_0`$ and whether the scaling behaviors agree with the analytic results; (3) to study the dependence of $`t_0`$ on the exponents $`\mu `$, $`\nu `$, and $`\lambda `$ for the runaway growth cases with $`\nu 1`$ and $`\lambda >1`$; and (4) to investigate whether the behaviors of the numerical solutions found in Paper I for the cases with $`\nu >1`$ are valid in general. In Section 2 we describe the rate coefficient and the initial conditions used in this survey. In Section 3 we provide a brief summary of the numerical methods developed in Paper I for solving the coagulation equation and additional information on the accuracy of the numerical results. The results are presented in Section 4, and the conclusions are summarized in Section 5.

2. RATE COEFFICIENT AND INITIAL CONDITIONS

In this paper we consider a rate coefficient of the form

$$A_{ij}=\frac{1}{2}\left(i^\mu j^\nu +i^\nu j^\mu \right)$$
$`(3)`$
with $`\mu \nu `$. Note that $`A_{ai,aj}=a^{\mu +\nu }A_{ij}`$ and $`A_{ij}i^\mu j^\nu `$ for $`ij`$, consistent with the definitions of the exponents $`\mu `$, $`\nu `$, and $`\lambda =\mu +\nu `$ in Section 1. Since the rate coefficient in equation (3) contains the exponents $`\mu `$ and $`\nu `$ as parameters explicitly, we can survey the entire $`(\mu ,\nu )`$ space by varying $`\mu `$ and $`\nu `$. This rate coefficient includes $`A_{ij}=(ij)^{\lambda /2}`$ (for $`\mu =\nu =\lambda /2`$) and $`A_{ij}=(i^\nu +j^\nu )/2`$ (for $`\mu =0`$), which have been used to model polymerization (e.g., Hendriks et al 1983), and the cases $`(\mu ,\nu )=(1,4/3)`$ and $`(2,2)`$, which have been used to model planetesimal accumulation and stellar merger (Malyshkin and Goodman 2001). It also has the nice property that it includes the three cases with exact analytic solutions to the coagulation equation: $`A_{ij}=1`$, $`(i+j)/2`$, and $`ij`$ for $`(\mu ,\nu )=(0,0)`$, $`(0,1)`$, and $`(1,1)`$, respectively.

We have obtained numerical solutions to the discrete form of the coagulation equation for the cases shown in figure 1. The ranges of $`\mu `$ and $`\nu `$ considered contain most of the values encountered in practical applications (but usually for other forms of $`A_{ij}`$). The numerical solutions were computed for the monodisperse initial conditions with $`n_0`$ particles of mass $`m_1`$, i.e., $`n_k(t=0)=n_0\delta _{k1}`$, where $`\delta _{k1}`$ is the Kronecker delta. Hereafter, we adopt units such that $`n_0=1`$, $`m_1=1`$, and $`A_{11}=1`$. With this set of units, $`m_k=k`$ and time is in units of $`1/(n_0A_{11})`$, the timescale for every particle of mass $`m_1`$ to merge with another particle of mass $`m_1`$.

3. NUMERICAL METHODS

In Paper I we have developed a numerical code that can yield accurate solutions to the discrete form of the coagulation equation (equation ) with a reasonable number of numerical mass bins. A detailed description of the code can be found in Paper I. In this section we provide a brief summary of the algorithm and its numerical parameters. We also provide additional information on the accuracy of the numerical results.

Our numerical code uses a combination of linearly and logarithmically spaced numerical mass bins. The first $`N_{bd}`$ numerical mass bins are linearly spaced with $`\stackrel{~}{m}_k=k`$. They have boundaries $`\stackrel{~}{m}_{k\pm 1/2}=k\pm 1/2`$ and width $`\mathrm{\Delta }\stackrel{~}{m}_k\stackrel{~}{m}_{k+1/2}\stackrel{~}{m}_{k1/2}=1`$. The next $`N_{bd}\times N_{\mathrm{dec}}`$ numerical mass bins are logarithmically spaced, with $`N_{bd}`$ bins per decade of mass; thus the $`k`$th mass is $`\stackrel{~}{m}_k=(\stackrel{~}{m}_{k+1/2}+\stackrel{~}{m}_{k1/2})/2`$, with $`\stackrel{~}{m}_{k+1/2}/\stackrel{~}{m}_{k1/2}=10^{1/N_{bd}}`$. There are in total $`N_{\mathrm{max}}=(N_{\mathrm{dec}}+1)N_{bd}`$ mass bins, and the mass of the most massive particles in the computational grid, $`m_{\mathrm{max}}`$, is approximately $`N_{bd}10^{N_{\mathrm{dec}}}`$. Initially, the number of “active” bins $`N_{\mathrm{bin}}N_{\mathrm{max}}`$. At the end of each time step, $`N_{\mathrm{bin}}`$ is increased (if necessary) to include all bins with $`N_k>n_{\mathrm{min}}`$, where $`N_k`$ is the total number of particles in bin $`k`$ and $`n_{\mathrm{min}}`$ is a numerical parameter; it is also increased if $`N_{N_{\mathrm{bin}}+1}`$ becomes comparable to the power-law extrapolation from $`N_{N_{\mathrm{bin}}}`$. Before $`N_{\mathrm{bin}}`$ reaches $`N_{\mathrm{max}}`$, $`N_{\mathrm{max}}`$ (or equivalently $`m_{\mathrm{max}}`$) has no effects on the numerical results. The numerical parameter $`n_{\mathrm{min}}`$ is specified in units of $`n_0`$ and, e.g., $`n_{\mathrm{min}}=10^{30}`$ in units of $`n_0`$ is equivalent to $`n_{\mathrm{min}}=1`$ and $`n_0=10^{30}`$ in physical units. Thus the effect of the numerical parameter $`n_{\mathrm{min}}`$ is similar to not allowing fractionally occupied numerical mass bins to interact.

The fundamental quantity evolved by our numerical code is the total mass $`M_k`$ in bin $`k`$. During a time step, the code calculates for each combination of $`i`$ and $`j`$ (with $`ijN_{\mathrm{bin}}`$) the mass loss from bins $`i`$ and $`j`$ due to mergers between particles in those bins and distributes the total mass of the merger products among the mass bins. Thus the code conserves mass exactly. For $`ijN_{bd}`$, it is correct to assume that the merging particles have masses $`\stackrel{~}{m}_i`$ and $`\stackrel{~}{m}_j`$ and that the merger products have mass $`\stackrel{~}{m}_i+\stackrel{~}{m}_j`$. For $`ij`$ and $`j>N_{bd}`$, we assume that the particles in bin $`i`$ have mass $`\stackrel{~}{m}_i`$ (which is exact for $`iN_{bd}`$) and that the mass distribution within bin $`j`$ follows a power-law distribution:

$$\rho _j(m)=c_j(m/\stackrel{~}{m}_{j1/2})^{q_j}\mathrm{for}\stackrel{~}{m}_{j1/2}<m\stackrel{~}{m}_{j+1/2}$$
$`(4)`$
where $`\rho _j(m)dm`$ is the total mass of particles with mass between $`m`$ and $`m+dm`$. The merger products have masses between $`\stackrel{~}{m}_i+\stackrel{~}{m}_{j1/2}`$ and $`\stackrel{~}{m}_i+\stackrel{~}{m}_{j+1/2}`$, and they are either added to a single bin $`k`$ (if $`\stackrel{~}{m}_{k1/2}\stackrel{~}{m}_i+\stackrel{~}{m}_{j1/2}`$ and $`\stackrel{~}{m}_i+\stackrel{~}{m}_{j+1/2}\stackrel{~}{m}_{k+1/2}`$) or distributed between bins $`k`$ and $`k+1`$. In equation (4), the power-law index $`q_j`$ is obtained from the masses in the adjacent bins,

$$q_j=\frac{\mathrm{log}\left(\frac{M_{j+1}}{\mathrm{\Delta }\stackrel{~}{m}_{j+1}}/\frac{M_{j1}}{\mathrm{\Delta }\stackrel{~}{m}_{j1}}\right)}{\mathrm{log}\left(\stackrel{~}{m}_{j+1}/\stackrel{~}{m}_{j1}\right)}$$
$`(5)`$
and the normalization constant $`c_j`$ from the constraint

$$_{\stackrel{~}{m}_{j1/2}}^{\stackrel{~}{m}_{j+1/2}}𝑑m\rho _j(m)=M_j.$$
$`(6)`$
Our numerical code uses the second-order Runge-Kutta method with a variable time step. The time step is continuously adjusted so that the fractional change of each $`M_k`$ per time step is less than $`\delta _M`$ and the mass loss from bin $`k`$ does not exceed $`M_k`$.

In Paper I we have compared in detail the numerical solutions from our code to the exact analytic solutions for $`A_{ij}=1`$, $`(i+j)/2`$, and $`ij`$, with the last case at $`t<t_0`$ only. The accuracy of the numerical solutions is extremely insensitive to $`\delta _M`$ and $`n_{\mathrm{min}}`$ and improves rapidly with increasing $`N_{bd}`$. Hereafter, unless stated otherwise, the numerical results were obtained using $`\delta _M=5\%`$, $`n_{\mathrm{min}}=10^{30}`$ (in units of $`n_0`$), and $`N_{bd}=40`$.

We report here several additional tests of our code. Ziff (1980) has constructed three forms of rate coefficients, with a parameter $`\gamma `$, for which a single moment $`_\gamma (t)=_{k=1}^{\mathrm{}}m_k^\gamma n_k(t)`$ of the mass spectrum $`n_k(t)`$ can be calculated analytically. We have obtained numerical solutions for a few of these rate coefficients (including both orderly and runaway growth cases) and have confirmed that the numerical results for the moment $`_\gamma (t)`$ agree with the analytic results, with accuracy similar to what was found for the three cases with exact analytic solutions.

In Section 4.3 we shall be interested in extending the calculations for some of the runaway growth cases with $`\nu 1`$ and $`\lambda >1`$ to $`t>t_0`$. Therefore, another test that we have performed is to extend the comparison for the case $`A_{ij}=ij`$ to $`t>t_0`$. For $`A_{ij}=ij`$, the evolution of the mass spectrum at $`t>t_0`$ depends on whether or not the runaway particle (gel) and the other particles (sol) interact. The (unmodified) coagulation equation is valid if there is no sol-gel interaction, and it has an exact analytic solution with $`n_k(t)t^1k^{5/2}`$ at large $`k`$ for all $`t>t_0`$ (Leyvraz and Tschudi 1981, Ziff et al 1983). (The coagulation equation can be modified to take into account sol-gel interaction, and an exact analytic solution also exists for this modified coagulation equation; Ziff et al 1983.) Since our code does not take into account sol-gel interaction and does not allow merger products with masses greater than $`m_{\mathrm{max}}`$ to interact, we expect the numerical solution at $`t>t_0`$ to agree with the analytic solution to the unmodified coagulation equation, except for $`m_km_{\mathrm{max}}`$. We have integrated the case $`A_{ij}=ij`$ up to $`t=1.25`$ and have found that the numerical solution at $`t>t_0=1`$ is in excellent agreement with the analytic solution for $`m_k\mathrm{}<0.01m_{\mathrm{max}}`$.

A quantity that will be discussed extensively in Section 4 is the logarithmic slope $`d\mathrm{ln}n/d\mathrm{ln}m`$ of the mass spectrum. For the numerical results, we use

$$d\mathrm{ln}n/d\mathrm{ln}m(\stackrel{~}{m}_k)=q_k1$$
$`(7)`$
where $`q_k`$ is defined in equation (5). Equation (7) is consistent with the power-law approximation used by the code since under this approximation the number distribution of particles within the numerical mass bin $`k`$ is $`n_k(m)=\rho _k(m)/mm^{q_k1}`$.

To determine the accuracy of the numerical results for $`d\mathrm{ln}n/d\mathrm{ln}m`$, we compare the numerical and analytic results for the three cases with exact analytic solutions. The analytic solutions are of the form $`nm^\tau \mathrm{exp}[b(t)m]`$ or $`d\mathrm{ln}n/d\mathrm{ln}m=\tau b(t)m`$ for $`m1`$, where $`\tau =0`$, $`3/2`$, and $`5/2`$ for $`A_{ij}=1`$, $`(i+j)/2`$, and $`ij`$, respectively. Thus $`nm^\tau `$ or $`d\mathrm{ln}n/d\mathrm{ln}m=\tau `$ for $`1mm_{}`$ when the characteristic mass $`m_{}(t)`$ defined in equation (9) is large. In figure 2(a) we show for each case the numerical and analytic $`d\mathrm{ln}n/d\mathrm{ln}m`$ at $`m\mathrm{}<m_{}(t)`$ at a given time. As we noted in Paper I, there is a small lag in the evolution of the numerical solutions for $`A_{ij}=(i+j)/2`$ and $`ij`$. Therefore, in these cases, the numerical results are compared to the analytic results at a slightly earlier time. There are fluctuations in the numerical results in the first decade of the logarithmically spaced mass bins ($`1\mathrm{}<m/N_{bd}\mathrm{}<10`$) due to the discreteness of the mass bins, but the fluctuations are $`\mathrm{}<0.015`$. The numerical results are much smoother and much more accurate outside this mass range. We can see from figure 2(a) that $`\tau `$ can be determined from the numerical results at $`1mm_{}`$ to better than $`\pm 0.001`$.

Figure 2(b) is similar to figure 2(a), but it shows the mass range $`m\mathrm{}>m_{}(t)`$. The differences between the analytic results, which decrease linearly with mass, and the numerical results are small, but there is a small curvature in the numerical results, and the numerical results become increasingly higher than the analytic results with increasing mass. (This is consistent with the observation in Paper I that the numerical solutions show a slightly slower exponential decay at the high-mass end of the mass spectrum.) As a result, if we fit the numerical results near $`d\mathrm{ln}n/d\mathrm{ln}m=20`$ (or $`30`$) to a straight line $`d\mathrm{ln}n/d\mathrm{ln}m=\theta bm`$, the resulting values for $`\theta `$ are greater than the correct values (which are $`\tau `$ as given above) by $`0.13`$$`0.15`$ (or $`0.29`$$`0.36`$). This is the accuracy to which we can check whether a numerical solution is consistent with $`d\mathrm{ln}n/d\mathrm{ln}m=\theta bm`$ and a given $`\theta `$ at $`mm_{}`$.

In Paper I we have discussed our numerical code in the context of numerical codes in the astrophysics literature. Other recent numerical codes include those by Krivitsky (1995), Hill and Ng (1996), and Tzivion et al (1999). As we mentioned in Section 1, the numerical code used by Krivitsky does not conserve mass and would have difficulty following the evolution of the mass spectrum accurately for a long time. The numerical code described by Hill and Ng conserves mass but uses a relatively simple algorithm for distributing merger products. We had in fact tried a similar algorithm for distributing merger products (Quinlan and Shapiro 1989) before we developed the algorithm based on the power-law approximation and had found that the high-mass end of the mass spectrum converges very slowly with increasing grid resolution ($`N_{bd}`$) if the rate coefficient increases steeply with the mass of the particles (i.e., if $`\nu `$ and/or $`\lambda `$ is large). Tzivion et al have developed a mass-conserving numerical code that evolves separately the total number ($`N_k`$) and mass ($`M_k`$) of particles in a numerical mass bin $`k`$. The numerical solutions obtained using this code appear to converge rapidly with increasing grid resolution for the case $`A_{m,m^{}}m+m^{}`$, but there was no demonstration that this is also true for rate coefficients with steeper mass dependence.

4. RESULTS

In this section we present the numerical solutions to the discrete form of the coagulation equation for the rate coefficient and initial conditions described in Section 2 (see also figure 1). For the cases with $`\nu 1`$ (Sections 4.1–4.3), whenever possible, we compare the properties of the numerical solutions to the predictions from self-similar analysis. Hereafter, unless otherwise stated, the self-similar analysis results cited can be found in van Dongen and Ernst (1985a, 1988).

4.1. Orderly Growth Cases with $`\nu 1`$ and $`\lambda <1`$

An example of the numerical results for the mass spectrum evolution for $`\nu 1`$ and $`\lambda <1`$ is shown in figure 3.3 In figure 3 and all subsequent figures, the numerical mass spectrum plotted is $`n(\stackrel{~}{m}_k)=N_k`$ for $`kN_{bd}`$ and $`n(\stackrel{~}{m}_k)=n_k(\stackrel{~}{m}_k)`$ for $`k>N_{bd}`$. In this and all other cases with $`\nu 1`$ and $`\lambda <1`$, the mass spectrum evolves in an orderly fashion. For these cases, we stopped the numerical integrations when the asymptotic behaviors of the solutions were clear and before $`N_{\mathrm{bin}}`$ reached $`N_{\mathrm{max}}`$.

For orderly growth with $`\nu 1`$ and $`\lambda <1`$, self-similar analysis predicts that self-similar solutions have the form

$$n_k(t)=m_{}(t)^2\phi [m_k/m_{}(t)]$$
$`(8)`$
where $`\phi (x)`$ is a scaling function, and the characteristic mass $`m_{}(t)`$ scales as $`t^{1/(1\lambda )}`$. Different definitions of $`m_{}(t)`$, which correspond to different scales for $`x=m_k/m_{}(t)`$ and $`\phi (x)`$, can be used. We adopt

$$m_{}(t)=_3(t)/_2(t)$$
$`(9)`$
where $`_{\mathrm{}}(t)_{k=1}^{\mathrm{}}m_k^{\mathrm{}}n_k(t)`$ is the $`\mathrm{}`$th moment of the mass spectrum. This choice of $`m_{}`$ is convenient because it can also be used for runaway growth with $`\nu 1`$ and $`\lambda >1`$ (Section 4.3).

In all cases, the numerical solution tends toward a self-similar solution of the form (8) at large $`t`$. This is illustrated in figure 4, where we plot the numerical solution at three different times in the form of $`\mathrm{log}m^2n`$ as a function of $`\mathrm{log}m/m_{}`$ for the case shown in figure 3. In the scaling limit (8), $`m^2n=x^2\phi (x)`$. We have evaluated the exponent $`z(t_i)=\mathrm{ln}[m_{}(t_{i+1})/m_{}(t_{i1})]/\mathrm{ln}(t_{i+1}/t_{i1})`$ from the numerical results for $`m_{}(t)`$ at output times $`t_{i1}`$, $`t_i`$, and $`t_{i+1}`$ for all cases with $`\nu 1`$ and $`\lambda <1`$. In all but two cases, the exponent $`z`$ at large $`t`$ agrees with $`1/(1\lambda )`$ to better than one part in $`1.5\times 10^3`$. For the two cases with $`(\mu ,\nu )=(1/6,2/3)`$ and $`(1/3,1/2)`$, the agreement between the numerical results at the end of the numerical runs, $`z=5.985`$ and $`5.963`$, and the analytic result, $`1/(1\lambda )=6`$, is slightly worse, but the numerically determined exponents are still slowly increasing with time at the end of the numerical runs, indicating that the numerical results have not completely reached the asymptotic regime.

Figure 5 shows the numerical results for the scaling function $`\phi (x)`$ for all cases with $`\nu 1`$ and $`\lambda <1`$. The location of the peak of $`x^2\phi `$ is not very sensitive to $`\mu `$ or $`\nu `$, and it is at $`x=0.33`$$`0.71`$. Self-similar analysis predicts that the scaling function $`\phi (x)`$ decays exponentially at large $`x`$. For $`\nu <1`$, $`\phi (x)x^\theta \mathrm{exp}(bx)`$ or $`d\mathrm{ln}\phi /d\mathrm{ln}x=\theta bx`$, with $`\theta =\lambda `$, at large $`x`$. The detailed behavior of $`\phi (x)`$ at large $`x`$ for $`\nu =1`$ depends on the specific form of $`A_{ij}`$. For $`A_{ij}`$ in equation (3) with $`\nu =1`$, the large-$`x`$ behavior of $`\phi (x)`$ is similar to that for $`\nu <1`$, but $`\theta =(\mu +3)/2`$ if $`1\mu <0`$ (van Dongen 1987b; see also Ernst et al 1984). In all cases, the numerical $`\phi (x)`$ decays exponentially at large $`x`$, and $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ at large $`x`$ is consistent with the analytic result (see Section 3 for a discussion of the accuracy of the numerical results at large $`x`$).

The behavior of the scaling function $`\phi (x)`$ at small $`x`$ is qualitatively different for $`\mu <0`$, $`\mu =0`$, and $`\mu >0`$. For the cases with $`\mu <0`$ (figure 5(a)), $`\phi (x)`$ also decays exponentially at small $`x`$, because light particles are rapidly accreted by heavy particles. Self-similar analysis predicts that $`\phi (x)x^a\mathrm{exp}(bx^\mu /\mu )`$ or $`d\mathrm{ln}\phi /d\mathrm{ln}x=a+bx^\mu `$ at small $`x`$, where $`a`$ and $`b`$ are constants that depend on the specific form of $`A_{ij}`$. For $`A_{ij}`$ in equation (3), $`a=2`$ if $`\nu >0`$ and $`a=1`$ if $`\nu =0`$. In figure 6, we plot the numerical results for $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ as a function of $`x^\mu `$ for the cases with $`\mu =1/6`$ to show that the numerical results indeed have the form $`d\mathrm{ln}\phi /d\mathrm{ln}x=a+bx^\mu `$ at large $`x^\mu `$ (or small $`x`$). Furthermore, in all cases with $`\mu <0`$, the value of $`a`$ from least-squares fit is in agreement with the analytic result to better than $`\pm 0.004`$.

For the cases with $`\mu =0`$, i.e., $`A_{ij}=(i^\nu +j^\nu )/2`$, the scaling function shows a power-law behavior $`\phi (x)x^\tau `$ at small $`x`$ (figure 5(b)). The numerical results for the exponent $`\tau `$ are $`0.000`$, $`1.001`$, $`1.033`$, $`1.109`$, $`1.216`$, $`1.347`$, and $`1.500`$ for $`\nu =0`$, $`1/6`$, $`1/3`$, $`1/2`$, $`2/3`$, $`5/6`$, and $`1`$, respectively. Self-similar analysis gives $`\tau =2p_\lambda /w`$, where $`p_\lambda `$ and $`w`$ depend on the specific form of $`A_{ij}`$. van Dongen and Ernst (1985b) have used this expression to derive analytic lower and upper bounds on $`\tau `$ for $`A_{ij}=(i^\nu +j^\nu )/2`$. Our numerical results are consistent with these bounds. Note that the exponent $`\tau `$ is discontinuous at $`\nu =0`$: $`\tau `$ decreases from $`1.5`$ at $`\nu =1`$ to $`1`$ as $`\nu 0`$, but $`\tau =0`$ at $`\nu =0`$ (recall that the $`\nu =0`$ and $`\nu =1`$ cases have exact analytic solutions). In contrast, for $`A_{ij}(i+j)^\nu `$ (which also has $`\mu =0`$), the exponent $`\tau `$ decreases smoothly from $`1.5`$ at $`\nu =1`$ to $`0`$ at $`\nu =0`$ (van Dongen and Ernst 1985c, Krivitsky 1995).

For the cases with $`\mu >0`$ (figure 5(c)), the scaling function $`\phi (x)`$ oscillates around a power law at small $`x`$, with the fractional amplitude of the oscillation decreasing as $`x0`$. This damped oscillatory approach to a power law is shown more clearly in figure 7, where we plot the logarithmic slope $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ as a function of $`\mathrm{log}x`$ for the cases with $`\mu =1/6`$. In all cases, the logarithmic slope tends to a constant value as $`x0`$, and the asymptotic value is consistent with the leading small-$`x`$ behavior predicted by self-similar analysis: $`\phi (x)x^{(1+\lambda )}`$ or $`d\mathrm{ln}\phi /d\mathrm{ln}x=(1+\lambda )`$. The oscillation appears to be periodic in the variable $`\mathrm{ln}x`$, but it is difficult to determine this accurately because of the limited number of cycles seen in our numerical results.

The damped oscillatory behavior at small $`x`$ for $`\mu >0`$ and $`\lambda <1`$ was not known previously. In their self-similar analysis, van Dongen and Ernst (1988) were unable to find higher-order corrections to the leading small-$`x`$ behavior of $`\phi (x)`$ for $`\mu >0`$ and $`\lambda <1`$, and they raised the possibility that physically acceptable self-similar solutions may not exist. As we have just shown, there are indeed physically acceptable self-similar solutions, and they are reached from monodisperse initial conditions. Based on our numerical results, we suggest that the leading small-$`x`$ behavior and the first correction could be of the form $`\phi (x)x^{(1+\lambda )}[1+f(x)\mathrm{cos}(B\mathrm{ln}x+C)]`$, where $`f(x)`$ is an increasing function of $`x`$, possibly $`Ax^\alpha `$ with $`\alpha >0`$. The failure to find the first correction in the self-similar analysis is probably due to the unusual form of this correction.

As we mentioned in Section 1, Krivitsky (1995) has obtained numerical solutions to the continuous form of the coagulation equation for $`A_{m,m^{}}(mm^{})^{\lambda /2}`$ with $`\lambda 1`$. Krivitsky concluded that the numerical solutions are self-similar at large time but that the asymptotic behavior at small mass is not a power law. The latter conclusion is different from ours and, we believe, incorrect. We can understand why Krivitsky reached this conclusion by examining figure 5(b) of Krivitsky (1995), where the numerical results for $`d\mathrm{ln}n/d\mathrm{ln}m`$ are shown for the case $`\lambda =0.4`$. By the last time shown, the numerical results have converged to a self-similar form for $`x\mathrm{}>10^4`$, and the logarithmic slope is indeed decreasing with decreasing $`x`$ over the range $`10^4\mathrm{}<x\mathrm{}<10^1`$. However, as we can see from our figure 7, over this range in $`x`$, $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ is in fact decreasing from a maximum to a minimum in its oscillatory approach to a constant value. Therefore, the incorrect conclusion was reached because Krivitsky did not evolve the numerical solutions for a sufficiently long time to see the true asymptotic behavior at small mass.

4.2. Orderly Growth Cases with $`\nu 1`$ and $`\lambda =1`$

On the borderline $`\lambda =1`$ and $`\nu 1`$, the numerical mass spectrum evolves in an orderly fashion, but with the characteristic mass $`m_{}(t)`$ increasing exponentially with time. For these cases, we set $`N_{\mathrm{dec}}=19`$ and stopped the numerical integrations as soon as $`N_{\mathrm{bin}}`$ reached $`N_{\mathrm{max}}`$, i.e., when the mass spectra extended over $`20`$ decades in mass.

We distinguish the cases with $`\mu =0`$ and $`\mu >0`$. As we discussed above, the case $`\mu =0`$ is $`A_{ij}=(i+j)/2`$ with exact analytic solution, and the numerical solution for this case is in excellent agreement with the analytic solution (see Section 3 and figures 2 and 5(b)). The mass spectrum tends toward a self-similar solution of the form (8) at large $`t`$, with $`m_{}(t)e^t`$ (see figure 8), and the scaling function $`\phi (x)x^{3/2}`$ at small $`x`$ and $`x^\theta \mathrm{exp}(bx)`$ with $`\theta =(\mu +3)/2=3/2`$ at large $`x`$.

For $`\mu >0`$, van Dongen and Ernst (1988) have derived a modified self-similar solution:

$$n_k(t)=(m_{}^2\mathrm{ln}m_{})^1\phi (m_k/m_{})$$
$`(10)`$
where $`(\mathrm{ln}m_{})^2=a+bt`$ and $`a`$ and $`b`$ are constants. The scaling function $`\phi (x)`$ is predicted to scale as $`x^2`$ at small $`x`$ and $`x^1\mathrm{exp}(bx)`$ at large $`x`$.

The numerical results for $`m_{}(t)`$ for the three cases with $`\mu >0`$ (and also the case $`\mu =0`$) are shown in figure 8. In each case, we have fitted the numerical results at the last two output times to $`\mathrm{ln}m_{}=a+bt`$ and $`(\mathrm{ln}m_{})^2=a+bt`$, and they are shown as dotted and solid lines, respectively. For $`(\mu ,\nu )=(1/2,1/2)`$, $`(\mathrm{ln}m_{})^2=a+bt`$ provides a reasonably good fit to the numerical results at large $`t`$. For $`(\mu ,\nu )=(1/3,2/3)`$ and, in particular, $`(1/6,5/6)`$, the numerical results at large $`t`$ show deviations from $`(\mathrm{ln}m_{})^2=a+bt`$. The deviations in the last two cases are probably due to the numerical results not having completely reached the asymptotic regime, but we cannot rule out the possibility that the asymptotic behavior is different from the analytic prediction.

If a numerical solution approaches the self-similar solution (10), we expect $`m^2n\mathrm{ln}m_{}x^2\phi (x)`$ and $`d\mathrm{ln}n/d\mathrm{ln}md\mathrm{ln}\phi /d\mathrm{ln}x`$. In figure 9 we show $`d\mathrm{ln}n/d\mathrm{ln}m`$ as a function of $`\mathrm{log}m/m_{}`$ at three different times for all cases with $`\mu >0`$ and $`\lambda =1`$. For $`(\mu ,\nu )=(1/2,1/2)`$, $`d\mathrm{ln}n/d\mathrm{ln}m`$ has converged to $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ at $`x=m/m_{}\mathrm{}>10^5`$, but the convergence at $`x\mathrm{}<10^5`$ is very slow and is not complete by the end of the numerical run. The range of $`x`$ over which $`d\mathrm{ln}n/d\mathrm{ln}m`$ has converged by the end of the numerical run is wider for smaller $`\mu `$. (A similar analysis of $`m^2n\mathrm{ln}m_{}`$ reveals a small increase in the normalization of $`m^2n\mathrm{ln}m_{}`$ with time in the range of $`x`$ where $`d\mathrm{ln}n/d\mathrm{ln}m`$ has converged. This increase is more pronounced for smaller $`\mu `$ and is probably due to $`m_{}`$ not having completely reached the asymptotic regime.) For $`(\mu ,\nu )=(1/6,5/6)`$, it is reasonably clear that $`d\mathrm{ln}\phi /d\mathrm{ln}x2`$ in the small-$`x`$ limit, consistent with the analytic prediction. For $`(\mu ,\nu )=(1/3,2/3)`$ and $`(1/2,1/2)`$, the small-$`x`$ behaviors are less certain because of the slow convergence at small $`x`$, but they also appear to be consistent with the analytic prediction. Finally, in all three cases, the large-$`x`$ behavior of the numerical $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ is consistent with the analytic prediction that $`d\mathrm{ln}\phi /d\mathrm{ln}x=1bx`$.

4.3. Runaway Growth Cases with $`\nu 1`$ and $`\lambda >1`$

In all cases with $`\nu 1`$ and $`\lambda >1`$, the numerical mass spectrum develops in a finite time $`t_0`$ a power-law tail, $`n_kk^\tau `$, at large mass that violates mass conservation, and runaway growth is expected to start at $`t_{\mathrm{crit}}=t_0`$ in the limit $`n_0\mathrm{}`$. For most cases with $`\nu 1`$ and $`\lambda >1`$, we stopped the numerical integrations as soon as $`N_{\mathrm{bin}}`$ reached $`N_{\mathrm{max}}`$ and the mass spectra extended over $`20`$ decades in mass; so the numerical solutions approach very close to but do not exceed $`t_0`$. To study the transition at $`t=t_0`$, we extended the integrations for the cases $`(\mu ,\nu )=(1/3,1)`$ and $`(2/3,2/3)`$ to $`t>t_0`$. Figure 10 shows the numerical results for the mass spectrum evolution for the case $`(\mu ,\nu )=(1/3,1)`$ at $`tt_0`$ (solid lines) and $`t>t_0`$ (dashed lines).

For runaway growth with $`\nu 1`$ and $`\lambda >1`$, self-similar analysis predicts that self-similar solutions close to but before $`t_0`$ have the form

$$n_k(t)=m_{}(t)^\tau \phi [m_k/m_{}(t)]$$
$`(11)`$
where the scaling function $`\phi (x)x^\tau `$ at small $`x`$, the characteristic mass $`m_{}(t)`$ diverges as $`(t_0t)^{1/\sigma }`$, and $`\sigma =\lambda +1\tau `$. With $`m_{}(t)`$ diverging at $`t_0`$, the self-similar solution (11) has $`n_k(t_0)k^\tau `$ at large $`k`$. In all cases, the numerical solution tends toward a self-similar solution of the form (11) as $`tt_0`$, and the numerical $`\phi (x)`$ is indeed a power law at small $`x`$. This is illustrated in figure 11, where we plot the numerical solution at three different times ($`<t_0`$) in the form of $`\mathrm{log}(m^2nm_{}^{\tau 2})`$ as a function of $`\mathrm{log}m/m_{}`$ for the case shown in figure 10, with the exponent $`\tau `$ determined from the numerical solution itself (see table 1 and discussion below). In the scaling limit (11), $`m^2nm_{}^{\tau 2}=x^2\phi (x)`$.

Despite the change in the form of the self-similar solution, the analytic predictions for the large-$`x`$ behavior of $`\phi (x)`$ are similar to those for orderly growth in Sections 4.1–4.2: $`\phi (x)x^\theta \mathrm{exp}(bx)`$ or $`d\mathrm{ln}\phi /d\mathrm{ln}x=\theta bx`$, where $`\theta =\lambda `$ if $`\nu <1`$, $`\theta =(\mu +3)/2`$ if $`\nu =1`$ and $`0<\mu <1`$, and $`\theta =5/2`$ if $`\nu =\mu =1`$. In all cases, the numerical $`\phi (x)`$ decays exponentially at large $`x`$, and $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ at large $`x`$ is consistent with the analytic result.

The analytic predictions for the exponents $`\tau `$ and $`\sigma `$ are $`\tau =(\lambda +3)/2`$ and $`\sigma =\lambda +1\tau =(\lambda 1)/2`$ (Hendriks et al 1983, van Dongen and Ernst 1985a, 1988). For comparison, we have determined $`\tau `$ and $`\sigma `$ from the numerical results for $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ and $`m_{}(t)`$, respectively. In figure 12 we show the numerical results for $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ for the cases with $`\lambda =4/3`$. The numerical results clearly converge to constant values at small $`x`$. The constant asymptotic values are consistent with $`\phi (x)x^\tau `$ or $`d\mathrm{ln}\phi /d\mathrm{ln}x=\tau `$ at small $`x`$ and directly yield the values of $`\tau `$. It is also clear from figure 12 that the asymptotic values and hence $`\tau `$ for these cases with the same $`\lambda `$($`=4/3`$) are different from each other and from the analytic prediction that $`d\mathrm{ln}\phi /d\mathrm{ln}x=\tau =(\lambda +3)/2=13/6`$ (dashed line in figure 12). The numerical results for the exponent $`\tau `$ for all cases with $`\nu 1`$ and $`\lambda >1`$ are listed in table 1. The only case where the exponent $`\tau `$ agrees with the analytic prediction $`(\lambda +3)/2`$ is the case $`\nu =\mu =1`$ (i.e., the case $`A_{ij}=ij`$ with exact analytic solution). In all other cases, the exponent $`\tau `$ is less than $`(\lambda +3)/2`$. For a given $`\lambda `$, the deviation of $`\tau `$ from $`(\lambda +3)/2`$ is largest for $`\nu =1`$ and smallest for $`\nu =\mu `$.

To determine the exponent $`\sigma `$, we fit the numerical results for $`m_{}(t)`$ at output times $`t_{i1}`$, $`t_i`$, and $`t_{i+1}`$ to $`m_{}(t)=C(t_0^{}t)^{1/\sigma ^{}}`$ to obtain $`C(t_i)`$, $`t_0^{}(t_i)`$, and $`\sigma ^{}(t_i)`$. In most cases, $`\sigma ^{}(t)`$ has converged to a constant value by the end of the numerical run and directly yields $`\sigma `$. In the remaining cases, we extrapolated $`\sigma ^{}(t)`$ to the limit $`m_{}(t)\mathrm{}`$ to obtain $`\sigma `$, but the difference between $`\sigma ^{}(t)`$ at the end of the numerical run and $`\sigma `$ is $`0.002`$. The numerical results for the exponent $`\sigma `$ are listed in table 1, together with $`\lambda +1\tau `$. Except for the case $`\nu =\mu =1`$, the exponent $`\sigma `$ is greater than the analytic prediction $`(\lambda 1)/2`$. For a given $`\lambda `$, the deviation of $`\sigma `$ from $`(\lambda 1)/2`$ is largest for $`\nu =1`$ and smallest for $`\nu =\mu `$. We note that the numerical results for $`\sigma `$ and $`\tau `$ are consistent with one another in that they satisfy the relation $`\sigma =\lambda +1\tau `$ for self-similar solutions of the form (11) to within $`\pm 0.001`$.

The procedure described above for determining the exponent $`\sigma `$ also yields the time $`t_0`$. The numerical results for $`t_0`$, in units of $`1/(n_0A_{11})`$, are listed in table 1. Since we expect $`t_01/[(\lambda 1)n_0A_{11}]`$ (see Paper I), it is convenient to parameterize $`t_0`$ as $`t_0=K/[(\lambda 1)n_0A_{11}]`$. The parameter $`K`$ is shown in figure 13. We find that $`K=1`$$`2`$ if $`\lambda \mathrm{}>1.1`$ and that, for a given $`\lambda `$, $`K`$ is smallest for $`\nu =1`$ and largest for $`\nu =\mu `$. For $`\nu =1`$, the parameter $`K`$ shows a maximum at $`\lambda 1.3`$. For $`\nu =\mu `$, $`K`$ increases monotonically with decreasing $`\lambda `$, but it is unclear whether $`K`$ approaches a constant value or diverges as $`\lambda 1`$. Finally, we note that the numerical results for $`t_0`$ are consistent with the bound $`t_01/[(\lambda 1)n_0A_{11}]`$ for $`A_{ij}`$ in equation (3) and monodisperse initial conditions and the stronger bound $`t_01/[(2^{\lambda 1}1)n_0A_{11}]`$ for $`\nu =\mu `$, derived analytically by Hendriks et al (1983).

We have found that the exponents $`\tau `$ and $`\sigma `$ are in general different from the analytic predictions. Since the analytic prediction for $`\sigma `$ follows from that for $`\tau `$ and the relation $`\sigma =\lambda +1\tau `$ for self-similar solutions of the form (11), and the numerical results for $`\sigma `$ and $`\tau `$ satisfy this relation, we have essentially a discrepancy in the exponent $`\tau `$. Let us examine the derivation of the analytic prediction for $`\tau `$, which can be summarized as follows (see Hendriks et al 1983, van Dongen and Ernst 1985a, 1988 for details). The mass flux from particles of mass $`m_im_k`$ to particles of mass $`m_i>m_k`$ is

$$F_k(t)=\underset{i=1}{\overset{k}{}}\underset{j=k+1i}{\overset{\mathrm{}}{}}m_iA_{ij}n_i(t)n_j(t).$$
$`(12)`$
It is assumed that solutions to the coagulation equation at $`t>t_0`$ violate mass conservation by having a non-zero and finite mass flux to the infinite-mass bin (i.e., in the limit $`k\mathrm{}`$), which is possible only if the mass spectrum at $`t>t_0`$ is of the form $`n_k(t)k^\tau ^{}`$ at large $`k`$. With $`n_k(t)k^\tau ^{}`$ at large $`k`$, the mass flux $`F_k(t)k^{\lambda +32\tau ^{}}`$ at large $`k`$. Thus $`\tau ^{}=(\lambda +3)/2`$ if the mass flux $`F_k(t)`$ is required to be non-zero and finite in the limit $`k\mathrm{}`$. Since the self-similar solution (11) has $`n_k(t_0)k^\tau `$ at large $`k`$, $`\tau =\tau ^{}=(\lambda +3)/2`$ if we assume that the large-$`k`$ behavior at $`t>t_0`$ is also valid at $`t=t_0`$. It should be noted that the arguments leading to $`\tau =(\lambda +3)/2`$ are not rigorous. In particular, as van Dongen and Ernst (1988) pointed out, one cannot exclude the possibility that the mass flux diverges at one instant of time, i.e., $`t_0`$. We have found that $`\tau <(\lambda +3)/2`$ in general. This implies that the mass flux $`F_k`$ at $`t_0`$, which is $`k^{\lambda +32\tau }`$ at large $`k`$, diverges as $`k\mathrm{}`$. Thus the numerical solutions violate mass conservation at $`t_0`$ by having a diverging mass flux to the infinite-mass bin.

Since the mass flux cannot diverge for all times $`t>t_0`$, do the solutions at $`t>t_0`$ have the form $`n_k(t)k^{(\lambda +3)/2}`$ at large $`k`$ for non-zero and finite mass flux to the infinite-mass bin? If so, how can this large-$`k`$ behavior at $`t>t_0`$ be reconciled with that at $`t=t_0`$? To answer these questions, we have extended the numerical integrations for the cases $`(\mu ,\nu )=(1/3,1)`$ and $`(2/3,2/3)`$ to $`t>t_0`$. The results for the mass spectrum evolution at $`t>t_0`$ for the case $`(\mu ,\nu )=(1/3,1)`$ are shown as dashed lines in figure 10. As in the test case $`A_{ij}=ij`$ discussed in Section 3, the numerical solution at $`m_km_{\mathrm{max}}`$ is affected by the finite value of the maximum particle mass, $`m_{\mathrm{max}}`$, in the computational grid. Otherwise, the mass spectrum at $`t>t_0`$ is indeed of the form $`n_k(t)k^{(\lambda +3)/2}`$ at large $`k`$, with the value of the exponent confirmed by an analysis of the logarithmic slope. Note, however, that the range of $`k`$ where $`n_k(t)k^{(\lambda +3)/2}`$ shrinks toward infinity as $`tt_0`$. The numerical solution for the case $`(\mu ,\nu )=(2/3,2/3)`$ shows the same large-$`k`$ behaviors at $`t>t_0`$. Therefore, we conclude that the transition at $`t=t_0`$ is accomplished as follows. As $`tt_0`$, the solution tends toward a self-similar solution of the form (11), with $`n_k(t)k^\tau `$ for $`1m_km_{}(t)`$, $`\tau <(\lambda +3)/2`$ in general, and $`m_{}(t)\mathrm{}`$. As $`tt_0`$, $`n_k(t)k^{(\lambda +3)/2}`$ for $`m_km_{}^{}(t)`$ and $`m_{}^{}(t)\mathrm{}`$.

4.4. Runaway Growth Cases with $`\nu >1`$

Examples of the numerical results for the mass spectrum evolution with $`n_{\mathrm{min}}=10^{30}`$ (solid lines) and $`10^{35}`$ (dotted lines) for $`\nu >1`$ are shown in the lower panels of figure 14. In all cases with $`\nu =2`$ and $`3/2`$ and the three cases with $`\nu =7/6`$ and $`\mu 0`$ (see, e.g., figure 14(a)), the behaviors of the numerical solutions for both $`n_{\mathrm{min}}=10^{30}`$ and $`10^{35}`$ are qualitatively similar to those found in Paper I for other forms of $`A_{ij}`$ with $`\nu >1`$. In the early stages, the mass spectrum appears to decay exponentially at large mass. After some evolution, the mass spectrum shows an exponential drop in an intermediate mass range and a power-law tail of the form $`n_kk^\nu `$ (or $`nm^\nu `$) at large mass. For the three cases with $`\nu =7/6`$ and $`\mu 1/3`$, the numerical solutions with $`n_{\mathrm{min}}=10^{30}`$ do not develop the $`m^\nu `$ tail when $`N_{\mathrm{bin}}`$ reaches $`N_{\mathrm{max}}`$ (see, e.g., figure 14(b)), and their behaviors are qualitatively similar to those for the runaway growth cases with $`\nu 1`$ and $`\lambda >1`$ (Section 4.3; figure 10). When $`n_{\mathrm{min}}`$ is decreased to $`10^{35}`$, the $`\mu =1/3`$ case does develop an $`m^\nu `$ tail when $`N_{\mathrm{bin}}`$ reaches $`N_{\mathrm{max}}`$ (figure 14(b)), but the other two cases do not. It is not feasible to compute solutions for much smaller $`n_{\mathrm{min}}`$, but we strongly suspect that the remaining two cases would develop the $`m^\nu `$ tail for sufficiently small $`n_{\mathrm{min}}`$.

The fact that the tail at large mass is of the form $`nm^\nu `$ or $`d\mathrm{ln}n/d\mathrm{ln}m=\nu `$ is illustrated in the upper panels of figure 14, where we plot the numerical results for $`d\mathrm{ln}n/d\mathrm{ln}m`$ at the specified times, just after the formation of the tail for the runs with $`n_{\mathrm{min}}=10^{35}`$. As found in Paper I, the time at which the power-law tail develops decreases slowly toward zero as $`n_{\mathrm{min}}`$ decreases (lower panels of figure 14). We do not repeat the demonstrations here, but it was shown in Paper I that numerical solutions with different maximum particle mass, $`m_{\mathrm{max}}`$, included in the computational grid are identical in the overlapping mass range and that the power-law tail simply extends to larger particle mass $`m`$ when $`m_{\mathrm{max}}`$ is increased. (It was also shown that the numerical solutions converge by $`N_{bd}=40`$ and $`\delta _M=5\%`$.) Therefore, in the limit $`n_{\mathrm{min}}0`$ and $`m_{\mathrm{max}}\mathrm{}`$, the numerical solutions for all cases with $`\nu >1`$ (with the possible exception of the cases with $`\nu =7/6`$ and $`\mu >1/3`$) should develop in an infinitesimal time power-law tails of the form $`nm^\nu `$ that extend to arbitrarily large $`m_{\mathrm{max}}`$. However, this mass spectrum is not self-consistent because the power-law tail implies a mass flux and, if $`1<\nu 2`$, a cumulative mass that diverge with $`m_{\mathrm{max}}`$ (see footnote 2 and Paper I). Thus, as in Paper I, the numerical results strongly suggest that there are no self-consistent solutions to the coagulation equation at any time if $`\nu >1`$. Furthermore, since the formation of the power-law tail in the coagulation equation solution with finite $`n_{\mathrm{min}}`$ probably corresponds to the onset of runaway growth in Monte Carlo simulations with finite $`n_0`$ (see Paper I) and decreasing $`n_{\mathrm{min}}`$ is equivalent to increasing $`n_0`$ (see Section 3), the time $`t_{\mathrm{crit}}`$ for the onset of runaway growth for $`\nu >1`$ should decrease slowly toward zero with increasing $`n_0`$.

5. CONCLUSIONS

We have conducted a systematic survey of numerical solutions to the coagulation equation (1) for a rate coefficient of the form $`A_{ij}(i^\mu j^\nu +i^\nu j^\mu )`$ and monodisperse initial conditions. The numerical results confirmed that there are three classes of rate coefficients with qualitatively different solutions to the coagulation equation.

For $`\nu 1`$ and $`\lambda =\mu +\nu <1`$ (Section 4.1), the numerical solution evolves in an orderly fashion and tends toward a self-similar solution of the form (8) at large time $`t`$. The time dependence of the characteristic mass $`m_{}(t)`$ at large $`t`$ and the behaviors of the scaling function $`\phi (x)`$ at small and large $`x=m_k/m_{}(t)`$ agree with the analytic predictions of van Dongen and Ernst (1985a, 1988). In particular, for the subset of cases with $`\mu >0`$, we disagreed with the earlier numerical study by Krivitsky (1995) and found that $`\phi (x)`$ does in fact approach the analytically predicted power-law behavior $`\phi (x)x^{(1+\lambda )}`$ at small $`x`$, but in a damped oscillatory fashion that was not known previously (figure 7). For $`\mu =0`$, we determined the exponent $`\tau `$ of the power-law behavior $`\phi (x)x^\tau `$ at small $`x`$.

On the borderline $`\nu 1`$ and $`\lambda =1`$ (Section 4.2), the numerical solution evolves in an orderly fashion. For $`\mu =0`$, i.e., $`A_{ij}i+j`$, the numerical solution is in excellent agreement with the exact analytic solution and tends toward a self-similar solution of the form (8) at large $`t`$. For $`\mu >0`$, the numerical solution appears to tend toward a self-similar solution of the form (10), but we had limited success in comparing the behaviors of $`m_{}(t)`$ and $`\phi (x)`$ at small $`x`$ to the analytic predictions because the convergence to the self-similar solution is very slow.

For $`\nu 1`$ and $`\lambda >1`$ (Section 4.3), the numerical mass spectrum $`n_k`$ develops in a finite time $`t_0`$ a power-law tail, $`n_kk^\tau `$, at large $`k`$ that violates mass conservation, and runaway growth/gelation is expected to start at $`t_{\mathrm{crit}}=t_0`$ in the limit that the initial number of particles $`n_0\mathrm{}`$. As $`tt_0`$, the numerical solution tends toward a self-similar solution of the form (11), with $`\phi (x)x^\tau `$ at small $`x`$ and $`m_{}(t)`$ diverging as $`(t_0t)^{1/\sigma }`$. The exponent $`\tau `$ is in general less than the analytic prediction $`(\lambda +3)/2`$, and the exponent $`\sigma `$ greater than the analytic prediction $`(\lambda 1)/2`$, but they satisfy the relation $`\sigma =\lambda +1\tau `$ for self-similar solutions of the form (11) (table 1). We studied the dependence of $`t_0`$ on the exponents $`\mu `$, $`\nu `$, and $`\lambda `$ and found that $`t_0=K/[(\lambda 1)n_0A_{11}]`$ and $`K=1`$$`2`$ if $`\lambda \mathrm{}>1.1`$ (table 1; figure 13). At $`t>t_0`$, $`n_kk^{(\lambda +3)/2}`$ for $`m_km_{}^{}(t)`$, with $`m_{}^{}(t)\mathrm{}`$ as $`tt_0`$.

For $`\nu >1`$ (Section 4.4), the behaviors of the numerical solution are qualitatively similar to those found in Paper I: the numerical mass spectrum develops a power-law tail of the form $`n_kk^\nu `$ at large $`k`$ that is not self-consistent, and the time at which the power-law tail develops decreases toward zero as the numerical parameter $`n_{\mathrm{min}}`$ decreases. The numerical results strongly suggest that there are no self-consistent solutions to the coagulation equation at any time and that runaway growth/gelation is instantaneous in the limit $`n_0\mathrm{}`$. They also indicate that the time $`t_{\mathrm{crit}}`$, in units of $`1/(n_0A_{11})`$, for the onset of runaway growth decreases slowly toward zero with increasing $`n_0`$, consistent with recent Monte Carlo simulation results (Malyshkin and Goodman 2001).

The results presented in this paper and Paper I suggest several problems for future investigations. First, as we pointed out in Section 4.1, the failure to find the first correction to the leading small-$`x`$ behavior $`\phi (x)x^{(1+\lambda )}`$ for the orderly growth cases with $`\mu >0`$ and $`\lambda <1`$ in previous self-similar analysis is probably due to the unusual, damped oscillatory form of this correction. Given the information provided by the numerical results, it may now be possible to derive the first correction analytically. We suggested that the first correction could, e.g., be of the form $`x^{(1+\lambda )}f(x)\mathrm{cos}(B\mathrm{ln}x+C)`$, where $`f(x)`$ is an increasing function of $`x`$. Second, we have found that the exponent $`\tau `$ for the runaway growth cases with $`\nu 1`$ and $`\lambda >1`$ is in general different from existing analytic prediction. It is important to investigate whether analytic prediction matching the numerical results could be derived. In this connection, it would be useful to obtain numerical solutions for other forms of $`A_{ij}`$ and determine whether the exponent depends on the specific form of $`A_{ij}`$. Finally, it should be emphasized that rate coefficients with $`\nu >1`$ do arise, and are of great interest, in astrophysics (see, e.g., Lee 1993, 2000, and references therein). For astrophysical applications, we are interested in systems with finite $`n_0`$ and interactions with the runaway particle. Thus, for $`\nu >1`$, a detailed comparison of numerical solutions with finite $`n_{\mathrm{min}}`$ and Monte Carlo simulations with finite $`n_0`$ should be conducted to test the correspondence between them, and the question of whether the coagulation equation can be modified to take into account the interactions between the runaway particle and the other particles should also be investigated.

ACKNOWLEDGEMENTS

I thank Martin Duncan, Stan Peale, and the referees for their helpful comments on the manuscript. This work was supported in part by the NASA PG&G Program under Grant NAG5-3646.

Table 1. Runaway growth cases with $`\nu 1`$ and $`\lambda >1`$.

| $`\lambda `$ | $`\mu `$ | $`\nu `$ | $`\tau `$ | $`\sigma `$ | $`\lambda +1\tau `$ | $`n_0A_{11}t_0`$ |
| --- | --- | --- | --- | --- | --- | --- |
|  |  |  |  |  |  |  |
| $`7/6`$ | $`1/6`$ | $`1`$ | $`2.012`$ | $`0.154`$ | $`0.155`$ | $`7.136`$ |
| $`7/6`$ | $`1/3`$ | $`5/6`$ | $`2.038`$ | $`0.128`$ | $`0.129`$ | $`9.146`$ |
| $`7/6`$ | $`1/2`$ | $`2/3`$ | $`2.054`$ | $`0.112`$ | $`0.113`$ | $`10.633`$ |
| $`7/6`$ | $`7/12`$ | $`7/12`$ | $`2.057`$ | $`0.110`$ | $`0.110`$ | $`10.855`$ |
| $`4/3`$ | $`1/3`$ | $`1`$ | $`2.076`$ | $`0.257`$ | $`0.257`$ | $`3.741`$ |
| $`4/3`$ | $`1/2`$ | $`5/6`$ | $`2.103`$ | $`0.230`$ | $`0.230`$ | $`4.284`$ |
| $`4/3`$ | $`2/3`$ | $`2/3`$ | $`2.112`$ | $`0.221`$ | $`0.221`$ | $`4.492`$ |
| $`3/2`$ | $`1/2`$ | $`1`$ | $`2.166`$ | $`0.334`$ | $`0.334`$ | $`2.451`$ |
| $`3/2`$ | $`2/3`$ | $`5/6`$ | $`2.184`$ | $`0.316`$ | $`0.316`$ | $`2.639`$ |
| $`3/2`$ | $`3/4`$ | $`3/4`$ | $`2.186`$ | $`0.314`$ | $`0.314`$ | $`2.664`$ |
| $`5/3`$ | $`2/3`$ | $`1`$ | $`2.269`$ | $`0.398`$ | $`0.398`$ | $`1.750`$ |
| $`5/3`$ | $`5/6`$ | $`5/6`$ | $`2.276`$ | $`0.390`$ | $`0.391`$ | $`1.808`$ |
| $`11/6`$ | $`5/6`$ | $`1`$ | $`2.380`$ | $`0.453`$ | $`0.453`$ | $`1.307`$ |
| $`11/6`$ | $`11/12`$ | $`11/12`$ | $`2.381`$ | $`0.451`$ | $`0.452`$ | $`1.317`$ |
| $`2`$ | $`1`$ | $`1`$ | $`2.500`$ | $`0.500`$ | $`0.500`$ | $`1.000`$ |

REFERENCES

Bak T A and Heilmann O J 1994 Post-gelation solutions to Smoluchowski’s coagulation equation J. Phys. A: Math. Gen. 27 4203–4209

Drake R L 1972 A general mathematical survey of the coagulation equation Topics in Current Aerosol Research (Part 2) ed G M Hidy and J R Brock (Oxford: Pergamon) pp 201–376

Ernst M H 1986 Kinetics of clustering in irreversible aggregation Fractals in Physics ed L Pietronero and E Tosatti (Amsterdam: North-Holland) pp 289–302

Ernst M H, Hendriks E M and Leyvraz F 1984 Smoluchowski’s equation and the $`\theta `$-exponent for branched polymers J. Phys. A: Math. Gen. 17 2137–2144

Hendriks E M, Ernst M H and Ziff R M 1983 Coagulation equations with gelation J. Stat. Phys. 31 519–563

Hill P J and Ng K M 1996 New discretization procedure for the agglomeration equation AIChE J. 42 727–741

Jeon I 1998 Existence of gelling solutions for coagulation-fragmentation equations Commun. Math. Phys. 194 541–567

Jeon I 1999 Spouge’s conjecture on complete and instantaneous gelation J. Stat. Phys. 96 1049–1070

Jullien R and Botet R 1987 Aggregation and Fractal Aggregates (Singapore: World Scientific)

Krivitsky D S 1995 Numerical solution of the Smoluchowski kinetic equation and asymptotics of the distribution function J. Phys. A: Math. Gen. 28 2025–2039

Lee M H 1993 $`N`$-body evolution of dense clusters of compact stars Astrophys. J. 418 147–162

Lee M H 2000 On the validity of the coagulation equation and the nature of runaway growth Icarus 143 74–86

Leyvraz F and Tschudi H R 1981 Singularities in the kinetics of coagulation processes J. Phys. A: Math. Gen. 14 3389–3405

Malyshkin L and Goodman J 2001 The timescale of runaway stochastic coagulation Icarus 150 314–322

McLeod J B 1962 On a recurrence formula in differential equations Quart. J. Math. Oxford 13 283–284

Quinlan G D and Shapiro S L 1989 Dynamical evolution of dense clusters of compact stars Astrophys. J. 343 725–749

Spouge J L 1985 Monte Carlo results for random coagulation J. Colloid Interface Sci. 107 38–43

Trubnikov B A 1971 Solution of the coagulation equations in the case of a bilinear coefficient of adhesion of particles Sov. Phys. Dokl. 16 124–126

Tzivion S, Reisin T G and Levin Z 1999 A numerical solution of the kinetic collection equation using high spectral grid resolution: A proposed reference J. Comput. Phys. 148 527–544

van Dongen P G J 1987a On the possible occurrence of instantaneous gelation in Smoluchowski’s coagulation equation J. Phys. A: Math. Gen. 20 1889–1904

van Dongen P G J 1987b Solutions of Smoluchowski’s coagulation equation at large cluster sizes Physica A 145 15–66

van Dongen P G J and Ernst M H 1985a Dynamic scaling in the kinetics of clustering Phys. Rev. Lett. 54 1396–1399

van Dongen P G J and Ernst M H 1985b Comment on “Large-time behavior of the Smoluchowski equations of coagulation” Phys. Rev. A 32 670–672

van Dongen P G J and Ernst M H 1985c Cluster size distribution in irreversible aggregation at large time J. Phys. A: Math. Gen. 18 2779

van Dongen P G J and Ernst M H 1988 Scaling solutions of Smoluchowski’s coagulation equation J. Stat. Phys. 50 295–329

Wetherill G W 1990 Comparison of analytical and physical modeling of planetesimal accumulation Icarus 88 336–354

Ziff R M 1980 Kinetics of polymerization J. Stat. Phys. 23 241–263

Ziff R M, Ernst M H and Henriks E M 1983 Kinetics of gelation and universality J. Phys. A: Math. Gen. 16 2293–2320

Figure 1. The exponents $`\mu `$ and $`\nu `$ of the cases for which we have obtained numerical solutions to the discrete form of the coagulation equation with rate coefficient $`A_{ij}=(i^\mu j^\nu +i^\nu j^\mu )/2`$. The orderly growth cases with $`\nu 1`$ and $`\lambda <1`$ and on the borderline $`\nu 1`$ and $`\lambda =1`$ are discussed in Sections 4.1 and 4.2, respectively. The runaway growth cases with $`\nu 1`$ and $`\lambda >1`$ are discussed in Section 4.3, and those with $`\nu >1`$ are discussed in Section 4.4.

Figure 2. (a) Comparison of numerical (solid lines) and analytic (dotted lines) results for $`d\mathrm{ln}n/d\mathrm{ln}m`$ at $`m\mathrm{}<m_{}(t)`$. The lines, from top to bottom, are $`d\mathrm{ln}n/d\mathrm{ln}m+C`$ for $`A_{ij}=1`$ at $`t=10^8`$, $`A_{ij}=(i+j)/2`$ at $`t=18`$, and $`A_{ij}=ij`$ at $`t=1`$; the constant $`C=0`$, $`1.45`$, and $`2.4`$, respectively. For $`A_{ij}=(i+j)/2`$ and $`ij`$, there is a small lag in the evolution of the numerical solutions, and the numerical results are compared to the analytic results at a slightly earlier time $`(1ϵ)t`$, where $`ϵ=3.5\times 10^4`$ and $`1.1\times 10^4`$, respectively. (b) Same as (a), but at $`m\mathrm{}>m_{}(t)`$. The lines, in decreasing steepness, are $`d\mathrm{ln}n/d\mathrm{ln}m`$ for $`A_{ij}=1`$ at $`t=10^8`$, $`A_{ij}=(i+j)/2`$ at $`t=18`$, and $`A_{ij}=ij`$ at $`t=1`$.

Figure 3. Evolution of the mass spectrum for $`(\mu ,\nu )=(1/3,1/3)`$. We plot $`\mathrm{log}m^2n`$ as a function of $`\mathrm{log}m`$ since $`m^2n`$ is the total mass per unit logarithmic mass interval: $`m^2nd\mathrm{ln}m=mn𝑑m`$.

Figure 4. Approach to self-similar solution as $`t\mathrm{}`$ for $`(\mu ,\nu )=(1/3,1/3)`$. The numerical solution is plotted for $`t=10^3`$ (solid line), $`10^4`$ (dashed line), and $`10^5`$ (dotted line) in the form of $`\mathrm{log}m^2n`$ as a function of $`\mathrm{log}m/m_{}`$.

Figure 5. (a) Scaling function $`\phi (x)`$ for the cases with $`\mu <0`$ and $`\nu 1`$. The solid (dotted) lines in increasing width are for $`\mu =1/2`$ ($`1/6`$) and $`\nu =0,1/3,2/3,1`$. (b) $`\phi (x)`$ for the cases with $`\mu =0`$ and $`\nu 1`$. The solid lines in increasing width are for $`\nu =0,1/6,\mathrm{},1`$. (c) $`\phi (x)`$ for the cases with $`\mu >0`$ and $`\lambda <1`$. The solid lines in increasing width are for $`\mu =1/6`$ and $`\nu =1/6,1/3,1/2,2/3`$. The dotted lines, offset vertically by $`10`$, are for $`\mu =1/3`$ and $`\nu =1/3`$ and $`1/2`$.

Figure 6. Logarithmic slope $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ as a function of $`x^\mu `$ for the cases with $`\mu =1/6`$ and $`\nu =0,1/3,2/3,1`$ (solid lines from top to bottom). The dotted lines are the asymptotes approached by the numerical results at large $`x^\mu `$ (or small $`x`$).

Figure 7. Logarithmic slope $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ for the cases with $`\mu =1/6`$ and $`\nu =1/6,1/3,1/2,2/3`$ (solid lines from top to bottom). The dashed lines show the leading small-$`x`$ behavior predicted by self-similar analysis: $`d\mathrm{ln}\phi /d\mathrm{ln}x=(1+\lambda )`$.

Figure 8. Characteristic mass $`m_{}(t)`$ for the cases with $`\lambda =1`$ and $`\nu 1`$. The points, from left to right, are the numerical results for $`\mu =0,1/6,1/3`$, and $`1/2`$. The dotted and solid lines show $`\mathrm{ln}m_{}=a+bt`$ and $`(\mathrm{ln}m_{})^2=a+bt`$ fitted to numerical results at the last two output times.

Figure 9. Logarithmic slope $`d\mathrm{ln}n/d\mathrm{ln}m`$ for the cases with $`\mu >0`$ and $`\lambda =1`$. The top lines are the numerical results for $`(\mu ,\nu )=(1/6,5/6)`$ at $`t=47`$ (solid line), $`54`$ (dashed line), and $`61`$ (dotted line), offset vertically by $`C=1`$. The middle lines are for $`(\mu ,\nu )=(1/3,2/3)`$ at $`t=63`$, $`74`$, and $`87`$, with $`C=0.5`$, and the bottom lines are for $`(\mu ,\nu )=(1/2,1/2)`$ at $`t=71`$, $`86`$, and $`102`$, with $`C=0`$.

Figure 10. Evolution of the mass spectrum for $`(\mu ,\nu )=(1/3,1)`$ at $`tt_0`$ (solid lines) and $`t>t_0`$ (dashed lines).

Figure 11. Approach to self-similar solution as $`tt_0`$ for $`(\mu ,\nu )=(1/3,1)`$. The numerical solution is plotted for $`t=3.73`$ (solid line), $`3.74`$ (dashed line), and $`3.7409`$ (dotted line) in the form of $`\mathrm{log}(m^2nm_{}^{\tau 2})`$ as a function of $`\mathrm{log}m/m_{}`$.

Figure 12. Logarithmic slope $`d\mathrm{ln}\phi /d\mathrm{ln}x`$ for the cases with $`\lambda =4/3`$. The solid lines from top to bottom are the numerical results for $`(\mu ,\nu )=(1/3,1)`$, $`(1/2,5/6)`$, and $`(2/3,2/3)`$, and the dashed line is the analytic prediction for the small-$`x`$ behavior: $`d\mathrm{ln}\phi /d\mathrm{ln}x=\tau =(\lambda +3)/2=13/6`$.

Figure 13. Parameter $`K=(\lambda 1)n_0A_{11}t_0`$ for the cases with $`\nu 1`$ and $`\lambda >1`$ listed in table 1. The cases with $`\nu =1`$ and $`\nu =\mu `$ are represented by squares and triangles, respectively, and the remaining cases are represented by crosses.

Figure 14. Evolution of the mass spectrum (lower panels) and logarithmic slope $`d\mathrm{ln}n/d\mathrm{ln}m`$ at the specified time (upper panels) for $`\mu =1/3`$ and $`\nu `$ equal to (a) $`2`$ and (b) $`7/6`$. The numerical solutions with $`n_{\mathrm{min}}=10^{30}`$ (solid lines) and $`10^{35}`$ (dotted lines) are shown. In the upper panels, the dashed lines indicate $`d\mathrm{ln}n/d\mathrm{ln}m=\nu `$.


format: text

text:
# Atomic topology and radial distribution functions of a-SiNx alloys. Ab initio simulations

## I Introduction

The properties of amorphous silicon-nitrogen alloys have attracted a great deal of attention in the last decade; a-SiNx has electrical, optical and mechanical features useful in a variety of industrial applications and the strong covalency of their bonding makes them the prototype of covalent materials. Their optical gaps depend strongly on the nitrogen content $`x`$ for $`0x1.33`$ so they can be tuned to fit specific needs in solar cells; their total and partial radial distribution functions (RDFs) are practically unknown, except for the stoichiometric content. Some semiempirical studies have been done on their electronic structure, optical gaps and RDFs and a first-principles approach has been used on an ad hoc-generated amorphous structure. Therefore, any ab initio approach that adequately generates, describes and predicts features of a-SiNx may have a wider applicability to deal with other covalently-bonded amorphous solids; in particular, systems like a-GeNx should be amenable to our approach .

Recently  we carried out ab initio Harris-functional-based studies of the atomic and electronic structure of pure and hydrogenated amorphous silicon, using 64 silicon atom cells plus hydrogens that diffuse within the cells. We now apply these methods to amorphous silicon-nitrogen alloys to test the adequacy and the predictive powers of our approach for covalent materials. To the best of our knowledge, these are the first ab initio thermally generated amorphous networks where studies of the atomic topology of a-SiNx are carried out; the content $`x`$ is in the range $`0x1.29`$ where $`x=y/(64y)`$ and $`y`$ is the number of nitrogen atoms.

The experimental and theoretical activity prior to 1990 is well documented in a paper by Robertson where pertinent references can be found . In this work Robertson reports his semiempirical tight-binding calculations for the optical gaps and also reports several experimental gaps for hydrogenated and non-hydrogenated a-SiNx alloys. Bethe lattice calculations have been done by Martín-Moreno et al.  and by San-Fabián et al.  also using semiempirical parameters, whereas Ordejón and Ynduráin  in a very nice work do non-parameterized calculations of a-SiNx where the equilibrium positions of Si and N atoms in clusters are ported to the alloy network constructed ad hoc. They obtain a wealth of information including optical gaps, however tetrahedral coordination of the silicon atoms and threefold planar coordination of the nitrogen atoms is assumed with interatomic distances of 2.33 Å for Si-Si and 1.74-1.76 Å for Si-N. A general characteristic of these calculations/simulations is that gap states, when considered, are introduced either by hand, progressively replacing Si by N, or by the algorithms that generate the random networks, unlike the procedure reported here. Recent semiempirical classical simulations by de Brito et al.  produced total RDFs and average nearest-neighbor results that are the subject of comparison with ours.

On the experimental side, as early as 1976 Voskoboynikov et al.  studied some RDFs and optical gaps of silicon-rich silicon-nitrogen films as a function of the gas ratio. It was then observed that the gaps increase as a function of the nitrogen content; the films seemed to contain hydrogen and large clusters of silicon. However, reliable experimental RDFs are scarce  and, apparently, only total ones for the stoichiometric amorphous composition exist , where a decomposition of the second peak of the total RDF into its partial contributions is also carried out.

It is clear that the atomic topology also determines the electronic properties of the amorphous samples, and therefore any understanding of the RDFs and the atomic distribution in the random networks is relevant in the characterization of the electronic and optical properties of these materials. In what follows we report the generation of random networks for amorphous silicon nitrogen alloys that leads to RDFs in good agreement with what is known experimentally and predict atomic structures for a variety of nitrogen contents. Electronic gap states are found in these structures that shall be dealt with in a future publication.

## II Method

Our 26 amorphous samples of a-SiNx were generated with FastStructure , a DFT code based on the Harris functional produced by Molecular Simulations, Inc.. The optimization techniques use a fast force generator to allow simulated annealing/molecular dynamics studies with quantum force calculations . The LDA parameterization invoked is that due to Vosko, Wilk and Nusair . An all electron calculation is carried out, and a minimal basis set of atomic orbitals was chosen with a cutoff radius of 5 Å for the amorphization and 3 Å for the optimization. The physical masses of nitrogen and silicon are used throughout and this allows us to see realistic randomizing processes of the atoms during the amorphization of the supercell. Finally, the forces are calculated using rigorous formal derivatives of the expression for the energy in the Harris functional .

In order to test the adequacy of calculations carried out with FastStructure we used it to obtain the size of the crystalline cell of $`\beta `$-Si3N4 that minimizes the energy at the $`\mathrm{\Gamma }`$-point. Fig. 1 shows the results of such calculation. The experimental crystalline volumen is given by 145.920 Å3  whereas the calculated volumen is 146.797 Å3; a deviation of 0.6%. For this reason we feel cautiously optimistic about the use of FastStructure to generate random networks of silicon-nitrogen alloys.

Since it has become increasingly clear that quenching from a melt generates undesirable structures  we took a different path . We amorphisized the crystalline diamond structures with a total of 64 atoms ($`(64y)`$ silicons and $`y`$ substitutional nitrogens) in the cell by slowly heating it, linearly, from room temperature to just below the corresponding melting point for each $`x`$, and then slowly cooling it to 0 K. To determine the melting temperatures for each $`x`$ we linearly interpolated between the pure silicon value and the stoichiometric compound, $`x=4/3=1.33`$, and then remained below these temperatures (Table I). Since the time step was the same for all runs, 6 fs, and the melting temperatures increased with $`x`$, the heating/cooling rate varied from $`2.30\times 10^{15}`$K/s for pure silicon, to $`3.11\times 10^{15}`$K/s for $`x=1.29`$. The atoms were allowed to move freely within each cell whose volume was determined by the corresponding density and content, Table I. The densities were taken from the experimental results of Guraya et al. . Once this first stage was completed, we subjected each cell to annealing cycles with intermediate quenching processes. Finally the samples were energy-optimized to make sure the final structures would be those of a local energy minimum.

| TABLE I | | | |
| --- | --- | --- | --- |
| Contents, Melting Temperatures and Densities | | | |
| for a-SiNx | | | |
| Sample | $`x`$ | Melting | Density |
|  |  | Temp. (K) | (g/cc) |
| Si64N0 | 0.000 | 1680 | 2.329 |
| Si59N5 | 0.085 | 1747 | 2.435 |
| Si54N10 | 0.185 | 1814 | 2.512 |
| Si49N15 | 0.306 | 1881 | 2.600 |
| Si44N20 | 0.455 | 1948 | 2.694 |
| Si39N25 | 0.641 | 2015 | 2.803 |
| Si34N30 | 0.882 | 2082 | 2.931 |
| Si33N31 | 0.939 | 2095 | 2.957 |
| Si32N32 | 1.000 | 2108 | 2.988 |
| Si31N33 | 1.065 | 2122 | 3.017 |
| Si30N34 | 1.133 | 2136 | 3.048 |
| Si29N35 | 1.207 | 2149 | 3.081 |
| Si28N36 | 1.286 | 2162 | 3.115 |

It should be kept in mind that our objective is always to generate realistic amorphous structures of a-SiNx and not, in any way, to mimic the experimental processes used to produce these alloys.

## III Results and discussion

We performed two runs for each $`x`$ value and from $`x=0.882`$ on (a-Si34N30) the number of nitrogens was increased one at a time to be able to map the interesting processes that occur for these contents (percolation of the Si-Si bonds, widening of the optical gaps, etc). Once the atomic structures were obtained, we calculated the corresponding total and partial RDFs for each of the 26 runs and averaged them by corresponding pairs. Of those 13 averaged plots we are reporting here total and partial RDFs for pure silicon, for the almost stoichiometric sample, a-Si28N36, $`x=1.29`$, and two intermediate ones: a-Si44N20, $`x=0.46`$ and a-Si32N32, $`x=1.00`$. Fig 2a for the pure amorphous silicon sample shows our results and the upper and lower bounds of the available experimental data ; the agreement is good since our RDF falls within these bounds and the four experimental peaks are correctly reproduced by our simulations. Figs. 2b to 2d show the variation of the partial RDFs for Si-Si, Si-N and N-N as a function of content and their contribution to the total RDF. In Fig. 2d the composition of the second peak of the total RDF for the nearly a-Si3N4 can be observed and it agrees completely with experiment , since it is formed by the average second-neighbor $`2n`$ contributions of mainly the N-N and Si-Si partials and to a lesser extent by the Si-N partial. The third peak is essentially due to the Si-N partial with a small contribution from the N-N partial but, as far as we know, there are no experimental results for comparison.

Figs. 2 show that as the nitrogen content increases the first peak of the total RDF (1.85 Å), which is due to the Si-N average nearest-neighbor ($`nn`$) contributions, increases systematically and the $`nn`$ Si-Si peak (2.45 Å) decreases systematically. The third peak moves toward low $`r`$ as $`x`$ increases since the N-N contribution becomes more predominant at high content (3.25 Å to 2.95 Å). In our structures there are no $`nn`$ nitrogens since the content is below stoichiometry and nitrogens have a marked tendency to bind to silicons. For the nearly stoichiometric sample, $`x=1.29`$, the Si-Si $`nn`$ contribution to the total RDF has practically dissapeared and this implies that there is a nitrogen atom between every pair of silicons. This is borne out by the results presented in Figs. 3 where a study of the the average coordination numbers $`cn`$ in the 13 random networks is depicted. The following cutoff radii were used: Si-Si, 2.55 Å; N-N, 3.35 Å; and Si-N, 2.15 Å, which are the positions of the minima after the first peaks of the corresponding partials. Fig. 3a shows the results of our simulations where it can be seen that the N-N plot flattens for $`x1.1`$, the percolation threshold of Si-Si bonds ; the Si-Si $`nn`$ go from 4 to practically 0. The Si-N graph refers to the $`nn`$ nitrogens around the silicon atoms and varies from 0 to 4, whereas the N-Si refers to the $`nn`$ silicons around nitrogens and indicates that nitrogens immediately surround themselves with practically 3 Si, saturating its valence. The crossing of the Si-Si and Si-N plots at $`x0.7`$, is in agreement with experiment (Davis et al. ). There is a crossing of the Si-Si, N-Si and N-N plots at $`x0.3`$ and a crossing of N-Si and Si-N at $`x1.0`$ which have been observed experimentally for hydrogenated alloys by Guraya et al. , Fig. 3b. However, due to the presence of hydrogen a curvature appears in the $`cn`$ for Si-Si, Si-N and N-Si so, in order to compare our results to this experiment, we did the following. We carried out the sum of N-H plus N-Si from the experiment, the average total number of atoms that surround a N, N-$``$, and plotted it along with our N-Si; we also did the sum of the experimental Si-N, Si-Si and Si-H, the average total number of atoms that surround a Si, Si-$``$, and plotted that along with our sum of Si-N plus Si-Si. This is presented in Fig. 3c. It is clear that our predictions closely agree with the integrated experimental results and validate our approach. The discrepancies are most likely due to existing dangling and floating bonds.

Fig. 4 is the comparison of our ab initio results with the classical Monte Carlo simulations of de Brito et al.  where they used empirical potentials developed a la Tersoff for the interactions between Si and N. It is clear that although the positions of some peaks are reproduced in both simulations, the general behavior of the total RDFs only agree qualitatively. Comparison of each of these simulations with experiment is presented in Fig. 5 where agreements and discrepancies can be appreciated.

## IV Conclusions

We have devised an ab initio approach for a-SiNx ($`0x1.29`$) that generates radial distribution functions in agreement with existing experimental results and realistic atomic structures where the simulated average coordination numbers and the experimental ones coincide. Total RDFs agree very well with experiment, where available, and partial RDFs show that the Si-Si $`nn`$ peak disappears as nitrogen increases indicating a tendency to form 6-atom arrangements as the content $`x`$ approaches the stoichiometric value. Experiment shows that for a-Si3N4 Si and N form closed rings, Si-N-Si-N-Si-N, typical of the crystalline $`Si_3N_4`$ structures. The growth of the Si-N peak as nitrogen increases bears out this behavior. No nitrogen-nitrogen bonds, including molecular nitrogen, are observed in the final structures even though for $`x>1`$ the starting diamond structure does contain nitrogens next to one another. For $`x1.1`$ the effects of the percolation threshold of the Si-Si bonds is observed in the N-N $`2n`$. For $`x0.7`$ the Si-Si and Si-N neighbors are practically the same, as found experimentally. Also, Si-Si, N-Si an N-N are practically the same for $`x0.3`$ as are Si-N and N-Si for $`x1.0`$. The integrated experimental results and our simulation agree. The first prominent peak in the total RDF of the nearly stoichiometric sample is due to Si-N and an analysis of the second peak indicates that N-N, Si-Si and Si-N contribute to it, in agreement with experiment. The third peak is mainly due to the Si-N, with a small contribution from the N-N; no experimental results exist for comparison. Our approach, being ab initio, is of wider applicability than classical or semiempirical ones and should be relevant for the understanding of the physics of amorphous covalent materials.

###### Acknowledgements.

AAV thanks DGAPA-UNAM for financing the project IN101798 and a sabbattical stay at MSI-San Diego, USA, where this work was begun. FA thanks CONACyT for supporting his PhD studies. This work was done on an Origin 2000 computer provided by DGSCA, UNAM.

Figure1. Energy vs volume for crystalline beta-Si3N4 obtained with FastStructure and experimental values taken from Wild et al. (see text). The experimental and calculated volumes agree to within 0.6%.

Figure2. Total and partial RDFs for (a) pure silicon; the light curves are the upper and lower experimental bounds (Ref. ); (b) a-Si44N20, x=0.46; (c) a-Si32N32, x=1.00; (d) the nearly stoichiometric sample a-Si28N36, x=1.29.

Figure3. Average coordination numbers ¡cn¿ as a function of x. (a) Our results. (b) Experimental results for hydrogenated alloys from Guraya et al.. (c) Comparison of the integrated results (see text). Lines are drawn as guides to the eye.

Figure4. Comparison of our results and those obtained by de Brito et al. (Ref. ) for the total RDF. The agreement is at best qualitative although the position of some of the peaks coincide.

Figure5. (a) Comparison of our simulations and experiment (Ref. ) for the stoichiometric sample. (b) Comparison of de Brito’s simulations and experiment for the same sample.