\documentclass[twocolumn,aps,superscriptaddress,prl,floatfix,showpacs]{revtex4}
\usepackage{amsmath,amssymb}
\usepackage{graphicx}
\begin{document}
\title{Percolation with Multiple Giant Clusters}
\author{E.~Ben-Naim}\email{ebn@lanl.gov}
\affiliation{Theoretical Division and Center for Nonlinear
Studies, Los Alamos National Laboratory, Los Alamos, New Mexico
87545 USA}
\author{P.~L.~Krapivsky}\email{paulk@bu.edu}
\affiliation{Center for Polymer Studies and Department of Physics,
Boston University, Boston, Massachusetts 02215 USA}
\begin{abstract}
We study mean-field percolation with freezing. Specifically, we
consider cluster formation via two competing processes:
irreversible aggregation and freezing. We find that when the
freezing rate exceeds a certain threshold, the percolation
transition is suppressed. Below this threshold, the system
undergoes a series of percolation transitions with multiple giant
clusters (``gels'') formed. Giant clusters are not self-averaging
as their total number and their sizes fluctuate from realization
to realization. The size distribution $F_k$, of frozen clusters of
size $k$, has a universal tail, $F_k\sim k^{-3}$. We propose
freezing as a practical mechanism for controlling the gel size.
\end{abstract}

\pacs{82.70.Gg, 02.50.Ey, 05.40.-a}

\maketitle

Percolation was originally discovered in the context of
polymerization and gelation \cite{flory,stock} and since, it has
found numerous applications in physics \cite{ah}, geophysics
\cite{sah}, chemistry \cite{floryB}, and biology \cite{pg}. It
plays an important role in a vast array of natural and artificial
processes ranging from flow in porous media \cite{BH} and cloud
formation \cite{drake,sein} to evolution of random graphs
\cite{jlr,gg}, combinatorial optimization \cite{hart}, algorithmic
complexity \cite{bbckw}, amorphous computing, and DNA computing
using self-assembly \cite{win}.

We study mean-field percolation using the framework of
aggregation. An aggregation process typically begins with a huge
number of molecular units (``monomers'') that join irreversibly to
form clusters (``polymers'').  At some time, a giant cluster
(``gel'') containing a finite fraction of the monomers in the
system is born, and it grows to engulf the entire system. In this
classic percolation picture only a single gel forms, but in many
natural and artificial processes the system freezes into a
non-trivial final state with multiple gels or even micro-gels
\cite{ms,sv}. In this study, we show that aggregation with
freezing naturally leads to formation of multiple gels and that
freezing is a useful mechanism for controlling the gel size.

We analyze a process of aggregation with freezing where there are
two types of clusters: active and frozen. Active clusters join by
binary aggregation into larger active clusters.  The aggregation
rate is proportional to the product of the two cluster sizes
\cite{drake,fran,aal}; this is equivalent to the gelation model of
Flory and Stockmayer where a chemical bond between two monomers
joins their respective polymers \cite{flory,stock}. In parallel,
all active clusters including the gel may become frozen at a
size-independent constant rate $\alpha$.  Frozen clusters and
frozen gels are passive: they do not interact with other (active
or passive) clusters.

This process is studied using the rate equation approach. The density
$c_k(t)$ of active clusters of mass $k$ at time $t$ (that is, made up
from $k$ monomers) satisfies
\begin{equation}
\label{ck-eq}
\frac{d c_k}{dt}=\frac{1}{2}\sum_{i+j=k} (ic_i)(jc_j) -m\,kc_k-\alpha\,c_k,
\end{equation}
with $m(t)$ the mass density of active clusters.  The first two
terms on the right-hand side describe how the cluster size
distribution changes due to aggregation, and the last term
accounts for loss due to freezing. The quantity $m(t)$ equals the
mass density of {\em all} active clusters including possibly a
giant cluster. When there is no active giant cluster, $m(t)\equiv
M_1(t)$ with the moments defined via $M_n(t)\equiv \langle
k^n\rangle =\sum_{k\geq 1} k^nc_k(t)$.

\smallskip
\noindent{\it The gelation transition.} Initially, all clusters are
finite in size, so $m=M_1$.  The moments $M_n$ provide a useful probe
of the dynamics.  From the governing equation (\ref{ck-eq}), the
second moment of the size distribution $M_2$ obeys the closed equation
$dM_2/dt=M_2(M_2-\alpha)$ and thus,
\begin{equation}
\label{M2-sol}
M_2(t)=\alpha\,\left[\left(\frac{\alpha}{\alpha_c}-1\right)\,
e^{\alpha t}+1\right]^{-1}.
\end{equation}
There is a critical freezing rate $\alpha_c=M_2(0)$.  For fast
freezing, $\alpha\geq \alpha_c$, the second moment is always finite
indicating that clusters remain finite at all times. In this case,
there is no gelation. For slow freezing, $\alpha<\alpha_c$, there is a
finite time singularity indicating that an infinite cluster, the gel,
emerges in a finite time \cite{finite}. This gelation time is
\begin{equation}
\label{tg} t_g=-\frac{1}{\alpha}\ln \left(1-\frac{\alpha}{\alpha_c}\right).
\end{equation}
The gelation point marks two phases. Prior to the gelation point, the
system contains only finite clusters that undergo cluster-cluster
aggregation (``coagulation phase''). Past the gelation point, the gel
grows via cluster-gel aggregation (``gelation phase''). We analyze
these two phases in order.

\smallskip
\noindent{\it Coagulation phase.}  Coagulation occurs for
$\alpha\geq\alpha_c$ at all times or for $\alpha<\alpha_c$ when
$t<t_g$. From (\ref{ck-eq}), the mass density of active clusters
satisfies $dm/dt=-\alpha\, m$, and thus, ordinary exponential decay
occurs,
\begin{equation}
\label{M1-sol}
m(t)=m(0)e^{-\alpha t}\,.
\end{equation}
For concreteness, we consider the monodisperse initial conditions
$c_k(0)=\delta_{k,1}$. In this case $M_n(0)=1$ and consequently,
$\alpha_c=1$. The cluster size distribution is obtained using the
standard transformation \cite{emh,sr,kb}, $c_k=e^{-\alpha
t}\,C_k$, and the modified time variable $\tau=\int_0^t dt'
e^{-\alpha t'}$ or explicitly
\begin{equation}
\label{tau}
\tau=\frac{1-e^{-\alpha t}}{\alpha}.
\end{equation}
This time variable increases monotonically with the physical time and
reaches $\tau\to 1/\alpha$ as $t\to\infty$. With these
transformations, Eq.~(\ref{ck-eq}) reduces to the no-freezing case
\hbox{$d C_k/d\tau=\frac{1}{2}\sum_{i+j=k} (iC_i)(jC_j) -
kC_k$.}  From the well-known solution of this equation \cite{ziff,bk},
the cluster-size distribution is
\begin{equation}
\label{ck-pre}
c_k(t)=\frac{k^{k-2}}{k!}\,\,\tau^{k-1}\,e^{-k\tau-\alpha t}.
\end{equation}
Generally, the size distribution decays exponentially at large
sizes and the typical cluster size is finite. The gelation time
(\ref{tg}) is simply $\tau_g=1$. Of course, no gelation occurs
when $\alpha>1$ because $\tau<1/\alpha<1$. The gelation point is
marked by an algebraic divergence of the size distribution
$c_k\sim (1-\alpha)k^{-5/2}$ for large $k$. We note that the mass
density decreases linearly with the modified time, $m=1-\alpha
\tau$, and that $m(\tau_g)=1-\alpha$.

\smallskip
\noindent{\it Gelation phase.} Past the gelation transition, a giant
cluster containing a finite fraction of the mass in the system
forms. In addition to cluster-cluster aggregation, cluster-gel
aggregation takes place with the giant cluster growing at the expense
of finite clusters. In parallel, all clusters including the gel may
freeze.

Formally, the size distribution (\ref{ck-pre}) generalizes to
\begin{equation}
\label{ck-post}
c_k(t)=\frac{k^{k-2}}{k!}\,\,\tau^{k-1}\,e^{-ku -\alpha t}
\end{equation}
with $u(t)=\int_0^t dt'\,m(t')$ yet to be determined. Statistical
properties of the size distribution are derived from the
generating function $c(z,t)=\sum_{k\geq 1} kc_k(t)e^{kz}$ that
equals
\begin{equation}
\label{czt}
c(z,t)=\tau^{-1}e^{-\alpha t} G(z+\ln\tau -u)
\end{equation}
where $G(z)=\sum_{k\geq 1}\frac{k^{k-1}}{k!}\,\,e^{kz}$ is the ``tree''
function \cite{jklp}.

During the gelation phase, active clusters consist of finite clusters,
the ``sol'', with mass $s$, and the gel with mass $g$ with the total
mass density, $m=s+g$, and these masses are coupled via the evolution
equations
\begin{subequations}
\begin{align}
\frac{dm}{dt}&= - \alpha s,\label{mt-eq}\\
\frac{ds}{dt}&=-\frac{s(m-s)}{1-s\tau e^{\alpha t}}-\alpha\,s.\label{st-eq}
\end{align}
\end{subequations}
The equation governing the total mass reflects that as long as the gel
is active, loss is due to freezing of finite clusters and the equation
for the sol mass follows from $ds/dt=-g\,M_2-\alpha M_1$, obtained by
summing (\ref{ck-eq}). To convert this equation into the explicit
Eq.~(\ref{st-eq}), the first two moments are $M_1=s$ and $M_2=
c_z(z=0)=s/(1-s\tau e^{\alpha t})$ where the identity
\hbox{$G'(z)=G/(1-G)$} \cite{wilf} was employed.

\begin{figure}[t]
 \includegraphics*[width=0.4\textwidth]{fig1.eps}
 \caption{The total mass, the sol mass, and the gel mass versus
 modified time $\tau$ for $\alpha=1/2$. These curves hold as long as
 the gel remains active.}
\label{sol}
\end{figure}

Equations (\ref{mt-eq}) and (\ref{st-eq}) are subject to the initial
conditions $m(t_g)=s(t_g)=1-\alpha$. Once the masses are found, the
formal solution (\ref{ck-post}) becomes explicit.  Results of
numerical integration of Eqs.~(\ref{mt-eq}) and (\ref{st-eq}) are
shown in Fig.~\ref{sol}.

Clearly, the longer the gel remains active, the larger it grows, but
its size can not exceed $g_{\rm max}=\lim_{t\to\infty}g(t)$ with the
obvious bound $g_{\rm max} <1-\alpha$ (Fig.~\ref{pg}). Furthermore, an
involved perturbation analysis \cite{bk-long} of Eqs.~(\ref{mt-eq})
and (\ref{st-eq}) yields
\begin{equation}
\label{g-pert}
g_{\rm max}\to
\begin{cases}
1-\frac{\pi^2}{6}\alpha&\alpha \downarrow 0,\\
C(1-\alpha)^2&\alpha \uparrow 1;
\end{cases}
\end{equation}
with $C=1.303892$. The quadratic dependence on the freezing rate
implies that the emerging gel is very small when $\alpha\uparrow
\alpha_c$. Thus, micro-gels, that may be practically indistinguishable
from large clusters, emerge.  We conclude that freezing can be used to
control the gel size, as gels of arbitrarily small size can be
produced using freezing rates just below criticality.

\begin{figure}[t]
 \includegraphics*[width=0.4\textwidth]{fig2.eps}
 \caption{The maximal gel size $g_{\rm max}$ and the probability of
 producing multiple gels $P_{\rm mult}$ versus the freezing
 rate $\alpha$.}
\label{pg}
\end{figure}

\smallskip
\noindent{\it Multiple giant clusters.} At any time during the
gelation phase, the gel itself may freeze. This freezing process is
random: the gel lifetime $T$ is an exponentially distributed random
variable, $P(T)=\alpha e^{-\alpha T}$. Until the gel freezes, the
system evolves deterministically, so the mass of the frozen gel is
$g(t_g+T)$. When the gel freezes, the total active mass $m(t)$ is
discontinuous: it exhibits a downward jump (Fig.~\ref{mass}).  Given
that the duration of the gelation phase is governed by a random
process, the mass of the frozen gel is also random. It fluctuates from
realization to realization, i.e., it is not a self-averaging quantity.

When the gel freezes, the system re-enters the coagulation phase
because all remaining clusters are finite. The initial conditions are
dictated by the duration of the preceding gelation phase and are
therefore also stochastic. Nevertheless, once the initial state is
set, the evolution in the coagulation phase is deterministic. The gel
is frozen so it no longer affects the evolution and only
cluster-cluster aggregation occurs. Let us assume that the gel freezes
at time $t_f$. Reseting time to zero, the first and the second moments
are simply given by Eqs.~(\ref{M1-sol}) and (\ref{M2-sol}),
respectively, with $M_n(0)$ replaced by $M_n(t_f)$. A second gelation
occurs if the freezing rate is sufficiently small,
$\alpha<M_2(t_f)$. Otherwise, the system remains in the coagulation
phase forever.

Because $M_2$ diverges at the gelation point, there is a time window
past the gelation time where the second moment exceeds the freezing
rate. If the gel freezes during this window, another percolation
transition is bound to occur. Therefore, there is no upper bound to
the number of frozen gels that may possibly form. Quantitatively, the
probability that a successive gelation occurs is $P_{\rm
mult}=\int_0^{t_*-t_g} dT\, P(T)$ with $M_2(t_*)=\alpha$
(Fig.~\ref{pg}); it can be obtained using perturbation analysis
\cite{bk-long}
\begin{equation}
\label{pmult}
P_{\rm mult}\to
\begin{cases}
\alpha \ln \frac{1}{e\alpha}&\alpha \downarrow 0,\\
0.450851&\alpha \uparrow 1.
\end{cases}
\end{equation}
In general, this probability is substantial and remarkably, it is
discontinuous at $\alpha_c$.

\begin{figure}[t]
\includegraphics*[width=0.4\textwidth]{fig3.eps}
\caption{The mass density $m$ versus time $\tau$. Shown are results of
a Monte Carlo simulation with system size $N=10^5$ and freezing rate
$\alpha=0.1$. The system alternates between the coagulation phase and
the gelation phase. In the former phase the mass decreases linearly
according to (\ref{M1-sol}) such that depletion occurs at time
$\tau=1/\alpha$. In the latter phase, the active mass decreases slower
than linear according to (\ref{mt-eq}) and (\ref{st-eq}).  The
gelation phase ends when the gel freezes.}
\label{mass}
\end{figure}

\smallskip
\noindent{\it Cyclic dynamics.} The general picture is now clear: the
process starts and ends in coagulation and throughout the evolution,
the system alternates between coagulation and gelation.  Once the
initial conditions are set, the behavior throughout the coagulation
phase and throughout the gelation phase are both deterministic.  Each
gelation phase ends with freezing of the active gel. Since the
duration of the gelation phase is random, the size of the giant
clusters, and the number of giant clusters are both random variables.
Generically, the system exhibits a series of percolation transitions,
each producing a frozen gel, so that overall, multiple gels are
produced.  The random freezing process governs the number of
percolation transitions as well as the size of the frozen gels.

Monte Carlo simulation of the stochastic process of aggregation
with freezing confirms this picture (Fig.~\ref{mass}). In the
simulations, we keep track the total aggregation rate
$R_a=N(M_1^2-M_2)/2$ and the total freezing rate $R_f=\alpha
NM_0$, where $N$ is the number of particles. Aggregation occurs
with probability $R_a/(R_a+R_f)$, and freezing occurs with the
complementary probability. A cluster is chosen for aggregation
with probability proportional to its size. Time is augmented by
$\Delta t=1/(R_a+R_f)$ after each aggregation or freezing event.

\smallskip
\noindent{\it Frozen clusters.}  The density of frozen clusters $F_k$
is found from $dF_k/dt=\alpha c_k$. By integrating Eq.~(\ref{ck-pre}),
we find the density $F_k(t_g)$ of clusters produced during the first
coagulation phase ($t_g\equiv\infty$ for $\alpha>1$)
\begin{equation}
\label{Pk-final}
F_k(t_g)=
\begin{cases}
\frac{\alpha}{k^2\cdot k!}\,\,\gamma(k,k)&\alpha\leq 1,\\
\frac{\alpha}{k^2\cdot k!}\,\,\gamma(k,k/\alpha)&\alpha\geq 1;
\end{cases}
\end{equation}
where $\gamma(n,x)=\int_0^x dy\,y^{n-1}\,e^{-y}$ is the incomplete
gamma function. When $\alpha\geq 1$, this quantity equals the final
density of frozen clusters, $F_k(\infty)=F_k(t_g)$, and at large
sizes, the behavior is as follows
\begin{equation}
F_k(\infty)\simeq
\begin{cases}
\frac{1}{2}\cdot k^{-3}                          &\alpha=1,\\
A(\alpha)\,k^{-7/2}\exp\left[-B(\alpha) k\right] &\alpha>1;
\end{cases}
\end{equation}
where $A=(2\pi)^{-1/2}\alpha^2/(\alpha-1)$ and
$B=\alpha^{-1}+\ln\alpha -1$.

\begin{figure}[t]
 \includegraphics*[width=0.4\textwidth]{fig4.eps}
 \caption{The size distribution of the frozen clusters for
   $\alpha=1/2$. The simulation results represent an average over
   $10^2$ independent realizations in a system of size $N=10^6$.}
\label{frozen}
\end{figure}

For the case $\alpha<1$, numerical simulations provide convincing
evidence that the tail behavior
\begin{equation}
F_k(\infty)\sim k^{-3}, \qquad \alpha<1
\end{equation}
is universal (Fig.~\ref{frozen}).  Indeed, Eq.~(\ref{Pk-final}) shows
that the first coagulation phase contributes $F_k(t_g)\simeq
\frac{\alpha}{2}k^{-3}$ to $F_k(\infty)$ and similarly, each
successive coagulation phase should contribute the very same $k^{-3}$
to the tail.  Because large clusters quickly merge into the gel, the
gelation phase yields a negligible contribution, and we conclude that
freezing leads to an additional non-trivial critical exponent
$\gamma=3$.

In summary, freezing results in a rich and interesting phenomenology
that includes multiple percolation transitions and multiple gels,
microgels that may be indistinguishable from finite clusters, cyclic
evolution with the system constantly alternating between coagulation
and gelation, lack of self-averaging, and an additional nontrivial
critical exponent.

Freezing also provides a practical mechanism for controlling
gelation. It may be used to engineer micro-gels of desired size by
implementing variable freezing rates. Slow freezing followed by
rapid freezing can be used to produce gels of prescribed size,
while near-critical freezing produces micro-gels of arbitrarily
small size.

There are a number of interesting potential generalizations of the
present work. One may consider situations with different freezing
mechanisms \cite{kb}, particularly for finite and infinite
clusters. We also solved the Stockmayer version \cite{stock,zs}
where $m=M_1$ always; in this case, there is no breakdown of
self-averaging and the final density of frozen clusters mimics the
critical behavior of active clusters, $F_k\sim k^{-5/2}$.  Another
natural generalization is percolation in finite dimensions where
the final distribution of frozen clusters should be an interesting
probe of the dynamics.

\acknowledgments
We acknowledge US DOE grant W-7405-ENG-36 (EBN) for support of this
work.

\begin{thebibliography}{99}

\bibitem{flory}  P.~J.~Flory,
                 J. Amer. Chem. Soc. {\bf 63}, 3083--3090 (1941).

\bibitem{stock}  W.~H.~Stockmayer, J. Chem. Phys. {\bf 11}, 45--55 (1943).

\bibitem{ah}     D.~Stauffer and A.~Aharony,
                 {\em Introduction to Percolation Theory}
                 (Taylor \& Francis,  London, 1992).

%\bibitem{bh1} A.~Bunde and S.~Havlin (editors), {\em Percolation and
%                 Disordered Systems: Theory and Applications}, Physica
%                 A {\bf 366}, (1998).

\bibitem{sah}    M.~Sahimi,
                 {\em Flow and Transport in Porous Media and Fractured Rock}
                 (VCH, Boston, 1995).

\bibitem{floryB} P.~J.~Flory,
                 {\em Principles of Polymer Chemistry}
                 (Cornell University Press, Ithaca, 1953).

\bibitem{pg}     P.~Grassberger,
                 Math. Biosci. {\bf 63}, 157 (1983).

\bibitem{BH}     S.~R.~Broadbent and J.~M.~Hammersley,
                 Proc. Cambridge Phil. Soc. {\bf 53}, 629 (1957).

\bibitem{drake}  R. L. Drake, in {\em Topics in Current Aerosol Researches},
                 eds. G. M. Hidy and J. R. Brock (Pergamon, New York,
                 1972), pp. 201.

\bibitem{sein}   J. H. Seinfeld and S. N. Pandis,
                 {\em Atmospheric Chemistry and Physics}
                 (John Wiley \& Sons,  New York, 1998).

\bibitem{jlr}    S.~Janson, T.~\L uczak, and A.~Rucinski,
                 {\em Random Graphs} (John Wiley \& Sons, New York, 2000).

\bibitem{gg}     G.~Grimmett,
                 {\em Percolation},
         (Springer, Berlin, 1999).

\bibitem{hart}   S.~Cocco, L.~Ein-Dor, and R.~Monasson, in
                 {\em New Optimization Algorithms in Physics},
                 eds. A.~K.~Hartmann and H.~Rieger
                 (Wiley--VCH, Weinheim, 2004).

\bibitem{bbckw}  B.~Bollob\'as, C.~Borgs, J.~T.~Chayes, J.~H.~Kim,
                 and D.~B.~Wilson, Rand.\ Struct.\ Alg. {\bf 18}, 201 (2001).

\bibitem{win}    E.~Winfree, F.~R.~Liu,  L.~A.~Wenzler, and N.~C.~Seeman,
                 Nature {\bf 394}, 539--544 (1998).

\bibitem{ms}     M.~J.~Murray and M.~J.~Snowden,
                 Adv. Coll. Int. Sci. {\bf 54}, 73 (1995).

\bibitem{sv}     B.~R.~Saunders and B.~Vincent,
                 Adv. Coll. Int. Sci. {\bf 80}, 1 (1999).

\bibitem{fran}   F.~Leyvraz, Phys. Rep. {\bf 383}, 95 (2003).

\bibitem{aal}    A.~A.~Lushnikov,
                 Phys. Rev. Lett. {\bf 93}, 198302 (2004).


\bibitem{finite} In a finite system, the gels contains a finite
                 fraction of the system mass.

\bibitem{emh}    E.~M.~Hendricks, J. Phys. A {\bf 22}, 99 (1984).

\bibitem{sr}     P.~Singh and G.~J.~Rodgers, J. Phys. A {\bf 29}, 437 (1996).

\bibitem{kb}     P.~L.~Krapivsky and E.~Ben-Naim,
                 J. Phys. A {\bf 33}, 5465 (2000).

\bibitem{ziff}   R. M. Ziff, E. M. Hendriks, and M. H. Ernst,
                 Phys. Rev. Lett. {\bf 49}, 593 (1982).

\bibitem{bk}     E.~Ben-Naim and P.~L.~Krapivsky, J. Phys. A {\bf 37}, L189
                 (2004); cond-mat/0408620.

\bibitem{jklp}   S.~Janson, D.~E.~Knuth, T.~\L uczak, and B.~Pittel,
                 Rand.\ Struct.\ Alg. {\bf 3}, 233 (1993).

\bibitem{wilf}   H.~S.~Wilf, {\em generatingfunctionology} (Academic
                 Press, Boston, 1990).

\bibitem{bk-long}     E.~Ben-Naim and P.~L.~Krapivsky, preprint.

\bibitem{zs}    R.~M.~Ziff and G.~Stell, J. Chem. Phys. {\bf 73}, 3492 (1980).

\end{thebibliography}

\end{document}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
