pp. X–XX TWO-SPECIES PARTICLE AGGREGATION AND STABILITY OF CO-DIMENSION ONE SOLUTIONS.

Systems of pairwise-interacting particles model a cornucopia of physical systems, from insect swarms and bacterial colonies to nanoparticle self-assembly. We study a continuum model with densities supported on co-dimension one curves for two-species particle interaction in $\mathbb{R}^2$, and apply linear stability analysis of concentric ring steady states to characterize the steady state patterns and instabilities which form. Conditions for linear well-posedness are determined and these results are compared to simulations of the discrete particle dynamics, showing predictive power of the linear theory. Some intriguing steady state patterns are shown through numerical examples.

1. Introduction.The collective behavior of interacting particle systems gives rise to emergent phenomena in physics, biology, chemistry, and other disciplines.Models of pairwise-interacting agents find applications in the biological contexts of locust swarms [3,36,35], animal flocks [8,23,29], and bacterial colonies [37].These mathematical approaches to swarming have also inspired algorithms for cooperative control of robotic vehicles [25].More questions for nonlocal particle systems arise in physical chemistry: the self-assembly of nanoparticles [6,17] and arrangement of ions into spheres [26,27] are just two examples.In the physical contexts of granular gasses [30] and molecular dynamics simulations of matter [16], particle systems also have a central role.
All of the above models, however multifaceted, share the same footing.Some number of particles interact with each other pairwise such that any two particles will repel each other when they are close and attract when they are far; typically, this attractive force disappears at very long distances.These interactions can generate rich steady states relevant to the models in which they arise.

ALAN MACKEY, THEODORE KOLOKOLNIKOV, AND ANDREA BERTOZZI
Consider the case when the forces arise due to a pairwise interaction energy where x i denotes the position in R d of the i th particle and P (r) is the potential energy between two particles.P (r) is usually a function with a unique minimum such that the force on one particle due to another, F (r) := −P (r), enjoys the repulsive-attractive properties mentioned above.In this framework, a steady state pattern can be understood as a minimizer of E.
We call the potential E confining if its minimizing configurations x 1 , . . ., x N stay contained inside a compact set as N → ∞.The question of whether or not a given function P will result in a confining potential has been addressed in terms of the notion of H-stability in statistical mechanics; see [10].For confining potentials, particles in ground states may reside in space-filling, co-dimension zero configurations or concentrate on co-dimension one manifolds.The question of which occurs is answered in [2] and the problem of characterizing ground states (or steady states) has been discussed in [19,20,33,39], and elsewhere.Applications where both codimension zero and one solutions are of importance include bacterial colony growth under stress, point vortex theory, and the Thomson problem [1,4,5,37].
It is a natural extension of the above work to consider the analogous problems for two particle species; i.e., when more than one type of particle is present in the interactions.Two-species models are relevant for the phenomena observed in [27], where two types of macroions in solution will self-recognize and assemble into hollow spherical structures.This self-recognition of particle species is a robust phenomenon observed in many of the numerical experiments considered in this paper.Two-species models also find application in large scale pedestrian movement [31], and the well-posedness of said models has been considered in [9]; a general treatment of well-posedness for the two-species problem is given in [15].Other applications include opinion formation in groups consisting of ordinary individuals and strong leaders [11] and two-species group consensus [12].Two-species bacterial aggregation driven by chemotaxis and diffusion is another area of active research, where [22] employs a two-species model for localized vortex formation in bacterial colonies.Global existence and finite time blowup are considered in [7] and [13], [40] treats the n-species problem, and [18] and [34] discuss the stability of uniform density and homogenous steady states.
Our numerical experiments have revealed phenomena which did not appear in the single species problem.Particle species either mix or segregate based on the relative strengths of the inter-species and intra-species forces, and occasionally settle in domains with irregular boundaries including cusps.Asymmetric steady states (which represent local minimizers of the potential) can be observed, and nontrivial structures form when particles are H-stable.These and other features indicate a substantial increase in complexity of the two-species problem over the single species problem.
In this paper, our objective is to characterize steady states formed in the twospecies aggregation problem in the absence of diffusion.Inspired by physical [26] and numerical experiments exhibiting steady states supported on or near co-dimension one manifolds, we consider the linear stability of such states.In particular, we study the continuum limit of the problem and the stability of a steady state consisting of two concentric rings of constant density.The theory put forth accurately predicts instabilities observed in numerical experiments and the breakup of ring solutions into fully two-dimensional patterns.The arguments presented here are restricted to the two-dimensional problem, but adapting the theory of [39] could generalize the results to higher dimensions.
2. Problem Description.Consider two species of particles, type I and type II, which occupy positions x 1 (t), . . ., x N1 (t), y 1 (t), . . ., y N2 (t) in R 2 .Steady state patterns are minimizers of the pairwise interaction energy , is simply a change of variables to simplify the calculations.
We are interested in gradient flow equations associated with (1): for i = 1, . . ., N 1 , and for j = 1, . . ., N 2 .In the above, g i (s) := − dVi ds (s), i = 1, 2, 3, give the interaction forces due to the potentials.The right-hand sides of (2) and (3) have been divided by N 1 as a simple rescaling of time, and in the second line the parameter µ := N 2 /N 1 has been introduced to facilitate the continuum limit considered next.One can think of the gradient flow either as an approximation to overdamped second order physical dynamics or simply as a means to identify minimizers of the energy.The gradient flow consists of a nonlocal PDE system of coupled advection equations similar to the Birkhoff-Rott equation for vortex sheets (c.f.[28,32]), and linear stability is reduced to a sequence of eigenvalue problems.Criteria for the stability of each element in a basis of perturbations, and for linear well-posedness of the concentric ring solution, are derived.Numerical examples are presented, which demonstrate strong agreement with the theory put forth.
In this work we consider the following potentials, which have all been considered in the literature for the single species problem [10,20,24,39]: the Morse potential and smoothed step discontinuity forces 2s with steady states pictured in figure 1. Combinations of all three of the above types (and others) are also plausible; for example, g 1 could arise from a power law, g 2 from the Morse potential, and g 3 from the tanh force.See figure 2. 3. The Continuum Limit.For the two-species case, we will say that an energy E such as (1) is confining if its minimizing configurations x 1 , . . ., x N1 , y 1 , . . ., y N2 stay contained inside a compact set as N 1 and N 2 → ∞.Under the assumption that the energy E is confining, the mass-normalized configurations of discrete particles approach continuum densities ρ 1 and ρ 2 .As we are interested in the stability of concentric ring solutions, we seek densities which are supported along one-dimensional curves Γ 1 (t) = Φ 1 (α, t) and Γ 2 (t) = Φ 2 (α, t) (parameterized by α ∈ D ⊂ R) which evolve with velocity fields v 1 and v 2 ; that is, The velocity fields v 1 and v 2 are determined from the respective densities by the continuum limits of equations ( 2) and (3): where we must assume that N 2 /N 1 → µ as N 1 , N 2 → ∞.To determine the dynamics of ρ 1 , ρ 2 , Φ 1 , and Φ 2 completely, we impose conservation of mass: which is implicit in the particle formulation of the problem.Equations ( 4), (5), and (6) specify a nonlocal coupled advection system determining Φ 1 and Φ 2 , from which ρ 1 and ρ 2 will be recovered later; see, for example, [32] and [39].
It is worth pointing out here that the ρ i are densities of measures singular with respect to Lebesgue measure on R 2 .Therefore, they should solve (6) weakly, or in the sense of distributions.We will assume there exist f i locally integrable on R such that for Borel sets where ρ s i admits the natural interpretation of the density along the surface One can now integrate by parts from (6) to define what it means for ρ i to be a solution: where the boundary term drops out because ψ is compactly supported.It follows that f (α, t) ≡ f (α, 0) =: f 0 (α).
We also rewrite (5) in terms of integrals along Γ 1 and Γ 2 : Appealing to (4) then yields and The two equations above determine Φ i .With these in hand, all that is left is to determine ρ i ; for this, ), the equations (5) evaluated for v 1 at x i and v 2 at y i reproduce (2) and (3) exactly, up to the time scaling introduced.
Before we proceed to the next section and linearize around the concentric ring steady state, it is worthwhile to consider the existence of such a steady state.If Φ 1 and Φ 2 parameterize concentric circles of radii R 1 and R 2 , we can take D = [−π, π) and Φ i (s) = Θ(s)R i e 1 (as in [39]) with With no motion in time, equations ( 7) and ( 8) give The (constant) densities and radii should satisfy ρ s i R i = m i , where m i is the total mass of particles of type i.In the discrete case, N i = m i and so we can say , and the above equations can be rewritten Appealing to some basic geometry, The first component of each integral cancels because it is odd on (−π, π), and so we are left with which determine R 1 and R 2 .So long as (9a) and (9b) have solutions, the concentric ring steady state exists.Of course, the integrands of (9a) and (9b) must be in and the g i have no singularities away from the origin.Assuming also that and observe that (9a) arises as ∂F ∂R1 = 0 and (9b) as ∂F ∂R2 = 0. Then if g 1 , g 2 , and g 3 are continuous (except perhaps at the origin-for commonly encountered potentials such as the Morse or Lennard-Jones potentials, this is the case) (9a) and (9b) will be satisfied at a maximum or minimum of F .To show, then, that solutions R 1 and R 2 exist, it suffices to show that F attains a minimum for some R 1 , R 2 > 0.
A general proof of this fact is difficult because the potentials V i will vary, and in some cases a concentric ring solution may not exist.However, for all cases pursued numerically, (9a) and (9b) have solutions R 1 , R 2 which do give rise to a steady state solution to ( 7) and (8).

4.
Linearization & Eigenvalue Problem.Recall that the rings have been parameterized as Φ i (s) = Θ(s)R i e 1 (where Θ is a rotation matrix).Consider now a small perturbation of each ring in the form Linearizing (7) and (8), then canceling the factor of e λt appearing in each term gives To simplify calculations below, define M = M (s, s ) := Θ −1 (s)Θ(s ), and and because Θ and M are unitary, Finally, In terms of these quantities, multiplying the linearized equations by Θ −1 and collecting terms multiplied by i (s) and i (s ) leaves R 1 sin(s − s) and u ⊗ v denotes the matrix with i, j entry u i v j .
Note that all the above matrices have even, periodic entries along the diagonals and odd, periodic entries off.With this in mind, consider (10) rewritten as (the superscripts are used to distinguish matrices, not as powers) where the diagonal entries of M 1 , M 2 , and M 3 are even and periodic, and off-diagonal entries are odd and periodic.It follows that is a constant diagonal matrix times 1 (s).For the other two terms, we hope for similar results to yield an eigenvalue problem in i and λ.Using the ansatz similar to that of [20], we compute the terms above involving M 2 and M 3 : The first entry is a linear combination of where a is a diagonal matrix.The third term will be similar and will give a diagonal matrix multiple of 2 , so that the equation for 1 becomes (a and b are matrix-valued functions) and the equation for x 2 sin(ns) .
Comparing coefficients of cos(ns) and sin(ns) results in an eigenvalue problem for x 1 , x 2 , y 1 , and y 2 : M 1 , . . ., M 6 are computed as in (12a-12d) and are shown below; here, u i , v, and w are functions of θ.M 1 and M 5 are diagonal and do not depend on n.
The first two matrices determine the stability of the type I particle ring with respect to frequency n perturbations, with the type II ring remaining fixed: The final two matrices M 5 and M 6 are analogous to M 1 and M 2 , except they determine the stability of the type II ring: after integrating by parts and using (9b), and 5. Linear Well-posedness.Here we consider the limit of the eigenvalue problem (13) as n → ∞.The goal is linear well-posedness; that is, to determine when the eigenvalues λ(n) of ( 13) satisfy λ(n) < 0 as n → ∞.That λ(n) → 0 as n → ∞ follows immediately from the Riemann-Lebesgue lemma; the requirement that the eigenvalues approach zero from below is important because it demonstrates that all but finitely many modes are stable.Intuitively, if modes of arbitrarily high frequency are unstable, the co-dimension one curve will break apart and the density will form a fully two-dimensional pattern.
Theorem 5.1 (Linear well-posedness).Assume that the forces have power series representations with a 0 , b 0 > 0 and c 0 = 0, valid in some neighborhood of the origin, where p 0 < p 1 < . . .etc. Define α = M 1 11 and β = M 5 11 .If R 1 = R 2 , then the concentric ring solution to (7,8) is linearly well-posed if and only if If R 1 = R 2 , g 1 = g 2 , and µ = 1, the concentric ring solution to (7, 8) is linearly well-posed if and only if (15a) holds and either r 0 is a nonnegative integer or r 0 > p 0 .
Before moving on to the proof, a remark is in order.Theorem 5.1 as stated does not cover the cases when R 1 = R 2 but g 1 = g 2 or µ = 1.However, (9a) and (9b) point out that unless µ = 1 and g 1 = g 2 , it is very unlikely that R 1 = R 2 ; for two arbitrary potentials g 1 , g 2 and ratio µ, it is a measure-zero type event that the radii equations would have such solutions.
Proof of Theorem 5.1.The analysis relies primarily on asymptotic expressions for the integrals occurring in M 1 , . . ., M 6 , which necessitates the assumptions of (14).Substituting ( 14) into the formulas for M 1 , . . ., M 6 leaves an eigenvalue problem where each entry of E is a potentially infinite sum of integrals.However, showing that E has negative eigenvalues is equivalent to showing its leading minors alternate sign, and it is easy to see that in each entry of E, only those terms which decay most slowly will affect the eigenvalues in the limit.
In practice, the values α = M 1 22 and β = M 2 22 must be evaluated analytically or numerically, because they are independent of n.Other entries of E may be evaluated asymptotically, and considering one of these entries gives an idea of how to proceed: where we used the trig identity By a similar use similar identities, it turns out that all entries of E reduce to linear combinations of one integral: and we are interested in the behavior of I for fixed c and p as n → ∞.
For c = 1 and p > −1/2 (which is necessary for the integrals to converge), an explicit formula with asymptotics is available from [39]: where C(p) > 0 is a positive constant depending on p.This asymptotic form may also be arrived at by stationary phase analysis.For c > 1, it can be shown readily via integration by parts that I decays faster than any polynomial: for any integer k, there exists a constant C(c, p, k) such that For M 2 and M 6 , ( 16) gives the relevant rates of decay.M 3 and M 4 are more complicated and the analysis breaks into cases.
Case I: R 1 = R 2 .In this case (17) shows that M 3 and M 4 approach zero faster than any of the other entries of E, and (13) asymptotically decouples into two quasi-single species problems Before moving on, it is worth pointing out that the decoupling has a pleasant physical interpretation.M 3 represents the effects of perturbations of the type II particle ring on the type I particle ring, and vice versa for M 4 .The fact that these vanish from (13) as n → ∞ means that, when R 1 = R 2 , high frequency perturbations of the ring of type I particles have no effect on the eventual stability of the ring of type II (and vice versa).
For each of the problems in (18) to have negative eigenvalues, it is necessary and sufficient that We turn first to M 2 : substituting in the assumed power series representation of g 1 and discarding all but the most slowly decaying terms, and The final entry is treated by integration by parts: The analogous work for M 6 looks almost exactly the same, except with q 0 replacing p 0 .With those expansions in hand, one can asymptotically compute the terms appearing in (19): and so we need only require that α < 0 and M 2 22 < 0. The asymptotic expression for the latter is negative so long as sin(πp 0 ) is, and this leads to (15a).The problem for M 5 + M 6 is nearly identical, and yields (15b) in exactly the same way.It is worth noting here that the criteria for linear well-posedness are very similar to those for the single-species case explored in [20] and [39].
Case II: R 1 = R 2 =: R. As mentioned earlier, this is very unlikely unless g 1 = g 2 and µ = 1; so, we will assume that is the case.Then M 1 = M 5 , M 2 = M 6 , and M 3 = M 4 , so E simplifies; however, the rate of decay of M 3 is not as fast now and so it must be taken into account.We determine when E has negative eigenvalues by checking when its leading minors alternate sign.Asymptotics for M 3 are necessary, but M 3 has the same form as M 2 with g 3 replacing g 1 and R replacing R 1 : The first minor of E is then (M 1 + M 2 (n)) 11 → α as n → ∞, so we must require that α < 0 .
The third minor begins to include terms from the cross-particle interaction force g 3 , and works out to be (to leading order in n) where C 1 and C 2 are some positive constants with respect to n.We have already required for the first and second minors that α < 0 and M 2 22 < 0, which is why the first term is negative and the second is positive.So we must require that r 0 is a nonnegative integer or 2p 0 + 1 < 4r 0 + 4.
The fourth minor works out to be (again to leading order) where C 1 and C 2 again denote positive constants.So we must require that r 0 is a nonnegative integer or r 0 > p 0 .
These restrictions also imply that the third minor is negative.This last minor gives the final requirement for the double ring solution to be linearly well-posed and yields the criterion in the theorem for the R 1 = R 2 case.Right: Separated particle ring.Forces g 1 and g 2 are the same as on the left, but g 3 = 1.01g 1 .
6. Numerical Examples.All numerical solutions here and in figure 1 were computed using a simple forward Euler scheme with an adaptive time step chosen as large as possible while requiring that the energy of the system (1) decreases at each time step.A scheme with higher order accuracy is not necessary, since we only seek a minimizer of the energy (1).Alternatively, choosing a time step based on an estimate of the local truncation error (as in [21]) is also efficient and yields the same results.All initial conditions are taken to be independently, uniformly distributed on a square.The theoretical predictions agree very well with numerical observations.Figure 3 shows an example in which the two particle species may mix or segregate based on the relative strengths of the inter-species and intra-species interactions.The destabilization of the alternating ring structure and appearance of mode two instability agrees exactly with the theory presented above.Figure 4 shows examples of mode five and mode one instabilities.The alternating, symmetric mode five arises from interaction forces defined in terms of The mode one instability is due to the interaction forces defined as Coupling effects of the rings on each other can be seen in figure 5, which shows coupling between type I particles (white) with a mode three instability and type II particles (black) with mode five.The interaction forces are defined in terms of (20) as with K = 0, 1, and 4, designed in [38] to exhibit pure mode 3 and 5 instabilities.The sequential disappearance of the instabilities in figure 5 corresponds to the eigenvalues of those modes becoming negative (i.e.eigenvalues from ( 13)), confirmed numerically.
Even when the concentric ring solution is not linearly well-posed, the stability or instability of low frequency modes is still borne out in the ground state.Figure 6 shows the occurrence of a mode two instability in such a scenario, which is correctly predicted by the linear stability theory even though the linearized equation is illposed.When multiple modes become unstable, it is possible that one, several, or all of the unstable modes will appear in the ground state.The question of those which do, and to what degree, is determined by the particular nonlinearity of the problem and is outside the scope of the linear theory.For the single species case, this was pointed out in [38] and the problem of which modes appear is still open.At this time, even less is known about the two-species problem.Weakly nonlinear analysis, considered in [33], may prove useful for this purpose and is one of several considerations for future work.There are several factors which limit the effectiveness of the linear theory when the stable ground state is far from the concentric ring solution.When R 1 = R 2 , as is the case in many examples with g 1 = g 2 , the inter-species interaction g 3 must be o(s −1/2 ) as s → 0 for the integrals in 13 to exist.This condition is not met in any of the cases of figure 1, and so the theory may not be applied.It is possible to replace g 1 (s) by g 1 ((1 + )s), and if is sufficiently small then the observed steady state is qualitatively indistinguishable from the unperturbed version but R 1 and R 2 are no longer equal.The theory does apply to this perturbed problem, but another difficulty arises: most or all modes are unstable.
This type of total instability also occurs in both panels of figure 2. While the instability of all low modes is certainly consistent with the observed steady states, Mode two becomes unstable for the power law potential.Left: g 1 (s) = g 2 (s) = s −0.15 , g 3 (s) = −s 0.15 .Right: it is largely uninformative.In the majority of cases when a mode 1 instability appears, all low frequency modes are unstable as well and the stable steady state is sometimes asymmetric (with the gradient flow coming to rest at a local minimum).Nevertheless, the linear stability theory may still provide some insight.In the right hand panel of figure 2, each low mode has one unstable eigenvector except for mode 3, which has two; however, only one and not a linear combination of both appears in the steady state.Why one and not the other or a combination of both appears is impossible to determine by the linear theory, similar to the problem encountered in [38].
That the theory works very well in the cases of figures 3, 4, 5, and 6 but not for figures 1 and 2 is not surprising-when the stable steady state is far from the concentric ring solution, the linear theory is less likely to apply.Two-particle minimizers also occasionally break symmetry (in the sense that the steady states for the type I and II particles differ by more than a simple reflection or rotation), pictured in the first panel of figure 2. Particle densities from figure 1 are asymmetric for certain parameters (see the sixth panel) and may be supported on domains with irregular geometry, including cusps.Some simulations show dependence on initial conditions, which did not manifest for the single species problem.Of course, it is not guaranteed that the gradient flow (2), (3) reach a global minimizer of the potential (1).This minor dependence on initial conditions seems to indicate that the energy landscape for the two-species problem is more complex, and that the gradient flow occasionally comes to rest at local minima.
Simple structures directly relevant to self-assembly, such as alternating particle rings or chains, may also be observed; see figures 3 and 7.The self-recognition and separation into two rings replicates the phenomena observed in [27], and occurs over a very long time scale (t ≈ 1.4 × 10 5 ) relative to that of the formation of the rings (t ≈ 10 2 ). 7. Discussion.Two-species particle aggregation systems have a rich solution structure, including densities with rings, spots, and radial or bilateral N -fold symmetry The first few panels show that the particles seem to form effective dipoles because the inter-species repulsion length scale is so small.When the number of particles increases, the confining nature of the potentials causes them to pack closer together and chains form, as in panel 5.As N increases further, the particles begin to form a two-dimensional lattice structure (panel 6).
which often concentrate on or near co-dimension one surfaces.Considering a continuum limit as the number of particles tends to infinity results in a PDE system formulation similar to that of the vortex sheet problem [28,32].Linear stability analysis successfully characterizes the steady states which form, verified numerically in section 6, and linear well-posedness of the PDE system is considered in section 5.
Future work could address the three and higher dimensional versions of the problem, weakly nonlinear analysis of bifurcations from rings to other steady states, the second-order problem, the n-species problem, and the inverse problem of constructing potentials with prescribed instabilities or patterns [38].In addition, the two-species system allows for the unique possibility of nontrivial H-stable ground states which are outside the scope of the co-dimension one analysis here; c.f. [24].
Figure 8 shows a power law example which completely leaves a co-dimension one manifold, and seems to exhibit an effective phase separation or surface tension arising from the nonlocal interactions.As the inter-species repulsion singularity becomes weaker than the intra-species repulsion singularity, black and white particles go from self-segregating to mixing.When the inter-species repulsion is substantially weaker, exhibited in the far right panel of figure 8, a regular alternating lattice structure emerges in large portions of the steady state.The formation of lattices in the single species problem is not unfamiliar; see [14] for similar phenomena in the single-species problem.The first two panels of figure 8 correspond to local, not global minimizers of the potential 1, and a steady state consisting of a straight line interface between the black and white particles has slightly lower energy.The steady states in the right two panels, however, have lower energy than the separated black/white half disks.The analysis carried out in this paper applies to solutions supported along one-dimensional curves, and as such does not apply to the phenomena in figure 8; however, the phase separation effects could be observed on a co-dimension one surface in R 3 .

Figure 5 .
Figure 5. Modes three and five stabilize each other as crossparticle attraction increases.

1 VFigure 7 .
Figure 7. Alternating particle chains arising from the Morse potential, numbers of particles N 1 = N 2 = N with N =8,20, 80, 200, 400, 800 from top left to bottom right.The first few panels show that the particles seem to form effective dipoles because the inter-species repulsion length scale is so small.When the number of particles increases, the confining nature of the potentials causes them to pack closer together and chains form, as in panel 5.As N increases further, the particles begin to form a two-dimensional lattice structure (panel 6).