Zoology of a non-local cross-diffusion model for two species

We study a non-local two species cross-interaction model with cross-diffusion. We propose a positivity preserving finite volume scheme based on the numerical method introduced in Ref. [15] and explore this new model numerically in terms of its long-time behaviours. Using the so gained insights, we compute analytical stationary states and travelling pulse solutions for a particular model in the case of attractive-attractive/attractive-repulsive cross-interactions. We show that, as the strength of the cross-diffusivity decreases, there is a transition from adjacent solutions to completely segregated densities, and we compute the threshold analytically for attractive-repulsive cross-interactions. Other bifurcating stationary states with various coexistence components of the support are analysed in the attractive-attractive case. We find a strong agreement between the numerically and the analytically computed steady states in these particular cases, whose main qualitative features are also present for more general potentials.


Introduction.
Multi-agent systems in nature oftentimes exhibit emergent behaviour, i.e. the formation of patterns in the absence of a leader or external stimuli such as light or food sources. The most prominent examples of these phenomena are probably fish schools, flocking birds, and herding sheep, reaching across scales from tiny bacteria to huge mammals. While self-interaction models for one particular species have been extensively studied, cf. Refs. [26,23,20,37,46] and references therein, there has been a growing interest in understanding and modelling interspecific interactions, i.e. the interaction among different types of species. One way to derive macroscopic models from microscopic dynamics consists in taking suitable scaling limits as the number of individuals goes to infinity. Minimal models for collective behaviour include attraction and/or repulsion between individuals as the main source of interaction, see [46,19,20,35] and the references therein. Attraction and repulsion are normally introduced through effective pairwise potentials whose strength and scaling properties determine the limiting continuum equations, see [39,9,8,16]. Usually localised strong repulsion gives rise to non-linear diffusion like those in porous medium type models [39], while longrange attraction remains non-local in the final macroscopic equation, see [16] and the references therein. In this paper we propose a finite-volume scheme to study two-species systems of the form ∂ t ρ = ∇ · ρ∇ W 11 ⋆ ρ + W 12 ⋆ η + ǫ(ρ + η) , (1a) with given initial data ρ(x, 0) = ρ 0 (x), and η(x, 0) = η 0 (x). (1c) Here, ρ, η are two unknown mass densities, W 11 , W 22 are self-interaction potentials (or intraspecific interaction potentials), W 12 , W 21 are cross-interaction potentials (or interspecific interaction), and ǫ > 0 is the coefficient of the cross-diffusivity. The non-linear diffusion term of porous medium type can be considered as a mechanism to include volume exclusion in cell chemotaxis [32,40,12], since it corresponds to very concentrated repulsion between all individuals. This model can also be easily understood as a natural extension of the well-known aggregation equation (cf. [37,46,3,18] ) to two species including a cross-diffusion term. Common interaction potentials for the one species case include power laws W (x) = |x| p /p, as for instance in the case of granular media models, cf. [2,47]. Another choice is a combination of power laws of the form W (x) = |x| a /a − |x| b /b, for −N < b < a where N is the space dimension. These potentials, featuring short-range repulsion and long-range attraction, are typically chosen in the context of swarming models, cf. [36,1,28,29,4,21,17]. Other typical choices include characteristic functions of sets (like spheres) or Morse potentials W (x) = −c a exp(−|x|/l a ) + c r exp(−|x|/l r ), or their regularised versions W p (x) = −c a exp(−|x| p /l a ) + c r exp(−|x| p /l r ), where c a , c r and l a , l r denote the interaction strength and radius of the attractive (resp. repulsive) part and p ≥ 2, cf. [26,22,21]. These potentials display a decaying interaction strength, e.g. accounting for biological limitations of visual, acoustic or olfactory sense. The asymptotic behaviour of solutions to one single equation where the repulsion is modelled by non-linear diffusion and the attraction by non-local forces has also received lots of attention in terms of qualitative properties, stationary states and metastability, see [11,15,27,13,14] and the references therein. Systems without cross-diffusion, ǫ = 0, were proposed in [24] as the formal mean-field limit of the following ODE systeṁ For symmetrisable systems, i.e. systems such that there exists some positive constant α > 0 with W 12 = αW 21 , they show the system can be assigned an interaction energy functional. As a result, the system admits a gradient flow structure and variational schemes can be applied to ensure existence of solutions, cf. [24,33]. However, in many contexts such a condition is too exclusive in the sense that lots of applications exhibit a lack of symmetry in the interactions between different species. In order to treat the system for general, and possibly different, cross-interactions W 12 , W 21 , they modify the well-known variational scheme to prove convergence even in the absence of gradient flow structure. These systems without cross-diffusion appear in modelling cell adhesion in mathematical biology with applications in zebrafish patterning and tumour growth models, see [30,25,41,48] for instance. In this paper we extend their cross-interaction model by a cross-diffusion term which is used to take into account the population pressure, i.e. the tendency of individuals to avoid areas of high population density. As cross-diffusion we choose the form introduced by Gurtin and Pipkin in their seminal paper [31]. Although their work is antedated by results of mathematicians and biologists interested in density segregation effects of biological evolution equations, cf. [44,43] and references therein, the particularity about their population pressure model is the occurrence of strict segregation of densities under certain circumstances, cf. [31,6,7]. This cross-diffusion term has been the basis to incorporate volume exclusions in models for e.g. tumour growth [5] or cell adhesion [38]. Hence, our model is of particular interest from a modelling point of view taking into account non-local interactions between the same species and different species as well as the urge of both species to avoid clustering. We discover a rich asymptotic behaviour including phenomena such as segregation of densities, regions of coexistence, travelling pulses -all of which are observed in biological contexts, cf. [42,45]. Existence of segregated stationary states under certain assumptions on the interaction potentials for small cross-diffusivity has been very recently obtained in [10]. Here we show that it is in fact possible to find explicit stationary states and travelling pulses for certain singular not necessarily decaying interaction potentials showing coexistence and segregation of densities. The rest of this paper is organised as follows: in Section 2 we discuss the basic properties of the system (1) in one dimension, in Section 3 we propose our numerical scheme which is used in Section 4 to explore the model and its long-time behaviour numerically. These insights are used to make reasonable assumptions on the support of the asymptotic solutions in order to derive analytic expressions for their shape and give a first classification of the zoology of the different stationary states. Finally we discuss in Section 5 how generic these phenomena are for different potentials and we draw the final conclusions of this work in Section 6.

2.
A non-local cross-diffusion model for two species. Throughout this paper we consider system (1) in one spatial dimension. Then the model reads for some initial data ρ(x, 0) = ρ 0 (x), and η(x, 0) = η 0 (x), and radially symmetric potentials W ij , for i, j = 1, 2. We can obtain some apriori estimates on solutions by using the following energy We note that for W ij ∈ W 2,∞ (R), along any solution (ρ, η) of system (2), there holds In the case of W ii = C ii x 2 /2 and W ij = ±C ij |x| for i = j with non-negative constants C ij , the estimate is also true, since and similarly for the terms in (2b), as long as ρ, η ∈ L ∞ (0, T ; L ∞ (R)). Thus the terms implying that ρ + η ∈ L 2 (0, T ; H 1 (R)). We deduce that the sum of both species remains continuous for almost all positive times -a property we will make use of later. Now, let us introduce our notion of steady states.
Proof. Clearly, the characterisation is sufficient, since the velocity field vanishes in each connected component of their supports if there exist constants c 1 , c 2 such that Eqs. (3) are satisfied. Conversely, if there holds we note that ρ, η, ∂ x (ρ + η) ∈ L 2 (R) by the definition of steady state, and therefore the right-hand sides are distributional derivatives of L 1 functions. By a well-known result (cf. e.g. [34], Lemma 1.2.1.), we deduce that there exist constants K 1 , K 2 ∈ R such that K 1 = ρ∂ x (W 11 ⋆ ρ + W 12 ⋆ η + ǫ(ρ + η)), Due to the integrabilty properties of the right-hand sides above, we infer that K 1 = K 2 = 0, and thus in the interior of any connected component of the supports of ρ and η, we obtain that there exist constants c 1 , c 2 ∈ R such that using the same argument as above.
Note that the assumption on the interiors of the supports of the species is purely technical and due to the regularity assumptions on our definition of stationary states. This avoids pathological cases such as functions supported on a fat Cantor set.
3. Numerical scheme. In order to solve system (2), we introduce a finite volume scheme based on Ref. [15]. The problem is posed on the domain Ω : and uniform size ∆x := x i+1/2 − x i−1/2 . Finally, the time interval [0, T ] is discretised by t n = n∆t, for n = 0, . . . , ⌈T /∆t⌉. We define the discretised initial data via We integrate system (2) over the test cell [t n , t n+1 ] × C i to obtain whereF n i+1/2 ,Ḡ n i+1/2 denote the flux on the boundary of cell C i , i.e.
Then the finite volume scheme for the cell averages ρ n i and η n i reads where we approximate the fluxes on the boundary, Eqs.(4), by the numerical fluxes using (·) + := max(·, 0) and (·) − := min(·, 0) to denote the positive part and the negative part, respectively. The velocity is discretised by centred differences: Here we have set where W l−k ij = W ij (x l − x k ), for i, j = 1, 2. This scheme has proven very robust for one species, and under a (more restrictive) CFL condition we can also prove the following result.
Proposition 3.1 (Non-negativity preservation). Consider system (2) with initial data ρ 0 , η 0 ≥ 0. Then for all n ∈ N the cell averages obtained by the finite volume method (5) satisfy ρ n i , η n i ≥ 0, granted that the following CFL condition is satisfied .
Proof. Let us assume that ρ n i , η n i ≥ 0, and we need to show that then ρ n+1 We can rearrange the terms so that Clearly, all terms in the second line are non-negative. The first line is non-negative if the CFL condition is satisfied. Application of the same procedure to η n+1 i yields the statement.
4. Numerical study. In this section we study system (2) numerically with emphasis on its long time behaviour. Throughout this chapter we use the self-interaction potentials and the cross-interaction potentials for the interspecific attractive-attractive and attractive-repulsive case, respectively. This choice of potentials allows us to compute steady states of system (2) explicitly. We find a wide range of different behaviours and properties, including segregation phenomena, for different cross-diffusivities and cross-interactions. Notice that the system is translationally invariant and therefore, if for symmetry considerations we can show that the centres of mass of both species in a stationary state are fixed and equal to some particular value, we can suppose that value to be zero without loss of generality. From numerical simulations we observe that steady states are compactly supported which motivates this ansatz when computing the profiles analytically. This is also due to the non-linear diffusion of porous medium type in the volume exclusion term. This chapter is subdivided into two sections addressing the mutually attractive case and the attractive-repulsive case, respectively. 4.1. Attractive-attractive case. Let us begin with the case of attractive interaction between both species, i.e. W 12 = W 21 = |x|. Upon exploring the system numerically, we find a vast variety of stationary patterns, including both symmetric and non-symmetric profiles whose occurrence and stability depends on the crossdiffusivity. In fact, the coefficient ǫ of the cross-diffusivity plays a crucial role in the bifurcations of these profiles, and will be discussed in the next section. Then, we study the system as the cross-diffusivity tends to zero and the stability of the steady states -a matter that seems closely intertwined with the bifurcations.

Steady states and behavioural bifurcation.
We begin by introducing the two types of symmetric steady states observed in the attractive-attractive case. Motivated by numerical simulations, we assume that the stationary distributions are compactly supported, i.e., is then only inhabited by the first species, but not η. Upon rearranging Eq. (3), we obtain The two non-local terms W 11 ⋆ ρ and W 12 ⋆ η can be computed individually. First the self-interaction terms becomes where are the mass and the first two moments of ρ, respectively. Then the cross-interaction term becomes where m 2 , M 2 denote the mass and the centre of mass of the second species. Due to symmetry and translational invariance of the solution, both M 1 and M 2 can be taken as zero without loss of generality. Upon substitution of the non-local terms in (7) and (8), Eq. (6) is simplified into , respectively. Using ρ(±c) = 0 at the boundary (where ρ + η vanishes identically), we get Finally, let us consider the interval [−b, b] where both species coexist. Again, ρ satisfies where the cross-interaction term W 12 ⋆ η can be further reduced, according to Notice that all terms on the right side are twice differentiable. Therefore from (10), ρ + η is twice differentiable in (−b, b), and upon differentiating Eq. (10) twice we obtain and similarly from the second equation in (3) The system of equations (11) and (12) can be solved by first introducing the decoupled system for u := ρ + η and v := ρ − η, giving by Thus, the solutions ρ and η are obtained as In fact, due to symmetry there holdsû 1 = 0, and Eqs. (13,14) can be simplified to Hence the symmetric steady states are determined uniquely by three parameters,û 2 , b and c, which are governed by algebraic equations. Since η is only supported on [−b, b], the condition for the total mass of η becomes From Eqs. (9,15), the condition for the total mass of ρ becomes Whenû 2 is eliminated, Eq. (16) provides a relation between b and c, i.e., Finally, consider the continuity of the sum of the densities ρ+η at Therefore b and c are in the zero locus of Eqs. (17,18) that are numerically solved, cf. Figure 1(a). Then the shape of the steady state is given by two parabola profiles on the parts only inhabited by the first species and cosine profiles where both species coexist: Figure 1(b) shows an excellent agreement between numerical and analytical steady states. Let us remark that Eq. (17) implies b = c in the case of m 1 = m 2 . As a consequence both species completely overlap and the profile is just that of a cosine, cf. Figure 2. Numerical simulations show that the Batman profiles are the only symmetric stationary distribution in a certain range of cross-diffusivities, namely (0, ǫ (1) ]. For ǫ ∈ (ǫ (1) , ǫ (2) ], a new family of profiles (called the second kind ) emerges coexisting with the Batman profiles in this range, cf. Figure 3. Finally, for ǫ > ǫ (2) only profiles of the second kind prevail. Since the steady states are a state of balance    between diffusion and attractive interactions, the second kind of profiles can be seen as states in which the attractive force is not strong enough to ensure the formation of a single group for η as observed in the Batman profiles. Similarly to the Batman profiles, we may determine parameters and their governing equations for profiles of second kind. In the symmetric case, using (3) the profiles are given by , and p is the fraction of mass in the corners of η, cf. Figure 3, (areas filled in red). Similarly, It is apparent that there are five unknowns b, c, d for the support, B for the amplitude in regions of coexistence, and p for the mass fraction. Correspondingly, we find four conditions in order to determine all parameters but p: for the mass of ρ and the continuity of the sum σ = ρ + η at x = c and x = b. Since p parameterises a family of solutions and describes both branches (as envelope) of the bifurcation diagram, cf. Figure 4, we are interested in finding the conditions leading to p min (ǫ), p max (ǫ) in the diagram, Figure 4. In order to determine the bifurcation diagram we run simulations with two different types of initial data -on the one hand we start the system with supp(η) ⊂ supp(ρ), on the other hand we initialise the system such that η is supported around ρ, cf. first row of Figure 5. The second row shows the stationary distribution asymptotically achieved with the respective initial data. We note that the mass fraction of η in the corners is different for both simulations albeit having used the same cross-diffusivity. The mass fraction in the left graph corresponds to p = p min and the mass fraction in the right graph to p = p max , respectively. Now we want to give conditions determining the envelopes p min (ǫ), p max (ǫ) of Figure 4. Let us impose non-negativity of η at x = b, i.e. η(b) ≥ 0. This is a reasonable assumption which is also reflected in the numerical simulations, cf. Figure 6(a). The figure shows steady states corresponding to the left initial data in Figure 5 as ǫ increases. While we observe a discontinuity of η at x = b for small ǫ, there is a critical value where η(b) = 0, for all ǫ > ǫ (1) . For the upper envelope we impose that the velocity field u 2 is non-negative at x = c since otherwise any small perturbation will render the stationary state unstable, i.e. mass would get transported into the interior, cf. Figure 6(b). These two conditions describe both envelopes in Figure 4.   In all four graphs the masses are m 1 = 0.1, m 2 = 0.6, and the cross-diffusivity is ǫ = 1.7. The first row depicts two different initial data -one (left) where η is included in ρ, and one (right) where η surrounds ρ. In the second row we present the corresponding steady states. Albeit having a similar make-up, they differ in their respective mass fraction of the corner, p. The left graph gives the minimal mass fraction p min while the right graph gives the maximal, pmax, respectively, cf. Vanishing diffusion regime. In this section we study the case of Batman profiles as ǫ → 0. Recall the two equations for b and c,  When ǫ is small, both b and c are O( √ ǫ), suggesting b = ǫ 1/2 b 0 + ǫ 1/2 b 1 + ǫb 2 + · · · , and c = ǫ 1/2 c 0 + ǫ 1/2 c 1 + ǫc 2 + · · · .
Upon substitution of the asymptotic expansions into Eq. (19) and (20), the leading order coefficients b 0 and c 0 satisfy Notice that both densities in the Batman profiles will converge to a Dirac Delta at zero with the respective masses while keeping their shape with this described asymptotic scaling for their supports. Asymmetric profiles. So far we only discussed symmetric steady states. However, there is an equally rich variety of non-symmetric stationary states, cf. Figure  7 and Figure 8. In Figure 7 we display the cases where the support only consists of     metric profiles for 0 < ǫ < ǫ (1) independent of the masses m 1 and m 2 . Only for larger cross-diffusivities, ǫ > ǫ (1) , asymmetric profiles can be observed. Moreover, there is a whole family of asymmetric profiles as can be seen in Figure 8. This is similar to the case of symmetric stationary states, parameterised by the mass fraction p.
Stability of steady states and symmetrising effect. Let us now discuss the numerical stability of the symmetric steady states. Here the bifurcation point ǫ (1) plays an important role, for the system exhibits a symmetrising effect whenever the crossdiffusivity lies below the critical one, in the sense that there is only one symmetric steady state attracting any initial data.
We fixed ǫ ∈ (0, ǫ (1) ) and chose ρ 0 = 2m 1 ½ [−0.5,0] and η 0 = 2m 2 ½ [0,0.5] for all combinations of masses of the form (m 1 , m 2 ) = 0.1 · (i, j) for i, j = 1 . . . 10. In all cases we observe that there is only one attractor, namely the Batman profile of the form given in Figure 1(b) and Figure 2(b) in the case m 1 = m 2 , respectively. For ǫ > ǫ (1) the system is not symmetrising anymore and small perturbations lead to different stationary states. This can be seen if p is varied in [p min , p max ], for it leads to different states. A similar argument holds for the asymmetric states, by shifting mass from one corner into the other, cf. Figure 8.

4.2.
Attractive-repulsive case. In this section we present the attractive -repulsive case, i.e. W 12 = |x| = −W 21 . Then the steady states have segregated densities, as asserted by the following proposition. Proof. Suppose the interior of a connected componente of supp(ρ) ∩ supp(η) is not empty. We know that both species satisfies Eqs. (3) in that connected component: Similar arguments as above imply that the interaction terms are twice differentiable in this interval, thus we differentiate twice and get 0 = m 1 + 2η + ǫ(ρ + η) ′′ and 0 = m 2 − 2ρ + ǫ(ρ + η) ′′ .

Steady states.
This section is dedicated to studying the steady states of the system with attractive-repulsive cross-interactions. Due to numerical simulations and the previous proposition we make the following assumption on the support where a < b ≤ −c < c ≤ d < e are some real numbers. Using Eqs. (3) we proceed similar as above, cf. Eqs. (7,8), to obtain for shape of the second species on the left part of the support and for the right part, respectively. Similar as above, we can see that the interaction terms are twice differentiable, therefore differentiating Eq. (3) in the support of ρ twice yields 0 = m 1 + ǫ(ρ) ′′ , and thus and analogously Concerning the first species, the parameters β, γ are determined by the continuity condition Eq. (22) and we obtain We can see that there are six unknowns, namely a, b, c, d, e, and M 2 with a total of five conditions: by imposing half of the mass of η to each side of ρ.

Case of strict segregation. Let us start by discussing the case
Then the condition on the mass yields We can solve η L (b) = 0 for a, Since half of the mass is located to the left of the first species, we get where we used Eq. (23a). Similarly, we solve η R (d) = 0 for e to obtain Using this expression we compute So we have determined c, b, d depending only on the masses and the second order moments of the second species, M 2 . We can substitute the values into Eqs. (23a, 23c) to determine a and e.
Critical ǫ and maximal M 2 . We are interested in a condition determining as to when segregation of species occurs. In fact there is a critical value of the crossdiffusivity, ǫ c , such that there only exist adjacent steady states for ǫ > ǫ c . For 0 < ǫ < ǫ c strictly segregated steady states occur if |M 2 | < M 2,max , where M 2 = 0 corresponds to the symmetric case. Figure 9 displays this behaviour. Let us derive an expression for ǫ c and M 2,max . For a fixed ǫ we may compute M 2,max . We begin with the case c = −b. We can solve for the critical M 2 , i.e.
Similarly, we can solve equation c = d for M 2 , which gives    (c) ǫ = 1/2. If ǫ = ǫ c both species touch at the points {−c, c} or are partially adjacent. If ǫ < ǫ c but we choose M 2 outside of the aforementioned range we observe steady states consisting of (partially) adjacent bumps. Figure 10 displays the steady states in the symmetric case, i.e. M 2 = 0, for attractiverepulsive cross-interactions. We observe a transition of behaviour for different values of ǫ, ranging from strictly segregated states to completely adjacent states. The numerical results agree perfectly with the results obtained analytically.
Vanishing diffusion regime. As we have seen in Figure 9, there is an ǫ c such that the steady states parameterised by M 2 ∈ [−M 2,max , M 2,max ] are segregated. In this section we consider the case of vanishing cross-diffusion. We can assume that ǫ < ǫ c and This is indeed a measure solution of system. To see this let us consideṙ Since we are looking for a steady state we observė We assume without loss of generality that Fixing X = 0 we get Y 2 = Y 1 + 2m 1 /m 2 and Y 1 ∈ [−2 m1 m2 , 0]. This is exactly the solution of the system as ǫ → 0, cf. Eq. (26).
Stability of steady states. Here we want to discuss the stability of the stationary states of the attractive-repulsive system. In general, the stationary states are not stable as small perturbations may lead to a completely different stationary state. It becomes clear in Figure 9, that perturbing η by shifting it to either side leads to a completely different stationary state. Although this is an arbitrarily small perturbation in any L p -norm, the translated profile is another stationary state. This is why these profiles are not stable. The same argument holds for symmetric stationary states. However, they are stable under symmetric perturbations since any symmetric initial data is attracted by the symmetric profile. Characterising fully the basin of attraction for each stationary state seems difficult. For perturbations shifting mass from η L to η R (or vice versa) there is no stationary state but the profile is then attracted by a travelling pulse solution.

Travelling pulses.
In addition to the convergence to steady states we observe travelling pulse solutions in the case of attractive-repulsive cross-interactions. There are two types of travelling pulses -those consisting of two bumps and those consisting of three. In our numerical study we do not observe more than three bumps, even in the case of exponentially decaying potentials. There are however metastable states where more bumps exist but after a sufficiently long time the collapse into two or three.
Two pulses. In order to compute these profiles, we assume [−a, a] denotes the initial support of u = η(0) and therefore [−a − x 0 , a − x 0 ] the initial support of ρ(0). We transform the system into co-moving coordinates, z = x − vt, and obtain the following conditions for the pulse profiles similarly to Eqs. (3). A computation similar to Eqs. (7,8), leads to the explicit form of the pulse on [−a, a] for some constantc 1 . Since u(z) is a parabola with roots ±a, u is symmetric. As a consequence we obtain M = v − m. By definition of M = zu(z)dz = 0, whence v = m. Hence the shape is given by Then the following consideration determines the boundary of the support, a, If there were travelling pulse solutions of this form they would satisfy the same equations as above. Then, The continuity of the sum suggests that u 1 (0) = u 2 (0) implies m = v. But then We solve this expression for a > 0 and find a = 3 √ 12ǫ. A comparison of the support of the adjacent solutions and the support of segregated solutions, cf. Eq. (28), shows that the adjacent solutions in fact only touch. Figure 11 shows the formation of two travelling pulses. We start with two indicator functions as initial data and let the system evolve. At about time t ≈ 2 we observe a fully established pulse profile. We let the system evolve further and compare the We transform to co-moving coordinates, z = x−vt, and obtain the following conditions for the profile whence we obtain Here Similarly, the profiles of the second species are given by and Again, we use the fact that the sum of both densities has to be continuous, i.e.
as well as Eqs. (30) hold. We consider the case of strictly segregated solutions first, i.e. b < −c, and c < d. Since then ρ(±c) = 0, we may deduce from Eq. (29) that v = m R − m L for the speed of propagation and for the shape of the first species (c is determined by the mass condition, Eq. (31)). Furthermore we obtain in terms of b. Similarly, we can get an expression for e in terms of d, i.e.
Using the expression for a, we obtain Note that Eqs. (32) completely determine the support and the profiles of the pulses. Figure 12 shows the formation of a triple pulse solution. We choose characteristic functions as initial data (dotted). The mass on the left is m L = 1/3 and, respectively, m R = 2/3 on the right. After some time the pulse profile is established. We compare the system (blue and red) at time t = 9 and time t = 24 with the analytical expression derived above (black). The figure displays a great agreement between our numerical result and the analytical. Once the profile is fully established it moves to the right at a constant speed. The numerical velocity is given by ∆x/∆t = 5/15 = 1/3. This is in perfect agreement with the analytically obtained results, i.e. v = m R − m L = 2/3 − 1/3 = 1/3. At this stage, let us draw our attention to two special cases. Remark 3 (Maximal M 2 ). Let us get back to the general case. We study the interval of M 2 . Assuming ǫ fixed, b = −c yields On the other hand, c = d gives where v = m R − m L , as above. It is worthwhile noting that in the case m L = m R both M 2,max and M 2,min coincide with Eqs. (24,25) for the stationary state.
Parallel to the consideration for (partially) adjacent steady states of the attractiverepulsive system we also find the existence of adjacent travelling pulse solutions.

5.
Generality. This section is dedicated to the study of more general or realistic potentials to understand whether the behaviours observed above are specific to our interaction potentials. Different cross-interaction and self-interaction potentials will be investigated. Even though analytic expressions for the steady states and travelling pulses seem no longer avaiablable, the behaviours are indeed generic and, in fact, even richer than the above particular model.

Different cross-interactions.
Let us begin by considering different crossinteraction potentials. We regard two types of potentials -power-laws and Morse-like potentials decaying at infinity, i.e. where p ∈ {1/2, 1, 3/2}. This choice of potentials is motivated as they are similar to the Newtonian cross-interaction. In both cases, we observe a very similar behaviour both in the mutually attractive case and the attractive-repulsive case, respectively. Figure 13 displays the Batman profile  for different cross-interaction potentials. In all simulations the same initial data, mass, and cross-diffusivity were used. Each steady state features the salient characteristics observed in the case W cr = |x|, i.e. a region of coexistence surrounded by regions inhabited by only one species. From the steady states we can also infer another information, namely, second type profiles exist and the point of bifurcation depends on the potential, for only Figure 13(d) exhibits a profile of second type. Similarly, we observe a symmetrising effect for small cross-diffusivities and asymmetric profiles.

Different self-interactions.
Here, we keep the cross-interaction potentials fixed as W 12 = |x| = ±W 21 and consider different self-interaction potentials of the form W (x) = |x| p /p, for p ∈ {3/2, 2, 4}. In each case we observe a very similar behaviour. We obtain the same variety including both Batman profiles and the profiles of second type. Again we observe that the system is symmetrising, however for a different ǫ (1) . In the attractive-repulsive case as well we observe the characteristic profiles and the formation of pulses.
6. Conclusions. In this paper we introduced a system of two interacting species with cross-diffusion. We used a positivity-preserving finite-volume scheme in order to study the system numerically. For a specific choice of potentials, the steady states can be constructed with parameters governed by algebraic equations. These numerically simulated and the analytically constructed stationary states and travelling pulses were found to agree with each other. Using the same scheme the model was explored for related potentials and the behaviours observed for the specific potentials turned out to be generic, when the cross-interaction potentials or the self-interaction potentials were exchanged. While this paper gives a first insight as to what qualitative properties can be expected from models taking the general form (1), there is still a lot of analytical work to be done. First and foremost, it is still an open problem to show existence of solutions to the systems. The formal gradient flow structure is lost when the crossinteraction potentials W 12 and W 21 are not proportional to each other, and the main problem is to find the right estimates for individual species since we only control the gradient of the sum of the densities.