Kinetic theory and numerical simulations of two-species coagulation

In this work we study the stochastic process of two-species coagulation. This process consists in the aggregation dynamics taking place in a ring. Particles and clusters of particles are set in this ring and they can move either clockwise or counterclockwise. They have a probability to aggregate forming larger clusters when they collide with another particle or cluster. We study the stochastic process both analytically and numerically. Analytically, we derive a kinetic theory which approximately describes the process dynamics. One of our strongest assumptions in this respect is the so called well-stirred limit, that allows neglecting the appearance of spatial coordinates in the theory, so this becomes effectively reduced to a zeroth dimensional model. We determine the long time behavior of such a model, making emphasis in one special case in which it displays self-similar solutions. In particular these calculations answer the question of how the system gets ordered, with all particles and clusters moving in the same direction, in the long time. We compare our analytical results with direct numerical simulations of the stochastic process and both corroborate its predictions and check its limitations. In particular, we numerically confirm the ordering dynamics predicted by the kinetic theory and explore properties of the realizations of the stochastic process which are not accessible to our theoretical approach.

1. Introduction. The theoretical study of coagulation and their kinetic description is of broad interest because of its vast applicability in diverse topics such as aerosols [28], polymerization [35,30], Ostwald ripening [21,8], galaxies and stars clustering [29], and population biology [25] among many others. In this work we propose a generalization of the stochastic process of coagulation. Herein we will consider the coagulation process among two different species: the aggregation will take place only when one element of one of the species interacts with an element of the other species. In particular, we will place particles and clusters of particles in a ring, where they will move with constant speed, either clockwise or counterclockwise. When two clusters (or two particles or one particle and one cluster) meet they have the chance to aggregate and form a cluster containing all particles involved in the collision. The direction of motion of the newborn cluster is chosen following certain probabilistic rules. We are interested in the properties of the realizations of such a stochastic process, and in particular in their long time behavior. Our main theoretical technique will be the use of kinetic equations, an approach we have outlined in [13]. This is of course just a possible extension of the theory of coagulation. We have designed it getting inspiration from self-organizing systems and in particular from collective organism behavior. Let us note that this is a field that has been studied using a broad range of different theoretical techniques [31,11,3,5,10]. Another field which has inspired ourselves is the study of the dynamics of opinion formation and spreading [12,32,6], which is represented for instance by the classical voter model [7,17,22]. As a final influence, we mention that clustering has been previously studied in population dynamics models [16] including swarming systems [18], and coagulation equations have been used in both swarming [26] and opinion formation models [27]. Despite of its simplicity, the two-species coagulation model could be related to some of these systems.
A particular system that has influenced the current developments is the collective motion of locusts. The experiment performed in [4] revealed that locusts marching on a (quasi one dimensional) ring presented a coherent collective motion for high densities; low densities were characterized by a random behavior of the individuals and intermediate densities showed coherent displacements alternating with sudden changes of direction. The models that have been introduced to describe this experiment assume that the organisms behave like interacting particles [4,34,14]. Related interacting particle models have been used to describe the collective behavior of many different organisms and analyzing the mathematical properties of such models has been a very active research area [11,5,10]. The two-species coagulation model could be thought of as a particular limit of some of these models or as a simplified version of them which still retains some desirable features.
The goal of the current work is not to describe the detailed behavior of any specific system. Instead, we will explore the mathematical properties of a stochastic process which has been designed by borrowing inspiration from different self-organizing systems. Therefore the focus will be on mathematical tractability. We will determine under which conditions consensus is reached and what form it adopts. In order to do this we introduce a kinetic theory that approximately describes the stochastic process. One of our main assumptions is to consider that the system is well-stirred, so the explicit dependence on the spatial coordinates can be dropped. This approximation is very common in other fields like in the study of reaction-diffusion systems or Boltzmann equations. The stochastic process, which will be just qualitatively described in this introduction, will be precisely introduced in section 4. The kinetic theory we will consider in order to theoretically describe this process, equation (1) below, will be heuristically introduced rather than derived as a proper limit of the stochastic process. Our present approach is an aposterioristic one: we will compare the theoretical predictions of our kinetic equation with direct numerical simulations of the stochastic process in order to check the validity of our theory.
In the following section we describe the kinetic theory that approximately describes the stochastic process under study in the well-stirred limit. We build our progress here based on our previous developments in [13]. The basic analysis of the resulting kinetic equations is carried out in this section too, and the results obtained are interpreted in the context of the underlying stochastic process. Part of the phenomenology of the ordering dynamics is already unveiled at this point. We devote the subsequent sections to a more detailed mathematical analysis of this kinetic theory. In particular, we study the advent of self-similarity of the solutions to the kinetic equations. We interpret this self-similar behavior as the formation of a giant cluster composed by all the particles in the system.
We postpone to the last section of the paper the direct numerical simulations of the stochastic process. In this last section we both check the predictions of our kinetic theory and explore the stochastic process beyond the kinetic level. Kinetic approximations neglect many sources of fluctuations and thus numerical simulations are required in order to describe many properties of the individual realizations of the stochastic process which are not reflected at the kinetic level. It is also important to remark here that the direct numerical simulations of the stochastic process do not assume the well-stirred limit. So in these simulations the spatial distribution of particles and clusters is explicitly taken into account. Let us mention again that the only reason why the well-stirred limit is considered in the theoretical part of this work is in order to allow for mathematical tractability, so there is no reason in imposing such a constraint in the numerical part. Furthermore, the numerical results justify the theoretical assumption under appropriate hypothesis on the model parameters.
Our analysis will be limited to the one-dimensional spatial situation with periodic boundary conditions (the dynamics is taking place in a circumference). Our present approach will be based on coagulation equations. This kinetic description is in general invalid for one-dimensional systems because it is known that spatial correlations do propagate in this dimensionality. So we assume a collision takes place when two clusters meet with a very small probability. The probability should be so small that all the particles travel the whole system several times before one collision happens on average. This way the system becomes well-stirred, and so we can neglect spatial correlations and treat the system as if it were zero dimensional, what allows mathematical tractability. The rate at which the collisions take place can be trivially absorbed by means of a rescaling of time, so we do not include it explicitly in the formulation of the kinetic equations. It will however be explicitly taken into account when we perform direct numerical simulations of the stochastic process in section 4.
We assume that the particles move in clusters of individuals and f + ( , t) (f − ( , t)) represents the number of clusters of size moving clockwise (counterclockwise) at time t. We note these functions should in principle only take integer values, however the kinetic approximation implicitly averages over many realizations, so they will be allowed to take any non-negative real value. Clusters are modified when they collide with other clusters, in such a way that the probability where Ψ(k, j; m, ) is the collision kernel: It states the probability with which a collision among clusters with k and j particles occurs and yields clusters with m and particles. One of our basic assumptions is its symmetry under reflections Ψ(k, j; m, ) = Ψ(j, k; , m). It is worth remarking that equation (1) is not the exact description of the coagulation process. It is in fact an approximation, but we will not derive it from the stochatic process. Our approach will be phenomenological: we postulate the validity of such a kinetic description and check it a posteriori by means of numerical simulations. Let us however briefly comment on that it is reasonable postulating this kinetic theory. The right hand side in the upper line of (1) is the gain term: it represents the possibility of formation of clusters of size from the collision of a cluster of size k and another one of size j; the sum takes into account all possible sizes j and k. The second line of this equation is the loss term: it takes into account the chance of disappearance of a cluster of size by its coagulation with a cluster of any other size.
Our analytical progress on this equation will be built by means of the introduction of the generating functions If we take the derivative these functions with respect to time we get As already mentioned, the focus of this paper is on mathematical tractability, and so we will only consider kernels such that this expression becomes a closed system of differential equations for the generating functions. Our two choices are described in the following.
1.1. Random kernel. In the present paper we will concentrate on two integrable kernels. Our first choice is Ψ(k, j, , m) = 0, otherwise.
Note that instead of writing down the Kronecker deltas explicitly in Eq. (1) one could include them into the definition of the collision kernel. In this case, instead of Eqs. (4) and (5), we would have Ψ( , m, k, j) = δ +m,k+j (δ k,0 + δ j,0 )/2.
This kernel expresses the following process: when two clusters collide a single cluster merges and it contains all particles of both. The newborn cluster selects its direction of motion, that could be either clockwise or counterclockwise, with equal probability. We will refer to (4) as the "random kernel". Substituting this kernel in (3) we find for the generating functions the following differential system The number of clusters traveling in either direction is given by Lemma 1.1. System (6) and (7) admits the conserved quantity Proof. The conservation law follows from substituting z = 1 in (6) and (7) and substracting both equations. Denoting N + (t) − N − (t) = C 0 and substituting this quantity back into the differential system we get This is a Bernoulli equation that can be integrated to yield The conclusion follows immediately. This is a linear ordinary differential equation that can be integrated to yield The second conclusion follows immediately after taking the long time limit in this formula.
These results indicate the system always becomes ordered in the long time. In fact, the number of remaining clusters equals the (absolute value) of the difference of the number of clusters traveling in either direction initially. So the direction of motion of the final ordered state is prescribed by the initial condition. This result should be interpreted probabilistically. The quantities N ± (t) and M ± (t) do not represent the dynamics of single realizations. They are instead the result of averaging these quantities over many realizations of the stochastic process. So they represent what happens on average. If the quantity C 0 is macroscopic, i. e. it is of the same order of magnitude of the initial total number of clusters, then the result will be observed in most of the realizations, up to some errors.
The case of the symmetric initial condition N + (0) = N − (0) is specially interesting. If this condition holds initially, then it holds for all times: N + (t) = N − (t) ∀t > 0. Furthermore, if we refer to these quantities just as N (t) (precisely because they are equal), we have These results are immediate consequences of the proof of Lemma 1.1. Equation (14) means that clusters in both directions disappear progressively in the long time. Note that equilibrium is approached exponentially fast when C 0 = 0 but the approach becomes algebraic when C 0 = 0. This case is special because it is the only one, for the random kernel, for which the kinetic theory does not predict order. However, if we go beyond the kinetic level and think about the underlying stochastic process, order should be achieved. This is so because states in which there exist clusters moving in both directions are active: collitions are possible and therefore the system could evolve through collisions towards any other possible state. On the other hand, ordered states are absorbing: once the system gets into one of this, collisions are no longer possible, and therefore there is no possible escape from them. Note the situation is reminiscent of that of systems undergoing absorbing state phase transitions, and in particular of the case of two symmetric absorbing barriers [15,33]. We devote much of the remainder of this paper to the study of this case. In section 2 we study the kinetic theory corresponding to this case. We show that the long time asymptotics of the distribution function adopts a self-similar form in this case.
In particular we have that where Φ is the self-similar profile, in a sense that is made totally precise in theorems 2.6 and 2.8, which are the main results of that section. This explains the fact that the number of clusters approaches zero as time evolves: clusters of smaller sizes progressively disappear as bigger and bigger clusters are formed. In section 4 we study direct numerical simulations of the aggregation process precisely in this case. We confirm the theoretical prediction finding that it is quite probable that order is achieved by means of the formation of one giant cluster containing all particles in the system. This fact reflects at the microscopic level the formation of the self-similar form for long times in the solution to the kinetic equation. So order is finally achieved in this case too, but it has a different nature than in the previous situation. Order has the form of a single giant cluster, not of several clusters; the direction of motion of the giant cluster is not prescribed at the initial time, but it is chosen at random with the same probability among the two different possibilities; and finally order is achieved after a longer transient in this case.
which describe collisions among clusters with k and j particles. After the collision a single cluster composed of all involved particles merges. The resulting cluster has a probability k/(k + j) of traveling in the same direction the cluster with k particles was traveling, and a probability j/(k + j) of traveling in the direction the cluster with j particles was traveling. We will refer to this as the "majority kernel". Our approach will be the same as in the other case and starts with equation (3). For the generating functions we find where the number of clusters N ± (t) = F ± (1, t). Denoting M ± (t) = ∂ z F ± (1, t) the number of particles, we find so the initial number of particles traveling in either direction does not change over time. This result apparently means the system could never become ordered. However, the actual meaning is that the number of particles moving in either direction, averaged over many realizations of the stochastic process, is constant. As we have already said the kinetic description implicitly assumes this average. In other words, in the dynamics dictated by this kernel the individual collisions do not conserve momentum, but momentum is conserved on average. In any case, individual realizations should become ordered, because as in the former case, order is an absorbing state. In particular, the situation is analogous to the previous one: we have a system with two symmetric absorbing barriers again. Result (18) could indicate the advent of dynamic self-similarity, just like in the former case, in the long time limit. Although for the majority kernel we did not develop the theory to the same extent as in the former case, we have built some progress which is detailed in section 3. In this section we show (computational) evidence of the existence of self-similar asymptotic behavior, although we do not clarify the conditions under which it develops.
2. Self-similar asymptotics for the random kernel. The goal of the present section is to carry out the rigorous analysis of the asymptotic behavior of our kinetic theory in the case that the random kernel is considered. We will divide the analysis in several sections. In section 2.1 we set the preliminary results for the analysis of a simplified situation in which we assume that all moments of both initial conditions 2.1. Preliminaries. As we have already mentioned, we will now concentrate on the case of a symmetric initial condition (N + (0) = N − (0)). Our goal is to derive its self-similar asymptotic behavior. Before moving to the general case, we start with a simplified situation: we assume that all moments of both initial condition are identical. This translates in considering equations (6) and (7) with initial conditions The probabilistic approach employed in the previous section implies the following structure of the initial condition: a n z n , a n ≥ 0.
We may assume N (0) = ∞ n=0 a n = 1 without loss of generality; in particular this implies F 0 (z) is holomorphic on the unit disc in C as a function of z. From now on we will consider z as complex variable taking its values on the unit disc of C and a real 0 ≤ t < ∞.
We start proving the following Lemma 2.1. Consider the system of differential equations (6)-(7) subject to the initial condition F + (z, 0) = F − (z, 0). Then we have F + (z, t) = F − (z, t) during the lapse of existence of the solutions.
Proof. Subtracting both equations one finds Using This equation can be integrated to yield so we conclude.
Corollary 1. The function F obeys the following differential equation: .

TWO-SPECIES COAGULATION 9
Proof. Choosing z = 1 we obtain and so Substituting this value of N (t) in the equation for F we reduce this second equation to an ODE of Bernoulli type. The solution of the initial value problem for F can be therefore found by means of the change of variables We find And the desired solution for F (z, t) follows from F (z, t) = 1/H(z, t) whenever H(z, t) = 0. Since |F 0 (z)| ≤ 1 for |z| ≤ 1, as can be deduced from the Taylor series in (19) by means of the triangle inequality, H(z, t) = 0 for every z in the unit disc of C and 0 ≤ t < ∞, so we conclude.
Corollary 2. The result in Lemma 2.1 holds for all time t > 0.
We have already seen in the proof of Lemma 2.3 that |F 0 (z)| ≤ 1 for |z| ≤ 1, as can be inferred from the triangle inequality applied to the Taylor series (19). We also have the strict inequality if we consider the strict triangle inequality instead.
By assumption we know F 0 (1) = 1; on the other hand we can have F 0 (z) = 1 for |z| = 1 at some points z = 1. In particular we have the following result: Proposition 1. Let F 0 (z) be as explained above. Let d = gcd(n|{a n = 0}), n ∈ N + . Then F 0 (z) = 1 if and only if z d = 1.
Proof. Suppose that a n = 0 only for n = md, m = 1, 2, . . . , d > 1. Then It then follows that F 0 (z) = 1 for z d = 1. Now we prove the reciprocal. If a n = 0 except for a subsequence then where a md contains the non-vanishing subsequence and d is the greatest common divisor of all the positive coefficients a n : d = gcd(n|{a n = 0}), n = 1, 2, · · · . In order to have F 0 (z) = 1 we need otherwise the strict triangular inequality applied to the Taylor series implies |F 0 (z)| < 1. So we have, by the definition of greatest common divisor, the existence of a pair of relative primes m 1 , m 2 > 0 such that z m1d = 1 and z m2d = 1. This automatically implies (z d ) xm1+ym2 = 1, for any two integers x and y. A direct consequence of Bézout lemma is the existence of pairs of integers (x, y) satisfying the diophantine equation xm 1 + ym 2 = 1; this leads to the desired conclusion z d = 1.

Remark 1.
The previous result immediately implies that the expression F 0 (z) − 1 has a family of zeros at the d−roots of unity.
This property in turn assures the existence of a "characteristic wavelength" in the following sense: If the initial condition is composed only by clusters having a number of particles which is multiple of some integer then only clusters with a number of particles multiple of the same integer will be generated.
The computations in the following section will assume that z = 1 is the only root of F 0 (z) = 1. If there were more roots the computations would be analogous.
and for a suitable constant C (the symbol C will be used to denote a generic constant whose value may change from line to line). On In the following we will obtain the asymptotic behavior of the solutions in this case. Note that an immediate corollary of Lemma 2.1 is that f + ( , t) = f − ( , t) ∀ t ≥ 0; therefore we will use the following Definition 2.4. We denote f ( , t) ≡ f + ( , t) = f − ( , t) whenever this equality holds.

2.2.
Formal asymptotics for the symmetric case. In this section we will carry out a formal analysis of the asymptotic behavior of f ( , t); these results will be rigourously justified in the following sections. Using the Cauchy formulas for the Taylor coefficients we find where the line integration is carried out in the counterclockwise direction. Then we have It is easy to see that the integral is well defined for all finite t.
We start our analysis with the following assumption on the initial data for an arbitrary positive constant C. This assumption is purely technical: we will use it in the first instance to simplify the analysis, but we will substitute it for a more general one further below. In this case we have F 0 (z) is holomorphic in |z| < 1 + δ, δ > 0 small enough, and 0 It is interesting to examine the integrand pole situated next to z = 1. Such a pole is situated at the root of If t = ∞ we know that the root is simple since 0 < F 0 (1) and z ∞ = 1. The implicit function theorem assures the existence of one root for t large enough and one has the approximation 2 Deforming contours and applying the residue theorem we obtain where δ is chosen small enough in order to avoid additional singularities in the integrand. By employing L'Hôpital method This gives the self-similar behavior of the solution when t → ∞, since In particular, in the region of order t one obtains So we obtain the self-similar structure 1 for long t, large , and both sharing the same magnitude. We note that 2F 0 (1) = M + (0) + M − (0) =: 0 , the total number of particles in the system (note this is a conserved quantity as the collisions we consider are mass conserving).
On the other hand, it is interesting to observe that in the region where is of order one the resulting integral term in the f ( , t) formula yields an order 1 t 2 term depending on the values of F 0 in regions where z is not in the neighborhood of one. The contribution to the mass from such a term approaches zero when → ∞, since this contribution is relevant only if is of order one. With the hypothesis of analyticity made we have that such a term is bounded by e −C t 2 , except for a multiplicative constant, when → ∞. Anyway this term induces a sort of "boundary layer" in the sense that the region where is of order one cannot be described using the self-similar function. Rigorous convergence results will be shown further below, where we will precisely state how the solution of the coagulation equation converges to the self-similar function. Also, both contributions to the dynamics, self-similar function and boundary layer, will be clearly identified in our numerical simulations in section 4.
2.3. The asymmetric case. Now we move to calculate the solution of the random kernel model when N we have, from equations (6), (7) and (8) The second equation may be integrated to give . (7) we get This equation of Bernoulli type can be explicitly solved making the substitution Equation (25) can be integrated to give In particular, the solution F − to (24) has the simple expression Our goal is to prove rigorous estimates for the asymptotic temporal behavior of the distribution functions f ± ( , t). The analysis will be split in the following sections. We start by proving the necessary estimates in the complex plane for the corresponding generating functions. We proceed analogously to what we did in the previous case of a totally symmetric initial condition. This is, we start looking for poles in the integrand of the integral expression for f − ( , t).
The main objective of this section is to prove that the exact solution for F − (z, t) is given by First we reformulate the problem in purely analytic terms. We can rescale the factor N 0 It is convenient to rewrite the solution of F − to check the condition that avoids singularities in the present case. We have Therefore the singularities arise at the zeroes of This function can have poles at the zeroes of f 0 (z) . At such points F − (z, t) = 0. Notice that the zeroes of g 0 are not problematic due to the cancellations of the zeroes of at the numerator and denominator. . This expression has no poles neither for t ∈ [0, ∞) and |z| ≤ 1 nor for t = ∞ and |z| < 1.
Proof. Expression (26) can be bounded using the fact |F − (z, t)| ≤ F − (1, t) for |z| ≤ 1; this inequality is obtained using the triangular inequality on the power series which defines this quantity, exactly as we did it in the previous section.
We have also shown the coagulation equation implies that the number of particles N − (t) is bounded and decreases as 1 t for large t. Then This shows no poles are present in this expression for finite t and |z| ≤ 1. Now we move to the case t = ∞. In this case our problem reduces to prove that the following function does not have any zero in the disk |z| < 1. We proceed by contradiction. We assume there is a zero for this expression at z 0 such that |z 0 | < 1. For long t we have the approximation Due to Rouche's Theorem, expressions (28) and (29) should have the same number of zeros inside the disk |z| = 1 for a sufficiently large t, provided F − (z, 0) has no zeros on |z| = 1. If this number is nonzero this implies that F − (z, t) blows up for a finite value of t at the interior of the disk |z| = 1 but this contradicts the estimate (27), so we conclude in this case. If F − (z 1 , 0) = 0 for some z 1 on |z| = 1 then we know this zero is isolated because the function is holomorphic. We may know apply Rouche's Theorem on the contour |z| = 1 − ; we can make arbitrarily small so we have |z | < 1 − for any possible pole z in the open disc |z| < 1, and by enlarging the time t we may keep under control the O(t −1 ) term in (29). This leads to the desired contradiction in the general case. Remark 2. Notice that there is at least a zero of (28) at the point z = 1.
Remark 3. In order to obtain convenient representation formulas for the selfsimilar asymptotics of the coagulation equation under study it is necessary to locate the remaining zeros of (28) in {|z| = 1}\{1}. As in the case of the totally symmetric initial conditions it is possible to have more than one zero at the unit disk |z| = 1.
For this to happen, by Proposition 1, we have the following necessary and sufficient condition In this case there are several zeros z d = 1. Physically the resulting solution will have gaps of length d. In summary, there is a complete analogy with what happened in the case of the totally symmetric initial condition, see the discussion in section 2.1.
In order to find the roots of the equation we reformulate the problem in terms of the function after assuming that g 0 (z) = 0 for every z that solves the equation. We have also used the fact that both f 0 and g 0 are bounded and the equality holds for f 0 = 0 only if g 0 = 0 and vice versa. We note the points z such that g 0 (z) = 0 are not important for our present purposes since at these points the equality holds only if f 0 (z) = 1, what in turn implies h 0 (z) = 1. From now on we drop the subindex "0" in order to simplify the notation of these functions.
Lemma 2.5. Let f (z) = ∞ n=1 a n z n and h(z) = ∞ n=1 b n z n be two holomorphic functions defined for z ∈ C on the closed disc |z| ≤ 1. Assume a n , b n ∈ R + ∪ {0} and f (1) = h(1) = 1. Then a complex number z 0 such that |z 0 | ≤ 1 is a solution to the equation Proof. The direct implication is obvious. So we will concentrate on proving the inverse implication in the following.
We can cast the equation into the following form using the fact that both f and h are bounded and the original equality holds for f = 0 only if h = 0 and vice versa.
Note that both f and h fulfill the inequality |f | ≤ 1, |h| ≤ 1, for |z| ≤ 1, so in particular these functions map the closed disc |z| ≤ 1 onto itself. The complex function w → we −w , where {w ∈ C | |w| = 1} , has winding number 1. This can be seen by noting that w = e iθ , θ ∈ [0, 2π), and so we have and so the phase θ − sin(θ) clearly reveals that this function winding number is 1. By invoking the Argument Principle we may conclude that the equation has at most one solution for w in the disc |w| < 1 for some fixed w 0 . Now let us focus on the |w| = 1 case. In this case we see that the phase θ − sin(θ) ∈ [0, 2π) and it is strictly increasing when θ ∈ [0, 2π). This tells us that the mapping w → we −w is univalued in |w| = 1 as well, yielding us the desired conclusion.
Remark 4. The same conclusion could be derived using the properties of a well known special function. Consider now the equation that defines the tree function T = T (z) [19]. Alternatively we have where W = W(z) is the Lambert Omega function [9]. The tree function is multivalued (as the Lambert function is), but we know, by using the properties of the Lambert function proved in [9], that for |T | ≤ 1, T (z) maps bijectively onto the z−plane, what guarantees the desired result. 2.4. Formal asymptotics for the asymmetric case. We now proceed to calculate the self-similar asymptotic form of the solution to the problem with identical initial conditions F + (z, 0) = F − (z, 0). We will assume /d ∈ N because otherwise f ( , t) ≡ 0. We already found the following explicit formula for the solution f ( , t) = 1 2πi where F 0 (z) is the initial condition of the generating function. We know that this function can be written as F 0 (z) = Q(z d ), where d is some positive integer and Q(z d ) is a function such that Q(1) = 1 and Q(w) = 1 for w ∈ C, |w| ≤ 1, and w = 1. In this case we have and changing variables z d = ζ we obtain f ( , t) = 1 2πid There is a pole in the integrand next to ζ = 1: + higher order terms, so we may write f ( , t) = 1 2πid .
The limit can be calculated with L'Hôpital rule and so we have In the limit t max{2/Q (1), 2}, d, and /t finite we find where we have used We thus see that the self-similar long-time solution 1 t 2 Φ( t ) is uniform in d; the number d is just a measure of the size of the clusters when the self-similar state is reached, what happens when d. Apparently, the time it takes to reach the self-similar state should depend on the gap index d: we have shown this time fulfills t max{4d/ 0 , 2}, and so, the bigger the gaps are, the longer the transient to the self-similar state would be. Note however that this condition trivially becomes t 2 because the gap index could be at most twice the number of particles. In consequence the length of the transient towards self-similarity is independent of d.
Now we move to calculate the self-similar asymptotics for the case of not identically distributed initial conditions. We assume N (0) = 1, and recall the solution for the generating function: Performing the change of variables F − (z, 0) = P (z d ) and G 0 (z) = Q(z d ) just like in the previous case we find where δ is small enough and is the location of the pole. Herein we have implicitly assumed that /d ∈ N; otherwise f − ( , t) ≡ 0 due to the orthogonality of the basis of plane waves. The limit can be calculated employing L'Hôpital rule .
For long times we find Putting all together and in the limit d, t max{ 2 P (1)+Q (1)/2 , 2} = 2 (while keeping /t constant) we obtain the self-similar structure This result is completely analogous to that of the previous case. Taking where the discrete Dirac comb (built as a sum of Kronecker deltas) explicitly signals the values of which correspond to a nonzero f ± ( , t). We emphasize this result is identical to the one obtained in the previous case, so this shows the asymptotic selfsimilar state is independent of the initial distributions of clusters f ± ( , 0) except for their zeroth mode (i. e. the number of clusters). Also as in the previous case, the cluster sizes in the self-similar regime fulfill d, and the transient time fulfills t 2. Note that the total number of particles is as always conserved.

2.5.
Convergence. In this section we prove rigorous convergence results to the self-similar profile formally calculated in the previous section An analogous convergence result for the case examined in section 2.2 is an immediate corollary of the present result, so we will not explicitely consider it here.
Theorem 2.6. There exists a suitable constant C such that for all 1 ≤ p < ∞ in the long time limit t → ∞, provided all moments of the initial condition are bounded.
Proof. We found in the previous section This equality holds as a consequence of the finiteness we imposed on the different moments of the initial condition, i. e.
Note this is equivalent to requiring the following bound for the initial condition for any C, δ > 0, see [20]. This assures that the transformed functions F ± (z, t) are analytic in the open disc |z| < 1 + δ, and so the contour deformation of the above complex integrals makes sense (note that this requirement is sharp). If there are no more poles in the integrand for |z| < 1 + δ apart from the aforementioned roots of unity, z d = 1, we have is the self-similar form. This result implies pointwise convergence of the solution f ± ( , t) to the self-similar form t −2 Φ( /t) uniformly in . Furthermore, the decay in is strong enough to have the estimate stated in the theorem.
Remark 5. Note that the integral can be bounded by a constant C if the integrand contains no poles. This is not necessarily the case for |ζ| = 1 + δ/2. But taking into account that Q(ζ) and P (ζ) are holomorphic in |z| < 1 + δ implies that the denominator (which is not zero on |ζ| ≤ 1, ζ = 1) is meromorphic, and thus it has a finite number of zeros in the open disk |ζ| < 1 + δ. A direct consequence of this fact is that we can choose some 0 < δ 1 < δ small enough such that the integral along the line |ζ| = 1 + δ 1 has no poles and it encloses an area in which ζ = 1 is the only pole.
2.6. The Sharp Case. In this section we are interested in obtaining asymptotic convergence results with the weakest possible hypothesis on f ± ( , 0). In particular, we will assume the finiteness of the first moment of the initial condition only: Note this is equivalent to have initial conditions of the form except for slow variation functions, for 1 < α ≤ 2. No further assumptions will be imposed in this section. As in the previous section we will only rigorously prove the formal asymptotics in section 2.4, as the corresponding result for section 2.2 is a direct corolary of this proof. We start proving the following technical result: Proof. Let us remind that for all z ∈ B 1 (0). Using this fact together with the analyticity of Q(z d ) in B 1 (0) the result follows.

Now we list the main result of this section
Theorem 2.8. The following convergence property holds true in the long time limit t → ∞, provided the first moment of the initial condition is bounded.
Proof. Consider our solution Expanding around ζ = 1 we find and a(t) = t(4 + t) 2(2 + t) 2 . Substituting we find and the second summand the rest g (t) = 1 2πid where we have implicitly assumed that /d is an integer (otherwise the integral vanishes) and where and δ > 0 is large enough so that the pole at ζ t lies in the area enclosed by |z| = 1+δ.
Computing the limit at the residue by means of l'Hôpital rule we find as in the previous cases, and the integral may be bounded as in the last section as well 1 2πid for long times and uniformly in . Note that we have made explicit the assumption of /d being an integer by means of the introduction of the Dirac comb. Let us now return to the rest g (t). We have the following bounds on it (the constant C may change from line to line): S(ζ, t). Integrating by parts we have To bound these integrals we take into account where 0 > 0 is a small enough constant and this last quantity, together with λ(t), admits also a constant upper bound. These considerations lead to Using that S 1 (ζ, t) = o(|ζ − 1|) and ∂ ζ S 1 (ζ, t) = o(1) we find where δ 1 > 0 is small enough. On the other hand we have C t 2 |ζ−1|≥δ1 because ∂ ζ S 1 (ζ, t) and S 1 (ζ, t) are bounded on {|ζ| = 1} ∩ {ζ = 1}, as found in section 2.3, where we have shown that H(ζ, t) is free from zeros in {|ζ| ≤ 1}∩{ζ = 1} uniformly in t. As a consequence we have As these estimates hold identically for f + ( , t) the final result follows. 3. Self-similar asymptotics for the majority kernel.
3.1. Formal calculations concerning self-similarity of the solution. In this section we concentrate in the model given by Eq. (15). Let us remind that the key result in this case was the conservation of the number of particles travelling in either direction given by Eq. (18). This phenomenology is reminiscent of the one given by the random kernel for the symmetric initial condition, for which order is not found at the kinetic level either. So, given the last section results, one would expect self-similarity of the solution in the case of the majority kernel too, but now independently of the initial condition. Following this reasoning, and given the scaling found for the solutions corresponding to the random kernel in the last section, we formally propose the following scaling behavior for the solutions F ± (z, t) to Eqs. (16)- (17). The scaling functions would obey the system where ξ = (z − 1)t is the self-similar variable. Note that making the substitution W + = ϕ + − ϕ + (0) and W − = ϕ − − ϕ − (0) we arrive at the following system of ordinary differential equations (ODEs) Changing variables ξ = −e τ we arrive at the four-dimensional differential system By changing variables again ψ ± = W ± + 1 we find the new differential system subject to the initial conditions ψ ± (−∞) = 1, H ± (−∞) = 0, and which long time behavior is H ± (+∞) = 0, ψ ± (+∞) = 1 − ϕ ± (0).
We will interpret this differential problem as a dynamical system. In this dynamical system one finds two invariant hyperplanes {H ± = 0} and one invariant plane {ψ + = ψ − , H + = H − }. All the fixed points of this system sit in the plane {H + = H − = 0}; in fact, this is a degenerated plane with all its points being fixed for the dynamical system. This system admits the first integral of motion E = ψ + ψ − − H + − H − , and it takes the value E = 1 for the initial conditions we are considering. All this implies that the orbits we are interested in, that are attracted by the plane {H + = H − = 0} in both the positive and negative infinite time limits, approach the hyperbola ψ + ψ − = 1 when τ → ±∞.
In the case ϕ ± (0) = 2 we can find an analytic expression for the heteroclinic connection. In this case ψ + = ψ − ≡ ψ and H + = H − ≡ H, and so We find the solution where τ 0 is an arbitrary real constant. In the original variables the solution reads where ξ 0 = −e τ0 is an arbitrary negative constant.
In the general case we have not been able to obtain the profiles ψ ± analytically, but we can anyway recover the original variables to find or alternatively where ξ 0 is an arbitrary negative constant. In the following section we will examine in a more rigorous way the existence of these formal expressions as well as their significance.

Continuum Dynamics and Laplace Transform.
First of all we note that the self-similar form we have calculated in the last section corresponds to the solution of the continuum version of Eq. (15), which reads In this section we will focus on this continuum version rather than on the discrete dynamics. We note that both dynamics, continuum and discrete, should behave analogously in the case of clusters composed by large numbers of particles, and in particular this number must much larger than d.
The Laplace transform of the solution to Eq. (44) iŝ which is defined for Re(z) > 0. The Laplace transformed version of this equation According to the formal developments in the previous section, we assume the following self-similar form of the solution in Laplace spacê where ξ = zt is the self-similar variable. This corresponds to the following selfsimilar form in real space that is reminiscent of the one found in the analysis of the random kernel. Note the self-similar profiles Φ ± can be recovered from the inverse Laplace transform where γ ∈ R is large enough and we have used the Bromwich integral formula [1].
Repeating the calculations of the previous section, employing the same notation, we reduce the problem to studying the differential system subject to the conditions and H ± (+∞) = 0, ψ ± (+∞) = 1 − ϕ ± (0).
As we have shown in the previous section this system encodes the form of the selfsimilar profiles Φ ± . This differential system, as a problem concerning real functions of real variables, was already studied in the previous section.
In order to proceed with the inverse Laplace transform we need to extend the solution to this differential system to the complex plane. The idea is to perform an analytic continuation by considering the solution of Eqs. (47a, 47b, 47c, 47d), whose trajectories numerically define the self-similar profiles ϕ ± , to C. In order to do so one needs to be sure that this system is free of singularities in a suitable region of the complex plane. In particular, the inverse Laplace transform formula (46) makes sense for those solutions free of singularities for Re(ξ) > 0. We have checked this is actually true by numerically integrating the differential system. The solution for the "+" fields for a particular initial condition is represented in Figs. 1, 2, 3 and 4; the solution of the "−" fields is completely analogous to this one. The initial condition was chosen to be the fixed point which serves as the α−limit of the heteroclinic connection in real variables plus a small positive perturbation in both Im(H + ) and Im(H − ). If both perturbations are of negative sign then the solution behaves analogously, but if they are of opposite sign then the solution grows unboundedly. By using a family of initial conditions, we have numerically checked that in fact there are no singularities in the whole strip |Im(τ )| ≤ π/2, Re(τ ) ∈ (−∞, ∞). This implies that the functions ϕ ± , considered as functions of ξ (considered also as a complex variable) are analytic in the half-plane Re(ξ) > 0 and therefore the function Φ ± can be obtained by means of the inversion formula for the Laplace transform, Eq. (46). This numerically shows the existence of scaling solutions to coagulation Eq. (44) obeying the self-similar scaling (45). On the other hand we note have clarified under which conditions these solutions are selected. It is reasonable to expect that for sufficiently symmetric initial conditions they will  indeed be selected. However it is not so clear that the same will happen for very asymmetric initial conditions. We leave this question as an open problem.
As a final remark let us mention that it is possible to analytically show that the functions ϕ ± (considered as functions of the variable τ , which is now considered to be complex) can be analytically extended to the region |Re(τ )| > L, |Im(τ )| ≤ π/2 for some real L sufficiently large. This is so because near the fixed points one can analytically continue the system solution by means of a linear stability analysis. In this case we have where the φ ± (τ ) are small functions. In this case we can solve the system to find We can be sure that the solution to the linearized system is free of singularities, and thus it can be used to safely analytically continue the solution into a neighborhood of τ = +∞ and − ∞ in the complex plane.

4.1.
Details of the simulation. In order to test the applicability of some of the theoretical results, we have performed extensive numerical simulations of a particular aggregation process. In our model, we consider initially N (t = 0) ≡ N 0 clusters, each one with only one particle. They are randomly and uniformly placed on the interval [0, 1] (we use periodic boundary conditions). Once the cluster locations have been set, we assign to them velocities +1 or −1. We denote by p 0 the initial fraction of clusters that have velocity +1 and, consequently, 1 − p 0 is the initial fraction with velocity −1. Given p 0 we have considered two different ways in the setting of the initial condition: (i) We run over the initial N 0 clusters and, for each one of them, we draw a random number u extracted from a uniform distribution in the (0, 1) interval. If u < p 0 we assign the cluster the velocity +1, otherwise it is given the velocity −1. In this case the number of clusters N + (t = 0) that initially have the velocity +1 is a random variable that follows a binomial distribution with average value N + (t = 0) = p 0 N 0 and variance σ Whatever the initial condition, the dynamical process evolves in the same way. Clusters move with their assigned velocity and, when two clusters collide, they coagulate with probability p (all results in this section take p = 0.1) such that a larger cluster containing all particles of both clusters is formed and the total number of clusters is reduced by one. The new aggregated cluster moves right or left randomly with probability 1/2. When two clusters do not coagulate, they simply pass through each other keeping their velocities. The process is repeated until no more collisions are possible. This could happen because only one cluster N = 1 (containing all particles) remains, or because all remaining clusters move with the same speed, either +1 or −1. The process is repeated M times (M realizations) starting with different random locations and velocities of the N 0 clusters. We denote by N (t) the random variable that counts the number of remaining clusters at time t.

Number and distribution of final clusters.
Let us define the random variable t ∞ as the time it takes a particular realization to reach a steady state in which no further evolution is possible. The remaining number of clusters at this time is denoted by N ∞ = N (t > t ∞ ). In Table 1 we present the average and the standard deviation of both quantities as a function of the initial number of clusters N 0 . We also show in the table the probability that N ∞ = 1, i.e. that there remains a single cluster containing all N 0 particles at the end of the run.
We have also computed from the simulations the probability distribution of the number N ∞ of remaining clusters, f (N ∞ ). It turns out, see evidence in Fig. 5, that the dependence of f (N ∞ ) on the intitial number of particles N 0 can be described by the scaling law:  Table 1. Average and standard deviation of the final number of clusters N ∞ and the time t ∞ it takes to reach a state where no further evolution is possible, as a function of the initial number of clusters N 0 . These results are for the case in which the number of initial clusters with velocity +1 follows a binomial distribution with a probability p 0 = 1/2. Averages are over M = 10 5 realizations for N 0 = 10, 10 2 , 10 3 , 10 4 and M = 4 × 10 4 realizations for N 0 = 10 5 . The last column is the probability that there is a single cluster at the end of the run.  Table 1: N 0 = 10, empty rhombi, N 0 = 10 2 , triangles, N 0 = 10 3 , squares, N 0 = 10 4 , circles and N 0 = 10 5 , filled rhombi. According to the scaling law discussed in the text we plot the distribution function for the vari- The solid line is the fit to the half-Gaussian distribution, as given by Eq. (49).
Moreover, G + (x) can be fitted to a half-Gaussian distribution:  Table 1 as a function of the initial number of clusters N 0 , as well as the analytical expressions in (50) (dotted lines). Note that the theoretical expressions agree very well with the numerical results, specially for large N 0 . and σ = √ 2. This implies a mean value and standard deviation: Both expressions are in good agreement with the numerical results of the simulation, specially for large N 0 , see Fig. 6. As indicated in Table 1, the time needed to reach equilibrium decreases with N 0 . This a priori surprising result indicates that the more particles there are initially, the faster the final state is reached. Obviously, as the inital density of particles decreases with increasing N 0 , particles are initially closer on average as N 0 increases and more collisions and coagulations are produced in the initial stages, so explaining this paradoxical result. The data, see Fig. 7, seem to suggest a power-law relationship of the form t ∞ ∼ N −1/4 0 . However, as for large N 0 the fluctuations are larger than the mean value, σ[t ∞ ] > t ∞ , the average value itself does not make much sense as a representative time to reach the asymptotic regime.
Finally, we have computed the probability distribution f (n) of the sizes of clusters at equilibrium. This is defined as the probability that a particle belongs to a cluster of size n. More precisely, f (n) is computed as the average (over realizations) number of clusters of size n multiplied by n and divided by N 0 , such that the normalization The results are plotted in Figs. 8 for N 0 = 10 and N 0 = 10 5 . Intermediate values of N 0 showing a similar behavior. It is not clear from the data whether a scaling law valid for all values of N 0 exists for these functions. Note the existence of the maximum at n = N 0 indicating that the most probable outcome is that a particle belongs to the single cluster containing all particles. This probability can be computed from Eq. (48) taking N ∞ = 1. In the limit of large N 0 the argument of the half-Gaussian distribution is x = 1 · N −1/2 0 → 0 implying that G + (0) = 1/ √ π, so that the value at the maximum is N −1/2 0 / √ π. As it can be seen in Table 1 and in Fig. 9, this result agrees well with the simulation results.

4.2.2.
Time dependence of fluctuations. We have computed the mean value Z(t) and the fluctuations σ 2 [Z(t)] of the variable Z(t) = N + (t) − N − (t), the difference between the number of clusters that have velocity +1, N + (t) and those that have velocity −1, N − (t). The data (not shown) indicates that Z(t) can be well approximated by a Gaussian distribution of mean Z(t) and r. m. s. σ[Z(t)]. We have found that, in accordance with the theoretical results, Z(t) is constant with time whereas the fluctuations can be decomposed as σ 2 [Z(t)] = σ 2 [Z(t = 0)] +σ 2 [Z(t)], beingσ 2 [Z(t)] independent on whether the initial condition for the number of clusters with velocity +1 was a fixed number [p 0 N 0 ] (σ[Z(t = 0)] = 0) or it followed a binomial distribution (σ 2 [Z(t = 0)] = 4p 0 (1 − p 0 )N 0 ). In the steady state, the Figure 8. Probability distribution f (n) of the sizes of clusters at equilibrium for N 0 = 10 (empty symbols) and N 0 = 10 5 (filled symbols, for clarity, in this case the vertical axis has been rescaled by an arbitrary factor of 100). Note the maximum at n = N 0 .
This expression has been checked in the case p 0 = 1/2 by plottingσ[Z(t)]N −1/2 0 vs √ N 0 t, see bottom panel of Fig.11, with a value a ∼ 0.25.

4.3.
Connection with theoretical results. It is obvious that the results of the simulations go much further beyond those of the kinetic theory. Yet, we have found agreement in two theoretical predictions: 1) the conservation in time of Z(t) and 2) the trend towards the formation of a cluster containing all particles in the system for p 0 = 1/2. It is difficult to compare many more predictions apart from these two, because the kinetic approach neglects many sources of fluctuations by its very nature, while the simulations retain all of them. Anyway, we have found result (53) especially interesting and, although a direct comparison with the theory is not possible, we will offer a possible theoretical explanation of it using quantities that can theoretically accessed.
Eq. (53) can be expanded around t = 0 to find This suggests the random variable Z(t), in the limit t → 0 + , performs a random walk with zero mean and diffusion constant a 2 N 2 0 / ≈ N 0 /16. Then the diffusive behavior is modified due to saturation effects as the system approaches the ordered state giving rise to the behavior described in (53). Now we will calculate a couple of quantities that might result of interest to interpret the value of the above found diffusion constant. From now on the condition N + (0) = N − (0) on the initial condition is imposed and all the clusters are initially of size = 1 (same conditions as in the simulations). First of all we consider the mean cluster size of all clusters traveling in the + direction (identical formulas hold for the other direction), which is given by the formula The second quantity of interest is the average size of the cluster the particles belong to, again for the + direction. It is given by the formula In both cases the calculations are performed following the techniques introduced in section 1. The second result means the particles are progressively located in clusters of bigger size, and per unit of time this size increases in N 0 /2 in average. This factor describes the average number of particles traveling in the opposite direction, and correspondingly it is the number of possible collisions for a given particle. The first result indicates how the mean cluster size grows in time. Per unit time this size increases by a factor (1/2) × (N 0 /2). Again we find that the mean cluster size increases with a velocity proportional to the number of possible collisions, but this time the proportionality factor 1/2 signals that after every collision one cluster disappears and the newborn cluster chooses its direction of motion randomly. In view of these results it is tempting to interpret the value of the diffusion constant (N 0 /4) × (N 0 /4) as the square of the rate at which the mean cluster size increases.