Synchronising and Non-synchronising Dynamics for a Two-species Aggregation Model

This paper deals with analysis and numerical simulations of a one-dimensional two-species hyperbolic aggregation model. This model is formed by a system of transport equations with nonlocal velocities, which describes the aggregate dynamics of a two-species population in interaction appearing for instance in bacterial chemotaxis. Blow-up of classical solutions occurs in finite time. This raises the question to define measure-valued solutions for this system. To this aim, we use the duality method developed for transport equations with discontinuous velocity to prove the existence and uniqueness of measure-valued solutions. The proof relies on a stability result. In addition, this approach allows to study the hyperbolic limit of a kinetic chemotaxis model. Moreover, we propose a finite volume numerical scheme whose convergence towards measure-valued solutions is proved. It allows for numerical simulations capturing the behaviour after blow up. Finally, numerical simulations illustrate the complex dynamics of aggregates until the formation of a single aggregate: after blow-up of classical solutions, aggregates of different species are synchronising or nonsynchronising when collide, that is move together or separately, depending on the parameters of the model and masses of species involved.


Introduction
Aggregation phenomena for a population of individuals interacting through an interaction potential are usually modelled by the so-called aggregation equation which is a nonlocal nonlinear conservation equation.This equation governs the dynamics of the density of individuals subject to an interaction potential K.In this work, we are interested in the case where the population consists of two species which respond to the interaction potential in different ways.In the one-dimensional case, the system of equations writes: (1.1) where θ α , χ α for α = 1, 2 are positive constants.
In this work, we are interested in the case where the interaction potential K in (1.1) is pointy i.e. satisfies the following assumptions: The aggregation equation arises in several applications in biology and physics.In fact, it is encountered in the modelling of cells which move in response to chemical cues.The velocity of cells a(ρ) depending on the distribution of nearby cells represents the gradient of the chemical substance which triggers the motion.Cells gather and form accumulations near regions more exposed to oxygen as observed in [20,24].We can also describe the movement of pedestrians using the aggregation equation as in [16] where the velocity of pedestrians is influenced by the distribution of neighbours.This equation can also be applied to model opinion formation (see [25]) where interactions between different opinions can be expressed by a convolution with the kernel K.
From the mathematical point of view, it is known that solutions to the aggregation equation with a pointy potential blow up in finite time (see e.g [11,5,2]).Then global-in-time existence for weak measure solutions has been investigated.In [5], existence of weak solutions for single species model has been obtained as a gradient flow.This technique has been extended to the two-species model at hand in [11].Another approach of defining weak solution for such kind of model has been proposed in [18,17] for the single species case.In this approach, the aggregation equation is seen as a transport equation with a discontinuous velocity a(ρ).Then solutions in the sense of duality have been defined for the aggregation equation.
Duality solutions has been introduced for linear transport equations with discontinuous velocity in the one-dimensional space in [3].Then it has been adapted to the study of nonlinear transport equations in [4,18,17].In [18,17], the authors use this notion of duality solutions for the one-species aggregation equation.Such solutions are constructed by approximating the problem with particles, i.e. looking for a solution given by a finite sum of Dirac delta functions.Particles attract themselves through the interacting potential K, when two particles collide, they stick to form a bigger particle.
In this work, we extend this approach to the two species case.To do so, we need to modify the strategy to the problem at hand.Indeed, collisions between particles of different species are more complex: particles can move together or separately after collision.This synchronising or non-synchronising dynamics implies several difficulties for the treatment of the dynamics of particles.In fact, particles of different species can not stick when they collide.Then an approximate problem is constructed by considering the transport equation with the a regularized velocity.Then measure valued solutions are constructed by using a stability result.
An important advantage of this approach is that it allows to prove convergence of finite volume schemes.Numerical simulations of the aggregation equation for the one-species case, which corresponds to the particular case of (1.1) when setting ρ 2 = 0, have been investigated by several authors.In [8] the authors propose a finite volume method consistent with the gradient flow structure of the equation, but no convergence result has been obtained.In [9], a Lagrangian method is proposed (see also the recent work [7]).For the dynamics after blow up, a finite volume scheme which converges to the theoretical solution is proposed in [19,6].In the two-species case, the behaviour is more complex since the interaction between the two species can occur and they may synchronise or not i.e. move together or separately depending on the parameters of the models and the masses of species.A numerical scheme illustrating this interesting synchronising or non-synchronising dynamics is provided in Section 6.In addition, a theoretical result on the convergence of the numerical approximation obtained with our numerical scheme towards the duality solution is given.Such complex interactions phenomena have been observed experimentally in [13].
System (1.1) can be derived from a hyperbolic limit of a kinetic chemotaxis model.In the case of two-velocities and in one space dimension, the kinetic chemotaxis model is given by ( where f ε α (x, v, t) stands for the distribution function of α-th species at time t, position x and velocity v, S ε (t, x) is the concentration of the chemical substance, T α [S](v, v ′ ) is the tumbling kernel from direction v ∈ V to direction v ′ ∈ V and ε > 0 is a small parameter.This tumbling kernel being affected by the gradient of the chemoattractant, is chosen as in [12] (1.3) where ψ α is a positive constant called the natural tumbling kernel and χ α is the chemosensitivity to the chemical S.This kinetic model for chemotaxis has been introduced in [21] to model the run-and-tumble process.Existence of solutions to this two species kinetic system has been studied in [14].
Summing and substracting equations (1.2) with respect to v = ±1 for f ε α yields (1.4) where ). Taking formally the limit ε → 0 in (1.5), we deduce that J ε α ⇀ χ α ∂ x S 0 ρ 0 α in the sense of distributions.Injecting in (1.4), we deduce formally that ρ 0 α satisfies the limiting equation: , where S 0 satisfies the elliptic equation: . This latter equation can be solved explicitly on R and S 0 is given by (1.7) Then we recover system (1.1).This formal computation can be made rigorous.The rigorous derivation of (1.6) from system (1.2) will be proved in this work.
The paper is organized as follows.We first recall some basic notations and notions about the duality solutions and state our main results.Section 3 is devoted to the derivation of the macroscopic velocity used to define properly the product a(ρ)ρ α and duality solutions.Existence and uniqueness of duality solutions are proved in Section 4, as well as its equivalence to gradient flow solutions.The convergence of the kinetic model (1.2) as ε → 0 towards the aggregation model (1.6)-(1.7) is shown in Section 5. Finally, a numerical scheme that captures the synchronising and non-synchronising behaviour of the aggregate equation is studied in Section 6, as well as several numerical examples showing the synchronising and non-synchronising dynamics.

Notations and main results
2.1.Notations.We will make use of the following notations.Let T > 0, we denote is the space of continuous functions that vanish at infinity.• M loc (R) is the set of local Borel measures, M b (R) those whose total variation is finite: is the space of time-continuous bounded Borel measures endowed with the weak topology.
We notice that if K satisfies (H2) and (H4), we have by taking y = −x in (H4) and using (H2) that, ∀ x ∈ R, ∂ x K(x)x ≤ λx 2 .We deduce that (H4) holds for ∂ x K i.e.: (2.1) We recall a compactness result on If there exists a sequence of bounded measures µ n ∈ M b (R) such that their total variations |µ n | (R) are uniformly bounded, then there exists a subsequence of µ n that converges weakly to µ in M b (R).
2.2.Duality solutions.For the sake of completeness, we recall the notion of duality solutions which has been introduced in [3] for one dimensional linear scalar conservation law with discontinuous coefficients.Let us then consider the linear conservation equation: with T > 0. We assume weak regularity of the velocity field b ∈ L ∞ (]0, T [×R) and b satisfies the so-called one-sided Lipschitz (OSL) condition: In order to define duality solutions, we introduce the related backward problem (2.4) We define the set of exceptional solutions E as follows The following result shows existence and weak stability for duality solutions provided that the velocity field satisfied the one-sided-Lipschitz condition.
Consider a sequence ρ n ∈ S M of duality solutions to , where ρ ∈ S M is the duality solution to Moreover, bn ρ n ⇀ bρ weakly in M b (]0, T [×R).
2.3.Main results.Up to a rescaling, we can assume without loss of generality that the total mass of each species is normalized to 1. Then we will work in the space of probabilities for densities.
The first theorem states the existence and uniqueness of duality solutions for system (1.1) and its equivalence with the gradient flow solution considered in [11].
Definition 2.4.(Duality solutions for system (1.1))We say that (ρ in the sense of duality on (t 1 , t 2 ) and â(ρ) = ∂ x K * ρ a.e.We emphasize that the final datum for (2.5) should be t 2 instead of T .
In our second main result, we prove the convergence of the kinetic model (1.2) towards the aggregation model.Theorem 2.6 (Hydrodynamical limit of the kinetic model).Assume that χ α (θ where ρ α is the unique duality solution of (1.6) and S = K * (θ 1 ρ 1 + θ 2 ρ 2 ) given in Theorem 2.5.
The condition χ α (θ 1 + θ 2 ) < 1 in the previous theorem is needed to guarantee that the tumbling kernel T α [S] defined in (1.3) is positive.
To conclude this Section on the main results, we emphasize that, a finite volume scheme to simulate (2.5) is proposed in Section 6 and its convergence towards duality solutions is stated in Theorem 6.3.

Macroscopic velocity
In this section, we find the representative â of a for which existence and uniqueness of duality solutions hold.To this end, we consider the similar system of transport equations to (1.1) associated to the velocity a n which converges to a. Next, the limit of the product a n (ρ n α )ρ n β for α, β = 1, 2 is computed.
3.1.Regularisation.We build a sequence (a n ) n∈N which converges to a by considering the sequence of regularised kernels be the sequence of regular kernels defined by ) in the sense of distributions.In addition, if we take x = − 1 n and y = 1 n in (H4), we have that n ] in the sense of distributions.Finally, we obtain that ∂ xx K n ≤ λ in the sense of distributions.
In the rest of the paper, the notation ∂ x K n will refer to the regularised kernels of Lemma 3.1.Given ∂ x K n , the velocity a n is defined similarly to (2.6) as In the following lemma, we show that if ρ n α and ρ n β admit weak limits ρ α and ρ β respectively in S M , then the limit of the product a n (ρ n α )ρ n β is â(ρ α )ρ β .Contrary to [22] where the twodimensional case is considered, this limiting measure does not charge the diagonal.
Proof.Before starting the proof of the lemma, we first introduce some notations which simplify the computations Step 1: Convergence almost everywhere in time of A n (t).
Since ∂ x K n (0) = 0, we have where I n (t) and II n (t) are defined by From the definition of ∂ x K n in Lemma 3.1, it follows that The estimate on ∂ x K n L ∞ in Lemma 3.1 and (H3) imply that with µ n and E n defined in (3.2).Let ε > 0. Since the set E n converges to the empty set, there exists N ∈ N such that ∀n ≥ N , For all n ≥ N , we observe that E n ⊂ E N , we have From the weak convergence of ρ n α , α = 1, 2, we note that the sequence µ n (t, •) converges weakly to µ(t, •).Since the total variation of µ n (t, •) is constant in n, the tight convergence is achieved.Then, there exists Hence, for all n ≥ N ′ , we get We deduce that I n (t) −→ 0.
Next, we show that II n (t) tends to zero.
where R is an integer which will be fixed later.From the construction of ∂ x K R in Lemma 3.1, we get Therefore, one has (3.4), by the same token as previously, there exists N such that for all n ≥ N , ) is a continuous function that vanishes on the diagonal (x, x) and we have The tight convergence of µ n to µ implies that there exists N ′′ > 0 such that for all n ≥ N ′′ Therefore for all n ≥ max{N ′ , N ′′ }, one has This implies that II n (t) converges to 0. Combining (3.5) and (3.6), we deduce that for almost every t ∈ [0, T ], A n converges to 0.
Step 2: Lebesgue's dominated convergence theorem For all t ∈ [0, T ], we have that Since A n converges almost everywhere to 0, T 0 A n (t)dt converges to zero from Lebesgue's dominated convergence theorem.
Proof.For x, y ∈ R, we compute: Using the λ-concavity of K, we deduce from (2.1) Since K n is also λ concave from the proof of Lemma 3.1, we get the one-sided Lipschitz estimate on a n by the same token as for a.

Existence and uniqueness of duality solutions
4.1.Proof of the existence of duality solutions in Theorem 2.5.The proof is divided into several steps.First, we construct an approximate problem for which the existence of duality solutions holds.Then, we pass to the limit in the approximate problem to get the existence of duality solutions thanks to the weak stability of duality solutions stated in Theorem 2.3 and recover Equation (2.5) from Lemma 3.2.Finally, we recover the bound on the second order moment.
Step 1: Existence of duality solutions for the approximate problem The macroscopic velocity a is replaced by an approximation a n defined in (3.1) and the following system is considered: m an approximation of ∂ x K n obtained by a convolution with a molifier.The solution ρ n,m α to the following equation is investigated.

dy).
Applying Theorem 1.1 in [10] gives the existence of solutions α is a duality solution.In addition, for φ ∈ C ∞ c (R) we have for α = 1, 2 the following estimate: Then, Using (4.3) and the density of C ∞ c (R) in C 0 (R), we deduce that ρ n,m α ∈ S M .Moreover, the equicontinuity of ρ n,m α in S M follows from (4.3) and the density of Theorem gives the existence of a subsequence in m of ρ n,m α which converges to a limit named ρ n α in S M .We pass to the limit when m tends to infinity in Equation (4.2) and obtain that ρ n α satisfies (4.1).
Step 2 : Extraction of a convergent subsequence of ρ n α and existence of duality solutions.
As above, there exists a subsequence of ρ n α in S M such that ρ n α ⇀ ρ α weakly in S M , for α = 1, 2. Let us find the equation satisfied by ρ α in the distributional sense.Let φ be in C ∞ c ([0, T ] × R).Since ρ n α satisfies (4.1) in the distributional sense, we have Using Lemma 3.2, we can pass to the limit in the latter equation and obtain, Thus ρ α satisfies (2.5) in the sense of distributions.From Proposition 3.3, the macroscopic velocity a n (ρ n ) satisfies an uniform OSL condition.Then, by weak stability of duality solutions in (see Theorem 2.3 (3)), we deduce that in the sense of duality in ]0, T [×R.
Step 3 : Finite second order moment.
From Equation (2.5), we deduce that the first and second moments satisfy in the sense of distributions Since ρ ini α ∈ P 2 (R) and â(ρ) is bounded from ((H2)), we deduce that the first two moments of ρ α are finite, then ρ α (t) ∈ P 2 (R) for t > 0.
Remark 4.1.If we define the weighted center of mass of the system x c as follows: We remark from straightforward computation that d dt x c = 0. Then the weighted center of mass is conserved for this system.4.2.Proof of the uniqueness of duality solutions in Theorem 2.5.Uniqueness relies on a stability estimate in Wasserstein distance, which is the metric endowed in P 2 (R).This Wasserstein distance d W is defined by (see e.g.[26,27]) , where Γ(µ, ν) is the set of measures on R 2 × R 2 with marginals µ and ν, i.e., ξ(y 1 )γ(dy 0 , dy 1 ) = R ξ(y 1 )ν(dy 1 ) .
The Wasserstein distance d W takes a more pratical form in the one-dimensional setting.Indeed, in one space dimension, we have (see e.g [23,26]) where F −1 ν and F −1 µ are the generalised inverse of cumulative distributions of ν and µ, defined by This Wasserstein distance can be extended to the product space P 2 (R) × P 2 (R).In the case at hand, we define W 2 (ν, µ) by µα are the generalised inverse of cumulative distributions of ν α and µ α for α = 1, 2, respectively.Using W 2 we prove a contraction inequality between duality solutions of (1.1).
Then W 2 (µ, ν) defined in (4.4) is bounded and satisfies the estimate: For the sake of clarity in the proof, we denote We also omit the argument t in notations F −1 α (t, x) and G −1 α (t, x).Computing the derivative of W 2 (µ, ν) 2 with respect to time, Straightforward and standard computations give that From the definition of â in (2.6), we get Setting z = F −1 1 (y) in the first integral and z = F −1 2 (y) in the second one yields Similarly, we get Using the oddness of ∂ x K, we can symmetrise the terms in the right-hand side of Similar computations can be carried out for Applying inequality (2.1) to (4.5) and using Young's inequality yields By definition of W 2 (4.4), we conclude Then the result follows from Gronwall's Lemma.
Proof of uniqueness.From Proposition 4.2, it is clear that if µ ini = ν ini , then µ = ν.We deduce uniqueness of duality solution in Theorem 2.5.

Convergence for the kinetic model
The convergence of the kinetic model (1.2) towards the aggregation model is analysed in this section.
Proof of Theorem 2.6.From the assumption χ α (θ 1 + θ 2 ) < 1 for α = 1, 2, we obtain that T α [S] defined in (1.3) is positive.Since T α [S] is a bounded and Lipschitz continuous function, we get the global in time existence of solutions to (1.2) and we have that To prove the convergence result stated in Theorem 2.6, we consider the zeroth and first order moments of the distribution f ε α (x, v, t) introduced previously.
), for α = 1, 2. From (1.2), these moments satisfy the following equations (5.1) From the first equation of (5.1), we deduce that ∀t ∈ [0, T ], , using the same token as in the proof of the existence, there exists ρ α ∈ S M such that ρ ε α ⇀ ρ α weakly in S M , for α = 1, 2. From the second equation of (5.1), we have We have that R ε converges weakly to zero in the sense of distributions.From Lemma 3.2, one obtains We conclude that 2 )ρ α in the sense of distributions.Passing to the limit in the first equation of (5.1), we deduce that ρ α satisfies (2.5) in the sense of distributions.We use uniqueness of duality solutions to conclude the proof.

Numerical simulations
This section is devoted to the numerical simulation of system (2.5).We provide a numerical scheme which preserves basic properties of the system such as positivity, conservation of mass for each species and conservation of the weighted center of mass.Moreover, we prove the convergence of the numerical approximation towards the duality solution defined in Theorem 2.5.
6.1.Numerical scheme and properties.Let us consider a cartesian grid of time step ∆t and space step ∆x.We denote x j = j∆x, j ∈ Z, t n = n∆t, n ∈ N.An approximation of ρ α (t n , x j ) denoted ρ n α,j is computed by using a finite volume approach where the flux F n α,j−1/2 is given by the flux vector splitting method (see [15]).Assuming that (ρ n α,j ) are known at time t n , we compute ρ n+1 α,j by the scheme: where (•) + := max{(•), 0}, (•) − := min{(•), 0} are respectively the positive and negative part of (•).Then we reconstruct where δ x j is the Dirac delta function at x j = j∆x.We first verify that this scheme allows the conservation of the mass and of the weighted center of mass.
Proposition 6.1.Let us consider (ρ ini 1 , ρ ini 2 ) ∈ P 2 (R) 2 such that for α = 1, 2, ρ ini α = j∈Z ρ 0 α,j δ x j .We assume that for n ∈ N * , (ρ n α,j ) j,n are given by the numerical scheme (6.1).Then the conservation of the mass of each species and of the weighted center of mass hold: Proof.Identity (6.2) can be obtained directly by summing over j ∈ Z the first equation in (6.1).
We now show (6.3).Multiplying by x j the first equation in (6.1) and summing over j ∈ Z, one gets Using a discrete integration by parts, one gets Finally, we get From (6.1), we have that Swapping indices i and j and using the oddness of Then (6.3) follows.
Lemma 6.2.Let (ρ ini 1 , ρ ini 2 ) be in P 2 (R) 2 such that ρ ini α = j∈Z ρ 0 α,j δ x j with j∈Z ρ 0 α,j = 1 and ρ α,j ≥ 0 for α = 1, 2. Assuming that for n ∈ N * , (ρ n α,j ) j,n are given by the numerical scheme (6.1).If the following CFL condition holds Then for all n ∈ N, ρ n α,j ≥ 0 and we have sup j,n Proof.This result is proved by induction.Let us assume that at time n, for all j ∈ Z, ρ n α,j is positive and sup j ân . From (6.1), it follows that (6.5) Using the condition (6.4) and the fact that sup j ân j ≤ ∂ x K L ∞ (θ 1 + θ 2 ), we get that ∆t ∆x ân j < 1. Therefore ρ n+1 α,j is positive as a linear combinaison of positive numbers.Then, recalling the expression of ân j given in (6.1), using the fact that ρ n+1 α,j , j ∈ Z, are positive and the conservation of the mass, (6.2) yields 6.2.Convergence of the numerical solution to the theoretical solution.In this part, we prove that the numerical scheme given in (6.1) converges to the duality solution obtained in Theorem 2.5.Theorem 6.3 (Convergence of the numerical scheme).Let T > 0, ∆x > 0 and ∆t > 0 such that (6.4) is satisfied and denote N t = T ∆t .Let ρ ini α ∈ P 2 (R), we define where (ρ n α,j ) j,n computed by (6.1).Then, we have where ρ α is the duality unique solution of Theorem 2.5 with initial data ρ ini α .
Proof of Theorem 6.3.For the initial data, it is clear that when ∆x → 0, we have ρ α,∆x (t = 0) ⇀ ρ ini α weakly.From Lemma 6.2, we get that for all j ∈ Z, n ∈ N, values of ρ n α,j are positive.
Step 1: Extraction of a convergent subsequence Equation (6.2) implies that the total variation of ρ α,∆x is fixed and independant of ∆x.
Therefore, there exists a subsequence of ρ α,∆x that converges weakly to Step 2: Modified equation satisfied by From the definition of ρ α,∆x in (6.6), we have Here and below we use < •, • > to denote the dual product in the sense of distributions.Discrete integration by parts yields where we use the notation φ n j := φ(x j , t n ).Using (6.5) and applying transformations to indices yields Taylor expansions gives the existence of ζ j in (x j , x j+1 ) and ζj in (x j−1 , x j ) such that Putting together, one obtains where R 1 α (∆x, ∆t) is given by From (6.6) and the definition of â in (2.6), we have where ân j are defined in (6.1).We get that From the Taylor expansion of ∂ x φ(x j , t): where R 2 α (∆x, ∆t) is defined as follows: The modified equation satisfied by ρ α,∆x in the distributional sense writes: From Lemma 6.2, we deduce that the terms R 1 α and R 2 α satisfy the estimates: , where C stands for a nonnegative constant.Passing to the limit and using the technical Lemma 3.2, we conclude that the limit ρ α satisfies (2.5) in the distributional sense with the expression (2.6) for the velocity.By uniqueness result in Theorem 2.5, we deduce that ρ α is the unique duality solution of (1.1).6.3.Dynamics of aggregates and numerical simulations.In this part, we carry out simulations of Equation (2.5) obtained thanks to scheme (6.1).Before numerically simulating the hydrodynamic behavior of the chemotaxis model, we first clarify the aggregate dynamics of this model, especially on the synchronising dynamics between aggregates of different species.
For the sake of simplicity, we choose θ 1 = θ 2 = 1 and K = 1 2 e −|x| in (2.5), which corresponds to the particular case of bacterial chemotaxis (see (1.7)).To illustrate the synchronising dynamics of the aggregates for (2.5), we consider the initial data given by sums of aggregates and look for a solution in the form We denote by u 1 and u 2 antiderivatives of ρ 1 and ρ 2 , respectively.Then the equation (2.5) reads (6.7) in the sense of distributions.Direct computation shows that Injecting these expressions into equation (6.7), the positions x k and y k satisfy the system of ODEs We recover the same system for particle solutions as in DiFrancesco and Fagioli [11] for two species.See also similar aggregate dynamics for single species in [5,18].In the case of one single species, the system of ODEs is determinant before any collision of aggregates, and after each collision, one can always 'restart' the particle system till final collapse of all aggregates.However, the case of collisions between particles of different species is more complex, since it does not necessarily imply whether the particles of different species will synchronise or not after colliding.In fact, as observed in the following simulations, both 'synchronising'(colliding particles of different species staying together) and 'non-synchronising' cases can occur, and the transitions between the synchronising types may happen, depending on the weighted attraction of other aggregates acting on them.
For illustration, we assume that two points of different species collide at a time t 0 .For instance, take x k (t 0 ) = y k (t 0 ) for some k, then at this time t 0 we have Note that γ k characterises external weighted attraction on ν k and µ k , depending on chemosensitivities, distances to other aggregates and the masses of all other aggregates.
Thus if χ 1 = χ 2 the velocity of species 1 and 2 is not the same at this time t 0 .However, with the special case at hand, K(x) = 1 2 e −|x| , we have We deduce that when x k < y k and x k → y k we have Obviously, in this case, particles µ k and ν k stay together if (y k − x k ) ′ (t) ≤ 0. On the other hand, when y k < x k and y k → x k we have In this case, particles µ k and ν k stay together when (x k − y k ) ′ (t) ≤ 0. Finally, to keep x k (t) = y k (t) for t ≥ t 0 , we need the condition (6.9) where γ k (t) is defined in (6.8).This relation characterises the weighted attraction of other aggregates acting on them.If the external weighted attraction on ν k and µ k (the left hand side of (6.9)) is small, they will stay together.When the external weighted attraction is big, the attraction between ν k and µ k is relatively weak and they will move separately, the one with bigger motility will move faster than the other.
We call (6.9) the synchronising condition for µ k and ν k .Similarly, we can get the synchronising condition for any µ i and ν j , ∀i, j.If more than two aggregates collide simultaneously, we can simply replace them by two aggregates of each species, each aggregate accumulating the total mass of each species.
In conclusion, according to the dynamics defined above, we can see that the initial aggregates will collapse such that they eventually form a single aggregate of the two species.The final aggregate can not separate, which is similar but illustrate more complex behaviour as one species case discussed in [18].Now we give some numerical examples showing "synchronising", "non-synchronising", transitions between "synchronising" and "non-synchronising" dynamical behaviours for the hydrodynamic model (2.5).