Discrete and continuum links to a nonlinear coupled transport problem of interacting populations

We are interested in exploring interacting particle systems that can be seen as microscopic models for a particular structure of coupled transport flux arising when different populations are jointly evolving. The scenarios we have in mind are inspired by the dynamics of pedestrian flows in open spaces and are intimately connected to cross-diffusion and thermo-diffusion problems holding a variational structure. The tools we use include a suitable structure of the relative entropy controlling TV-norms, the construction of Lyapunov functionals and particular closed-form solutions to nonlinear transport equations, a hydrodynamics limiting procedure due to Philipowski, as well as the construction of numerical approximates to both the continuum limit problem in 2D and to the original interacting particle systems.


Introduction
The starting point of the results presented in this paper is the following question 1 : Can one design a system of interacting particles that converges in some suitable limit to the following system of nonlinearly coupled system of transport equations: with initial conditions u(0, x) = u 0 (x) and v(0, x) = v 0 (x) (x ∈ R)? Here u and v refer to mass concentrations of some chemical species which are participating in a non-competitive manner in a joint transport process. The background of the question (and interest of M. Mimura) is connected to the role pheromones play in influencing the aggregation phenomenon, one of the main survival mechanisms in insects, birds and animal colonies; we refer the reader to [6] for more on this context. It is worth noting the coupled structure of the transport fluxes resembles situations arising in cross-diffusion and thermo-diffusion. Compare [3] for the thermodynamical foundations of cross-and thermo-diffusion and [9] for a nice paper illustrating the role of cross-diffusion mechanisms towards pattern formation in chemical systems. Our own interest in this framework targets at the fundamental understanding of well-observed optimal self-organization behaviours (e.g lane formation in counter-flows) exhibited by the motion of pedestrian flows (cf. e.g. [7] and references cited therein).
Interestingly, due to the symmetry in the structure of the equations, the system (1a) -(1b) admits a direct interpretation from the porous media theory point of view, which later turns out to be very useful in understanding mathematically the particle system origin of this transport problem.
We assume that u and v denote two populations (of pedestrians, ants, chemical species, etc.) that like to travel together. Think, for instance, of a pair of large families of individuals that wish to reach perhaps a common destination or target, under the basic assumption that besides some kind of social pairwise repulsion and adherence to the same drift there are no other interactions in the crowd made of the two populations. This basic situation can be modelled as a system of continuity equations ∂ t u + div(uw) = 0, where w is the common drift to which the two populations adhere. The velocity vector w is assumed now to comply with Darcy's law In (2), K µ ∈ (0, ∞) denotes the permeability coefficient (usually a tensor for a heterogeneous region) and p is the total (social) pressure in the system. Now, making the ansatz on the structure of the pressure and then summing up the above continuity equations, we obtain the system (1a)-(1b), where for simplicity we take K ≡ 1.
The paper is organised as follows. In Section 2, we provide some basic analytic understanding of (1a)-(1b) by transforming the system to an equivalent one, showing the local well-posedness, constructing a special class of solutions and proving the preservation of relative entropy and the consequences this has on the large-time behaviour of the system. In Section 3, we introduce a stochastically interacting many-particle system to approximate (1a)-(1b). Finally, Section 4 presents numerical illustrations of the particle system, indicating numerical evidence on the expected convergence.

Analytical results
In this section, we provide a couple of analytical results on the continuum model. We first transform the system (1a)-(1b) to an equivalent one. Using this transformation, we ensure in a straightforward way the local existence of classical solutions. In addition, we construct a special class of solutions and show remarkable properties of these solutions, especially concerning the preservation of the relative entropy.

An equivalent system
Defining w := u + v, we see that w solves the following porous media-like equation: We transform the system (1a)-(1b) posed for (u, v) into the following system for (w, u): Conversely, suppose that (w, u) satisfies the system (4a)-(4b). Then (u, v), where v = w − u, satisfies the original system (1a)-(1b). Therefore, the two systems are equivalent.
The transformation has two advantages. First, the new system (4a)-(4b) is only one-sided coupled in the sense that one can solve (4a) independently to obtain w, and then substitute to find u from (4b) with w given. Second, (4a) is the famous Boussinesq's equation of groundwater flow, while (4b) is the standard continuity equation. Both equations have been studied extensively and have a rich literature. Therefore, we can apply existing methods and techniques to handle them from the mathematical analysis point of view.

A general solution to the continuity equation by the method of characteristics
Let V be a given velocity field and f 0 : R → R be a given function. We first seek solutions for the following general continuity equation We consider the following ordinary differential equation (ODE): The solution of this ODE is X(t) = F (x, t). Conversely, we also can regard x as a function of X(t), i.e., x = G(X(t), t), where G : R × R (y, t) → G(y, t) ∈ R and G(y, 0) = y. Lemma 2.1 (Solving the continuity equation, see e.g. [1]). The function solves the continuity equation (5).

Classical solutions
The first result of this paper refers to the local existence of classical solutions of (4a)-(4b). Let T > 0 be sufficiently large but fixed and let (w 0 , u 0 ) be given. We say that the couple (w, u), where w, u : [0, T ] × R → R is a classical solution to the system (4a)-(4b) if w, u ∈ C 2,1 ([0, T ], R) and satisfy (4a)-(4b).
Theorem 2.2. Suppose that w 0 and u 0 are continuous functions in R with for some ε > 0 and all x ∈ R. There exists T * ∈ (0, T ) such that the system (4a)-(4b) has a classical solution in Proof. This theorem is a direct consequence of [10, Theorem 3.1] for the (global) existence of the Boussinesq's solution and of the Peano's theorem for the local existence of the characteristic trajectory.

A special class of solutions
Due to the particular structure of the system (4a)-(4b), namely (4a) being the Boussinesq's equation and (4b) being the continuity equation, we are able to construct a special class of solutions. We consider a solution profile of quadratic functions for w and then find u accordingly. The idea of the former has been used before, see for instance [13].
Step 1. We rewrite (4a) as We consider solutions to (8) of the form By substituting this form in (8), we obtain the following system which finally leads to Therefore, Step 2. Substituting (11) back into (4b), we obtain the following continuity equations in terms of u.
We now apply Lemma 2.1 to solve this equation. The ODE (6) becomes which gives Therefore, we obtain Concluding, we have obtained a special solution to the system (4a)-(4b) as follows for some a, b ∈ R. If one considers non-negative solutions, then one should take the positive parts of these expressions.

Preservation of the relative entropy and consequences
We observe that our original system is symmetric in the sense that if we swap u and v in (1a)-(1b) then the system remains unchanged. Therefore, if the initial data u 0 and v 0 are equal, then it is expected that u and v will be equal at any later time, which is a necessary condition for uniqueness. Two mathematical questions naturally arise at this point: (i) How to prove equality of u and v rigorously?
(ii) If u 0 and v 0 are not equal, can we still quantify the distance between u(t) and v(t) in terms of the initial data?
In this section, we provide affirmative answers to these questions using the concept of relative entropy and the total variation metric. We generalise the results of this section (and of the previous one) to a more general system in Section 2.6. It will become clear that structure of the system matches nicely with the concept of the relative entropy.
We now recall the definition of the relative entropy, the total variation metric, as well as a relationship between the twos. We refer the reader to the survey paper [11] for more information.
Let f (x) dx and g(x) dx be two probability densities on R. The relative entropy of f with respect to g is defined by The total variation distance between f (x) dx and g(x) dx is defined as Note that the relative entropy is always non-negative and it is equal to 0 if and only if f = g. Although it is not a distance (it satisfies neither the triangle inequality nor the symmetry condition), it is a useful quantity to measure the difference between two probability measures and has been used extensively in the literature. In addition, it also provides an upper-bound for the total variation distance T V (f, g) by Pinsker's inequality, see for instance [11, Theorem 1.1], as Theorem 2.3. Suppose that u, v are classical solutions to the system (1a)-(1b) that decay sufficiently fast at infinity. Then the function t → H(u(t)||v(t)) is constant, i.e., for any 0 < t < T * it holds H(u(t)||v(t)) = H(u 0 ||v 0 ).
Proof. We calculate the time-derivative of t → H(u(t)||v(t)) as follows (the time variable t is dropped in the right-hand side for simplicity of notation) Note that (14) follows due to ∂ ∂t u = 0 which is a consequence of conservation of mass. In (15) we have used integration by parts where the boundary terms vanish due to the assumption on the decay of the solution.
Corollary 2.4. For any 0 < t < T * , it holds that Proof. This inequality is direct consequence of the Pinsker's inequality (12) and Theorem 2.3. Heuristically, note that if u 0 = v 0 and (u, v) satisfies the system (1a)-(1b), then so does (v, u). To guarantee the uniqueness, it follows that u = v = 1 2 w.

Theorem 2.3 offers a much stronger result.
Proof. This is a direct consequence of Theorem 2.3 (or Corollary 2.4). Since u 0 = v 0 , we have H(u 0 ||v 0 ) = 0. Then it follows from Theorem 2.3 that H(u(t)||v(t)) = 0 for all t > 0, which in turn implies that u(t) = v(t) = 1 2 w(t) for all t > 0.

Generalisations
It is worth noting that Theorem 2.2 and Theorem 2.3 can be extended to a more general system of the form The transformed system for (w, u), where w = (u + v), now becomes For instance, if f (z) = z m−1 for some m > 1, then the equation for w be- Remark 2.7. We note that the common relation that makes the relative entropies in both Theorem 2.3 and Theorem 2.6 vanish is Tracing back this relation in the calculations, this property appears because of the combination of three ingredients: the formula of the relative entropy, the symmetry of the system, and the formulas of the continuity equations. The last two properties together form the structure of the system.
i) The continuity equations provide that ii) The symmetry of (1a) and (1b) means that X = Y (so that if we swap u and v, the system is unchanged).
In other words, we find that the relative entropy is constant essentially due to the structure of the system.

Particle system approach
In this section, we introduce a many-particle system that includes coupled weakly interacting stochastic differential equations. We formally show that the empirical measures associated to this system converge to solutions of the original system (1a)-(1b). The rigorous proof will be given in a separate paper.
We consider the following particle system: are independent standard Wiener processes, {V ε } ε≥0 are a sequence of smooth functions which are chosen later on. Note that the system in (18) can be seen as a generalisation of the many-particle system arising in [8] to our model of coupled interactions of two species. Remark also that in [5,4], the authors studied similar systems but in the absence of the stochastic noise.
We define the following empirical measures We now formally derive the system (1a)-(1b) in two steps: Step 1: Hydrodynamic limit, as n tends to infinity: where (u ε , v ε ) solves a system which depends on V ε and with some viscous terms.
Step 2: Viscosity limit, as ε tends to 0: where (u, v) solves the original system.
The derivation explains the choice of scalings occurring in the many-particle system. Now, we perform the first step.
Step 1 (Hydrodynamic limit): Let f be a sufficiently smooth function. By definition (19) of the empirical measure u n,ε t , we have Using Itô's lemma and definition of the empirical measures in (19), we derive

Numerical simulations
In this section, we illustrate numerically in 2D the solution of (18) for specific initial data and explore numerically to which extent the continuum model (1a)-(1b) can be approximated based on (18).

Continuum system
We naturally extend the one-dimensional model to two dimensions using the following notations: Denote by T * the final observation time. The populations with boundary conditions These boundary conditions ensure the conservation of mass in the system, which is also preserved by the suitable finite-volume scheme. To simulate this system, we use finite-volume discretisation on an equidistant grid where the fluxes adhere to a flux limiter. To approximate u and v, we use a first-order upwind discretisation. Fluxes are approximated with a second-order central-difference approximation.
The semi-discrete system of ODE's is non-stiff. We use an explicit integration method to maintain the positivity of the solution and acquire and to put no constraints on the Jacobi matrix. Together with the finite-volume space discretisation, this allows for discontinuous initial data. We use a four-stage Runge-Kutta integration scheme to perform the time integration.

Multi-particle system
Recall the multi-particle system formulation (18). We simulate the system in the same domain as (20). As a potential function V ε , we use where r ∈ R represents the interparticle distance, c is the interaction range parameter and ε > 0, modelling a repulsive effect for r > 0. This potential formulation is consistent with the potential function description from [8]. In addition the stochasticity allows for modelling the diffusive behaviour present in the continuum system.
Given a particle configuration at time t, we use the Euler-Maruyama method (a stochastic variant of the Euler time-integration method) to compute the configuration in t + ∆t. The positions in time step t k are updated with: Here, W i are samples of a standard normal distribution. This term emerges from the distribution of the standard Wiener process: W t ∼ N (0, t). We preserve the conservation of mass by implementing reflective boundaries. With these boundaries, we mimic the zero-flux boundaries in the continuum system.
We compute the density by approximating the empirical measure defined in (19). We smoothen the particle positions with a Gaussian kernel Φ h . This allows us to compare the multi-particle system to its continuum counterpart. This empirical measure approximation µ h (t) for particles X 1 , . . . , X N is defined as µ h (t) =

Discussion
The simulations point out a qualitative agreement with the analytical results. Finding the exact relation between the number of particles N and interaction parameter ε is challenging. This is due to how density is measured in the multiparticle system (by a finite-radius approximation of the Dirac distribution) and the hidden scaling restrictions that exist on how N and 1/ε go to infinity. This is illustrated by the following experiments.
For the related problem in [8] we observe the convergence rate (23). We believe that (23) holds in our context as well. From this we induce the condition that 1/ε, N → ∞ under the condition N grows much faster than 1/ε.
Our numerical experiments indicate that if N and 1/ε increase such that this condition is not respected, the time for the system to reach an equilibrium grows to infinity.
We analyse the density after final time T * for a varying set of parameters. We define the residual norm r i of the particle system as a discrete variant of (23) by performing a sequence of n simulations to find observed densities µ i h (T * ) and measuring the L 2 -norm distance between the densities of simulation i − 1 and i, Figure 1 depicts the residual defined in (24) for a sequence of simulations with N = 8192 fixed and ε → 0. This figure illustrates the transition from systems that converge towards an equilibrium (for ε < 2 −9 ) to stationary systems (for ε > 2 −9 ).
For ε = (1/N ) α and small α we observe the convergence in density profiles. Figure 1: L 2 -norm of density residual of particle system after t = T for subsequent simulations, for fixed N and ε → 0. Finally, the size of the smoothing length h also plays a significant role in representing the interpolated density. The finite range of the Dirac interpolation implies that some mass is lost at the boundaries of the domain. This effect is visible when comparing the density profiles at the boundaries of the domain. Otherwise, a larger smoothing length increases the convergence rate and decreases the distance to the macroscopic density profile.
Further research is required to find an appropriate measure in which experiments converge to the expected macroscopic limit inside Ω as well as a correct relation between N and ε.