Sorting Phenomena in a Mathematical Model For Two Mutually Attracting/Repelling Species

TL;DR

This paper investigates a mathematical model for two interacting species with attraction and repulsion, demonstrating phase separation and spatial sorting, with proofs of stationary solutions' existence, shape, and conditions for uniqueness.

Contribution

It extends existing models to include different intra- and inter-species attractions, proving phase separation and stationary solution properties using advanced mathematical techniques.

Findings

Support of stationary solutions has zero Lebesgue measure.

Support of the sum of densities is simply connected.

Existence and shape of segregated solutions are established.

Abstract

Macroscopic models for systems involving diffusion, short-range repulsion, and long-range attraction have been studied extensively in the last decades. In this paper we extend the analysis to a system for two species interacting with each other according to different inner- and intra-species attractions. Under suitable conditions on this self- and crosswise attraction an interesting effect can be observed, namely phase separation into neighbouring regions, each of which contains only one of the species. We prove that the intersection of the support of the stationary solutions of the continuum model for the two species has zero Lebesgue measure, while the support of the sum of the two densities is simply connected. Preliminary results indicate the existence of phase separation, i.e. spatial sorting of the different species. A detailed analysis in one spatial dimension follows. The existence and shape of segregated stationary solutions is shown via the Krein-Rutman theorem. Moreover, for small repulsion/nonlinear diffusion, also uniqueness of these stationary states is proved.

Equations

$${\partial }_t\rho =\mathrm{div}\left[\rho \nabla \left(a(\rho )-W\ast \rho \right)\right],$$

$${\partial }_t\rho =\mathrm{div}\left[\mu (\rho ,v)\nabla \rho -\chi (\rho ,v)\rho \nabla v\right],\tau {\partial }_{t}v=\eta \Delta v+k(\rho ,v).$$

$$\mathcal{F}[\rho ]=\int A(\rho )dx-\frac{1}{2}\int \rho W\ast \rho dx\text{mbox}whereA(\rho )={\int }_0^{\rho }a(\xi )d\xi ,$$

$$\mathcal{E}[{\rho }_1,{\rho }_2]=\int f({\rho }_1,{\rho }_2)dx-\frac{1}{2}\int {\rho }_1{S}_1\ast {\rho }_{1}dx-\frac{1}{2}\int {\rho }_2{S}_2\ast {\rho }_{2}dx-\int {\rho }_{1}K\ast {\rho }_{2}dx,$$

$$\begin{cases}\partial_{t}\rho_{1}=\mbox{div}\big{[}\rho_{1}\nabla\big{(}f_{\rho_{1}}(\rho_{1},\rho_{2})-S_{1}\ast\rho_{1}-K\ast\rho_{2}\big{)}\big{]}&\\ \partial_{t}\rho_{2}=\mbox{div}\big{[}\rho_{2}\nabla\big{(}f_{\rho_{2}}(\rho_{1},\rho_{2})-S_{2}\ast\rho_{2}-K\ast\rho_{1}\big{)}\big{]}\,.&\end{cases}$$

$$f({\rho }_1,{\rho }_2)=ϵ({\rho }_1+{\rho }_2{)}^2/2,ϵ>0,$$

$$\begin{cases}\partial_{t}\rho_{1}=\mbox{div}\big{[}{\epsilon}\rho_{1}\nabla(\rho_{1}+\rho_{2})-\rho_{1}\nabla S_{1}\ast\rho_{1}-\rho_{1}\nabla K\ast\rho_{2}\big{]}\\ \partial_{t}\rho_{2}=\mbox{div}\big{[}{\epsilon}\rho_{2}\nabla(\rho_{1}+\rho_{2})-\rho_{2}\nabla S_{2}\ast\rho_{2}-\rho_{2}\nabla K\ast\rho_{1}\big{]}\end{cases}.$$

$${\partial }_t\rho ={\partial }_x\left[\mu (\rho ){\partial }_x\rho -\gamma (\cdot ,\rho )\right],$$

$$\Delta v+K{e}^{2v}-1=0,$$

$${\rho }^n\mapsto {\rho }^{n+1}={\mathrm{argmin}}_{\rho }\left(F[\rho ]+\frac{1}{2\Delta t}{d}_2^2({\rho }^n,\rho )\right),$$

$$\mathcal{G}[\rho]=\frac{1}{2}\int\rho\big{(}\epsilon\rho+W\ast\rho\big{)}dx,$$

$${\rho }_t=[\rho (ϵ\rho +W\ast \rho {)}_x{]}_x,$$

$${\partial }_{t}U=\text{mbox}div\left(D(U)\nabla U\right)-\text{mbox}div\left(\left(\begin{array}{cc}{\rho }_1& 0\\ 0& {\rho }_2\end{array}\right)\nabla \left(\begin{array}{c}{S}_1\ast {\rho }_1+K\ast {\rho }_2\\ S_2\ast {\rho }_2+K\ast {\rho }_1\end{array}\right)\right),$$

$$D(U)=\left(\begin{array}{cc}{\rho }_1{f}_{{\rho }_1,{\rho }_1}& {\rho }_1{f}_{{\rho }_1,{\rho }_2}\\ {\rho }_2{f}_{{\rho }_1,{\rho }_2}& {\rho }_2{f}_{{\rho }_2,{\rho }_2}\end{array}\right).$$

$$\begin{cases}\partial_{t}\rho_{1}=\mbox{div}\big{(}{\epsilon}\rho_{1}\nabla(\rho_{1}+\rho_{2})-\rho_{1}\nabla V_{1}\big{)},\\ \partial_{t}\rho_{2}=\mbox{div}\big{(}{\epsilon}\rho_{2}\nabla(\rho_{1}+\rho_{2})-\rho_{2}\nabla V_{2}\big{)}\,,\end{cases}$$

$$\text{mbox}supp({\rho }_1^{\infty })\cap \text{mbox}supp({\rho }_2^{\infty })\subseteq \{\nabla V_1=\nabla V_2\}.$$

$$ϵ\nabla ({\rho }_1^{\infty }+{\rho }_2^{\infty })(x)=\nabla V_1(x)=\nabla V_2(x),$$

$$\text{mbox}Let{V}_1=S_1\ast {\rho }_1+K\ast {\rho }_2,V_2=S_2\ast {\rho }_2+K\ast {\rho }_1.$$

$$w:={\rho }_1+{\rho }_2,\zeta :=\frac{{\rho }_1-{\rho }_2}{w},$$

$${\partial }_{t}w$$

$$w\partial_{t}\zeta=\left({\epsilon}w\nabla w-w\frac{1+\zeta}{2}\nabla V_{1}-w\frac{1-\zeta}{2}\nabla V_{2}\right)\cdot\nabla\zeta-\frac{1-\zeta^{2}}{2}\mbox{div}\big{(}w\nabla(V_{1}-V_{2})\big{)}.$$

$$\begin{cases}0=\big{(}{\epsilon}\rho_{1}(\rho_{1}+\rho_{2})_{x}-\rho_{1}S^{\prime}_{1}\ast\rho_{1}-\rho_{1}K^{\prime}\ast\rho_{2}\big{)}_{x}\\ 0=\big{(}{\epsilon}\rho_{2}(\rho_{1}+\rho_{2})_{x}-\rho_{2}S^{\prime}_{2}\ast\rho_{2}-\rho_{2}K^{\prime}\ast\rho_{1}\big{)}_{x},\end{cases}$$

$$\mathcal{F}[{\rho }_1,{\rho }_2]=\frac{ϵ}{2}\int ({\rho }_1+{\rho }_2{)}^{2}dx-\frac{1}{2}\int {\rho }_1{S}_1\ast {\rho }_{1}dx-\frac{1}{2}\int {\rho }_2{S}_2\ast {\rho }_{2}dx-\int {\rho }_{1}K\ast {\rho }_{2}dx.$$

$$\mathcal{M}=\left\{(\rho_{1},\rho_{2})\in(L^{2}({\mathbb{R}})\cap L^{1}({\mathbb{R}}))^{2}~{}\Big{|}~{}\rho_{i}\geq 0,\,\int_{\mathbb{R}}\rho_{i}~{}dx=m_{i}\,,\,i=1,2\right\},$$

$$\displaystyle 0=\int\rho_{1}\Big{(}{\epsilon}(\rho_{1}+\rho_{2})_{x}-S_{1}^{\prime}\ast\rho_{1}-K^{\prime}\ast\rho_{2}\Big{)}U_{x}\,dx,$$

$$\displaystyle 0=\int\rho_{2}\Big{(}{\epsilon}(\rho_{1}+\rho_{2})_{x}-S_{2}^{\prime}\ast\rho_{2}-K^{\prime}\ast\rho_{1}\Big{)}V_{x}\,dx,$$

$$\mathcal{F}[{\rho }_1,{\rho }_2]=\frac{ϵ}{2}{\int }_{\mathbb{R}}({\rho }_1+{\rho }_2{)}^{2}dx-\frac{1}{2}{\int }_{\mathbb{R}}{\rho }_1{S}_1\ast {\rho }_{2}dx-\frac{1}{2}{\int }_{\mathbb{R}}{\rho }_2{S}_2\ast {\rho }_{2}dx-{\int }_{\mathbb{R}}{\rho }_{1}K\ast {\rho }_{2}dx,$$

$$\frac{d}{d\rho_{1}}{\mathcal{F}}[\rho_{1},\rho_{2}](\mu)=\lim_{\delta\rightarrow 0}\frac{1}{\delta}\Big{(}{\mathcal{F}}[\rho_{1}+\delta\mu,\rho_{2}]-{\mathcal{F}}[\rho_{1},\rho_{2}]\Big{)},$$

$$\mathcal{F}[{\rho }_1+\delta \mu ,{\rho }_2]-\mathcal{F}[{\rho }_1,{\rho }_2]$$

$$-\frac{1}{2}{\int }_{\mathbb{R}}({\rho }_1+\delta \mu )S_1\ast ({\rho }_1+\delta \mu )dx+\frac{1}{2}{\int }_{\mathbb{R}}{\rho }_1{S}_1\ast {\rho }_{1}dx$$

$$-{\int }_{\mathbb{R}}({\rho }_1+\delta \mu )K\ast {\rho }_{2}dx+{\int }_{\mathbb{R}}{\rho }_{1}K\ast {\rho }_{2}dx$$

$$=\frac{ϵ}{2}{\int }_{\mathbb{R}}(\delta \mu {)}^2+2\delta \mu ({\rho }_1+{\rho }_2)dx-\frac{1}{2}{\int }_{\mathbb{R}}2\delta \mu S_1\ast {\rho }_{1}dx$$

Full text

Preliminary results indicate the existence of phase separation, i.e. spatial sorting of the different species. A detailed analysis is given in one spatial dimension. The existence and shape of segregated stationary solutions is shown via the Krein-Rutman theorem. Moreover, for small repulsion/nonlinear diffusion, also uniqueness of these stationary states is proved.

1. Introduction

The interplay between (nonlinear) diffusion and nonlocal attractive/repulsive interactions arises in a variety of contexts in the natural-, life-, and social sciences. As a non exhaustive list of examples, let us mention granular media physics [7, 19], astrophysics [13, 81], semiconductors [41, 81], chemotaxis [55, 72, 49, 2, 70, 42], ecology, animal swarming and aggregation [62, 66, 2, 18, 63, 77], alignment [71, 43], and opinion formation [76, 1]. One way of deriving such models from first order microscopic systems of (stochastic) interacting particles $x_{1},\ldots,x_{n}$ is the following. Each particle $x_{i}$ is driven by nonlinear forces due to short range repulsion, see [63] - respectively undergoes an independent Brownian motion. Further, it moves towards higher concentrations of particles of its own kind, respectively those of an external signal, [73]. At the macroscopic and continuum level, this set of rules results in a nonlocal partial differential equation

[TABLE]

respectively, a chemotaxis-type system

[TABLE]

For $\mu(\rho,v)=\rho a^{\prime}(\rho)$, $\chi(\rho,v)=\chi_{0}$, $\tau=0$, and $k(\rho,v)=\rho-\beta v$, thus $v=(\beta I-\Delta)^{-1}\rho$, system (2) is equivalent to (1), with $W$ being the Newtonian- or Bessel-potential ($\beta=0$ or $\beta>0$). Especially when $a(\rho)=\log(\rho)$, i.e. $\mu(\rho,v)=1$ in this case, then a prototype model for nonlocal aggregation phenomena results, namely the simplified, classical parabolic-elliptic Keller-Segel model [55], or in physics a model for gravitational self-interacting clusters [81].

In (1), $\rho=\rho(x,t)$ denotes the density of particles, $a:[0,\infty)\rightarrow[0,+\infty)$ is a $C^{1}$ monotone increasing function with $a^{\prime}(0)=0$, and $W(x)=\tilde{W}(|x|)$ is a $C^{1}$ potential with $\tilde{W}^{\prime}(r)<0$ for all $r>0$. The PDE (1) can be interpreted as the gradient flow of the functional

[TABLE]

w.r.t. the $d_{2}$ Wasserstein distance arising in optimal transport theory. Compare also the different Lyapunov functionals and energies used for (2) in [81, 42, 70], e.g. $A(\rho)=\rho\log\rho$. Roughly speaking, this means that the velocity $\tilde{v}=-\nabla\left(a(\rho)-W\ast\rho\right)$ in the continuity equation $\rho_{t}+\mathrm{div}\left(\rho\tilde{v}\right)=0$ can be interpreted as the sub-differential of $\mathcal{F}$ in the metric space $\mathcal{P}_{2}$ of probability measures with finite second moment and metric $d_{2}$, see [3] for more details.

Models of type (1), (2) with two or more species being involved are considered in the context of chemotaxis [40, 48, 82, 83], opinion formation [39], pedestrian dynamics [4], and population biology [30, 36]. A reasonable generalization of (3) for two species is

[TABLE]

where $f:[0,+\infty)\times[0,+\infty)\rightarrow{\mathbb{R}}$, $f\in C^{1}$, and $S_{1}$, $S_{2}$, $K$ are $C^{1}$ and radially decreasing potentials like $W$ above. The first term of $\mathcal{E}$ typically represents local repulsion, while the nonlocal terms model attractive forces. We call $S_{1}$ and $S_{2}$ self-interaction potentials and $K$ cross-interaction potential. For $(\mathcal{P}_{2}({\mathbb{R}}^{d}),d_{2})\times(\mathcal{P}_{2}({\mathbb{R}}^{d}),d_{2})$ with the natural product topology, the formal gradient flow of $E[\rho_{1},\rho_{2}]$ w.r.t. this metric structure is given by

[TABLE]

The functionals we investigate in this paper are of the form (4), with

[TABLE]

leading to

[TABLE]

Besides modelling local repulsion and global attraction between different types of particles (cf. the introduction and references in [22]) our particular motivation is to consider cell sorting due to differential attraction and the resulting pattern formation. Similar cells or species of equal size (consequently the terms $(\rho_{1}+\rho_{2})$ in (6)), but with different reactions to attraction forces undergo a reorganization process, where cells with stronger self-attraction finally sort into the center of the total cell population and those with weaker self-attraction to the outside. Differential attraction can either be by an external (chemical) signal or directly between the species. Such phenomena are observed in developmental processes. In the first case one may assume that $S_{1}$, $S_{2}$, and $K$ are multiples of the same kernel, modelling indirectly a chemo-attractant; like the fundamental solution of the elliptic equation for the chemo-attractant in the prototype Keller-Segel system, which can be rewritten into the single prototype aggregation equation. We are interested in suitable conditions for the self- and the cross-attraction, which result in the above mentioned sorting phenomenon. In [40] for differential attraction, i.e. different chemotactic sensitivities, of two species towards higher concentrations of one chemo-attractant, it was proved that if the solution for the more strongly attracted species blows up, than also the second one blows up at the same time. The amount of mass, which concentrates in the joint blowup is different though, and controlled by the system parameters. This last, formal result, hints towards a possible separation of the main amount of masses of the two species. The more strongly attracted species accumulates more mass in the blowup than the other one. Cell sorting due to differential adhesion/attraction was also discussed in [80, 58], where stochastic particle models and numerical studies of continuum models were considered. Based on experimental results, see e.g. [61], sorting due to differential chemotaxis and chemical spiral waves was modeled and simulated in [78, 68]. In [54] a rigorous analysis was given.

We are interested in steady states of (6) and minimizers of (4), hence let us first discuss some relevant results in the single species setting (2) and (1). Varying $\chi$ and $k$ in (2) plays a similar role as varying $W$ in (1) for the pattern formation properties of the respective solution. Though a one-to-one connection between all variants of aggregation equations and chemotaxis systems as in the case of the Newtonian- and Besselpotential for $(-W)$ and the prototype Keller-Segel system has not been established, both types of models and their dynamics are strongly intertwined.

A wealth of mathematical results on existence of global solutions, blowup phenomena and pattern formation exists for these models, see e.g. the summaries in [46, 47, 69]. The longtime behavior of the prototype models is especially interesting in two dimensions as a biological model for the peculiar chemotactic self-organisation of Dictyostelium discoideum [49], and for self-graviational collapse and star formation in astrophysics in dimension three [13]. In both cases, see also [66], blowup phenomena in the respective dimension, and point support of stationary Dirac-type solutions are of crucial relevance for the respective application. Non-trivial stationary solutions and their qualitative features have been analyzed in [72] for the general system (2) with Neumann boundary conditions, for $\chi(\rho,v)=\chi_{0}\mu(\rho,v)\rho\phi^{\prime}(v)$, and $\tau=1$. Especially the one dimensional steady states and their stability are well understood. These results should be compared with aggregation equation analoga. With $\rho=C\exp(\chi_{0}\phi(v))$ the steady state analysis in [72] reduces to analyzing the elliptic equation $\eta\Delta v+k(\psi(v,\lambda),v)=0$, where $\psi$ is the flux of $r^{\prime}(s)=\chi(r,s)/\mu(r,s)=\chi_{0}\rho\phi^{\prime}(v)$.

In [2] density dependent diffusion-drift equations for aggregation of the form

[TABLE]

are analyzed with e.g. $\mu(\rho)=m\rho^{m-1}$ and $\gamma(t,x,\rho)$ being proportional to $\rho$, respectively depending on a functional of the density distribution $\rho(t,\cdot)$. Thus integral terms are included, as discussed in [62, 66], and (7) can be compared with (1). The connections between chemotaxis systems and (7) are given, and a Lyapunov function is constructed, i.e. a functional decreasing along solutions as time increases. This provides a general theorem on global existence of solutions. Asymptotic convergence to steady states is proved, which can be non-trivial, for instance plateau-like.

Free energies in such models and in chemotaxis are based on an entropic term, e.g. $A(\rho)=\rho\log\rho-\rho+1$ (or just $A(\rho)=\rho\log\rho$) in (3). In [81] not only global minimizers of (3) are studied, but the whole solution set. Here $\rho$ is a steady solution of the governing PDE - the associated gradient flow - if and only if it is a critical point of the variational problem with constraint $\int\rho=M$. Due to the logarithmic term in the entropy we have $a(\rho)\sim\log\rho$ and inverting this relation to an entropy variable $\rho\sim\exp a(\rho)$ one can obtain simplified systems for stationary states. With this exponential transformation and the Moser-Trudinger inequality as given in [65, 64] the critical mass $8\pi$ for graviational collapse is deduced. Existence, uniqueness, stability and symmetry breaking of stationary solutions are proved via the precise connections between the free energy and the respective Vlasov-Fokker-Planck (chemotaxis-like) equation.

In [42] a Lyapunov functional - c.f. the energy functional in (3) - for the prototype chemotaxis-model in two dimensions is provided, and extended to more general equations in [70]. Also here the connection to the exponential transformation and the elliptic equation given in [65] is notified. Geometric criteria for not necessarily trivial stationary states are derived. Now any metric $ds^{2}$ on a two-dimensional sphere determines a Gauss curvature function $K$ satisfying the Gauss-Bonnet formula, ${\int\int}_{S^{2}}Kd\tilde{\mu}=4\pi$. Here $d\tilde{\mu}$ is the volume element of $S^{2}$. To characterize all $K$ belonging to metrices $ds^{2}$ and relating to the standard metric $ds_{0}^{2}$ via $ds^{2}=p\,ds_{0}^{2}$, with a positive function $p$ on the sphere, one has to determine $p=p(K)$ uniquely. In [65] it is proved that the transformation $p=\exp(v)$ reduces this question to solving

[TABLE]

on the sphere, which is done by a variational approach. This is a specific form of the elliptic problem analyzed in [72], where the more general transformation $\rho=C\exp(\chi_{0}\phi(v))$ was used for the steady states analysis of generalized chemotaxis systems. With $k(\psi(v,\lambda),v)=Ke^{2v}-1$, i.e. $\psi(v,\lambda)=Ke^{2v}$, respectively $\chi_{0}\phi=2v$ in [72], the above elliptic PDE results. Qualitative features of simplified chemotaxis systems with non-linear diffusion have been discussed in [74, 75]. See also further references therein.

Existence and uniqueness for (1) are discussed e.g. in [28], the last part of the book [3], and earlier via the related works on chemotaxis, c.f. [49, 45, 46, 47, 74, 75]. In [8] also singular potentials $W$ are considered, e. g. Coulomb potential as in electrodynamics, $W(x)=|x|$, $W(x)=|x|^{\alpha}$, $\alpha>2-d$. We also refer to [17, 14, 15] for further discussion.

For $a(\rho)=\log\rho$ in (3), without any nonlocal effects, but with the addition of a confinement external potential $V$, in [50] weak solutions to the linear Fokker-Planck equation in a gradient flow setting were derived by constructing the finite time-step minimizing movement

[TABLE]

and taking the limit $\Delta t\rightarrow 0$. This idea was applied to (1) in [28], and to a more general metric framework in [3].

The large time asymptotics of (1) depend on the competing effects of $\int A(\rho)$, which promotes particle spreading, such that the density $\rho$ stabilises at a constant state (which is zero if we consider ${\mathbb{R}}^{d}$), whereas the nonlocal term drives the particles towards aggregation. A Dirac delta results, which is located at the (preserved) center of mass of the system. To prove existence of a global minimum of $\mathcal{F}$ is therefore a challenging problem. In [5] this issue was tackled in detail by using Lions’ concentrated compactness technique.

The limiting case $a=0$ was analyzed in [33, 24], see also the references therein. The case $a(\rho)=\rho$ was derived as a nonlocal repulsive effect under the action of a potential $V_{\epsilon}$ which converges weakly to a Dirac delta as $\epsilon\rightarrow 0$. For a formal argument see [63], see also [67] for a rigorous derivation via interacting particle systems with short-repulsion (the repulsion range shrinking to zero as the number of particles goes to infinity). The results in [22] and [5] provide a clear picture of the behaviour of the corresponding functional

[TABLE]

with a given diffusion constant $\epsilon>0$, and a $W^{1,1}$ attractive potential $W\leq 0$. If $\epsilon\geq\|W\|_{L^{1}}$, then $\mathcal{G}$ is uniformly convex, and hence zero is its global minimizer. This suggests that in this case $\rho\rightarrow 0$ for large times. On the other hand, if $\epsilon<\|W\|_{L^{1}}$, a nontrivial global minimizer exists in the class $\int\rho\,dx=1$. This suggests existence of nontrivial steady state for large times, but the solution to the Cauchy problem could still decay to zero if a large enough variance is present initially. This problem is still open in the case of slow diffusion, see [6]. In one dimension a more refined analysis can be given. In [22] it was proved that a unique steady state with given mass and center of mass exists for

[TABLE]

provided that $\epsilon<\|W\|_{L^{1}}$ and $W$ is radially increasing, negative, and supported on the whole of ${\mathbb{R}}$. Such a steady state is also the unique global minimizer for $\mathcal{G}$ for fixed mass and fixed center of mass, it is compactly supported, symmetric, and in $W^{1,\infty}$, with a shape similar to a Barenblatt profile for the porous medium equation, see [79]. Such steady state is locally stable for large times, see the recent [38]. The result in [22] was partly extended to more general nonlinear diffusions in [23], where the existence of a unique diffusion constant for a given support of the steady state is proven. Both results rely on the Krein-Rutman theorem in order to characterise the steady states as eigenvectors of a certain nonlocal operator. Further generalizations and completions of some open questions in those papers were recently given in [53, 52]. One major issue solved in [22] is to prove that a one-to-one correspondence exists between the diffusion constant (eigenvalue) and the support of the steady state. The results in [22] improved a previous result in [20] about the existence and uniqueness of nontrivial steady states for small diffusion constant $\epsilon\ll 1$ performed via an implicit function theorem argument.

Now let us come back to the functional (4). The first term of $\mathcal{E}$ typically represents local repulsion. Apart from nonlocal cross-diffusion via the kernels $S_{i}$ also $f$, when containing mixed terms in $\rho_{1}$ and $\rho_{2}$, introduces cross diffusion. Indeed, (5) can be formally rewritten as

[TABLE]

with $U=(\rho_{1},\rho_{2})$,

[TABLE]

If $D(U)$ is symmetric and semi-positive definite, then the theory in [51] applies, see also the recent [34] in the context of Wasserstein gradient flows. However, in our situation $D(U)$ is never symmetric. This makes existence proofs for solutions of (5) in full generality a delicate problem. In most cases, $D(U)$ does not even have a semi-positive definite symmetric part.

A comprehensive existence and uniqueness theory for (5) with $f\equiv 0$ was given in [35], where convergence of the JKO scheme leads to existence of weak measure valued solutions for ‘mildly singular’ potentials. A suitable notion of displacement convexity for systems provides a uniqueness result. An implicit-explicit variation of the JKO scheme yields an existence result in [35], also for nonlocal terms which are not of gradient flow type, i.e. for non symmetric cross-interaction terms. Without the nonlocal interaction part, an existence theory with a direct application of the JKO scheme is given in [57] for a fully coupled system of two second-order parabolic degenerate equations arising as a thin film approximation to the Muskat problem. The peculiar form of the cross-diffusion in this case allows for an energy functional with good coercivity properties. In [57] the regularity needed in order to identify the suitable Euler-Lagrange equation in the variational scheme was obtained. Although not strictly related to our model, let us also mention the hybrid variational scheme generalizing the JKO scheme. This has been introduced in [16] for the parabolic-parabolic Keller-Segel model in ${\mathbb{R}}^{2}$ working in the product space ${\mathscr{P}}_{2}({\mathbb{R}}^{2})\times L^{2}({\mathbb{R}}^{2})$. Using a modified Wasserstein distance between vector-valued densities on ${\mathbb{R}}$, in [84] the variational structure of systems with degenerate diffusion and nonlinear reaction terms was investigated. A more general approach has been discussed in [59], where a convexity concept for reaction-diffusion systems was developed, that allows to analyze some important examples.

In system (6) the symmetric part of $D(U)$ is not semi-positive definite for all $\rho_{1},\rho_{2}\geq 0$ in general. A brute-force JKO approach would imply a good dissipation estimate only for $w=\rho_{1}+\rho_{2}$, but not for both $\rho_{1}$ and $\rho_{2}$. Thus a regularising effect for the sum $w$ is still possible. On the other hand, such a degenerate dissipation estimate does not prevent the formation of discontinuities for $\rho_{1}$ and $\rho_{2}$ separately. In order to partly validate this hypothesis, we focus on stationary states in one space dimension, and show that discontinuous stationary patterns may arise. More precisely, we prove that $\rho_{1}$ and $\rho_{2}$ can separate completely in the stationary state, i.e. they feature a jump discontinuity at exactly the same point, with the sum $w=\rho_{1}+\rho_{2}$ remaining smooth at that point. We call such a configuration a fully segregated steady state.

Segregation in multi-species systems with nonlinear diffusion terms like in (6) and possible reaction terms has been widely investigated in the literature. We refer to [11, 10, 12, 9] and references therein. It is well known that segregated initial data produce a unique segregated solution. The problem of existence of a mixed solution is partly open. Difficulties in our case occur due to the presence of the nonlocal terms, which typically do not allow for a local comparison principle and may produce aggregation phenomena. The emergence vs. non-emergence of segregated steady states in two species systems with nonlocal attraction has recently been treated in [31], where the repulsive effect of the nonlinear diffusion has been replaced by an upper bound for $\rho_{1}+\rho_{2}$. Then the minimisation problem for $\mathcal{E}$ with $f=0$ and $0\leq\rho_{1}+\rho_{2}\leq 1$ with all interaction potentials being multiples of a given function $K$ is analyzed. Sharp conditions on these multiplying factors are provided, which results in complete segregation of the two species. In [60] rearrangements for three already separated domains has been considered in a different context.

More specifically on system (6) (almost parallel to our result), the recent [27] shows how to construct explicit segregated and non-segregated steady states with power-law interaction potentials (with possible confinement effects) and produces numerical evidence that segregation may not occur for large $\epsilon$. At the same time, [56] shows that initially segregated initial data yield segregated solutions for all times for a system with the same (cross-)diffusion term of (6) and external potentials (under suitable assumptions on the latter). We also mention at this stage the result in [25], in which an existence theory for a reaction diffusion system with the same (cross-)diffusion of (6) has been proven via JKO scheme with $BV$ initial conditions, without any initial separation assumption.

Let us know briefly summarise our results. For the one-dimensional case, we first investigate conditions for existence or non-existence of non trivial stationary solutions. For small diffusion coefficient ${\epsilon}$ we show existence of segregated stationary states via the implicit function theorem. We extend the strategies in [21, 22] to the multi species case. We also relate stationary solutions to energy minimizers in a rigorous way and provide some results characterizing their structure. In particular we verify that for stationary solutions and energy minimizers, the sum $w$ is supported on a connected interval and we give a rigorous result on the segregation in the case of dominant self-attraction. In a case of weak self-attraction we can characterized the support of the minimizers by the same arguments as in [31].

For the interaction kernels in (6) we assume unless further noted:

• (A1)

$S_{1},S_{2},K\in C^{2}({\mathbb{R}})$,

• (A2)

$S_{1},S_{2},K$ are radially symmetric and decreasing w.r.t. the radial variable,

• (A3)

$S_{1},S_{2},K$ are nonnegative and have finite mass on ${\mathbb{R}}$.

The paper is organized as follows. In Section 2 we give some preliminary results on phase separation for (6). In Section 3 we analyze the relation of stationary solutions with critical points of the associated energy functional and provide sufficient conditions on $\epsilon$ and on the interaction kernels yielding non existence of non-trivial steady states. Section 4 is devoted to the existence analysis of stationary solutions involving phase-separation via the Krein-Rutman based approach given in [22]. Finally, in Section 5 existence and uniqueness results for stationary solutions are proved in case of small repulsion, corresponding to small nonlinear diffusion in (6). We complement our results with some numerical simulations in Section 6 showing segregated behavior as well as mixing and diffusion-dominated behaviour for large diffusions.

2. Segregation due to differential aggregation

First, we provide some preliminary results in arbitrary dimensions about pattern formation for model (6), and more specifically on the emergence or non-existence of segregation. Consider the canonical model

[TABLE]

where $V_{1}$ and $V_{2}$ are two given smooth external potentials. We can rewrite (6) in this form, with $V_{1}$ and $V_{2}$ being determined by convolutions of $\rho_{1}$ and $\rho_{2}$ with aggregation kernels, c.f. (11).

Proposition 2.1**.**

Let $V_{i}\in L^{1}_{loc}({\mathbb{R}}^{d})\cap W^{1,\infty}({\mathbb{R}}^{d})$ be given external potentials. If $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ is a $C^{1}$ (weak) stationary solution of (9), then we have

[TABLE]

Proof.

Let $x\in\mbox{supp}(\rho_{1}^{\infty})\cap\mbox{supp}(\rho_{2}^{\infty})$. Then the weak formulation of (9) implies

[TABLE]

from which the assertion follows. ∎

An immediate consequence of Proposition 2.1 is the following result, which deals with the special case of all the interaction kernels being multiples of the fundamental solution of the Laplace equation.

Proposition 2.2**.**

Let $K$ be the fundamental solution of the Laplace equation in ${\mathbb{R}}^{d}$, and $S_{1}=\sigma_{1}K$, $S_{2}=\sigma_{2}K$ with $\sigma_{1}\leq 1\leq\sigma_{2}$ and $\sigma_{1}\neq\sigma_{2}$. Then, every $C^{1}$ stationary solution $(\rho_{1},\rho_{2})$ of (6) is fully segregated, i. e. $\mbox{supp}(\rho_{1})\cap\mbox{supp}(\rho_{2})$ has empty interior.

Proof.

[TABLE]

Assume that there is a non-empty open set $\mathcal{O}\subset\mbox{supp}(\rho_{1})\cap\mbox{supp}(\rho_{2})$. Then Proposition 2.1 implies $\mathcal{O}\subset\{\nabla V_{1}=\nabla V_{2}\}$. Hence $V_{1}-V_{2}=c=const.$ on $\mathcal{O}$, i.e. $(S_{1}-K)*\rho_{1}-(S_{2}-K)*\rho_{2}=c$. With the assumptions on the kernels, we can apply the Laplace operator and obtain $(\sigma_{1}-1)\rho_{1}(x)+(1-\sigma_{2})\rho_{2}(x)=0$, for all $x\in\mathcal{O}$. The assumptions on $\sigma_{1}$ and $\sigma_{2}$ imply $\rho_{1}(x)=\rho_{2}(x)=0$, which is a contradiction. ∎

Proposition 2.2 shows segregation of steady states in a particular case. However, this effect occurs for a much wider class of aggregation kernels $S_{1},S_{2},K$, as we will argue and partially prove below.

First, take a closer look at the dynamics of segregation by using a transformation of variables similar to the one in [12], where strong reaction terms induced the segregation though. Let

[TABLE]

where the relative difference $\zeta$ is only considered on the support of the total density $w$. Adding both equations in (9) and using the notion of (11) yields

[TABLE]

Thus the dynamics of $w$ are governed by a porous medium equation with additional convective terms. The dynamics of $\zeta$ are obtained by subtracting the equations for $\rho_{1}$ and $\rho_{2}$ and then inserting the equation for $w$, which yields

[TABLE]

As in [12] the evolution of $\zeta$ is governed by a first-order equation, which gives particular insight into the dynamics of the system. The first of the two terms on the right-hand side of the above equation is along the flux of $w$, which determines the spatial position and shape of the solution. The crucial part for the segregation dynamics is the second term, a reaction term w.r.t. $\zeta$ and two fixed points $\zeta=\pm 1$, corresponding to segregation. Depending on the sign of $\mbox{div}\left(w\nabla(V_{1}-V_{2})\right)$ one of them is stable. Thus there is always some dynamics towards a segregated state driven by the differences in the attraction forces. Instead of pursuing a time-dependent analysis inspired by the above considerations, we restrict our search for segregated states to the analysis of stationary solutions, and leave the stability vs. instability analysis of segregated patterns to future work.

3. Steady states vs energy minimization

We now explore the relation between one dimensional steady states of (6), namely solutions $(\rho_{1},\rho_{2})$ of

[TABLE]

and the minimisers of the energy functional

[TABLE]

Note, that $\rho_{1},\rho_{2}\in L^{2}({\mathbb{R}})$ is a natural setting for the minimization of ${\mathcal{F}}$, since assumption (A3) ensures that the nonlocal terms in ${\mathcal{F}}$ are finite. More precisely, we will look for minimizers within the set

[TABLE]

for given masses $m_{1}$, $m_{2}$. The functional setting for the weak solutions of (13) should involve some space derivative, since the cross-diffusion term in (13) is not of type $\Delta F(\rho_{1},\rho_{2})$ with some nonlinear vector field $F$. Therefore one cannot integrate by parts twice in the distributional formulation. For a weak formulation involving spatial derivatives of $\rho_{1}$ and $\rho_{2}$ and allowing for discontinuities of each species at the same time, we choose $BV({\mathbb{R}})$ as the functional setting for the steady states, whose sum $w=\rho_{1}+\rho_{2}$ is Lipschitz continuous. Such a choice is partly supported by the results in [25], in which the preservation of the $BV$ regularity is proven for a similar system.

Definition 3.1**.**

The pair $(\rho_{1},\rho_{2})\in\mathcal{M}\cap BV({\mathbb{R}})^{2}$ is a weak solution to (13) if $w=\rho_{1}+\rho_{2}\in\mathrm{Lip}({\mathbb{R}})$, and

[TABLE]

for arbitrary $U,V\in C^{1}_{c}({\mathbb{R}})$.

Proposition 3.1**.**

Let $\left(\rho_{1},\rho_{2}\right)\in\mathcal{M}\cap BV({\mathbb{R}})^{2}$ be a minimizer of

[TABLE]

such that $w=\rho_{1}+\rho_{2}\in\mathrm{Lip}({\mathbb{R}})$. Then $(\rho_{1},\rho_{2})$ is a weak solutions of (13) according to Definition 3.1. Moreover, every weak solution $(\rho_{1},\rho_{2})$ of (13) according to Definition 3.1 is a critical point of ${\mathcal{F}}$.

Proof.

Let $(\rho_{1},\rho_{2})\in\mathcal{M}\cap BV({\mathbb{R}})^{2}$ be a minimizer of ${\mathcal{F}}$. We calculate the first Gateaux derivative of ${\mathcal{F}}$

[TABLE]

along an arbitrary direction $\mu\in L^{2}({\mathbb{R}})$, such that $(\rho_{1}+\delta\mu,\rho_{2})\in\mathcal{M}$. We have

[TABLE]

Dividing by $\delta$ and taking the limit $\delta\rightarrow 0$ results in

[TABLE]

Let $U\in C^{1}_{c}({\mathbb{R}})$ be an arbitrary vector field and let $\mu_{\gamma}=\partial_{x}(\rho_{1,\gamma}\,U_{\gamma,x})$, where the subscript $\gamma$ denotes convolution with a standard (compactly supported) $C^{\infty}$-mollifier. Then

[TABLE]

is well defined, since $w$ is Lipschitz continuous. For $\gamma\searrow 0$ one obtains the first equation of (13) in weak form. The second equation follows similarly.

In order to prove the last assertion in the statement, we work by contradiction. Let $\rho_{1}$ and $\rho_{2}$ be as assumed and let a direction $\mu\in L^{2}$, such that $(\rho_{1}+\delta\mu,\rho_{2})\in\mathcal{M}$ for which, e.g.

[TABLE]

By a density argument, one finds a function $U\in C^{2}_{c}({\mathbb{R}})$ such that

[TABLE]

[TABLE]

which is well defined since $\rho_{1}+\rho_{2}$ is Lipschitz. Thus $(\rho_{1},\rho_{2})$ cannot be a steady state. ∎

Next, introduce a technical result that will be useful to investigate sufficient conditions for a steady state to be a local minimizer of $\mathcal{F}$ in the next subsection.

Lemma 3.1**.**

The second Gateaux derivatives for ${\mathcal{F}}$ on $(\rho_{1},\rho_{2})$

[TABLE]

are given by

[TABLE]

where $\mu,\nu\in L^{2}({\mathbb{R}})^{2}$ are arbitrary and such that $(\rho_{1}+\delta\mu,\rho_{2})$ and $(\rho_{1},\rho_{2}+\delta\nu)$ are in $\mathcal{M}$.

Proof.

Computing the upper left element of $H[\mu,\nu]$ gives

[TABLE]

Therefore, we obtain

[TABLE]

The other entries of $H[\mu,\nu]$ can be computed similarly. ∎

3.1. Non-Existence of Steady States

We now establish a simple necessary condition for the existence of non trivial steady states, based on the idea that ${\mathcal{F}}$ is strictly convex when the diffusion part is dominant. Thus zero is the unique global minimizer.

Lemma 3.2**.**

Assume that the Fourier transforms of the interaction kernels satisfy

[TABLE]

for all $\xi\in{\mathbb{R}}$. Then there does not exist any global minimiser $(\rho_{1},\rho_{2})\in\mathcal{M}$.

Proof.

Thanks to assumption (A3), there exists $C>0$ such that

[TABLE]

Applying the Fourier transform we get

[TABLE]

[TABLE]

Therefore, conditions (15) imply that the above quadratic form is positive definite, and hence ${\mathcal{F}}[\rho_{1},\rho_{2}]>0$ for all $\rho_{1},\rho_{2}\in L^{2}({\mathbb{R}})$. Now, assume by contradiction that there exists a global minimiser $(\rho_{1},\rho_{2})$ with the prescribed mass constraints. Then ${\mathcal{F}}[\rho_{1},\rho_{2}]>0$. We now rescale $(\rho_{1},\rho_{2})$ by a parameter $\lambda>0$ in such a way to preserve the total mass of both components, i.e.

[TABLE]

notice that $(\rho_{1,\lambda},\rho_{2,\lambda})\in\mathcal{M}$. We can use (16) as follows

[TABLE]

The latter right-hand side converges to zero as $\lambda\rightarrow+\infty$. Therefore, for a large enough $\lambda$ the value of the functional on $(\rho_{1,\lambda},\rho_{2,\lambda})$ can be made smaller than ${\mathcal{F}}[\rho_{1},\rho_{2}]$, thus contradicting the fact that $(\rho_{1},\rho_{2})$ is a nontrivial global minimiser. ∎

We now establish a reasonable necessary condition for the existence of a nontrivial steady state.

Theorem 3.1**.**

Under conditions (15), there exists no nonzero stationary solution for (13).

Proof.

Let $(\rho_{1},\rho_{2})\in\mathcal{M}$. We use Lemma 3.1 to compute the second derivative $H[\mu,\nu]$ of ${\mathcal{F}}[\rho_{1},\rho_{2}]$ in the direction $(\mu,\nu)$. If

[TABLE]

then $\det(H[\mu,\nu])>0$. Indeed, applying the Fourier transform on both sides of the inequality above, we obtain the equivalent condition

[TABLE]

Now from (15)$(ii)$ we deduce

[TABLE]

which implies $\det(H[\mu,\nu])>0$. From (15)$(i)$ one deduces that the second derivative of ${\mathcal{F}}$ in (arbitrary) direction $(\mu,\nu)$ is positive definite. Hence, ${\mathcal{F}}$ is strictly convex on $L^{2}({\mathbb{R}})\times L^{2}({\mathbb{R}})$. Now, assume $(\bar{\rho}_{1},\bar{\rho}_{2})$ is a non trivial steady state. Since ${\mathcal{F}}$ is strictly convex, then the critical point $(\bar{\rho}_{1},\bar{\rho}_{2})$ is also the unique global minimiser. This contradicts Lemma 3.2. ∎

Remark 3.1**.**

Conditions (15) generalise the diffusion-dominated condition established in [22] for the one species case. In that case, the only interaction kernel $S$ should satisfy $\|S\|_{L^{1}}>{\epsilon}$ in order to ensure the existence of a non trivial steady state, or equivalently the condition ${\epsilon}\geq\|S\|_{L^{1}}$ implies that no global minima exist with fixed positive mass. In the case with two species, the two conditions $\|S_{i}\|_{L^{1}}<{\epsilon}$ $i=1,2$, ensure that

[TABLE]

and therefore the trace condition in (15) is trivially satisfied. This condition has a similar interpretation to the one in [22] for the one species case, i. e. the diffusion part is stronger than the (self) attraction parts, and so the spreading behaviour of particles dominates in the large time dynamics. However, (15)$(ii)$ is a more specific feature of the two species case. In order to provide a heuristic interpretation of such condition, let us consider for simplicity Gaussian potentials of the form

[TABLE]

for some positive constants $A,B,\sigma,\lambda$. Applying the Fourier transform we see that condition (i) is equivalent to

[TABLE]

Condition (15)$(ii)$ is satisfied e.g. if

[TABLE]

i.e. if $K$ has a larger mass than $S$ and the variance of $K$ is much smaller compared to that of $S$. Roughly speaking, this means that cross interaction should be relevant only at very small distances between particles of different species. This is not surprising. Consider for instance two separated patterns for $\rho_{1}$ and $\rho_{2}$, and assume that condition (15)$(i)$ is satisfied. Then, the lack of cross diffusion (which only appears when both $\rho_{1}$ and $\rho_{2}$ are positive) pushes the two patterns towards ‘spreading’. The only effect that could prevent the diffusive behaviour to dominate is the cross-interaction. Now, assume that (15)$(ii)$ is satisfied, for instance in the form (18). Such a condition somehow ensures that the cross-interaction potential is not exerting any long-range confining effect on the particles to compensate the diffusive effect, because it only acts at small distances. When the two patterns touch each other, the cross-interaction potential will only interfere within the mixing area (still because of (15)$(ii)$ ), still not enough to produce a ‘global confinement’ and prevent the whole solution from decaying.**

The above remark pinpoints an important fact about the minimization problem for $\mathcal{F}$: when the two species are separated and far from each other, the only coupling mechanism between them is the cross-interaction term ruled by the potential $K$. Hence, if we impose the conditions ensuring that the two species would feature a nontrivial global minimum in absence of any coupling (${\epsilon}<A$ in the above example), we immediately see that the functional can be diminished by decreasing the distance between the two species. Hence, if the two species are separated at the global minimum state, common sense would suggest that their supports are glued together. Proving that this is actually what happens is the main goal of this paper, and is tackled in the next sections from different viewpoints.

3.2. Structure of Energy Minimizers

The previous interpretation of the supports of $\rho_{1}$ and $\rho_{2}$ staying separated (possibly gluing together at one edge of their boundaries) is supported by the next result, which refers to a specific case, namely $K$ being sufficiently smooth with unique maximum at zero, and $S_{i}=\sigma_{i}K$, $i=1,2$, with positive real numbers $\sigma_{i}$. Under these conditions, provided that $\sigma_{1}+\sigma_{2}>2$, we can prove that the supports of the two components must necessarily have an empty intersection at local minimisers of ${\mathcal{F}}$, which corresponds to segregation. The next theorem is true also in multiple dimensions, but for simplicity we only state it in one space dimension.

Theorem 3.2**.**

Assume that $K\in C({\mathbb{R}})$ with a unique maximum at zero and being strictly radially decreasing in a neighbourhood. Let $S_{i}=\sigma_{i}K$ for some $\sigma_{1},\sigma_{2}>0$ with $\sigma_{1}+\sigma_{2}>2$. Let $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ be a local minimizer of ${\mathcal{F}}$ on $\mathcal{M}$, and

[TABLE]

Then $\mathcal{S}$ has zero Lebesgue measure.

Proof.

Assume that $\mathcal{S}$ has positive Lebesgue measure, then we can choose a subset $\mathcal{O}$ of positive measure and some constant $C>0$ such that

[TABLE]

For a fixed $d>0$ (e.g. half the diameter of $\mathcal{O}$) we can choose $z_{1},z_{2}\in\mathcal{O}$ with distance greater or equal than $d$ and $0<\gamma<{d}/4$ sufficiently small such that the balls $B_{\gamma}(z_{i})$ of radius $\gamma$ around $z_{1}$ and $z_{2}$ intersected with $\mathcal{O}$ have positive Lebesgue measure and are disjoint. Define

[TABLE]

and $u_{2}=-u_{1}$. Since $u_{1},u_{2}\in L^{2}({\mathbb{R}})$ have mean zero, for $\delta>0$ sufficiently small, $\rho_{i}^{\infty}\pm\delta u_{i}$ is nonnegative and has the same mass as $\rho_{i}^{\infty}$. Hence $(\rho_{1}^{\infty}\pm\delta u_{1},\rho_{2}^{\infty}\pm\delta u_{2})$ is admissible for the minimization of ${\mathcal{F}}$ with $\delta<<1$. It is hence straightforward to see that at a minimizer the first variation of ${\mathcal{F}}$ in direction $(u_{1},u_{2})$ vanishes. With the homogeneity of ${\mathcal{F}}$, we find

[TABLE]

Now, due to $u_{1}+u_{2}=0$ we obtain

[TABLE]

And we have

[TABLE]

For $\gamma$ sufficiently small, the local strict radial decrease around zero implies that the infimum is larger than the supremum above (note that $2\gamma<d-2\gamma$) and thus, $\int_{{\mathbb{R}}^{d}}u_{1}K*u_{1}~{}dx$ is positive. Hence, ${\mathcal{F}}[u_{1},u_{2}]<0$ and $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ cannot be a minimizer of ${\mathcal{F}}$. ∎

In the case of one weaker interaction, e.g. $\sigma_{2}<1$ one can give an improved result on the support by following the proof of [31], where only the attractive part of the energy with an upper bound on the total density $\rho_{1}+\rho_{2}$ was considered. The idea is to write the energy as a functional of $\rho_{1}+\rho_{2}$ and $\rho_{2}$ and apply Riesz-Sobolev symmetrization to these functions in order to obtain states of lower energy. Since such a symmetrization leaves the $L^{2}$-norm of $\rho_{1}+\rho_{2}$ unchanged it can be applied directly to our case:

Theorem 3.3**.**

Assume that $K\in C({\mathbb{R}})$ is strictly radially decreasing. Let $\sigma_{1},\sigma_{2}>0$ with $\sigma_{1}+\sigma_{2}>2$ and $\sigma_{2}<1$. Let $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ be a global minimiser with compact support of ${\mathcal{F}}$ on $\mathcal{M}$ with $m_{1}=m_{2}$. Then $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ is radially symmetric around its center of mass $X_{0}$ and there exist $b>a>0$ such that

[TABLE]

Of course the phase separation does not mean that the supports are at positive distance, due to the attractive cross interactions we expect that the supports are glued together. This is indeed true in general for stationary solutions. The proof is the same as the corresponding result for a single species in [22, Lemma 4.1], we hence only provide a sketch of the proof.

Proposition 3.2**.**

Let $\sigma_{1},\sigma_{2}>0$, $K$ be strictly decreasing with respect to the radial variable and $\left(\rho_{1}^{\infty},\rho_{2}^{\infty}\right)\in\mathcal{M}\cap BV({\mathbb{R}})^{2}$ be a weak solution of (13) according to Definition 3.1. Then $supp(\rho_{1}^{\infty}+\rho_{2}^{\infty})$ is a connected interval.

Sketch of the Proof.    Let us just mention the main idea of the proof: Assume there exists an interval $[a,b]$ with $a<b$ not in the support of $\rho_{1}^{\infty}+\rho_{2}^{\infty}$, but such that the intersections of $(-\infty,a)$ and $(b,+\infty)$ with the support of $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ are both nonempty. Then we can construct a smooth velocity $V$ field equal to $+1$ on $(-\infty,a)$ and $-1$ on $(b,\infty)$. Consider then the infinitesimal change of the energy on the curve $(u_{1}(s,x),u_{2}(s,x))$ solving (locally) the Cauchy problem

[TABLE]

for $i=1,2$. Since $\frac{d}{ds}\mathcal{F}[u_{1}(s),\rho_{2}^{\infty}]|_{s=0}=\frac{d}{ds}\mathcal{F}[\rho_{1}^{\infty},u_{2}(s)]|_{s=0}=0$ in view of the fact that $\left(\rho_{1}^{\infty},\rho_{2}^{\infty}\right)$ is a stationary state according to Definition 3.1, due to the computation

[TABLE]

with a few manipulations we obtain

[TABLE]

Now, since all the involved kernels $S_{1},S_{2},K$ are even, with a few manipulations similar to those in [22, Lemma 4.1], we get that the support of $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ must be empty in at least one between $(-\infty,a)$ and $(b,+\infty)$, which is a contradiction. This proves that $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ cannot be a stationary state according to Definition 3.1. ∎

We can provide the same result for an arbitrary minimizer of the interaction energy with a similar proof, respectively the Lagrangian version:

Proposition 3.3**.**

Let $\sigma_{1},\sigma_{2}>0$, $K$ be strictly decreasing with respect to the radial variable and $\left(\rho_{1}^{\infty},\rho_{2}^{\infty}\right)\in\mathcal{M}$ be a minimizer of ${\mathcal{F}}$. Then $supp(\rho_{1}^{\infty}+\rho_{2}^{\infty})$ is a connected interval.

Proof.

Suppose that the difference between $supp(\rho_{1}^{\infty}+\rho_{2}^{\infty})$ and its convex hull contains an interval $[a,b]$. Let $0<d<\frac{b-a}{2}$ and define $(\tilde{\rho}_{1},\tilde{\rho}_{2})\in\mathcal{M}$ via

[TABLE]

Then we have

[TABLE]

while the attractive interaction terms in the energy are smaller for $(\tilde{\rho}_{1},\tilde{\rho}_{2})$. Hence, ${\mathcal{F}}[\tilde{\rho}_{1},\tilde{\rho}_{2}]<{\mathcal{F}}[\rho_{1}^{\infty},\rho_{2}^{\infty}]$, which is a contradiction. ∎

Despite the result of compact support for $\rho_{1}^{\infty}+\rho_{2}^{\infty}$, which holds for arbitrary positive interactions, we expect mixing in case of dominant cross-attraction, i.e. $\sigma_{1}<1$ and $\sigma_{2}<1$ similar to [31].

Note that Proposition 3.3 is not a special case of Proposition 3.2, since we cannot argue that $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ are Lipschitz-continuous, hence we cannot apply the result Proposition 3.1. In order to do so we will need to verify the Lipschitz property, for which Proposition 3.3 appears to be crucial. In this direction we now provide a result closing the gap between energy minimizers and stationary solutions. At least for $\sigma_{1}+\sigma_{2}>2$ and $\sigma_{2}<1$ the assumptions of the following theorem are satisfied:

Theorem 3.4**.**

Assume $K\in C^{1}({\mathbb{R}})$ to be strictly radially decreasing. Let $\sigma_{1},\sigma_{2}>0$ and let $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ be a global minimiser with compact support of ${\mathcal{F}}$ on $\mathcal{M}$. Moreover, assume that for some $M>1$ there exist $a_{1}<a_{2}<\ldots<a_{M}$ such that supp$(\rho_{1}^{\infty}+\rho_{2}^{\infty})=\cup_{i=1}^{M-1}[a_{i},a_{i+1}]$ and that for each $i=1,\ldots,M-1$ we have

[TABLE]

with $\{k,\ell\}=\{1,2\}$. Then $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ is Lipschitz-continuous.

Proof.

Let $i$ be such that $k=1$ and $\ell=2$ (the opposite case is analogous). Moreover, let $\varphi\in C({\mathbb{R}})$ be an arbitrary continuous function supported in a compact subset $(a_{i},a_{i+1})$. Then by standard arguments on the first variation we see that ${\mathcal{F}}[\varphi,\rho_{2}^{\infty}]=0$, which implies

[TABLE]

on $(a_{i},a_{i+1})$. The continuous differentiability of $K$ implies that $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ is Lipschitz continuous (with modulus independent of $i$) on $(a_{i},a_{i+1})$. Hence, it suffices to show that $\rho_{1}^{\infty}+\rho_{2}^{\infty}$ is continuous on the finite set of points $\{a_{i}\}$ in order to obtain Lipschitz continuity on ${\mathbb{R}}$. Due to the uniform Lipschitz continuity on the subintervals, the left- and right-sided limits exist for each subinterval.

First, let $\rho_{k}^{\infty}>0$ in $(a_{1},a_{2})$ and assume $\lim_{x\rightarrow a_{1}^{+}}\rho_{k}^{\infty}(x)>0$. For $\delta_{1},\delta_{2}>0$ being sufficiently small define

[TABLE]

Then it is straight-forward to show that

[TABLE]

which strictly smaller than ${\mathcal{F}}[\rho_{1}^{\infty},\rho_{2}^{\infty}]$ for $\delta_{1},\delta_{2}$ sufficiently small and hence a contradiction to $(\rho_{1}^{\infty},\rho_{2}^{\infty})$ being a minimiser. Thus,

[TABLE]

A completely analogous proof shows continuity at $x=a_{M}$. Finally consider the limit at $x=a_{i}$, $1<i<M$, if $\lim_{x\rightarrow a_{i}^{+}}\rho_{k}^{\infty}(x)=\lim_{x\rightarrow a_{i}^{+}}\rho_{\ell}^{\infty}(x)$. Define a similar perturbation in a neighbourhood $(a_{i}-\delta_{2},a_{i}+\delta_{2})$ increasing the smaller and decreasing the larger density. With the same kind of expansion as at $a_{1}$ one can obtain a state with lower energy, which is a contradiction. ∎

4. Existence of segregated states via Krein-Rutman theorem

Hinted by the results in the previous section, we now provide reasonable sufficient conditions for the existence of segregated stationary states for system (13). We remark that the system (13) is translation invariant. Therefore, we shall fix the center of mass to zero for simplicity.

Definition 4.1**.**

We call a stationary solution $(\rho_{1},\rho_{2})$ to (13) according to Definition 3.1 a symmetric segregated steady state if there exist $L_{1},L_{2}>0$ such that

[TABLE]

where $\rho_{1},\rho_{2}$ are even functions, and $C^{1}$-regular inside their respective supports, such that $w=\rho_{1}+\rho_{2}$ is monotone decreasing on $[0,L_{2}]$ with $w(L_{2})=0$.

Recall that Definition 3.1 prescribes $w(x):=\rho_{1}(x)+\rho_{2}(x)$ being Lipschitz continuous on ${\mathbb{R}}$.

Our strategy to prove existence of segregated steady states is a non-trivial extension of the strategy proposed in [22], and adapted in [23]. The main idea is to fix $L_{1}$ and $L_{2}$ and then look at the stationary equations for $\rho_{1}$ and $\rho_{2}$ as eigenvalue conditions for a suitable integral operator. The existence of the eigenvectors will be proven by using the strong version of the Krein-Rutman theorem.

Theorem 4.1** (Krein-Rutman).**

Let $X$ be a Banach space, $K\subset X$ be a solid cone, such that $\lambda K\subset K$ for all $\lambda\geq 0$ and $K$ has a nonempty interior $K^{o}$. Let $T$ be a compact linear operator on X, which is strongly positive with respect to $K$, i.e. $T[u]\in K^{o}$ if $u\in K\setminus\{0\}$. Then,

• (i)

the spectral radius $r(T)$ is strictly positive and $r(T)$ is a simple eigenvalue with an eigenvector $v\in K^{o}$. There is no other eigenvalue with a corresponding eigenvector $v\in K$.

• (ii)

$|\lambda|<r(T)$* for all other eigenvalues $\lambda\neq r(T)$.*

We shall work with the following class of kernels:

[TABLE]

which, together with our nonnegativity assumption (A3) imply in particular that all the kernels are supported on ${\mathbb{R}}$.

Moreover, having fixed $L_{1}<L_{2}$, we assume for all $x\in(0,L_{1})$ that

[TABLE]

[TABLE]

for all $x\in(L_{1},L_{2})$. Assumptions (22) and (23) are met for instance in the significant case $S_{i}=\sigma_{i}K$, $\sigma_{1}>1>\sigma_{2}$. We assume further that

[TABLE]

Let $\rho_{1}$ and $\rho_{2}$ be symmetric steady states, then we can rephrase (13) as

[TABLE]

for suitable constants $C_{1},C_{2}>0$. Assuming segregation, we rewrite (25) as

[TABLE]

Let $\bar{w}=\rho_{1}(L_{1})=\rho_{2}(L_{1})$. Let $p(x)=-\rho_{1}^{\prime}(x)$ be restricted to the interior of $I_{1}$ and $q(x)=-\rho_{2}^{\prime}(x)$ to the interior of $I_{2}$. By differentiating (26) we obtain

[TABLE]

where

[TABLE]

A symmetrization yields for $x>0$, that

[TABLE]

To simplify notation, define the non-negative functions, for $x,y>0$,

[TABLE]

Then we can rewrite system (27) as

[TABLE]

The result of this section is stated in the following Proposition.

Proposition 4.1**.**

Let assumptions (A1), (A3), (21), (22), (23), and (24) be satisfied. Let $0<L_{1}<L_{2}$ be fixed. Then, there exists a unique (up to mass normalisation) symmetric segregated steady state $\left(\rho_{1},\rho_{2}\right)$ to (13) with ${\epsilon}={\epsilon}(L_{1},L_{2})>0$.

Proof.

Consider the Banach space

[TABLE]

equipped with the $W^{1,\infty}$-norm for the first two components and with the standard one-dimensional Euclidean norm for the third component. We use the notation $P=(p,q;w)$ for all $P\in X_{L_{1},L_{2}}$, with $p\in C^{1}([0,L_{1}])$, $q\in C^{1}([L_{1},L_{2}])$, and $w\geq 0$. For a given $P\in X_{L_{1},L_{2}}$ define $T_{L_{1},L_{2}}[P]\in X_{L_{1},L_{2}}$ as

[TABLE]

[TABLE]

The operator $T_{L_{1},L_{2}}$ is compact on the Banach space $X_{L_{1},L_{2}}$ as all the involved kernels are $C^{2}$ on compact intervals, and hence by Arzelá’s theorem the image of the unit ball in $X_{L_{1},L_{2}}$ is pre-compact. It is easy to show that the set

[TABLE]

is a solid cone in $X_{L_{1},L_{2}}$ and that

[TABLE]

where $K_{L_{1},L_{2}}^{o}$ denotes the interior of $K_{L_{1},L_{2}}$. We show that $T$ is strongly positive on the solid cone. First observe that, due to (22) and (23), all components of $T_{L_{1},L_{2}}[(p,q;w)](x)$ are nonnegative if $(p,q;w)\in K_{L_{1},L_{2}}$ with $(p,q;w)\neq 0$. Moreover, setting $f(x)$ as the first component of $T_{L_{1},L_{2}}[(p,q;w)](x)$ as above, we get

[TABLE]

since $(p,q;w)\neq 0$. Hence, $T_{L_{1},L_{2}}[(p,q,w)]$ belongs to $H_{L_{1},L_{2}}$ and therefore the operator is strongly positive on the cone $K_{L_{1},L_{2}}$. The strong version of the Krein-Rutman theorem (Theorem 4.1) then gives the existence of an eigenvalue ${\epsilon}={\epsilon}(L_{1},L_{2})$ such that

[TABLE]

with eigenspace generated by a vector $(q,p,w)\in K_{L_{1},L_{2}}^{o}$. Now, let $\rho_{2}(x)=\int_{x}^{L_{2}}q(x)dx$. The second and third equation in the eigenvalue condition imply $\bar{w}=\rho_{2}(L_{1})$. Let $\rho_{1}(x)=\bar{w}+\int_{x}^{L_{1}}p(x)dx$. Clearly, $\rho_{1}(L_{1})=\rho_{2}(L_{1})$. Then, a computation similar to the one that led to the definition of the operator $T_{L_{1},L_{2}}$ after differentiating (26) implies that $\rho_{1}$ and $\rho_{2}$ are symmetric segregated steady states to (13). ∎

Remark 4.1**.**

The above result shows that a segregated steady state $(\rho_{1},\rho_{2})$ exists for arbitrary positive constants $L_{1}$ and $L_{2}$ such that $\mathrm{supp}(\rho_{1})=[-L_{1},L_{1}]$ and $\mathrm{supp}(\rho_{2})=[-L_{2},-L_{1}]\cup[L_{1},L_{2}]$ for some diffusion constant $\epsilon>0$ which depends on $L_{1}$ and $L_{2}$. Clearly, we obtain a one-parameter family of solutions upon multiplication of $(\rho_{1},\rho_{2})$ by an arbitrary positive constant. Such a constant is uniquely determined by the total mass $m=m_{1}+m_{2}$ with $m_{1}=\int_{\mathbb{R}}\rho_{1}(x)dx$ and $m_{2}=\int_{\mathbb{R}}\rho_{2}(x)dx$. The two values $m_{1}$ and $m_{2}$ are not determined explicitly here. We expect the value of $L_{1}$ to be determined uniquely once $L_{2}$ and the two masses $m_{1}$ and $m_{2}$ are fixed. Most importantly, this approach does not provide an explicit information on the value of the diffusion constant $\epsilon$. A similar situation occurs in the problem studied in [23] for the one species equation with general power-law diffusion. Reasoning as in [22], we expect that the diffusion constant $\epsilon$ obeys a monotone increasing law of the form $\epsilon=\overline{\epsilon}(L_{2})$. If such a conjecture is true, then for any given $\epsilon$ in the range of the map $\overline{\epsilon}$ and for any given $m_{1},m_{2}>0$ there exist unique $L_{1},L_{2}>0$ and a unique (up to $x$-translations) segregated state $(\rho_{1},\rho_{2})$ with $m_{1}=\int_{\mathbb{R}}\rho_{1}(x)dx$, $m_{2}=\int_{\mathbb{R}}\rho_{2}(x)dx$, $\mathrm{supp}(\rho_{1})=[-L_{1},L_{1}]$, and $\mathrm{supp}(\rho_{2})=[-L_{2},-L_{1}]\cup[L_{1},L_{2}]$. Such a conjecture will be addressed in full generality in a future study. In the next section we are able to prove that such a statement is true for a small enough diffusion coefficient.**

5. Existence and uniqueness of steady states for small diffusion

Here we prove existence and uniqueness of a symmetric segregated steady state for fixed masses, and a fixed small diffusion coefficient. We refer to Definition 4.1 to denote segregated states. Similar to [21], we formulate the problem in the pseudo-inverse formalism and then use an implicit function theorem argument. Throughout this section we shall require

[TABLE]

The main result of this section reads as follow.

Theorem 5.1** (Existence of segregated steady states for small diffusion).**

Assume that the interaction kernels fulfil the additional assumption (30). Then, there exists a constant $\epsilon_{0}$ such that for all $\epsilon\in(0,\epsilon_{0})$ the stationary equation (13) admits a unique solution in the sense of Definition 4.1 with fixed masses $m_{1}=\int\rho_{1}dx$ and $m_{2}=\int\rho_{2}dx$.

5.1. Formulation in the pseudo-inverse variable

Assume that $(\rho_{1},\rho_{2})$ is a segregated state with the structure as in Definition 4.1, and

[TABLE]

Hence, $w=\rho_{1}+\rho_{2}$ has unit mass and is supported on $[-L_{2},L_{2}]$. Let

[TABLE]

Let $u:[0,1)\rightarrow{\mathbb{R}}$ be the pseudo-inverse of $F$

[TABLE]

Set $u_{i}(z)=u(z)\mathbf{1}_{J_{i}}(z)$, $i=1,2$, with

[TABLE]

[TABLE]

Then (13) can be rewritten as

[TABLE]

Recall that $\rho_{1}$ and $\rho_{2}$ are symmetric, which implies that the pseudo-inverse $u$ satisfies $u(1-z)=-u(z)$. Moreover, since $w$ is strictly positive on $(-L_{2},L_{2})$, $u$ is strictly increasing on $(0,1)$ and Lipschitz continuous on the compact subintervals of $(0,1)$. Due to $w$ being zero at $z=0$ and $z=1$, $u$ is expected to have an infinite slope at the boundary.

Since $w=\rho_{1}+\rho_{2}$ and since $\rho_{1}$, $\rho_{2}$ have disjoint supports, we have

[TABLE]

The idea is to solve (31) for small ${\epsilon}$, hinted by the fact that the case ${\epsilon}=0$ has the unique solution $u_{1}\equiv u_{2}\equiv 0$, which corresponds to $\rho_{1}$ and $\rho_{2}$ being two Dirac’s deltas with masses $z_{1}$ and $1-z_{1}$ respectively. Similar to [20], we expect that the support of $w$ is small for small ${\epsilon}$. This suggests the linearisation $u_{i}=\delta v_{i}$, $i=1,2$, with $v_{1}$, $v_{2}$ being odd functions defined on $J_{1}$, $J_{2}$ respectively. A simple scaling argument suggests that $\delta={\epsilon}^{\frac{1}{3}}$, and therefore

[TABLE]

Multiplying the first equation by $\delta\partial_{z}v_{1}$ and the second one by $\delta\partial_{z}v_{2}$, we get

[TABLE]

Taking the primitives w.r.t. $z$, we obtain for $z\in J_{1}$, respectively $z\in J_{2}$ that

[TABLE]

[TABLE]

with integration constants $\alpha_{1}$, $\alpha_{2}$, which are obtained by substituting $z=1$ into (33)

[TABLE]

and imposing the continuity condition for $(\partial_{z}u)^{-1}$ in $\tilde{z}=(1+z_{1})/2$

[TABLE]

The symmetry requirements imply

[TABLE]

Now we multiply (32) by $\delta\partial_{z}v_{1}$, (33) by $\delta\partial_{z}v_{2}$, integrate w.r.t. $z$ and introduce $G_{i},H$ such that $G_{i}^{\prime}=S_{i}$, $i=1,2$, and $H^{\prime}=K$. Here $G_{1},G_{2},H$ can be chosen odd with $G_{1}(0)=G_{2}(0)=H(0)$. The integration constants can be recovered by prescribing

[TABLE]

After some manipulations we obtain

[TABLE]

as well as

[TABLE]

A solution $(v_{1},v_{2})$ to (5.1), (5.1) should be extended to $[0,1/2]$ to obtain an odd profile on the whole interval $[0,1]$.

Remark 5.1**.**

The functional system (5.1)-(5.1) has been obtained manipulating the original system (32)-(33) through integrations w.r.t. the independent variable $z$ and dividing by $\delta^{2}$. Therefore, seeking for a solution to (5.1)-(5.1) for $\delta=0$ somehow involves the computation of the first order term in the expansion of the r.h.s. of (32)-(33) with respect to $\delta^{2}$. This computation will be performed in the next subsection. In some sense, in terms of the linearisation $u_{i}=\delta v_{i}$, we are computing the first order term w.r.t. $\delta$ of $u_{i}$ near $\delta=0$.**

5.2. A functional equation

System (5.1), (5.1) can be viewed as a functional equation in the following form. Introduce the Banach space

[TABLE]

And consider the standard norm on $\Omega$

[TABLE]

For $\alpha>0$, consider the norm

[TABLE]

and set $\Omega_{\alpha}:=\left\{(v_{1},\mu,v_{2},\lambda)\in\Omega\,:\,\,\interleave(v_{1},\mu,v_{2},\lambda)\interleave_{\alpha}<+\infty\right\}$.

We now define the convex subset $U_{\alpha}\subset\Omega$

[TABLE]

In order to formulate (5.1), (5.1) only on $[1/2,1]$, we again use the symmetrised function

[TABLE]

for given $G:{\mathbb{R}}\rightarrow{\mathbb{R}}$. Moreover, we use the notation $\tilde{J}_{i}=J_{i}\cap[1/2,1]\,,\quad i=1,2$. Let $(v_{1},\mu,v_{2},\lambda)\in\Omega$, and $\delta>0$. Define

[TABLE]

[TABLE]

Here $F_{1}$ and $F_{2}$ also depend on the $z$-variable. For simplicity, we drop this dependence throughout the rest of the section. Substituting $z=\tilde{z}$ in (37) and $z=1$ in (38), we define

[TABLE]

by substituting $v_{1}(\tilde{z})=\mu$ and $v_{2}(1)=\lambda$ in both expressions.

Define the extension of $F_{1}$ and $F_{2}$ to $\delta=0$ by using Taylor expansion together with the symmetry properties of the involved kernels. Let

[TABLE]

Due to the assumptions on $S_{1}$, $S_{2}$, and $K$, the constants $S_{1}^{\prime\prime}(0)$, $S_{2}^{\prime\prime}(0)$ and $K^{\prime\prime}(0)$ are non-positive. After some tedious calculations we obtain the natural definition

[TABLE]

Our goal is to solve

[TABLE]

for $\delta>0$ small enough. First we will solve the case $\delta=0$ and then prove that a solution still exists when $\delta$ is close to zero. The solution for $\delta=0$ is already partially explicit from formulas (39), (40). We only need to determine $\mu$ and $\lambda$. To do so, we set the $m$ and $\ell$ components above equal to zero and obtain

[TABLE]

which should be solved under the constraint $\mu<\lambda$. The second equation in (41) can be rewritten as

[TABLE]

which describes a cubic hyperbola, asymptotic to the straight lines $\mu-\lambda=0$ and $2\lambda+\mu=0$ in the $(\mu,\lambda)$ plane. Consider the branch of such a curve in the region $0<\mu<\lambda$, intersecting the $\mu=0$ axis at $(\mu,\lambda)=(0,\lambda_{0})$ with

[TABLE]

Such a branch describes a monotone increasing function $\lambda=\bar{\lambda}(\mu)$ on $\mu>0$, asymptotic to $\lambda=\mu$ as $\mu\rightarrow+\infty$. Now, summing up the two equations in (41) we get the following additional condition on $\lambda$ and $\mu$:

[TABLE]

The function $\lambda=\tilde{\lambda}(\mu)$ is also monotone, and attains the value

[TABLE]

On the other hand, $\tilde{\lambda}(\mu)$ is asymptotic to the straight line $\lambda=((C_{2}-C_{1})/C_{2})^{1/3}\mu$ as $\mu\rightarrow+\infty$. Since $((C_{2}-C_{1})/C_{2})^{1/3}<1$, the two curves $\lambda=\bar{\lambda}(\mu)$ and $\lambda=\tilde{\lambda}(\mu)$ intersect at exactly one point in the region $0<\mu<\lambda$. Hence, $\mu$ and $\lambda$ are uniquely determined.

At this stage, the solution $(v_{1},v_{2})$ can be easily recovered as the pseudo-inverse variables associated to the densities $\tilde{\rho}_{1}$ and $\tilde{\rho}_{2}$

[TABLE]

which corresponds to two Barenblatt profiles centered at $x=0$ (with possibly different masses and supports) such that the resulting profile $w=\tilde{\rho_{1}}+\tilde{\rho_{2}}$ is continuous at $x=\pm\mu$. This proves that

[TABLE]

has a unique solution $(v_{1},\mu,v_{2},\lambda)\in\Omega$. From now on we call such solution $(v_{0,1},\mu_{0},v_{0,2},\lambda_{0})$.

Remark 5.2** (Inner and outer species).**

In the above computation we arbitrarily assigned $\rho_{1}$ to the role of the ‘inner’ species in the segregated state. All of the above procedure also works with $\rho_{2}$ as inner species, and therefore we should speak of two possible segregated states rather than just one. The only parameters characterizing each species with respect to this pattern are $S_{1}^{\prime\prime}(0)$ and $S_{2}^{\prime\prime}(0)$, which model self-attraction forces. These constants affect the constants $C_{1}$ and $C_{2}$ above, which are non-negative. When $C_{2}<C_{1}$, the self-attraction force of the first species is stronger than that of the second one, which suggests that $\rho_{1}$ will concentrate faster than $\rho_{2}$, thus producing a pattern in which the first species occupies the inner region whereas the second species occupies the outer region. On the other hand, a steady state with the reversed order can still be constructed under the condition $C_{2}<C_{1}$, but we conjecture it to be unstable. This is supported by numerical simulations shown in Section 6, in which we see that the inner species is the one featuring the ‘more concentrated’ Barenblatt-type profile. We observe that the slope of the function $\lambda=\tilde{\lambda}(\mu)$ above is negative if $C_{2}<C_{1}$ and positive if $C_{1}<C_{2}$. This implies that the values of $\mu$ and $\lambda$ solving (41) are smaller when $C_{2}<C_{1}$ compared to the case $C_{2}>C_{1}$. Hence, $\rho_{1}$ and $\rho_{2}$ have smaller support when $C_{2}<C_{1}$. We deduce that the ‘correct’ steady state is ‘more concentrated’ than the unstable one.**

Remark 5.3** (Non-symmetric segregation).**

As shown by the numerical simulations in Section 6, segregation may emerge via a non-symmetric structure, in which $\rho_{1}+\rho_{2}$ is supported on a connected interval and $\rho_{1}$ and $\rho_{2}$ feature exactly one jump discontinuity, see figure 2 below. This is typically the case for instance when the two species are initially separated. The mathematical proof of the existence of non-symmetric segregated steady states can be carried out similarly to the symmetric case. We omit the details and restrict to the symmetric case for simplicity.**

5.3. Solution via implicit function theorem for $\delta>0$

In this section we adapt the strategy of [21, Section 4] to our two-species problem. In order to solve problem (5.1), (5.1) for $\delta>0$ small enough, we analyse the operator

[TABLE]

Our first goal is to prove that $\mathcal{F}$ is Frechét continuously differentiable in a neighborhood of $(v_{1,0},\mu_{0},v_{2,0},\lambda_{0})$, the solution of $\mathcal{F}[(v_{1},\mu,v_{2},\lambda);0]=0$. By Taylor expansion in (5.1), (5.1) w.r.t. $\delta$ around $\delta=0$ one sees that $\mathcal{F}$ has a continuous first derivative with respect to $\delta$. We omit the details.

For fixed $\delta>0$, the Jacobian of $\mathcal{F}[\cdot;\delta]$ w.r.t. the first four variables has the block structure

[TABLE]

where partial derivatives with respect to $v_{1}$ and $v_{2}$ are meant to be Frechét derivatives. We now compute such terms. Consider perturbations $(w_{1},a,w_{2},b)\in\Omega_{1/2}$ such that

[TABLE]

We notice that (42) is satisfied if $\interleave(w_{1},a,w_{2},b)\interleave_{1/2}$ is small enough, since $v_{0,1}$ and $v_{0,2}$ have their gradient bounded from below by a positive constant. Upon extending $w_{1}$ and $w_{2}$ odd to $[1-\tilde{z},\tilde{z}]$ and $[0,1-\tilde{z}]\cup[\tilde{z},1]$ respectively, one has

[TABLE]

Now we compute the partial derivatives of $\mathcal{F}$. For simplicity, we drop the –0 indices to denote the $\delta=0$ state and avoid indicating the respective interval for the variable $z$.

Let $(w_{1},a,w_{2},b)\in\Omega_{1/2}$. By extending all involved functions to $[0,1/2]$ in the usual symmetric form, we get for $\delta>0$

[TABLE]

and in the limit $\delta\searrow 0$ we obtain

[TABLE]

This limit is so far just formal. The same holds for the $\delta\searrow 0$ limits computed below. However, we shall prove later that these are actually rigorous limits. Similarly,

[TABLE]

One can easily see that, in the $\delta\searrow 0$ limit one gets

[TABLE]

We now compute for $\delta>0$,

[TABLE]

and a Taylor expansion arguments shows

[TABLE]

Similarly,

[TABLE]

with the $\delta\searrow 0$ limit

[TABLE]

Turning to the second row of $D\mathcal{F}$ we get

[TABLE]

and a simple symmetry argument shows that

[TABLE]

Next we compute

[TABLE]

and a simple computation shows in the $\delta\searrow 0$ limit,

[TABLE]

Similar to the above, we have

[TABLE]

We conclude the second row of the Jacobian by computing

[TABLE]

with the $\delta\searrow 0$ limit being

[TABLE]

Now looking at the third row of $D\mathcal{F}$, we compute

[TABLE]

and for $\delta\searrow 0$ this shows

[TABLE]

Then, we have

[TABLE]

with the $\delta\searrow 0$ limit

[TABLE]

We continue computing the third row of $D\mathcal{F}$ with

[TABLE]

and the $\delta\searrow 0$ limit

[TABLE]

To conclude the third row of the Jacobian, we have

[TABLE]

and the $\delta\searrow 0$ limit is

[TABLE]

Finally, we analyse the last row of the Jacobian of $\mathcal{F}$.

[TABLE]

and the $\delta\searrow 0$ limit can be easily computed to be

[TABLE]

Then, we have

[TABLE]

with the $\delta\searrow 0$ limit being

[TABLE]

We then continue with

[TABLE]

and

[TABLE]

The last term is

[TABLE]

and the $\delta\searrow 0$ limit is

[TABLE]

The above computations show for $\delta$ small enough, that $D\mathcal{F}[(v_{0,1},\mu_{0},v_{0,2},\lambda_{0});\delta]$ is a bounded linear operator from $\Omega$ into itself, and that $D\mathcal{F}$ is continuous at $\delta=0$ in the operator norm. This is easily seen via Taylor expansion, using bounds on the $L^{\infty}$ norms, and symmetry properties.

Lemma 5.1**.**

$D\mathcal{F}[(v_{0,1},\mu_{0},v_{0,2},\lambda_{0});\delta]$* is a bounded linear operator from $\Omega_{1/2}$ to $\Omega_{1}$ for $\delta>0$ small enough.*

Proof.

Due to the structure of the spaces $\Omega_{\alpha}$, we only need to check the following. Define for $z\in J_{2}$ and $\delta\geq 0$,

[TABLE]

We need to prove that

[TABLE]

Consider all the above terms separately, and for notational reasons omit the dependence on $(v_{1},\mu,v_{2},\lambda)$. By Taylor expansion, and using simple symmetry properties, we get

[TABLE]

for some $\bar{x}(\zeta)$ and $\tilde{x}(\zeta)$ in the interval $[0,\lambda\delta]$. Now, it can be proved (cf. a similar argument in [21, Lemma 4.1]) that $(\lambda-v_{2}(z))/(1-z)^{1/2}$ is uniformly bounded in $[\tilde{z},1]$. This provides the desired estimate for the related term in (45). One can check that

[TABLE]

We now estimate

[TABLE]

and some tedious Taylor expansions imply

[TABLE]

This proves the assertion. ∎

Lemma 5.2**.**

$D\mathcal{F}[(v_{0,1},\mu_{0},v_{0,2},\lambda_{0});0]$* is a linear isomorphism between $\Omega_{1/2}$ and $\Omega_{1}$ for $\delta>0$ small enough.*

Proof.

We observe that the Jacobian of $\mathcal{F}$ at $\delta=0$ has the structure

[TABLE]

Given $(h_{1},\alpha,h_{2},\beta)\in\Omega_{1}$, we have to prove that

[TABLE]

admits a unique solution $(w_{1},a,w_{2},b)\in\Omega_{1/2}$ with

[TABLE]

for some $C>0$ independent of $(h_{1},\alpha,h_{2},\beta)$.

As a first step, we claim that $\frac{\partial F_{1}}{\partial v_{1}}$ is invertible as a map from $L^{\infty}$ to $L^{\infty}$ at $\delta=0$. To see this, we use (43). For $h_{1}(z)=\frac{\partial F_{1}}{\partial v_{1}}(w_{1})$ we get

[TABLE]

and the assertion follows. Therefore, the proof will be completed if we can show that the sub-matrix

[TABLE]

is an invertible operator in the components $(a,w_{2},b)$. First, we prove that

[TABLE]

which is equivalent to

[TABLE]

This is always satisfied since $\lambda>\mu$. Condition (47) implies that the linear system

[TABLE]

has a unique solution $(a,b)$. Now we only need to determine $w_{2}$. By subtracting the last two rows of the linear system (46), and by some simple manipulation, we obtain

[TABLE]

Since $(\lambda-v_{2}(z))/(1-z)$ is uniformly positive on $[\tilde{z},1]$ (cf. a similar proof in [21, Lemma 1.4]), we obtain the desired assertion, by dividing the above identity by $\sqrt{1-z}$. ∎

We are now ready to prove the main theorem of this section, Theorem 5.1, as well as one of the most important results in this paper.

Proof of Theorem 5.1.

The results in this section, in particular Lemma 5.1 and Lemma 5.2, together with the implicit function theorem on Banach spaces (see e. g. [32, Theorem 15.1]), imply that the functional equation $\mathcal{F}[(v_{1},\mu,v_{2},\lambda);\delta]=0$ has a solution for $\delta$ small enough. Here $\mathcal{F}$ is defined in Subsection 5.2 (see in particular (37)) and at the beginning of Subsection 5.3. Hence, we obtain a solution $(v_{1},v_{2})$ to (5.1)-(5.1). The computations in Subsection 5.1 imply that $(u_{1},u_{2})$ with $u_{i}=\delta v_{i}$, $i=1,2$, is a solution to (31) once we achieve enough regularity for the $v_{i}$. This follows easily from (5.1)-(5.1). Indeed, since $\delta$ is very small and since $S_{1}$, $S_{2}$, and $K_{0}$ are bounded away from zero, we can easily recover the $\partial_{z}v_{i}$ after differentiating (5.1)-(5.1) with respect to $z$. In particular, we find that the $\partial_{z}v_{i}$ are bounded away from zero. Hence, we can divide the resulting equations by $\partial_{z}v_{1}$ and $\partial_{z}v_{2}$ respectively. Similar to the above, we can once again differentiate (using the $C^{1}$ regularity of the $v_{i}$) and obtain (31) for $(u_{1},u_{2})$. The usual change of variable transforming pseudo inverse variables to densities shows that $(\rho_{1},\rho_{2})$ solves (13) where $\rho_{i}=\partial_{x}F_{i}$, $i=1,2$, and $F_{i}$ is the pseudo-inverse of $u_{i}$ for $i=1,2$. ∎

6. Numerical Simulations

Here we present some examples of sorting phenomena and mixing by solving (6) numerically in one space dimension. We use a particle method introduced for equations in gradient flow form in [29, 26], which is equivalent to a finite difference scheme for the pseudo-inverse equation, see also [44] and [37] for scalar conservation laws.

In this section we denote by $\rho$ and $\eta$ the two densities and by $u$ and $v$ the corresponding pseudo-inverse functions, that are solutions of the following system

[TABLE]

for $z\in\left[0,1\right]$. This is in order to avoid confusion w.r.t. the indices for discretization. Clearly, the masses of $\rho$ and $\eta$ are normalized to one. The main issue for the equations above is how to treat the cross-diffusion part numerically. We intentionally left the cross-diffusion terms above in the form of ‘external potentials’. Given $N\in\mathbb{N}$ we consider a partition of the interval $\left[0,1\right]$, $\left\{z_{i}\right\}_{i=1}^{N}$ and call $u(z_{i})=u_{i}$, $v(z_{i})=v_{i}$, and $m=\frac{1}{N}$. The discretization in space the reads

[TABLE]

for $i=1,...,N$. Integrating the above ODE system we get two families of particles $\left\{u_{i}(t)\right\}_{i=1}^{N}$, $\left\{v_{i}(t)\right\}_{i=1}^{N}$, for $t\in\left[0,T\right]$ and reconstruct the discrete densities as follows

[TABLE]

Given the two initial conditions $\rho^{0}$ and $\eta^{0}$ the initial positions for the ODE system are determined via the atomization

[TABLE]

for $i=2,...,N$. The cross-diffusion part is reconstructed at each time iteration using the discretized density. For the nonlocal part we choose

[TABLE]

where $K$ is a normalized Gaussian potential. For a diffusion coefficient $\epsilon=1$, we show segregation phenomena in Figures 1, 2 for two different choices of initial conditions. Two different types of segregation are possible. In Figure 1 the initial data for the two species are perfectly matching, this produces symmetric segregated states in the large time limit. In Figure 2 the two initial data are shifted, this produces a non symmetric segregation in which the two species form two adjacent patterns with connected support. Mixing is shown in Figures 3, 4 for the diffusion dominated regime, namely $\sigma_{1}+\sigma_{2}<2$ and diffusion coefficient $\epsilon>1$. Again, different situations may arise depending on the initial data. In the former case (perfectly overlapping initial conditions) the two species are almost entirely overlapping for large times, whereas in the latter case they overlap in a proper subset of the support of $\rho+\eta$. In all simulations we have $N=50$ and final time $T=2$.

Acknowledgments

MDF and SF were supported by the EU Erasmus Mundus MathMods programme, www.mathmods.eu, by the GNAMPA (Italian group of Analysis, Probability, and Applications) projects Geometric and qualitative properties of solutions to elliptic and parabolic equations, respectively Analisi e stabilità per modelli di equazioni alle derivate parziali nella matematica applicata. MDF and SF acknowledge hospitality by the Institut für Angewandte Mathematik: Analysis und Numerik, Westfälische Wilhelms-Universität Münster, which hosted them in 2012 when the first ideas of this work came out. This work was done while MB and AS were PIs of the Excellence Cluster Cells in Motion (CiM), funded by the DFG.