Convergence of a finite volume scheme for a system of interacting species with cross-diffusion

In this work we present the convergence of a positivity preserving semi-discrete finite volume scheme for a coupled system of two non-local partial differential equations with cross-diffusion. The key to proving the convergence result is to establish positivity in order to obtain a discrete energy estimate to obtain compactness. We numerically observe the convergence to reference solutions with a first order accuracy in space. Moreover we recover segregated stationary states in spite of the regularising effect of the self-diffusion. However, if the self-diffusion or the cross-diffusion is strong enough, mixing occurs while both densities remain continuous.

governing the evolution of two species ρ and η on an interval (a, b) ⊂ R for t ∈ [0, T ).
The system is equipped with nonnegative initial data ρ 0 , η 0 ∈ L 1 + (a, b) ∩ L ∞ + (a, b). We denote by m 1 the mass of ρ 0 and by m 2 the mass of η 0 , respectively, On the boundary x = a and b, we prescribe no-flux boundary conditions such that the total mass of each species is conserved with respect to time t ≥ 0. While the self-interaction potentials W 11 , W 22 ∈ C 2 b (R) model the interactions among individuals of the same species (also referred to as intraspecific interactions), the crossinteraction potentials W 12 , W 21 ∈ C 2 b (R) encode the interactions between individuals belonging to different species, i.e. interspecific interactions. Here C 2 b (R) denotes the set of twice continuously differentiable functions on R with bounded derivatives. Notice that the convolutions W i j ψ, with ψ a density function defined on [a, b], are defined by extending the density ψ by zero outside the interval [a, b]. The two positive parameters , ν > 0 determine the strengths of the self-diffusion and the cross-diffusion of both species, respectively. Nonlinear diffusion, be it self-diffusion or cross-diffusion, is biologically relevant. As a matter of fact, around the second half of the twentieth century biologists found that the dispersal rate of certain insects depends on the density itself, leading to the nonlinear diffusion terms we incorporated in the model, cf. [13,32,34,36,37]. At the same time we would like to stress that the self-diffusion terms are relevant for the convergence analysis below.
It is the interplay between the non-local interactions of both species and their individual and joint size-exclusion, modelled by the non-linear diffusion [4][5][6]10,12,40], that leads to a large variety of behaviours including complete phase separation or mixing of both densities in both stationary configurations and travelling pulses [11,19].
While their single species counterparts have been studied quite intensively [15,21,35,41] and references therein, related two-species models like the system of our interest, Eq. (1), have only recently gained considerable attention [11,17,19,24,25]. One of the most striking phenomena of these interaction models with cross-diffusion is the possibility of phase separation. Since the seminal papers [5,33] established segregation effects for the first time for the purely diffusive system corresponding to (1) for W i j ≡ 0, i, j ∈ {1, 2} and = 0, many generalisations were presented. This includes reaction-(cross-)diffusion systems [3,7,17] and references therein, and by adding non-local interactions [2,11,19,24] and references therein. Ref. [24] have established the existence of weak solutions to a class of non-local systems under a strong coercivity assumption on the cross-diffusion also satisfied by system (1).
Typical applications of these non-local models comprise many biological contexts such cell-cell adhesion [16,38,39], for instance, as well as tumour models [26,31], but also the formation of the characteristic stripe patterns of zebrafish can be modelled by these non-local models [42]. Systems of this kind are truly ubiquitous in nature and we remark that 'species' may not only refer to biological species but also to a much wider class of (possibly inanimate) agents such as planets, physical or chemical particles, just to name a few.
Since system (1) is in conservative form a finite volume scheme is a natural choice as a numerical method. This is owing to the fact that, by construction, finite volume schemes are locally conservative: due to the divergence theorem, the change in density on a test cell has to equal the sum of the in-flux and the out-flux of the same cell. There is a huge literature on finite volume schemes, first and foremost [28]. Therein, the authors give a detailed description of the construction of such methods and address convergence issues. Schemes similar to the one proposed in Sect. 2 have been studied in [9] in the case of nonlinear degenerate diffusion equations in any dimension. A similar scheme for a system of two coupled PDEs was proposed in [22]. Later, the authors in [14] generalised the scheme proposed in [9] including both local and non-local drifts. The scheme was then extended to two species in [19]. All the aforementioned schemes have in common that they preserve nonnegativity-a property that is also crucial for our analysis.
Before we define the finite volume scheme we shall present a formal energy estimate for the continuous system. The main difficulty in this paper is to establish positivity and reproducing the continuous energy estimate at the discrete level. The remainder of the introduction is dedicated to presenting the aforementioned energy estimate. Let us consider where the second equality holds due to the no-flux boundary conditions. Upon rearranging we get whence, upon adding both, we obtain where σ = ρ + η and denote the advective parts associated to ρ and η, respectively. The advective parts can be controlled by using the weighted Young inequality to get for some α > 0. In choosing 0 < α < we obtain where C ρ = W 11 ρ + W 12 η 2 L 2 and C η = W 22 η + W 21 ρ 2 L 2 . From the last line, Eq.
(2), we may deduce bounds on the gradient of each species as well as on their sum. As mentioned above the crucial ingredient for this estimate is the positivity of solutions.
The rest of this paper is organised as follows. In the subsequent section we present a semi-discrete finite volume approximation of system (1) and we present the main result, Theorem 2.4. Section 3 is dedicated to establishing positivity and to the derivation of a priori estimates. In Sect. 4 we obtain compactness, pass to the limit, and identify the limiting functions as weak solutions to system (1). We conclude the paper with a numerical exploration in Sect. 6. We study the numerical order of accuracy and discuss stationary states and phase segregation phenomena.

Numerical scheme and main result
In this section we introduce the semi-discrete finite volume scheme for system (1). To begin with, let us introduce our notion of weak solutions. and respectively, for any ϕ ∈ C ∞ c ([0, T ) × (a, b); R). Here we have set V k = −W k 1 ρ − W k 2 η, for k ∈ {1, 2}, and σ = ρ + η, as above.
Notice that the existence of weak solutions to system (1) will follow directly from the convergence of the numerical solution. Indeed, our analysis relies on a compactness argument which does not suppose a priori existence of solution to system (1).
To this end we first define the following space discretisation of the domain.

Definition 2.2 (Space discretisation)
To discretise space, we introduce the mesh where the control volumes are given by We assume that the measure of the control volumes are given by |C i | = Note that x 1/2 = a, and x N +1/2 = b, cf. (Fig. 1).
We also define x i = (x i+1/2 + x i−1/2 )/2 the centre of cell C i and set x i+1/2 = x i+1 − x i for i = 1, . . . , N − 1. We assume that the mesh is regular in the sense that there exists ξ ∈ (0, 1) such that for h : and, as a consequence, ξ h ≤ x i+1/2 ≤ h, as well.
On this mesh we shall now define the semi-discrete finite volume approximation of system (1). The discretised initial data are given by the cell averages of the continuous initial data, i.e.
for all i ∈ I . Throughout, we write ρ i (resp. η i ) to denote the approximations of the two densities on the i-th finite volume cell, C i . Next, we introduce the discrete versions of the cross-diffusion and the interaction terms. We set where for k, l = 1, 2, and for the cross-diffusion term, respectively. Then the scheme reads for i ∈ I . Here the numerical fluxes are given by for i = 1, . . . , N − 1, with the numerical no-flux boundary condition where we introduced the discrete gradient du i+1/2 as As usual, we use (z) ± to denote the positive (resp. negative) part of z, i.e.
At this stage, the numerical flux (9b) may look strange since • the cross-diffusion term is approximated as a convective term using that with σ = ρ + η and ∂σ ∂ x is considered as a velocity field. This treatment has already been used in [9] and allows to preserve the positivity of both discrete densities (ρ, η) (see Lemma 3.1), which is crucial for the convergence analysis.
• In this new formulation, the velocity field is split in two parts both treated by an upwind scheme. One part comes from the cross-diffusion part, and the second one comes from the non-local interaction fields. This splitting is crucial to recovering a consistent dissipative term for the discrete energy estimate corresponding to Eq. (2).

Definition 2.3 (Piecewise constant approximation)
For a given mesh T h we define the approximate solution to system (1) by for all (t, x) ∈ [0, T ] × C i , with i = 1, . . . , N . Moreover, we define the following approximations of the gradients Furthermore, in order to define dρ h and dη h on the whole interval (a, b) we set them to zero on (a, x 1 ) and (x N , b).
Notice that the discrete gradients (dρ h , dη h ) are piecewise constant just like (ρ h , η h ) however not on the same partition of the interval (a, b). In a similar fashion we define the piecewise constant interpolation of the discrete advection fields, i.e., . . , N − 1, and zero at the boundary.
We have set out all definitions necessary to formulate the convergence of the numerical scheme (9).

Theorem 2.4 (Convergence to a weak solution)
is a weak solution as in Definition 2.1. Furthermore we have ρ, η ∈ L 2 (0, T ; H 1 (a, b)); (iii) as a consequence system (1) has a weak solution.

A priori estimates
This section is dedicated to deriving a priori estimates for our system. In order to do so we require the positivity of approximate solutions and their conservation of mass, respectively. The following lemma guarantees these properties.

Lemma 3.1 (Existence of nonnegative solutions and conservation of mass)
Assume that the initial data (ρ 0 , η 0 ) are non-negative. Then there exists a unique nonnegative approximate solution (ρ h , η h ) h>0 to the scheme (9a)-(9c). Furthermore, the finite volume scheme conserves the initial mass of both densities.
Proof On the one hand we notice that the right-hand side of (9a)-(9b) is locally Lipschitz with respect to (ρ i , η i ) 1≤i≤N . Hence, we may apply the Cauchy-Lipschitz theorem to obtain a unique continuously differentiable local-in-time solution.
On the other hand to prove that this solution is global in time, we show the nonnegativity of the solution together with the conservation of mass and argue by contradiction.
On a given mesh, let some initial data, ρ i (0), η i (0) ≥ 0, be given for i = 1, . . . N . We rewrite the scheme in the following way. where Then let t ≥ 0 be the maximal time for all densities to remain nonnegative, i.e.
If t < ∞, then there exists a nonincreasing sequence (t k ) k∈N such that t k > t , t k → t as k → ∞ and there exists i k ∈ {1, . . . , N } verifying Since the index i k takes a finite number of integer values, we can extract a nonincreasing subsequence of (t k ) k∈N still labeled in the same manner such that there exists an index j 0 ∈ {1, . . . , N } and where t k → t , as k goes to infinity.
By the above computation, Eq. (10), we see that, if ρ j 0 +1 (t ) > 0 or respectively hence there exists τ > 0 such that for any t ∈ [t , t + τ ), we have ρ j 0 (t) > ρ j 0 (t ) = 0, which cannot occur since ρ j 0 is continuous and for t k > t , ρ j 0 (t k ) < 0 for any k ∈ N with t k → t when k goes to infinity.
= 0 then by uniqueness of the solution, we have that ρ j 0 ≡ 0 for t ≥ t , which contradicts again that ρ j 0 (t k ) < 0 for any k ∈ N large enough.
Finally we get the conservation of mass, by the no-flux condition. Analogously, the second species remains nonnegative and its mass is conserved as well. As a consequence of the control of the L 1 -norm of (ρ h , η h ) we can extend the local solution to a global, nonnegative solution.
Now, we are ready to study the evolution of the energy of the system on the semidiscrete level. The remaining part of this section is dedicated to proving the following lemma-an estimate similar to (2) for the semi-discrete scheme (9).

Lemma 3.2 (Energy control) Consider a solution of the semi-discrete scheme
where the constant C > 0 is given by Proof Upon using the scheme, Eq. (9a), we get due to the conservation of mass, ensured by Eq. (9c). By discrete integration by parts and the no-flux condition, Eq. (9c), we obtain where, in the last equality, we substituted the definition of the numerical flux, Eq. (9b). Let us defineρ for i ∈ {1, . . . , N − 1}, and note that thenρ i+1/2 ∈ [ρ i , ρ i+1 ] by concavity of the log.
Here, and throughout, we use the shorthand notation Reordering the terms, we obtain Thus, usingρ i+1/2 ∈ [ρ i , ρ i+1 ] and the monotonicity of log, we note that This is easy to see, for, if ρ i = ρ i+1 , we observe dlog ρ i+1/2 = 0 and Eqs. (14) hold with equality. In the case of ρ i < ρ i+1 we observe whence we infer the inequality. The same argument can be applied in order to obtain the second line of Eq. (14). Thus we may infer from Eq. (13) that Note that the definition ofρ i+1/2 in Eq. (12), is consistent with the case ρ i = ρ i+1 and there holdsρ Furthermore, we notice that where we employed Eq. (12). Hence we have A similar computation can be applied to the second species, which yields Upon adding up Eqs. (15) and (16), we obtain where R h is given by Finally we notice that for any α > 0, by Young's inequality. Observing, that for k = 1, 2, and by conservation of positivity and mass, it gives for k = 1, 2, Thus Eq. (18) becomes where C 2α is given in Eq. (11). Finally, substituting the latter estimate into Eq. (17), we obtain, upon using Eq. (8), for any solution (ρ i ) i∈I , (η i ) i∈I of the semi-discrete scheme (9). Hence choosing α = /2 concludes the proof.

Corollary 3.3 (A priori bounds)
Let (ρ i ) i∈I , (η i ) i∈I be solutions of the semi-discrete scheme (9). Then there exists a constant C > 0 such that Hence the discrete version of the classical entropy functional is bounded from below. Therefore, we integrate the inequality of Lemma 3.2 in time and get which proves the statement.
Thanks to a classical discrete Poincaré inequality [28, Lemma 3.7] and [8], we get uniform L 2 -estimates on the discrete approximation (ρ h , η h ) h>0 . Lemma 3.4 Let (ρ i ) i∈I , (η i ) i∈I be the numerical solutions obtained from scheme (9). Then there holds

Proof of Theorem 2.4
This section is dedicated to proving compactness of both species, the fluxes, and the regularising porous-medium type diffusion. Upon establishing the compactness result we identify the limits as weak solutions in the sense of Definition 2.1.
First by application of Lemma 3.1, we get existence and uniqueness of a nonnegative approximate solution (ρ h , η h ) to (9a)-(9b). Hence the first item of Theorem 2.4 is proven. Now let us investigate the asymptotic h → 0.

Strong compactness of approximate solutions
We shall now make use of the above estimates in order to obtain strong compactness of both species, (ρ h , η h ) in L 2 (Q T ). h>0 be the approximation to system (1) obtained by the semi-discrete scheme (9). Then there exist functions ρ, η ∈ L 2 (Q T ) such that Proof We invoke the compactness criterion by Aubin and Lions [23]. Accordingly, a set P ⊂ L 2 (0, T ; B) is relatively compact if P is bounded in L 2 (0, T ; X ) and the set of derivatives {∂ t ρ ρ ∈ P} is bounded in a third space L 1 (0, T ; Y ), whenever the involved Banach spaces satisfy X → → B → Y , i.e. the first embedding is compact and the second one continuous. For our purpose we choose X : = BV (a, b), B := L 2 (a, b), and Y := H −2 (a, b). The first embedding is indeed compact, e.g. Ref.
[1, Theorem 10.1.4] and the second one is continuous.
In the second step we show the time derivatives are bounded in L 1 (0, T ; H −2 (a, b)). To this end, let ϕ ∈ C ∞ c ((a, b)). Throughout, we write ·, · for ·, · H −2 ,H 2 for the dual pairing. Making use of the scheme, there holds having used the scheme, Eq. (9a). Next we set perform a discrete integration by parts and use the no-flux boundary conditions, Eq. (9c), to obtain Using the definition of the numerical flux, Eq. (9b), we get Let us begin with the self-diffusion part. Using the Cauchy-Schwarz inequality, we estimate the discrete gradient and ρ itself by Corollary 3.3 and Lemma 3.4 where we used that ϕ ∈ H 1 ⊂ L ∞ and the regularity of the mesh, ξ > 0, cf. Eq. (4).
Next, we address the cross-diffusion and non-local interactions terms using the same argument. For instance for the cross-diffusive part, we have where we used Corollary 3.3 and Lemma 3.4 again. The non-local interaction term is estimated in the same way, thus there holds a, b)), which concludes the proof.
From the latter result we can prove the convergence of the discrete advection field dV 1,h and dV 2,h defined as in Definition 2.3.
We define V k,h and V k as On the one hand from the strong convergence of (ρ h , η h ) to (ρ, η) in L 2 (0, T ; L 2 (a, b)) and the convolution product's properties, we obtain On the other hand, we have for any hence there exists a constant C > 0 such that Integrating over x ∈ [x i , x i+1 ) and summing over i ∈ {1, . . . , N − 1}, we get that Notice that ( (20) and (21) we get that dV k,h − V k L 2 (0,T ;L 2 (a,b)) goes to zero as h tends to zero.

Weak compactness for the discrete gradients
In the previous section we have established the strong L 2 -convergence of both species, (ρ h ) h>0 and (η h ) h>0 . However, in order to be able to pass to the limit in the crossdiffusion term ρ h (dρ h + dη h ) we need to establish weak convergence in the discrete gradients in L 2 . This is done in the following proposition.

Passing to the limit
We have now garnered all information necessary to prove Theorem 2.4. For brevity we shall only show the convergence result for ρ, as it follows for η similarly, using the same arguments. a, b)) be a test function. We introduce the following notations: On the other hand, we set and multiply the scheme, Eq. (9a), by the test function and integrate in time and space to get where When h tends to zero and from the strong convergence of (ρ h , η h ) h>0 to (ρ, η) in L 2 (Q T ), the strong convergence of (dV k,h ) h>0 to V k in L 2 (Q T ) and the weak convergence of the discrete gradient (dU h ) h>0 to − ∂σ ∂ x in L 2 (Q T ), it is easy to see that when h → 0. Therefore it suffices to prove that ε(h) → 0, as h goes to zero, which will be achieved by proving that

The self-diffusion part D h − D 1,h
On the one hand, after a simple integration we get Hence, we have and observing that we obtain, in conjunction with the Cauchy-Schwarz inequality and the a priori bounds established in Corollary 3.3, and Lemma 3.4, that in the virtue of the estimate Eq. (19).

The cross-diffusion part
Let us now treat the cross-diffusion part. This term is more complicated since it involves the piecewise constant functions ρ h and dU h , which are not defined on the same mesh. Thus, on the one hand we reformulate the discrete cross-diffusion term C 1,h as C 1,h = C 10,h + C 11,h with where a direct computation and the application of Corollary 3.3 and Lemma 3.4 yield On the other hand, the term C h can be rewritten as the term C h can be decomposed as C h = C 00,h + C 01,h with Similarly to (25), the first term C 00,h can be estimated as whereas the second term C 01,h is compared to C 11,h Using a second order Taylor expansion of ϕ at x i and x i+1 , it yields that hence we get from Corollary 3.3 and Lemma 3.4 that Gathering Eqs. (25), (26), and (27), we finally obtain that The advective part The evaluation of A h − A 1,h is along the same lines of the cross-diffusion terms C h − C 1,h since the latter is treated as an advective term. Hence, thanks to Lemma 4.2, we get that Finally by definition of ε(h) and using Eq. (23) together with Eqs. (24), (28), and (29), we obtain that is, ε(h) → 0, when h → 0, which proves that (ρ, η) is a weak solution to Eq. (1). This proves the second item of Theorem 2.4. Finally the last item concerning the existence of solutions to (1) is a direct consequence of the convergence.

A fully discrete implicit scheme
In this section we shall comment on a discrete-in-time version of the semi-discrete scheme (9). To this end we replace the time derivative in Eq. (9a) by simple forward differences and obtain the following implicit and fully-discrete scheme where t > 0. System (30) gives rise to two approximating sequences (ρ n i ) 1≤i≤N and (η n i ) 1≤i≤N , for 0 ≤ n ≤ M where M := T / t and the discrete time instances t n := n t; cf. Theorem 5.1. Here the numerical fluxes are given by for i = 1, . . . , N − 1, with the numerical no-flux boundary condition for n = 0, . . . , M. Recall that Similarly to Definition 2.3 we define the piecewise constant interpolation by . . , N , and n = 0, . . . , M. Moreover, we define the discrete approximation of the spatial gradients as . . , N −1 and n = 0, . . . , M. As above, we set the discrete gradients to zero on (a, x 1 ) and (x N , b). Furthermore, we define the discrete time derivative as for (t, x) ∈ [t n , t n+1 ) × C i , for i = 1, . . . , N and n = 0, . . . , M − 1.

Theorem 5.1 (Existence and uniqueness result)
Let ρ 0 i , η 0 i be nonnegative initial data with mass m 1 and m 2 , respectively, and assume the following time step restriction condition Then there exists a unique nonnegative solution (ρ n i , η n i ) to scheme (30), (30a) and (30b).
Proof We show existence first and prove uniqueness later. Suppose we are given (ρ n i ) 1≤i≤N and (η n i ) 1≤i≤N from some previous iteration. In order to construct the next iteration we shall employ Brouwer's fixed point theorem. It is easy to verify that the set is a convex and compact subset of R N × R N . Hence, we define the fixed point operator S : X → R N ×R N by setting (ρ , η ) = S(ρ, η) where (ρ , η ) are implicitly given as for any (ρ, η) ∈ X , where F and G denote the numerical fluxes for i = 1, . . . , N − 1 where dU is computed from (ρ, η), with the numerical no-flux boundary condition F 1/2 = F N +1/2 = 0, and G 1/2 = G N +1/2 = 0.
Notice that for any given (ρ, η) ∈ X , the viscosity terms in front of the discrete gradients involved in the definition of the fluxes F i+1/2 and G i+1/2 are indeed nonnegative, hence the couple (ρ , η ) is well defined since it corresponds to the unique solution of a classical fully implicit scheme in time with an upwind discretisation for the convective terms and a centred approximation for diffusive terms [29]. Moreover, since ρ n and η n are nonnegative and using the monotonicity of the numerical flux with respect to (ρ , η ), we prove that both densities ρ and η are also nonnegative. Furthermore, using the nonnegativity and the no-flux conditions, we get which yields that S(X ) ⊂ X . Finally, S is continuous as the composition of continuous functions. Thus, we may apply Brouwer's fixed point theorem to infer the existence of a fixed point, (ρ n+1 , η n+1 ). It now remains to show uniqueness of the fixed point.
To treat in a systematic way the boundary conditions and simplify the presentation, we define ghost values for (ρ, η) by setting for α ∈ {ρ, η}, and k ∈ {1, 2}, Then we consider two solutions (ρ,η) and (ρ, η) to (30). Setting h(x) := x 2 /2, s := ρ −ρ and r := η −η, we get after substituting the two solutions to (30), for i = 1, . . . , N , and a similar relation for (r i ) 1≤i≤N . Applying a Taylor expansion on h(ρ) = h(ρ) + h (ρ) s, withρ a convex combination of ρ andρ, we may write where A i , B i−1 , C i+1 are nonnegative coefficients given by Now, we multiply equation (34) by sign(s i ) and sum over i = 1, . . . , N , hence using that x → x ± is Lipschitz continuous and observing that and in a similar way, Gathering these latter inequalities and from (4), it gives that Finally, from the nonnegativity and the preservation of mass, we have hence under the condition (31), we conclude that s = r = 0 and the uniqueness follows.
It is worth to mention here that the condition (31) is not optimal since we only use the discrete L 1 -estimate on ρ and η to control the discrete gradient and the L ∞ -norm.

Sketch of the proof of Theorem 5.2
It is easily observed that the total mass is conserved due to the discrete no-flux boundary conditions, cf. (30b). Together with the nonnegativity we were able to prove a semi-discrete version of the energy estimate Eq. (2) which is at the heart of the convergence result. Similarly as above, we are able to prove a fully discrete version of the energy estimate which then reads with C > 0 as in Corollary 3.3. The inequality follows from the convexity of x(log x − 1) since Multiplying this expression by x i / t and summing over i = 1, . . . N and n = 0, . . . , M we obtain The right-hand side is then substituted by the scheme and simplified along the lines of the proof of Lemma 3.2. Since the computations are exactly the same we omit them here for brevity and only note that it is important to set to obtain the right sign in the numerical artefacts in Eq. (14), which now read Following the lines of the proof Lemma 3.2 and Corollary 3.3 we obtain the fully discrete a priori bounds for some constant C > 0, where we used 'd x ' to denote the discrete spatial gradient as before. Again, an application of Aubin-Lions Lemma provides relative compactness in the space L 2 ((0, T ) × (a, b)) of the piecewise constant interpolations ρ h , η h . As above, the discrete gradients are uniformly bounded in L 2 ((0, T ) × (a, b)) and their weak convergence is a consequence of the Banach-Alaoglu Theorem. Identifying the limits as a weak solution to system (1) is shown in the same way as above and we leave it as an exercise for the reader.
Remark 1 (Explicit scheme) We do not consider an explicit scheme here since its analysis is much more complicated due to the lack of uniform estimates. Indeed, an explicit scheme requires a CFL condition on the time step which appears in the stability analysis or the energy estimate provided in Lemma 3.2. In our case, it leads to where ρ n+1/2 belongs to the interval [ρ n , ρ n+1 ] or [ρ n+1 , ρ n ]. The control of this reminder term would require some lower bounds estimates on the density ρ (see for instance [30]).

Numerical examples and validation
In this section we perform some numerical simulations of system (1) using our scheme, Eqs. (9). In Sect. 6.1 we test our scheme by computing the error between the numerical simulation and a benchmark solution on a finer grid. Furthermore we determine the numerical convergence order. In Sect. 6.2 we compute the numerical stationary states of system (1) and discuss the implication of different cross-diffusivities and the selfdiffusivities, respectively. Let us note here, that in the case of no regularising porous-medium diffusion, i.e. = 0, and certain singular potentials the stationary states of system (1) are even known explicitly [19]. This allows us to compare the numerical solution directly to the analytical stationary state in Sect. 6.2.1.
Throughout the remainder of this section we apply the scheme (9) to system (1) using different self-diffusions, , and cross-diffusions, ν.

Error and numerical order of convergence
This section is dedicated to validating our main convergence result, Theorem 2.4. Due to the lack of explicit solutions we compute the numerical solution on a fine grid and consider it a benchmark solution. We then compute numerical approximations on coarser grids and study the error in order to obtain the numerical convergence order.
In all our simulations we use a fourth order Runge-Kutta scheme to solve the ordinary differential equations Eq. (9) with initial data Eq. (5). The discrepancy between the benchmark solution and numerical solutions on coarser grids is measured by the error Here ρ ex , η ex denote the benchmark solutions. We use this quantity to study the convergence of our scheme as the grid size decreases.

No non-local interactions
Let us begin with the purely diffusive system. We consider system (1) without any interactions, i.e. W i j ≡ 0, for i, j = 1, 2, and we choose ν = 0.5 and = 0.1.
In Fig. 2 we present the convergence result as we decrease the grid size. We computed a benchmark solution on a grid of x = 2 −10 on the time interval [0, 10] with t = 0.05. Figure 2a shows the convergence for symmetric initial data whereas Fig. 2b shows the convergence of the same system in case of asymmetric initial data. In both cases we overlay a line of slope one and we conclude that the numerical convergence is of order one.

Gaussian cross-interactions
Next, we add non-local self-interaction and cross-interactions. We choose smooth Gaussians with different variances. These potentials, like the related, more singular Morse potentials, are classical in mathematical biology since oftentimes the availability of sensory information such as sight, smell or hearing is spatially limited [18,20,27].

Fig. 2
In the purely diffusive system all interaction kernels are set to zero. Both graphs show the convergence to the benchmark solution. The triangular markers denote the discrepancy between the numerical solution and the benchmark solution. A line of slope one is superimposed for the ease of comparison. On the left we start the system with symmetric initial conditions, on the right we start with asymmetric initial data. In both cases the scheme has numerical convergence order 1 For the intraspecific interaction we use for the interspecific interactions.
We consider system (1) with the diffusive coefficients ν = 0.4 and ∈ {0.1, 0.5}, and we initialise the system with 9]. Here the constant c is such that ρ and η have unit mass. Figure 3 depicts the simulation with Gaussian kernels as interaction potentials. In Fig. 3a, b we present the error plots corresponding to ∈ {0.1, 0.5}. Again we observe convergence to the reference solution with a first order accuracy in space.

General behaviour of solutions and stationary states
In this section we aim to study the asymptotic behaviour of system (1) numerically. Let us begin by going back to the set up of Sect. 6.1.2. We note that the potentials were chosen in such a way that there is an attractive intraspecific force, and the crossinteractions are chosen as attractive-repulsive explaining the segregation observed in

Fig. 4
We choose Gaussian interaction kernels of different strengths and ranges for the self-interaction and the cross-interaction, respectively. The graphs correspond to the simulationed stationary states for self-diffusivities = 0.1, and = 0.5, respectively Fig. 4a. For larger self-diffusivity, = 0.5, we see that some mixing occurs. In the absence of any individual diffusion we would have expected adjacent species with a jump discontinuity at their shared boundary [19]. However, this phenomenon is no longer possible as we have a control on the gradients of each individual species, by Lemma 3.3, rendering jumps impossible thus explaining the continuous transition.
In the subsequent section we shall push our scheme even further by dropping the smoothness assumption on our potentials.

Case of singular potentials
In this section we go beyond the limit of what we could prove in this paper. On the one hand we consider more singular potentials and on the other hand we consider vanishing individual diffusion. We study system (1) for ∈ {0, 0.02, 0.04, 0.06, 0.09}, and (a) In the case ν = 0.05, = 0, we obtain a great agree-a ment of the numerically computed stationary states and the analytical stationary states described in [19].
(b) Adding individual diffusion may still lead to segregated stationary states. However both specie remain continuous as they mix.
(c) The case ν = 0.5, = 0 leads to adjacent stationary states. Again we see an excellent agreement of the numerical stationary states and the analytical ones [19].

(d)
The regularising effect of the individual diffusion, by Corollary 3.3, becomes apparent immediately. Instantaneously both species become continuous as they start to intermingle in a small region. This region grows as we keep increasing the individual diffusion.
for the self-interaction terms and for the cross-interactions. The system is posed on the domain [0, 5] with a grid size of x = 2 −8 . Note that the case of = 0 corresponds to the absence of individual diffusion, see Fig. 5a, c. By virtue of Corollary 3.3, it is the individual diffusion that regularises the stationary states, in the sense that we will not observe any discontinuities in either ρ or η. As we add individual diffusion we can see the immediate regularisation. While stationary states may still remain segregated, as it is shown in Fig. 5b, adjacent solutions are not possible anymore (see Fig. 5d).   In the case of attractive-repulsive interspecific interactions, i.e.
we expect both species to segregate [19]. We initialise the system with the following symmetric initial data as symmetric initial data are known to approach stationary states [19]. Here the constant c normalises the mass of ρ and η to one. Figure 5a, c show our scheme performs well even in regimes we are unable to show convergence due to the lack of regularity in the potentials as well as the lack of regularity due to the absence of the porous medium type self-diffusion. While the schemes developed in [14,19] are asymptotic preserving, their convergence to weak solutions of the respective equations could not be established. We reproduce the steady states of [19] that exhibit phase separation phenomena. Figure 6 displays the stationary state in the case ν = 0.09 and = 0. Even in the case of no regularising individual diffusion and with Newtonian cross-interactions we observe a numerical convergence order of one.
In the case of attractive-attractive cross-interactions, i.e. W 12 (x) = |x| = W 21 (x), we observe an interesting phenomenon. Even in the absence of the individual diffusion, i.e. = 0, some additional mixing occurs even though we expect sharp boundaries, see Fig. 7, due to numerical diffusion. This is in contrast to the finite volume schemes proposed in [14,19].  7 We choose m 1 = 0.6 and m 2 = 0.1 in order to be able to compare the stationary state with the explicit one given in [19]. We can see a strong resemblance between the numerical stationary state and the one obtained analytically. However there are some regimes of mixing due to numerical diffusion

Energy dissipation
It is known that system (1) has a formal gradient flow structure, cf. [24], whenever W 12 = W 21 . In this case, the evolution of system 1 is such that it decays the energy functional E(ρ, η) := 1 2 ρW 11 ρ + ηW 22 η dx + ρW 12 η dx Here, we present two examples, one corresponding to the potential W ii (x) = x 2 2 , and W i j (x) = |x|, for i, j = 1, 2 and i = j, cf. Fig. 8, the other one corresponding to for i, j = 1, 2 and i = j, cf. Fig. 9.
In the first case, we choose the initial data where c > 0 normalises the mass to one. Due to the long-range we observe an attraction of the two initial bumps until they meet. They begin to mix until they are completely merged. The graph in the last panel shows the decay of the associated energy to a constant one which corresponds to the energy of the steady state. It appears the energy  is dissipated at an exponential rate but an analytic result for systems, corresponding to that of a single equation, cf. [21], is not known to our knowledge. In the second case, for Gaussian potentials, we change the computational domain to (0, π), for convenience. We choose the initial data ρ(x) = sin(2x) 2 , and cos(2x) 2 , cf. Fig. 9. We observe the formation of nearly segregated clusters which, as the evolution continues, as begin to merge due to the nonlocal interaction. However, the short-range cross-interaction is working against this trend which explains that the evolution slows down just before the merging, a phenomenon which is also observed in meta-stability. After the 5 clusters have merged into three the profile stabilises which is reflected in the evolution of the energy, cf. graph in the panel. We still observe a decay, however, after a strong initial decrease the energy decays much slower for a while before going to the constant corresponding to the stationary state. The explanation lies in the increase of the internal energy. Initially, we have ρ + η ≡ 1 which is a minimiser of the internal energy. Due to the nonlocal interactions the system wants

Fig. 9
Evolution of mixed initial data, ρ = sin(2x) 2 and η = cos(2x) 2 on the domain (0, π). The interactions are linear combinations of Gaussians modelling self-attraction whereas the cross-interactions are short-range repulsive and long-range attractive. The associated energy decays abruptly at the merging and separation between aggregates. The slower decay of the energy before t ≈ 10 is due to the trade off between the local (internal) energy and the nonlocal interaction energy. After the rearrangement of the five initial clusters to only three, the energy stabilises to rearrange but at the cost of increasing the internal energy to allow for a decrease in the interaction energy. This beautifully portrays the interplay of local and nonlocal effects. Similar effects are known in the context of meta-stability.

Conclusion
In this paper we presented a finite volume scheme for a system of non-local partial differential equations with cross-diffusion. We were able to reproduce a continuous energy estimate on the discrete level for our scheme. These discrete estimates for the approximate solution are enough to get compactness results and we are able to identify the limit of the approximate solutions as a weak solution of the equation. We complement the analytical part with numerical simulations. These back up our convergence result and we are also able to apply the scheme in cases in which we cannot show convergence. To this end we pushed the scheme to regimes of singular potentials also lacking the regularising porous medium type self-diffusion terms. Comparing them with the explicit stationary states from [19] we conclude the scheme performs well even in regimes it was not designed for.