Luca Roatta
email: lucaroa@gmail.com

Discretization of space and time: mass-energy relation, accelerating expansion of the Universe, Hubble constant

Assuming that space and time can only have discrete values, we obtain the expression of the gravitational potential energy that at large distance coincides with the Newtonian. In very precise circumstances it coincides with the relativistic mass-energy relation: this shows that the Universe is a black hole in which all bodies are subjected to an acceleration toward the border of the Universe itself. Since the Universe is a black hole with a fixed radius, we can obtain the density of the Universe and the value of the Hubble constant.

Introduction

Let's assume, as work hypothesis, the existence of both discrete space and discrete time, namely spatial and temporal intervals not further divisible; this assumption leads to some interesting consequences. Here we find that the discrete gravitational potential energy is very similar to the Newtonian, with which coincides at large distances.

It is also shown that the expression of the gravitational potential energy of a body inside a black hole, far from its center and far from its border, corresponds to the relativistic mass-energy relation: it means that the Universe is a black hole. In the same circumstances the body undergoes a constant acceleration toward the border of the Universe itself.

So, if we suppose that neither space nor time are continuous, but that instead both are discrete, and following the terminology used in a previous document [START_REF] Roatta | Discretization of space and time in wave mechanics: the validity limit[END_REF] , we call l 0 the fundamental length and t 0 the fundamental time.

Gravitational potential energy

In a previous document [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] we obtained that in a discrete context the modified expression of the gravitational force exerted by a body A having mass M on a body B having mass m at the distance r>d from A is:

f (r )=-mc 2 r (r -l 0 ) d (r-d) 2 (1)
where

d= GM c 2 (2)
is the radius of the black hole of mass M.

If r»l 0 Eq. ( 1) becomes:

f (r )=-mc 2 d (r -d ) 2 (3)
In a continuous context, if we want to calculate the potential energy, we have to perform an integration; in a discrete context, however, the integral should be replaced by a summation.

But, if r»l 0 , we can perform the integration without problems, because the effects of the discretization are negligible. To give an idea of when this might be valid, it must be noted that, by virtue of the value obtained [START_REF] Roatta | Discretization of space and time: determining the values of minimum length and minimum time[END_REF] for l 0 , already at the level of the atomic nucleus this approximation is certainly valid, being the size of the latter about 20 orders of magnitude greater than l 0 . This makes perfectly permissible, even in the discrete context, the use of integrals and derivatives at astronomical scale.

From the mathematical point of view, the structure of Eq. ( 3) is equivalent to the Newtonian; from the physical point of view, the approximation introduced allows us to use the same reasonings that lead to determine the traditional potential energy, and then we have:

U (r )=-mc 2 d r -d (4)
We can immediately verify that at very large distances (r»d), Eq. ( 4) coincides with the Newtonian U(r)=-GMm/r.

Eq. ( 4) shows that, for r=d=GM/c 2 , U(r) is infinite and negative; this means that we must provide infinite energy to move the body away from this position: that is, it can not be removed.

So far we have considered the case where r>d.

If d>r, the reasoning have to change, because in the first case, in analogy with Newtonian physics, it has been imposed U(∞)=0 because the two bodies do not interact at infinite distance. Here this condition can no longer be imposed: r can not extend indefinitely, because of the limitation d>r. When d>r Eq. ( 1) becomes [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] :

f (r )=-mc 2 r (r +l 0 ) d ( d-r ) 2 (5)
and Eq. (3) becomes:

f (r )=-mc 2 d (d-r ) 2 (6)
Examining Eq. ( 5), we can note that the two bodies do not interact for r=0: in this case the condition to impose is U(0)=0.

We must make a clarification: strictly speaking we could not consider r=0, because we have imposed the condition r»l 0 ; however, we have seen that this approximation is already valid if the order of magnitude of r is the atomic nucleus, a very little distance that at astronomical scale can be surely approximated to zero. Under these conditions, we can again perform the integration instead of the summation, and the result for U(r) is an expression that coincides exactly with Eq. ( 4), which therefore can be used either for the case r>d either for the case d>r.

Relativistic mass-energy relation

If we consider the case d>r and we also impose the condition d»r, Eq. ( 4) becomes: U (r )=mc 2 (7) We can obtain this result only imposing strict conditions: r»l 0 and d»r; that is, if the body is far away from the center of a black hole and far away from its edge, but anyway inside.

We got a result apparently strange that tells us that the gravitational potential energy of a body can be reduced to its energy mc [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] ; but as well it tells us that this can happen only if the body is inside a black hole. This is a different interpretation of mc 2 with respect to the relativistic one.

Since the relation E=mc [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] , that in a discrete context coincides with the gravitational potential energy U(r)=mc [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] , has been widely demonstrated both theoretically and experimentally [START_REF] Gu | A Sensitive Test of Mass-Energy Relation[END_REF] , the conditions in which we live are the ones, the only ones, that allow the correctness of this relation; then we must necessarily conclude that we live in a black hole, albeit far away from the edge.

So we also have to necessarily conclude that the Universe is a black hole. This is not a completely new concept: some papers [5][6] investigate more or less deeply this possibility. However, in this discrete theory the fact that the Universe is a black hole emerges as a consequence of the correctness of the mass-energy relation.

Let's now calculate the radius that the universe should have if it were a black hole.

There are very different estimations of the mass of the Universe [START_REF] Valev | Estimations of total mass and energy of the universe[END_REF] : if we consider a mass of about α•10 53 kg, where α is of the order of some units, we obtain that the radius of the black hole, given by Eq. ( 2), is approximately 0.7•α•10 26 m, that is about a few tens of billion light years. This is a very rough assessment, but it is compatible with the estimate of the size of the observable Universe. This strengthens the assertion that we are really inside a black hole and that this black hole is the whole Universe.

Accelerating expansion of the Universe

If we consider the case d>r and we also impose the condition d»r, the same conditions imposed to obtain U(r)=mc 2 , Eq. ( 6) becomes:

f (r )=- mc 2 d (8)
and the acceleration of the body is:

a(r)=- c 2 d (9)
Again, we can obtain this result only imposing strict conditions: r»l 0 and d»r; that is, if the body is far away from the center of a black hole and far away from its edge, but anyway inside.

We have found [START_REF] Roatta | Discretization of space and time: a slight modification to the Newtonian gravitation which implies the existence of black holes[END_REF] that, when d>r, the acceleration a body undergoes is toward the edge of the black hole, that in this case is the Universe.

Substituting the value of d shown by Eq. (2) in Eq. ( 9) we obtain:

a( r)=- c 4 GM ( 10
)
where M is the mass of the Universe. So a body, in the same conditions that allow the correctness of the mass-energy relation, undergoes a constant acceleration toward the edge of the Universe whose numeric value is a(r)~1.2/α•10 -9 m/s 2 , where α is of the order of some units.

It is a very small amount, but not null, which tells us that, with the conditions r»l 0 and d»r, all bodies are moving towards the edge of the Universe with motion uniformly accelerated. The acceleration, however, greatly increases in proximity of the universe edge, as indicated by Eq. ( 6).

This means that the Universe is undergoing an accelerated expansion? It is not so: in this discrete theory, the Universe is a black hole with a well defined radius depending only on its mass and thus constant. Instead, the objects inside the black hole (the Universe) are attracted toward the border with accelerated motion. This is not the first time that the term c 4 /GM appears in a document [START_REF] Bruneton | Notes on several phenomenological laws of quantum gravity[END_REF] , but in this discrete theory indicates the constant acceleration a body undergoes when it is inside a black hole, far away from its edge.

Where are we?

Experimental verifications of the relation E=mc 2 indicate a very precise correspondence [START_REF] Bruneton | Notes on several phenomenological laws of quantum gravity[END_REF] ; from the point of view of this discrete theory it is equivalent to say that our position in the Universe is very close to its center: according to Eq. ( 4) the higher is the precision of the mass-energy relation, the more we are close to the center of the Universe.

The difference between Eq. ( 4) and mc 2 when d»r is:

Δ E=U (r )-mc 2 =mc 2 1 1- r d -mc 2 ∼mc 2 (1+ r d )-mc 2 ∼mc 2 r d (11) So the relative difference is Δ E mc 2 = U ( r)-mc 2 m c 2 = r d (12)
In a document [START_REF] Pritchard | [END_REF] , a group of physicists at MIT reported that the mass-energy relation has been verified with an accuracy of about 10 -6 . If we could assume that this very high accuracy coincides with the relative difference in Eq. ( 12), we should deduce that our actual position in the Universe is at r=10 -6 d from the center, that is at about 0.7•α•10 4 light years, where α is of the order of some units; it's a very little distance if compared with the radius of the Universe. It seems that the center of the Universe could be inside the Milky Way.

Measuring the mass-energy relation, we can determine our position in the Universe in terms of the ratio r/d. Measuring also the acceleration evidenced in the expansion of the Universe, and using both Eq. (4) and Eq. ( 6), we can obtain the values for r and d.

Density of the Universe

We have seen that in this discrete theory the universe is a black hole and its radius, given by Eq. ( 2), is constant. So also its density is constant:

ρ= M 4 3 π d 3 = 3 c 6 4 π G 3 M 2 (13)
In a continuous context, from the equations of General Relativity has been found that

ρ c = 3 H 0 2 8 π G (14)
where H 0 is the Hubble constant and ρ c is the critical density: if the universe has an higher density, its expansion will stop and it will collapse back on itself (closed universe); if the universe has a lower density, its expansion will never stop (open universe); if the universe has exactly the critical density, its expansion will stop after an infinite time (flat universe).

Recent measurements [START_REF] Valev | Estimations of total mass and energy of the universe[END_REF] using different methods have shown with high accuracy that the density of the universe is very close to the critical density, so the universe should be flat.

In this discrete theory there is no critical density, but only density, because the universe in not expanding nor contracting. The model of the continuous context that is the closest to the discrete one is the flat universe, in which the density of the universe is exactly the critical density.

Equating Eq. ( 13) and Eq. ( 14) we obtain:

H 0 =√ 2 c 3 GM (15)
that, for Eq. ( 2), can be written as

H 0 =√ 2 c d ( 16
)
There is still a quite high uncertainty on the measured value of H 0 , but if we use the widely accepted value of H 0 ~70 km s -1 Mpc -1 , we obtain d~6 Gpc~18 Gly.

From Eq. ( 15) we can also find the mass of the universe in terms of G and H 0 :

M =√ 2 c 3 GH 0 (17)
that indicates a value of M~3•10 53 kg. Eq. ( 17) is very close to the Hoyle-Carvalho formula, obtained [START_REF] Valev | Estimations of total mass and energy of the universe[END_REF] in a different way.

Conclusion

Assuming that space and time can only have discrete values, we obtained the expression of the gravitational potential energy, that under very strict conditions coincides with mc 2 . Under the same conditions a body undergoes a constant acceleration that tends to move it away from the center of the Universe. These conditions show that the Universe is a black hole and that our position is close to its center. We also obtained the constant density of the universe and the value of the Hubble constant.


T O Gallouët

M Laborde

L Monsaingeon

An unbalanced Optimal Transport splitting scheme for general advection-reaction-diffusion problems

In this paper, we show that unbalanced optimal transport provides a convenient framework to handle reaction and diffusion processes in a unified metric framework. We use a constructive method, alternating minimizing movements for the Wasserstein distance and for the Fisher-Rao distance, and prove existence of weak solutions for general scalar reaction-diffusion-advection equations. We extend the approach to systems of multiple interacting species, and also consider an application to a very degenerate diffusion problem involving a Gamma-limit. Moreover, some numerical simulations are included.

Introduction

Since the seminal works of Jordan-Kinderlehrer-Otto [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], it is well known that certain diffusion equations can be interpreted as gradient flows in the space of probability measures, endowed with the quadratic Wasserstein distance W. The well-known JKO scheme (a.k.a. minimizing movement), which is a natural implicit Euler scheme for such gradient flows, naturally leads to constructive proofs of existence for weak solutions to equations or systems with mass conservation such as, for instance, Fokker-Planck equations [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF], Porous Media Equations [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF], aggregation equation [START_REF] Carrillo | Global-in-time weak measure solutions and finite-time aggregation for nonlocal interaction equations[END_REF], double degenerate diffusion equations [START_REF] Otto | Double degenerate diffusion equations as steepest descent[END_REF], general degenerate parabolic equation [START_REF] Agueh | Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory[END_REF] etc. We refer to the classical textbooks of Ambrosio, Gigli and Savaré [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] and to the books of Villani [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF] for a detailed account of the theory and extended bibliography. Recently, this theory has been extended to study the evolution of interacting species with mass-conservation, see for examples [START_REF] Di | Measure solutions for non-local interaction PDEs with two species[END_REF][START_REF] Zinsl | Geodesically convex energies and confinement of solutions for a multicomponent system of nonlocal interaction equations[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF][START_REF] Kinderlehrer | A wasserstein gradient flow approach to poisson-nernst-planck equations[END_REF][START_REF] Carlier | A splitting method for nonlinear diffusions with nonlocal, nonpotential drifts[END_REF].

Nevertheless in biology, for example for diffusive prey-predator models, the conservation of mass may not hold, and the classical optimal transport theory does not apply. An unbalanced optimal transport theory was recently introduced simultaneously in [START_REF] Chizat | An interpolating distance between optimal transport and Fischer-Rao[END_REF][START_REF] Chizat | Unbalanced optimal transport: geometry and Kantorovich formulation[END_REF][START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF][START_REF] Liero | Optimal transport in competition with reaction: the Hellinger-Kantorovich distance and geodesic curves[END_REF], and the resulting Wasserstein-Fisher-Rao (WFR) metrics (also referred to as the Hellinger-Kantorovich distance HK) allows to compute distances between measures with variable masses while retaining a convenient Riemannian structure. See section 2 for the definition and a short discussions on this WFR metric. We also refer to [START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF][START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with dirichlet boundary conditions[END_REF] for earlier attempts to account for mass variations within the framework of optimal transport.

The WFR metrics can be seen as an inf-convolution between Wasserstein/transport and Fisher-Rao/reaction processes, and is therefore extremely convenient to control both in a unified metric setting. This allows to deal with non-conservative models of population dynamics, see e.g. [START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF][START_REF] Kondratyev | A fitness-driven crossdiffusion system from population dynamics as a gradient flow[END_REF]. In [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF], the first and third authors proposed a variant of the JKO scheme for WFR-gradient flows corresponding to some particular class of reaction-diffusion PDEs: roughly speaking, the reaction and diffusion were handled separately in two separate FR, W metrics, and then patched together using a particular uncoupling of the inf-convolution, namely WFR 2 ≈ W 2 + FR 2 in some sense (see [18, section 3] for a thorough discussion). However, the analysis was restricted to very particular structures for the PDE, corresponding to pure WFR gradient-flows.

In this work we aim at extending this splitting scheme in order to handle more general reactiondiffusion problems, not necessarily corresponding to gradient flows. Roughly speaking, the structure of our splitting scheme is the following: the transport/diffusion part of the PDE is treated by a single Wasserstein JKO step

ρ k W -------→ transport ρ k+1/2 ,
and the next Fisher-Rao JKO step

ρ k+1/2 FR ------→ reaction ρ k+1
handles the reaction part of the evolution. As already mentioned, the WFR metric will allow to suitable control both steps in a unified metric framework. We will first state a general convergence result for scalar reaction-diffusion equations, and then illustrate on a few particular examples how the general idea can be adapted to treat e.g. prey-predator systems or very degenerate Hele-Shaw diffusion problems. In this work we do not focus on optimal results and do not seek full generality, but rather wish to illustrate the efficiency of the general approach.

Another advantage of the splitting scheme is that is well adapted to existing Monge/Kantorovich/Wasserstein numerical solvers, and the Fisher-Rao step turns out to be a simple pointwise convex problem which can be implemented in a very simple way. See also [START_REF] Chizat | A tumor growth hele-shaw problem as a gradient flow[END_REF][START_REF] Chizat | Scaling algorithms for unbalanced transport problems[END_REF] for a more direct numerical approach by entropic regularization. Throughout the paper we will illustrate the theoretical results with a few numerical tests. All the numerical simulations were implemented with the augmented Lagrangian ALG2-JKO scheme from [START_REF] Benamou | An augmented lagrangian approach to wasserstein gradient flows and applications[END_REF] for the Wasserstein step, and we used a classical Newton algorithm for the Fisher-Rao step.

The paper is organized as follows. In section 2 we recall the basic definitions and useful properties of the Wasserstein-Fisher-Rao distance WFR. Section 3 contains the precise description of the splitting scheme and a detailed convergence analysis for a broad class of reaction-diffusion equations. In section 4 we present an extension to some prey-predator multicomponent systems with nonlocal interactions. In section 5 we extend the general result from section 3 to a very degenerate tumor growth model studied in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF], corresponding to a pure WFR gradient flow: we show that the splitting scheme captures fine properties of the model, particularly the Γ-convergence of discrete gradient flows as the degenerate diffusion parameter of Porous Medium type m → ∞. The last section 6 contains an extension to a tumor-growth model coupled with an evolution equation for the nutrients.

Preliminaries

Let us first fix some notations. Throughout the whole paper, Ω denotes a possibly unbounded convex subset of R d , Q T represents the product space [0, T ] × Ω, for T > 0, and we write M + = M + (Ω) for the set of nonnegative finite Radon measures on Ω. We say that a curve of measures t → ρ t ∈ C w ([0, 1]; M + ) is narrowly continuous if it is continuous with respect to the narrow convergence of measures, namely for the duality with C b (Ω) test-functions.

Definition 2.1. The Fisher-Rao distance between ρ 0 , ρ 1 ∈ M + is FR(ρ 0 , ρ 1 ) := min (ρt,rt)∈AFR[ρ0,ρ1] ˆ1 0 ˆΩ |r t | 2 dρ t (x)dt,
where the admissible set

A FR [ρ 0 , ρ 1 ] consists in curves [0, 1] t → (ρ t , r t ) ∈ M + × M such that t → ρ t ∈ C w ([0, 1]; M + ) is narrowly continuous with endpoints ρ t (0) = ρ 0 , ρ t (1) = ρ 1 , and
∂ t ρ t = ρ t r t
in the sense of distributions D ((0, 1) × Ω).

The Monge-Kantorovich-Wasserstein admits several equivalent definitions and formulations, and we refer e.g. to [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Villani | Optimal transport, volume 338 of Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF][START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF][START_REF] Santambrogio | Optimal Transport for Applied Mathematicians[END_REF] for a complete description. For our purpose we shall only need the dynamical Benamou-Brenier formula: Theorem 2.2 (Benamou-Brenier formula, [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF][START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF]). There holds

W 2 (ρ 0 , ρ 1 ) = min (ρ,v)∈AW[ρ0,ρ1] ˆ1 0 ˆΩ |v t | 2 dρ t dt, (2.1)
where the admissible set A W [ρ 0 , ρ 1 ] consists in curves (0, 1) t → (ρ t , v t ) ∈ M + × M(Ω; R d ) such that t → ρ t is narrowly continuous with endpoints ρ t (0) = ρ 0 , ρ t (1) = ρ 1 and solving the continuity equation

∂ t ρ t + div(ρ t v t ) = 0
in the sense of distributions D ((0, 1) × Ω).

According to the original definition in [START_REF] Chizat | An interpolating distance between optimal transport and Fischer-Rao[END_REF] we have

Definition 2.3. The Wasserstein-Fisher-Rao distance between ρ 0 , ρ 1 ∈ M + (Ω) is WFR 2 (ρ 0 , ρ 1 ) := inf (ρ,v,r)∈AWFR[ρ0,ρ1] ˆ1 0 ˆΩ(|v t (x)| 2 + |r t | 2 ) dρ t (x)dt, (2.2)
where the admissible set

A WFR [ρ 0 , ρ 1 ] is the set of curves t ∈ [0, 1] → (ρ t , v t , r t ) ∈ M + × M(Ω; R d ) × M such that t → ρ t ∈ C w ([0, 1], M + ) is narrowly continuous with endpoints ρ |t=0 = ρ 0 , ρ |t=1 = ρ 1
and solves the continuity equation with source

∂ t ρ t + div(ρ t v t ) = ρ t r t .
Comparing definition 2.3 with definition 2.1 and Theorem 2.2, this dynamical formulation à la Benamou-Brenier shows that the WFR distance can be viewed as an inf-convolution of the Wasserstein and Fisher-Rao distances W, FR. From [START_REF] Chizat | An interpolating distance between optimal transport and Fischer-Rao[END_REF][START_REF] Chizat | Unbalanced optimal transport: geometry and Kantorovich formulation[END_REF][START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF] the infimum in (2.2) is always a minimum, and the corresponding minimizing curves t → ρ t are of course constant-speed geodesics WFR(ρ t , ρ s ) = |t -s|WFR(ρ 0 , ρ 1 ). Then (M + , WFR) is a complete metric space, and WFR metrizes the narrow convergences of measures (see again [START_REF] Chizat | An interpolating distance between optimal transport and Fischer-Rao[END_REF][START_REF] Chizat | Unbalanced optimal transport: geometry and Kantorovich formulation[END_REF][START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF]). Interestingly, there are other possible formulations of the distance in terms of static unbalanced optimal transportation, primal-dual characterizations with relaxed marginals, lifting to probability measures on a cone over Ω, duality with subsolutions of Hamilton-Jacobi equations, and we refer to [START_REF] Chizat | An interpolating distance between optimal transport and Fischer-Rao[END_REF][START_REF] Chizat | Unbalanced optimal transport: geometry and Kantorovich formulation[END_REF][START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF][START_REF] Liero | Optimal transport in competition with reaction: the Hellinger-Kantorovich distance and geodesic curves[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF] for more details.

As a first useful interplay between the distances WFR, W, FR we have

Proposition 2.4 ([18]). Let ρ 0 , ρ 1 ∈ M + 2 such that |ρ 0 | = |ρ 1 |. Then WFR 2 (ρ 0 , ρ 1 ) W 2 (ρ 0 , ρ 1 ).
Similarly for all µ 0 , µ 1 ∈ M + (with possibly different masses) there holds

WFR 2 (µ 0 , µ 1 ) FR 2 (µ 0 , µ 1 ).
Finally, for all ν 0 , ν

1 ∈ M + 2 such that |ν 0 | = |ν 1 | and all ν ∈ M + , there holds WFR 2 (ν 0 , ν) 2(W 2 (ν 0 , ν 1 ) + FR 2 (ν 1 , ν)).
Moreover, we have the following link between the reaction and the velocity in (2.2), which was the original definition in [START_REF] Kondratyev | A new optimal transport distance on the space of finite radon measures[END_REF]: Proposition 2.5 ([18]). The definition (2.3) of the WFR distance can be restricted to the subclass of admissible paths (v t , r t ) = (∇u t , u t ) for potentials u t ∈ H 1 (dρ t ) and continuity equations

∂ t ρ t + div(ρ t ∇u t ) = ρ t u t .
This shows that (M + , WFR) can be endowed with the formal Riemannian structure constructed as follow: any two tangent vectors ξ 1 = ∂ t ρ 1 , ξ 2 = ∂ t ρ 2 can be uniquely identified with potentials u i by solving the elliptic equations

ξ i = -div(ρ∇u i ) + ρu i .
Then the Riemaniann tensor is naturally constructed on the H 1 (dρ) scalar product, i-e

g ρ (ξ 1 , ξ 2 ) := u 1 , u 2 H 1 (dρ) = ˆΩ(∇u 1 • ∇u 2 + u 1 u 2 )dρ.
This is purely formal, and we refer again to [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF] for discussions. Given a functional

F(ρ) := ˆΩ F (ρ) + ˆΩ ρV + 1 2 ˆΩ(K * ρ)ρ,
this Riemannian structure also allows to compute WFR gradients as

grad WFR F(ρ) = -div ρ∇ δF δρ + ρ δF δρ = grad W F(ρ) + grad FR F(ρ),
where δF δρ = F (ρ) + V + K * ρ denotes the Euclidean first variation of F with respect to ρ. In other words, the Riemannian tangent vector grad WFR F(ρ) is represented in the previous H 1 (dρ) duality by the scalar potential u = δF δρ .

An existence result for general parabolic equations

In this section, we propose to solve scalar parabolic equations of the form

   ∂ t ρ = div(ρ∇(F 1 (ρ) + V 1 )) -ρ(F 2 (ρ) + V 2 ) ρ| t=0 = ρ 0 ρ∇(F 1 (ρ) + V 1 )| ∂Ω • ν = 0 (3.1)
in a bounded domain Ω ⊂ R d with Neumann boundary condition and suitable initial conditions. Our goal is to extend to the case

F 1 = F 2 , V 1 = V 2
the method initially introduced in [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF] for variational WFR-gradient flows, i-e (3.1) with

F 1 = F 2 and V 1 = V 2 .
We assume for simplicity that F 1 : R → R is given by

F 1 (z) =    z log z -z (linear diffusion) or 1 m1-1 z m1 (Porous Media diffusion) , (3.2)
and F 2 : R → R is given by

F 2 (z) = 1 m 2 -1 z m2 , for some m 2 > 1. (3.3)
Note that we cannot take F 2 (z) = z log z -z because the Boltzmann entropy is not well behaved (neither regular nor convex) with respect to the Fisher-Rao metric in the reaction step, see [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF][START_REF] Liero | Optimal transport in competition with reaction: the Hellinger-Kantorovich distance and geodesic curves[END_REF][START_REF] Liero | Optimal Entropy-Transport problems and a new Hellinger-Kantorovich distance between positive measures[END_REF] for discussions. In addition, we assume that

V 1 ∈ W 1,∞ (Ω) and V 2 ∈ L ∞ (Ω).
We denote E 1 , E 2 : M + → R the energy functionals

E i (ρ) := F i (ρ) + V i (ρ),
where

F i (ρ) := ´Ω F i (ρ) if ρ L |Ω +∞ otherwise, and V i (ρ) := ˆΩ V i ρ.
Although more general statements with suitable structural assumptions could certainly be proved, we do not seek full generality here and choose to restrict from the beginning to the above simple (but nontrivial) setting for the sake of exposition.

Definition 3.1. A weak solution of (3.1) is a curve [0, +∞) t → ρ(t, •) ∈ L 1 + ∩ L ∞ (Ω) such that for all T < ∞ the pressure P 1 (ρ) := ρF 1 (ρ) -F 1 (ρ) satisfies ∇P 1 (ρ) ∈ L 2 ([0, T ] × Ω), and ˆ+∞ 0 ˆΩ(ρ∂ t φ -∇V 1 • ∇φρ -∇P 1 (ρ) • ∇φ -ρ(F 2 (ρ) + V 2 )φ) dx dt = -ˆΩ φ(0, x)ρ 0 (x) dx for every φ ∈ C ∞ c ([0, +∞) × R d ).
Note that the pressure P 1 is defined so that the diffusion term div(ρ∇F 1 (ρ)) = ∆P 1 (ρ), at least for smooth solutions.

The starting point of our analysis is that (3.1) can be written, at least formally as,

∂ t ρ = div(ρ∇(F 1 (ρ) + V 1 )) -ρ(F 2 (ρ) + V 2 ) ↔ ∂ t ρ = -grad W E 1 (ρ) -grad FR E 2 (ρ).
Our splitting scheme is a variant of that originally introduced in [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF], and can be viewed as an operator splitting method: each part of the PDE above is discretized (in time) in its own W, FR metric, and corresponds respectively to a W/transport/diffusion step and to a FR/reaction step. More precisely, let h > 0 be a small time step. Starting from the initial datum ρ 0 h := ρ 0 , we construct two recursive sequences (ρ k h ) k and (ρ

k+1/2 h ) k such that          ρ k+1/2 h ∈ argmin ρ∈M + ,|ρ|=|ρ k h | 1 2h W 2 (ρ, ρ k h ) + E 1 (ρ) , ρ k+1 h ∈ argmin ρ∈M + 1 2h FR 2 2 (ρ, ρ k+1/2 h ) + E 2 (ρ) . (3.4)
With our structural assumptions on F i , V i and arguing as in [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF], the direct method shows that this scheme is well-posed, i-e that each minimizing problem in (3.4) admits a unique minimizer. We construct next two piecewise-constant interpolating curves k+1/2 h , for all t ∈ (kh, (k + 1)h].

(3.5)

Our main results in this section is the constructive existence of weak solutions to (3.1):

Theorem 3.2. Assume that ρ 0 ∈ L 1 + ∩ L ∞ (Ω).
Then, up to a discrete subsequence (still denoted h → 0 and not relabeled here), ρ h and ρh converge strongly in L 1 ((0, T ) × Ω) to a weak solution ρ of (3.1).

Note that any uniqueness for (3.1) would imply convergence of the whole (continuous) sequence ρ h , ρh → ρ as h → 0, but for the sake of simplicity we shall not address this issue here.

The main technical obstacle in the proof of Theorem 3.2 is to retrieve compactness in time. For the classical minimizing scheme of any energy E on any metric space (X, d), suitable time compactness is usually retrieved in the form of the total-square distance estimate 1 2h k≥0

d 2 (x k , x k+1 ) E(x 0 ) -inf E.
This usually works because only one functional is involved, and E(x 0 ) -inf E is obtained as a telescopic sum of one-step energy dissipations E(x k+1 ) -E(x k ). Here each of our elementary step in (3.1) involves one of the W, FR metrics, and we will use the WFR distance to control both simultaneously: this strongly leverages the inf-convolution structure, the WFR distance being precisely built on a compromise between W/transport and FR/reaction. On the other hand we also have two different functionals E 1 , E 2 , and we will have to carefully estimate the dissipation of E 1 during the FR reaction step (driven by E 2 ) as well as the dissipation of E 2 during the W transport/diffusion step (driven by E 1 ).

We start by collecting one-step estimates, exploiting the optimality conditions for each elementary minimization procedure, and postpone the proof of Theorem 3.2 to the end of the section.

Optimality conditions and pointwise L ∞ estimates

The optimality conditions for the first Wasserstein step ρ k → ρ k+1/2 in (3.4) are by now classical, and can be written for example

-∇ϕ k+1/2 h h ρ k+1/2 h = ∇P 1 (ρ k+1/2 h ) + ρ k+1/2 h ∇V 1 a.e. ( 3.6)
Here ϕ k+1/2 h is an optimal (backward) Kantorovich potential from ρ

k+1/2 h to ρ k h . Lemma 3.3. For all k 0, ρ k+1/2 h L 1 = ρ k h L 1 (3.7)
and for all constant

C such that V 1 C, ρ k h (x) (F 1 ) -1 (C -V 1 (x)) a.e. ⇒ ρ k+1/2 h (x) (F 1 ) -1 (C -V 1 (x)) a.e. (3.8)
Proof. The Wasserstein step is mass conservative by construction, so the first part is obvious. The second part is a direct consequence of a generalization [36, lemma 2] of Otto's maximum principle [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF].

Remark 3.4. Note that if ρ k h M , we may take C = F 1 (M ) + V 1 L ∞ in (3.8). Formally, this corresponds to taking ρ(x) := (F 1 ) -1 (C -V 1 (x)) as a stationary Barenblatt supersolution for ∂ t ρ = div(ρ∇(F 1 (ρ) + V 1 )
) at the continuous level. In addition, if V 1 ≡ 0 we recover Otto's maximum principle [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF] in the form

ρ k+1/2 L ∞ ρ k L ∞ .
For the second Fisher-Rao reaction step, the optimality condition has been obtained in [18, section 4.2] in the form

ρ k+1 h -ρ k+1/2 h ρ k+1 h = - h 2 ρ k+1 h F 2 (ρ k+1 h ) + V 2 a.e. ( 3.9)
As a consequence we have Lemma 3.5. There is C ≡ C(V 2 ) > 0 such that for h h 0 (V 2 ) small enough we have

ρ k+1 h (x) (1 + Ch)ρ k+1/2 h (x) a.e., (3.10)
and for all M > 0 there is

c ≡ c(M, V 2 ) such that if ρ k+1/2 h ∞ M then (1 -ch)ρ k+1/2 h (x) ρ k+1 h (x) a.e. (3.11)
Note in particular that this immediately implies

supp ρ k+1 h = supp ρ k+1/2 h , (3.12)
which was to be expected since the reaction part ∂ t ρ = -ρ(F 2 (ρ) + V 2 ) of the PDE (3.1) preserves strict positivity.

Proof. We start with the upper bound: inside supp ρ k+1 h , (3.9) and F 2 0 give

ρ k+1 h (x) -ρ k+1/2 h (x) = -h ρ k+1 h (x)(F 2 (ρ k+1 h (x)) + V 2 (x)) -hV 2 (x) ρ k+1 h (x) h V 2 ∞ ρ k+1 h (x) whence ρ k+1 h (x) 1 1 -h V 2 ∞ ρ k+1/2 h (x).
Taking squares and using

1 (1 -h V 2 ∞ ) 2 = 1 + 2 V 2 L ∞ h + O(h 2 ) 1 + 3 V 2 L ∞ h
for small h gives the desired inequality.

For the lower bound (3.11), we first observe that since F 2 0 and from (3.10) we have

F 2 (ρ k+1 h ) F 2 ((1 + Ch)ρ k+1/2 h ) F 2 (2M ) if h is small enough. Then (3.9) gives inside supp ρ k+1 ρ k+1 h (x) -ρ k+1/2 h (x)) = -h ρ k+1 h (x)(F 2 (ρ k+1 h (x)) + V 2 (x)) -h(F 2 (2M ) + V 2 ∞ ) ρ k+1 h (x), hence ρ k+1 h (x) 1 (1 + h(F 2 (2M ) + V 2 ∞ )) 2 ρ k+1/2 h (x) (1 -ch)ρ k+1/2 h (x)
for small h.

Combining Lemma 3.3 and Lemma 3.5, we obtain at the continuous level Proposition 3.6. For all T > 0 there exist constants M T , M T such that for all t ∈ [0, T ],

ρ h (t) L 1 ∩L ∞ , ρh (t) L 1 ∩L ∞ M T and ρ h (t) -ρh (t) L 1 hM T uniformly in h 0.
Note from the second estimate that strong L 1 ((0, T ) × Ω) convergence of ρ h will immediately imply convergence of ρh to the same limit.

Proof. By induction combining (3.8) and (3.10), we obtain, for all t ∈ [0, T ], [36, lemma 2]. The L 1 bound is even easier: since the Wasserstein step is mass preserving, we can integrate (3.10) in space to get

ρ h (t) L ∞ , ρh (t) L ∞ C T , where C T is a constant depending on V 1 L ∞ , see
ρ k+1 h L 1 (1 + Ch) ρ k+1/2 h L 1 = (1 + Ch) ρ k+1 h L 1 .
For t T ⇔ k T /h the L 1 bounds immediately follow by induction, with (1+Ch) T /h e CT . and we conclude again by induction.

In order to compare now ρ h and ρh , we take advantage of the above upper bound to write ρ k+1/2 h M T as long as kh T . Taking c = c(M T ) in (3.11) and combining with (3.10), we have

-chρ k+1/2 h ρ k+1/2 h -ρ k+1 h Chρ k+1/2 h a.e.
Integrating in Ω we conclude that

ρ h (t) -ρh (t) 1 = ρ k+1 h -ρ k+1/2 h 1 h max{c, C} ρ k+1/2 h 1 h max{c, C}M T = hM T
and the proof is complete.

Energy dissipation

Our goal is here to estimate the crossed dissipation along each elementary W, FR step. Testing ρ = ρ k h in the first Wasserstein step in (3.4), we get as usual

1 2h W 2 (ρ k+1/2 h , ρ k h ) F 1 (ρ k h ) -F 1 (ρ k+1/2 h ) + ˆΩ V 1 (ρ k h -ρ k+1/2 h
).

(3.13) Since V 1 is Globally Lipschitz we can first use standard methods from [START_REF] Di | Measure solutions for non-local interaction PDEs with two species[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF]

to control ´Ω V 1 (ρ k h - ρ k+1/2 h ) in terms of W 2 (ρ k+1/2 h , ρ k h )
, and suitably reabsorb in the left-hand side to obtain

1 4h W 2 (ρ k+1/2 h , ρ k h ) F 1 (ρ k h ) -F 1 (ρ k+1/2 h ) + C T h. (3.14)
The dissipation of F 1 along the Fisher-Rao step is controlled as Proposition 3.7. For all T > 0 there exists a constant C T > 0 such that, for all k

0 and k ≤ T /h , F 1 (ρ k+1 h ) F 1 (ρ k+1/2 h ) + C T h. (3.15)
Proof. We first treat the case of

F 1 (z) = 1 m1-1 z m1 with m 1 > 1. Since F 1 is increasing, we use (3.10) to obtain F 1 (ρ k+1 h ) -F 1 (ρ k+1/2 h ) ((1 + Ch) m1 -1) m 1 -1 ˆΩ(ρ k+1/2 h ) m1 Ch ρ k+1/2 m1-1 L ∞ ρ k+1/2
L 1 , and we conclude from Proposition 3.6.

In the second case F 1 (z) = z log(z) -z, we have

F 1 (ρ k+1 h ) = ˆ{ρ k+1 h e -1 } ρ k+1 h log(ρ n+1 h ) + ˆ{ρ k+1 h e -1 } ρ k+1 h log(ρ k+1 h ) -ˆΩ ρ k+1 h .
Note from Proposition 3.6 that the z contribution in

F 1 (z) = z log z -z is immediately controlled by | ´ρk+1 h - ´ρk+1/2 h | ρ k+1 h -ρ k+1/2 h L 1
hM T , so we only have to estimate the z log z contribution. Since z → z log z is increasing on {z e -1 } and using (3.10), the second term in the right hand side becomes

ˆ{ρ k+1 h e -1 } ρ k+1 h log(ρ k+1 h ) ˆ{ρ k+1 h e -1 } (1 + Ch)ρ k+1/2 h log((1 + Ch)ρ k+1/2 h ) ˆ{ρ k+1 h e -1 } ρ k+1/2 h log(ρ k+1/2 h ) + Ch ˆ{ρ k+1 h e -1 } ρ k+1/2 h log(ρ k+1/2 h ) +(1 + Ch) ˆ{ρ k+1 h e -1 } ρ k+1/2 h log(1 + Ch) ˆ{ρ k+1 h e -1 } ρ k+1/2 h log(ρ k+1/2 h ) + C T h, where we used ρ k+1/2 h L 1
M T from Proposition 3.6 as well as log(1 + Ch) Ch in the last inequality. Using the same method with the bound from below (3.11) on {ρ k+1 h e -1 } (where z → z log z is now decreasing), we obtain similarly

ˆ{ρ k+1 h e -1 } ρ k+1 h log(ρ k+1 h ) ˆ{ρ k+1 h e -1 } ρ k+1/2 h log(ρ k+1/2 h ) + C T h.

Combining both inequalities gives

ˆΩ ρ k+1 h log(ρ k+1 h ) ˆΩ ρ k+1/2 h log(ρ k+1/2 h ) + C T h
and the proof is complete.

Summing (3.14) and (3.15) over k we obtain

1 2h N -1 k=0 W 2 (ρ k+1/2 h , ρ k h ) F 1 (ρ 0 ) -F 1 (ρ N h ) + C T , (3.16)
where N = T h .

In the above estimate we just controlled the dissipation of F 1 along the FR/reaction steps, and the goal is now to similarly estimate the dissipation of F 2 along the Wasserstein step. Testing ρ = ρ k+1/2 h in the second Fisher-Rao step in (3.4), we obtain

1 2h FR 2 (ρ k+1 h , ρ k+1/2 h ) F 2 (ρ k+1/2 h ) -F 2 (ρ k+1 h ) + ˆΩ V 2 (ρ k+1/2 h -ρ k+1 h ). (3.17) Since we assumed V 2 ∈ L ∞ (Ω) and because ρ h (t) = ρ k+1 h remains close to ρh (t) = ρ k+1/2 h
in L 1 uniformly in t, h by Proposition 3.6, we immediately control the potential part as

ˆΩ V 2 (ρ k+1/2 h -ρ k+1 h ) V 2 ∞ C T h. (3.18)
For the internal energy we argue exactly as in the proof Proposition 3.7 (for the Porous Media part, since we chose here F 2 (z) = 1 m2-1 z m2 ), and obtain

F 2 (ρ k+1/2 h ) -F 2 (ρ k+1 h ) C T h. (3.19)
Combining (3.17), (3.18) and (3.19), we immediately deduce that

1 2h N -1 k=0 FR 2 (ρ k+1/2 h , ρ k+1 h ) C T , (3.20)
where N = T h as before.

Finally, we recover an approximate compactness in time in the form Proposition 3.8. There exists a constant C T > 0 such that for all h small enough and k N = T /h ,

1 h N -1 k=0 WFR 2 (ρ k h , ρ k+1 h ) 4F 1 (ρ 0 ) + C T . (3.21)
Proof. Adding (3.16) and (3.20) gives

1 h N -1 k=0 W 2 (ρ k h , ρ k+1/2 h ) + FR 2 (ρ k+1/2 h , ρ k+1 h ) 2 F 1 (ρ 0 ) -F 1 (ρ N h ) + C T + 2C T 2F 1 (ρ 0 ) + C T , since in any case F 1 (z) = 1 m1-1 z m1 0 and F 1 (z) = z log z -z -1 is bounded from below on the bounded domain Ω, hence F 1 (ρ N h ) -C Ω uniformly. It then follows from Proposition 2.4 that W 2 (ρ k h , ρ k+1/2 h ) + FR 2 (ρ k+1/2 h , ρ k+1 h ) 1 2 WFR 2 ρ k h , ρ k+1 h
in the left-hand side, and the result immediately follows.

Estimates and convergences

From the total-square distance estimate (3.21) we recover as usual the approximate 1 2 -Hölder estimate

WFR(ρ h (t), ρ h (s)) + WFR(ρ h (t), ρh (s)) C T |t -s + h|
(WFR(ρ h (t), ρ(t)) + WFR(ρ h (t), ρ(t))) → 0.
In order to pass to the limit in the nonlinear terms, we first strengthen this WFR-convergence into a more tractable L 1 convergence. The first step is to retrieve compactness in space: Proposition 3.9. For all T > 0, ρ h and ρh satisfies

P 1 (ρ h ) L 2 ([0,T ];H 1 (Ω)) C T .
(3.24)

Proof. From (3.6) and the L 1 ∩ L ∞ bounds from Proposition 3.6 we see that

ˆΩ |∇P 1 (ρ k+1/2 h )| 2 1 2h 2 ˆΩ |∇ϕ k+1/2 h | 2 (ρ k+1/2 h ) 2 + 1 2 ˆΩ |∇V 1 | 2 (ρ k+1/2 h ) 2 C T 2h 2 ˆΩ |∇ϕ k+1/2 h | 2 ρ k+1/2 h + 1 2 ∇V 1 2 ∞ ˆΩ(ρ k+1/2 h ) 2 C T W 2 (ρ k+1/2 h , ρ k h ) h 2 + 1 since ϕ k+1/2 h
is the optimal (backward) Kantorovich potential from ρ k+1/2 h to ρ k h . Multiplying by h > 0, summing over k, and exploiting (3.16) gives

P 1 (ρ h ) 2 L 2 ([0,T ];H 1 (Ω)) N -1 k=0 h P 1 (ρ k+1/2 h ) 2 H 1 C T (F 1 (ρ 0 ) -F 1 (ρ N h ) + 1) C T ,
where we used as before F 1 (ρ N h ) -C Ω in the last inequality.

We are now in position of proving our main result:

Proof of Theorem 3.2. Exploiting (3.21) and (3.24), we can apply the extension of the Aubin-Lions lemma established by Rossi and Savaré in [START_REF] Rossi | Tightness, integral equicontinuity and compactness for evolution problems in Banach spaces[END_REF] to obtain that ρh converges to ρ strongly in L 1 (Q T ) (see [START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF]). By diagonal extraction if needed, we can assume that the convergence holds in L 1 (Q T ) for all fixed T > 0. Then by Proposition 3.6 we have

ρ h -ρ L 1 (Q T ) ρ h -ρh L 1 (Q T + ρh -ρ L 1 (Q T ) C T h + ρh -ρ L 1 (Q T ) → 0 hence ρ h → ρ as well. Moreover, since P 1 (ρ h ) is bounded in L 2 ((0, T ), H 1 (Ω)) we can assume that ∇P 1 (ρ h ) ∇P 1 (ρ) in L 2 ((0, T ), H 1 (Ω)
) for all T > 0. Exploiting the Euler-Lagrange equations (3.6)(3.9) and arguing exactly as in [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF]Theorem 4], it is easy to pass to the limit to conclude that

ˆΩ ρ(t 2 )ϕ -ρ(t 1 )ϕ = - ˆt2 t1 ˆΩ ∇P (ρ) • ∇ϕ + ρ∇V 1 • ∇ϕ -ρ(F 2 (ρ) + V 2 )ϕ for all 0 < t 1 < t 2 and ϕ ∈ C 1 b (Ω). Since ρ ∈ C([0, T ]; M + WFR )
takes the initial datum ρ(0) = ρ 0 and WFR metrizes the narrow convergence of measures, this is well-known to be equivalent to our weak formulation in Definition 3.1, and the proof is complete. Remark 3.10. In the above proofs one can check that Theorem 3.2 extends in fact to all C 1 nonlinearities F 2 such that F 2 C for some C ∈ R. Likewise, we stated and proved our main result in bounded domains for convenience: all the above arguments immediately extend to Ω = R d at least for F 1 (z) = 1 m1-1 z m1 0. The only place where we actually used the boundedness of Ω was in the proof of Proposition 3.8, when we bounded from below F 1 (ρ N h ) -C Ω in order to retrieve the total-square distance estimate. When Ω = R d and F 1 (z) = z log z -z a lower bound F 1 (ρ N h ) -C T still holds, but the proof requires a tedious control of the second moments m 2 (ρ) = ´Rd |x| 2 ρ hence we did not address this technical issue for the sake of brevity.

Application to systems

In this section we shall try to illustrate that the previous scheme is very tractable and allows to solve systems of the form

   ∂ t ρ 1 = div(ρ 1 ∇(F 1 (ρ 1 ) + V 1 [ρ 1 , ρ 2 ])) -ρ 1 (G 1 (ρ 1 ) + U 1 [ρ 1 , ρ 2 ]), ∂ t ρ 2 = div(ρ 2 ∇(F 2 (ρ 2 ) + V 2 [ρ 1 , ρ 2 ])) -ρ 2 (G 2 (ρ 2 ) + U 2 [ρ 1 , ρ 2 ]), ρ 1|t=0 = ρ 1,0 , ρ 2|t=0 = ρ 2,0 . (4.1)
For simplicity we assume again that Ω is a smooth, bounded subset of R d . Then the system (4.1) is endowed with Neumann boundary conditions,

ρ 1 ∇(F 1 (ρ 1 ) + V 1 [ρ 1 , ρ 2 ]) • ν = 0 and ρ 2 ∇(F 2 (ρ 2 ) + V 2 [ρ 1 , ρ 2 ]) • ν = 0 on R + × ∂Ω,
where ν is the outward unit normal to ∂Ω. In system of the form (4.1), we allow interactions between densities in the potential terms

V i [ρ 1 , ρ 2 ] and U i [ρ 1 , ρ 2 ].
In the mass-conservative case (without reaction terms), this system has already been studied in [START_REF] Di | Measure solutions for non-local interaction PDEs with two species[END_REF][START_REF] Laborde | On some non linear evolution systems which are perturbations of Wasserstein gradient flows[END_REF][START_REF] Carlier | A splitting method for nonlinear diffusions with nonlocal, nonpotential drifts[END_REF], using a semi-implicit JKO scheme introduced by Di Francesco and Fagioli, [START_REF] Di | Measure solutions for non-local interaction PDEs with two species[END_REF]. This section combines the splitting scheme introduced in the previous section and semi-implicit schemes both for the Wasserstein JKO step and for the Fisher-Rao JKO step.

For the ease of exposition we keep the same assumptions for F i and G i as in the previous section, i.e the diffusion terms F i satisfy (3.2) and the reaction terms G i satisfy (3.3). Moreover, since the potentials depend now on the densities ρ 1 and ρ 2 , we need stronger hypotheses: we assume that

V i : L 1 (Ω; R + ) 2 → C 1 (Ω) are continuous and verify, uniformly in ρ 1 , ρ 2 ∈ L 1 (Ω; R + ), V i [ρ 1 , ρ 2 ] W 1,∞(Ω) K(1 + ρ 1 L 1 (Ω) + ρ 2 L 1 (Ω) ), ∇(V i [ρ 1 , ρ 2 ]) -∇(V i [µ 1 , µ 2 ]) L ∞ (Ω) K( ρ 1 -µ 1 L 1 (Ω) + ρ 2 -µ 2 L 1 (Ω) ). (4.2)
The interacting potentials we have in mind are of the form

V i [ρ 1 , ρ 2 ] = K i,1 * ρ 1 + K i,2 * ρ 2 , where K i,1 , K i,2 ∈ W 1,∞ (
Ω) and then V i satisfies (4.2). For the reaction, we assume that the potentials U i are continuous from L 1 (Ω) 2 + to L 1 with moreover

U i [ρ 1 , ρ 2 ] -K, ∀ ρ 1 , ρ 2 ∈ L 1 (Ω; R + ) (4.3)
for some K ∈ R, and

U i [ρ 1 , ρ 2 ] L ∞ (Ω) K M , ∀ ρ 1 L 1 (Ω) , ρ 2 L 1 (Ω) M (4.4)
for some nondecreasimg function K M 0 of M . The examples we have in mind are of the form

U 1 [ρ 1 , ρ 2 ] = C 1 ρ 2 1 + ρ 1 , U 2 [ρ 1 , ρ 2 ] = -C 2 ρ 1 1 + ρ 1 for some constants C i ≥ 0, or nonlocal reactions U i [ρ 1 , ρ 2 ](x) = ˆΩ K i,1 (x, y)ρ 1 (y) dy + ˆΩ K i,2 (x, y)ρ 2 (y) dy
for some nonnegative kernels K i,j ∈ L 1 ∩L ∞ . Such reaction models appear for example in biological adaptive dynamics [START_REF] Perthame | Transport equations in biology[END_REF]. Definition 4.1. We say that (ρ 1 , ρ 2 ) : R + → L 1 + ∩ L ∞ + (Ω) is a weak solution of (4.1) if, for i ∈ {1, 2} and all T < +∞, the pressure

P i (ρ i ) := ρ i F i (ρ i ) -F i (ρ i ) satisfies ∇P i (ρ i ) ∈ L 2 ([0, T ] × Ω), and
ˆ+∞ 0 ˆΩ(ρ∂ t φ i -ρ i ∇V i [ρ 1 , ρ 2 ] • ∇φ i -∇P i (ρ i ) • ∇φ i -ρ i (G i (ρ i ) + U i [ρ 1 , ρ 2 ])φ i ) dx dt = -ˆΩ φ i (0, x)ρ i,0 (x) dx, (4.5) for all φ i ∈ C ∞ c ([0, +∞) × R d ).
Then, the following result holds,

Theorem 4.2. Assume that ρ 1,0 , ρ 2,0 ∈ L 1 ∩ L ∞ + (Ω) and that V i , U i satisfy (4.2)(4.
3)(4.4). Then (4.1) admits at least one weak solution.

Note that this result can be easily adapted to systems with an arbitrary number of species N 2, coupled by nonlocal terms

V i [ρ 1 , . . . , ρ N ] and U i [ρ 1 , . . . , ρ N ].
Remark 4.3. A refined analysis shows that our approach would allow to handle systems of the form

∂ t ρ 1 -div(ρ 1 ∇(F 1 (ρ 1 ) + V 1 )) = -ρ 1 h 1 (ρ 1 , ρ 2 ), ∂ t ρ 2 -div(ρ 2 ∇(F 2 (ρ 2 ) + V 2 )) = +ρ 2 h 2 (ρ 1 ),
where h 1 is a nonnegative continuous function and h 2 is a continuous functions. Indeed since h 1 ≥ 0 the reaction term is the first equation is nonpositive, hence ρ 1 (t) L ∞ (Ω) C T . Then it follows that -h 2 (ρ 1 ) satisfies assumptions (4.3) and (4.4). A classical example is

h 2 (ρ 1 ) = ρ α 1 and h 1 (ρ 1 , ρ 2 ) = ρ α-1 1 ρ 2
, where α 1, see for example [START_REF] Pierre | Global existence in reaction-diffusion systems with control of mass: a survey[END_REF] for more discussions.

As already mentioned, the proof of theorem 4.2 is based on a semi-implicit splitting scheme. More precisely, we construct four sequences ρ

k+1/2 1,h , ρ k+1 1,h , ρ k+1/2 2,h
, ρ k+1 2,h defined recursively as

           ρ k+1/2 i,h ∈ argmin ρ∈M + ,|ρ|=|ρ k i,h | 1 2h W 2 (ρ, ρ k i,h ) + F i (ρ) + V i (ρ|ρ k 1,h , ρ k 2,h ) ρ k+1 i,h ∈ argmin ρ∈M + 1 2h FR 2 (ρ, ρ k+1/2 i,h ) + G i (ρ) + U i (ρ|ρ k 1,h , ρ k 2,h ) , (4.6)
where the fully implicit terms

F i (ρ) := ´Ω F i (ρ) if ρ L |Ω +∞ otherwise and G i (ρ) := ´Ω G i (ρ) if ρ L |Ω +∞ otherwise ,
and the semi-implicit terms

V i (ρ|µ 1 , µ 2 ) := ˆΩ V i [µ 1 , µ 2 ]ρ and U i (ρ|µ 1 , µ 2 ) := ˆΩ U i [µ 1 , µ 2 ]ρ.
In the previous section, the proof of theorem 3.2 for scalar equations strongly leveraged the uniform L ∞ (Ω)-bounds on the discrete solutions. Here an additional difficulty arises due to the nonlocal terms ∇V i [ρ 1 , ρ 2 ] and U i [ρ 1 , ρ 2 ], which are a priori not uniformly bounded in L ∞ (Ω).

Using assumption (4.3) we will first obtain a uniform L 1 (Ω)-bound on ρ 1 , ρ 2 , and then extend proposition 3.6 to the system (4.1). This in turn will give a uniform W 1,∞ control on V i [ρ 1 , ρ 2 ] and L ∞ control on U i [ρ 1 , ρ 2 ] through our assumptions (4.2)-(4.3)-(4.4), which will finally allow to argue as in the previous section and give L ∞ control on ρ 1 , ρ 2 .

Numerical simulations for a diffusive prey-predator system are presented at the end of this section.

Properties of discrete solutions

Arguing as in the case of one equation, the optimality conditions for the Wasserstein step and for the Fisher-Rao step first give Lemma 4.4. For all k 0 and i ∈ {1, 2}, we have

ρ k+1/2 i,h L 1 = ρ k i,h L 1 . (4.7)
Moreover, there exists

C i ≡ C(U i ) > 0 (uniform in k) such that ρ k+1 i,h (x) (1 + C i h)ρ k+1/2 i,h (x) a.e. (4.8)
Proof. The first part is simply the mass conservation in the Wasserstein step, and the second part follows the lines of the proof of (3.10) in Lemma 3.5 using assumption (4.3).

As a direct consequence we have uniform control on the L 1 -norms:

Lemma 4.5. For all T > 0 there exist constants C T , C T > 0 such that, for all t ∈ [0, T ],

ρ i,h (t) L 1 , ρi,h (t) L 1 C T and V i [ρ 1,h (t), ρ 2,h (t)] W 1,∞ , V i [ρ 1,h (t), ρ2,h (t)] W 1,∞ C T . (4.9)
Proof. Integrating (4.8) and iterating with (4.7), we obtain for all t T and k T /h

ρ k+1 i,h L 1 (1 + C i h) ρ k i,h L 1 (1 + C i h) k ρ i,0 L 1 e CiT ρ i,0 L 1 .
Then (4.9) follows from our assumption (4.2) on the interactions.

Combining (4.8) and (4.9), we deduce Proposition 4.6. For all T > 0, there exists M T such that for all t ∈ [0, T ],

ρ i,h (t) L ∞ , ρi,h (t) L ∞ M T .
Then, there exists c i ≡ c(M T , U i ) ≥ 0, such that, for all k T /h and h h 0 (U 1 , U 2 ),

(

-c i h)ρ k+1/2 i,h 1
ρ k+1 i,h . In particular, there exist M T > 0 such that for all t ∈ [0, T ],

ρ i,h (t) -ρi,h (t) L 1 hM T .
Proof. The first L ∞ estimate can be found in [36, Lemma 2], and the rest of our statement can be proved exactly as in Lemma 3.5 and Proposition 3.6.

Estimates and convergences

Since we proved that V 1 [ρ 1,h , ρ 2,h ] and V 2 [ρ 1,h , ρ 2,h ] are bounded in L ∞ ([0, T ], W 1,∞ (Ω))
, we can argue exactly as in the previous section for the Wasserstein step and obtain

1 4h W 2 (ρ k+1/2 i,h , ρ k i,h ) F i (ρ k i,h ) -F i (ρ k+1/2 i,h ) + C T h, ( 4
) G i (ρ k+1/2 i,h ) -G i (ρ k+1 i,h ) + C T h. (4.11)
The dissipation of F i along the Fisher-Rao step is obtained in the same way as Proposition 3.7 and we omit the details: Proposition 4.7. For all T > 0 and i ∈ {1, 2}, there exist constants C T , C T > 0 such that, for all k 0 with hk T ,

F i (ρ k+1 i,h ) F i (ρ k+1/2 i,h ) + C T h, G i (ρ k+1/2 i,h ) G i (ρ k+1 i,h ) + C T h.
From (4.10) and (4.11) this immediately gives a telescopic sum

1 2h W 2 (ρ k i,h , ρ k+1/2 i,h ) + FR 2 (ρ k+1/2 h , ρ k h ) 2[F i (ρ k i,h ) -F i (ρ k+1 i,h )] + C T h
which in turn yields an approximate 1 2 -Hölder estimate (with respect to the WFR distance) as in Proposition 3.8. The rest of the proof of Theorem 4.2 is then identical to section 3 and we omit the details.

Numerical application: prey-predator systems

Our constructive scheme can be implemented numerically, by simply discretizing (4.6) in space. We use the augmented Lagrangian method ALG-JKO from [START_REF] Benamou | An augmented lagrangian approach to wasserstein gradient flows and applications[END_REF] to solve the Wasserstein step, and the Fisher-Rao step is just a convex pointwise minimization problem. Indeed, it is known [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF][START_REF] Lisini | Cahn-Hilliard and thin film equations with nonlinear mobility as gradient flows in weighted-Wasserstein metrics[END_REF] that

FR 2 (ρ, µ) = 4 √ ρ - √ µ 2 L 2
, hence the Fisher-Rao step in (4.6) is a mere convex pointwise minimization problem of the form: for all x ∈ Ω (and omitting all indexes ρ i,h ),

ρ k+1 (x) = argmin ρ≥0 4 √ ρ -ρ k+1/2 (x) 2 + 2hF (ρ) .
This is easily solved using any simple Newton procedure.

Figure [START_REF] Agueh | Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory[END_REF] shows the numerical solution of the following diffusive prey-predator system

∂ t ρ 1 -∆ρ 1 -div(ρ 1 ∇V 1 [ρ 1 , ρ 2 ]) = Aρ 1 (1 -ρ 1 ) -B ρ1ρ2 1+ρ1 , ∂ t ρ 2 -∆ρ 2 -div(ρ 2 ∇V 2 [ρ 1 , ρ 2 ]) = Bρ1ρ2 1+ρ1 -Cρ 2 ,
.

Here the ρ 1 species are preys and ρ 2 are predators, see for example [START_REF] Murray | Mathematical biology[END_REF], the parameters A = 10, C = 5, B = 70, and the interactions are chosen as

V 1 [ρ 1 , ρ 2 ] = |x| 2 * ρ 1 -|x| 2 * ρ 2 , V 2 [ρ 1 , ρ 2 ] = |x| 2 * ρ 1 + |x| 2 * ρ 2 .
In (4.1) this corresponds to

G 1 (ρ 1 ) = A ρ 2 1 2 , G 2 (ρ 2 ) = 0, U 1 [ρ 1 , ρ 2 ] = Bρ 2 1 + ρ 1 -A, U 2 [ρ 1 , ρ 2 ] = - Bρ 1 1 + ρ 1 + C.
Of course, U 1 and U 2 satisfy assumptions (4.3) and (4.4), and then Theorem 4.2 gives a solution of the prey-predator system. As before, we shall disregard the uniqueness issue for the sake of simplicity. Figure [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF] depicts the mass evolution of the prey and predator species: we observe the usual oscillations in time with phase opposition, a characteristic behaviour for Lotka-Volterra types of systems.

Application to a tumor growth model with very degenerate enery

In this section we take interest in the equation This equation is motivated by tumor growth models [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF][START_REF] Perthame | Traveling wave solution of the Hele-Shaw model of tumor growth with nutrient[END_REF] and exhibits a Hele-Shaw patch dynamics: if ρ 0 = χ Ω0 then the solution remains an indicator ρ(t) = χ Ω(t) and the boundary moves with normal velocity V = -∇p| ∂Ω(t) , see [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF] for a rigorous analysis in the framework of viscosity solutions.

       ∂ t ρ = div(ρ∇p) + ρ(1 -p), p 0 and p(1 -ρ) = 0 0 ρ 1, ρ |t=0 = ρ 0 . ( 5
At least formally, we remark that (5.1) is the Wasserstein-Fisher-Rao gradient flow of the singular functional

F(ρ) := F ∞ (ρ) -ˆΩ ρ,
where

F ∞ (ρ) := 0 if ρ 1 a.e, +∞ otherwise.
Indeed, the compatibility conditions p 0 and p(1 -ρ) = 0 in (5.1) really mean that the pressure p belongs to the subdifferential ∂F ∞ (ρ), and (5.1) thus reads as the gradient flow

∂ t ρ = div(ρ∇u) -ρu, u = p -1 ∈ -∂F(ρ).
However, this functional is too singular for the previous splitting scheme to correctly capture the very degenerate diffusion. Indeed, the naive and direct approach from section 3 would lead to

         ρ k+1/2 h ∈ argmin ρ 1, |ρ|=|ρ k h | 1 2h W 2 (ρ, ρ k h ) -´Ω ρ , ρ k+1 h ∈ argmin ρ 1 1 2h FR 2 (ρ, ρ k+1/2 h ) -´Ω ρ .
Since the Wasserstein step is mass-conservative by definition, the ´ρ term has no effect in the first step and the latter reads as "project ρ k h on {ρ 1} w.r.t to the W distance". Since the output of the reaction step ρ k+1 h 1, the Wasserstein step will never actually project anything, and the diffusion is completly shut down. As an example, it is easy to see that if the initial datum is an indicator ρ 0 = χ Ω0 then the above naive scheme leads to a stationary solution ρ k+1 h = ρ k+1/2 h = ρ 0 for all k 0, while the real solution should evolve according to the aforementioned Hele-Shaw dynamics ρ(t) = χ Ω(t) [START_REF] Alexander | Quasi-static evolution and congested crowd transport[END_REF][START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF]. One could otherwise try to write a semi-implicit scheme as follows: 1) keep the projection on {ρ 1} in the first Wasserstein step. As in [START_REF] Maury | Handling congestion in crowd motion modeling[END_REF] a pressure term p k+1/2 h appears as a Lagrange multiplier in the Wasserstein projection. 2) in the FR/reaction step, relax the constraint ρ 1 and minimize instead ρ k+1 ∈ argmin 1 2h FR 2 (ρ) + ´ρp k+1/2 -´ρ , and keep iterating. This seems to correctly capture the diffusion at least numerically speaking, but raises technical issues in the rigorous proof of convergence and most importantly destroys the variational structure at the discrete level (due to the fact that the reaction step becomes semi-explicit).

We shall use instead an approximation procedure, which preserves the variational structure at the discrete level: it is well-known that the Porous-Medium functional [START_REF] Braides | Γ-convergence for beginners[END_REF]. In the spirit of [START_REF] Sandier | Gamma-convergence of gradient flows with applications to ginzburg-landau[END_REF], one should therefore expect that the gradient flow ρ m of F m (ρ)-´ρ converges to the gradient flow ρ ∞ of the limiting functional F(ρ) = F ∞ (ρ) -´ρ. Implementing the splitting scheme for the regular energy functional F m (ρ) -´ρ gives a sequence ρ h,m , and we shall prove below that ρ h,m converges to a solution of the limiting gradient flow as m → ∞ and h → 0. However, it is known [START_REF] Fleißner | Gamma-convergence and relaxations for gradient flows in metric spaces: a minimizing movement approach[END_REF] that the limit depends in general on the interplay between the time-step h and the regularization parameter (m → ∞ here), and for technical reasons we shall enforce the condition mh → 0 as m → ∞ and h → 0.

F m (ρ) := ´Ω ρ m m-1 if ρ m ∈ L 1 (Ω) +∞ otherwise Γ-converges to F ∞ as m → ∞, see
Note that [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF] already contained a similar approximation m → ∞ but without exploiting the variational structure of the m-gradient flow, and our approach is thus different. The above gradient-flow structure was already noticed and fully exploited in the ongoing work [START_REF] Chizat | A tumor growth hele-shaw problem as a gradient flow[END_REF], where existence and uniqueness of weak solutions is proved and numerical simulations are performed needless of any splitting an using directly the WFR structure. Here we rather emphasize the fact that the splitting does capture delicate Γ-convergence phenomena.

In order to make this rigorous, we fix a time step h > 0 and construct two sequences (ρ

k+1/2 h,m ) k and (ρ k h,m ) k , with ρ 0 h,m = ρ 0 , defined recursively as            ρ k+1/2 h ∈ argmin ρ∈M + , |ρ|=|ρ k h | 1 2h W 2 (ρ, ρ k h,m ) + F m (ρ) -´Ω ρ , ρ k+1 h ∈ argmin ρ∈M + 1 2h FR 2 (ρ, ρ k+1/2 h ) + F m (ρ) -´Ω ρ . (5.2)
As is common in the classical theory of Porous Media Equations [START_REF] Luis | The porous medium equation: mathematical theory[END_REF], we define the pressure as the first variation

p m := F m (ρ) = m m -1 ρ m-1 .
We accordingly write

p k+1/2 h,m := m m -1 (ρ k+1/2 h,m ) m-1 and p k+1 h,m := m m -1 (ρ k+1 h,m ) m-1
for the discrete pressures. As in section 3 we denote by ρ h,m (t), p h,m (t) and ρh,m (t), ph,m (t) the piecewise constant interpolations of ρ k+1 h,m , p k+1 h,m and ρ

k+1/2 h,m , p k+1/2
h,m , respectively. Our main result is Theorem 5.1. Assume that ρ 0 ∈ BV (Ω), ρ 0 1, and mh → 0 as h → 0 and m → ∞. Then for all T > 0, ρ h,m , ρh,m both converge to some ρ strongly in L 1 ((0, T ) × Ω), the pressures p h,m , ph,m both converge to some p weakly in L 2 ((0, T ), H 1 (Ω)), and (ρ, p) is the unique weak solution of (5.1).

Since we have a WFR gradient-flow structure, uniqueness should formally follows from the -1 geodesic convexity of the driving functional E ∞ (ρ) -´Ω ρ with respect to the WFR distance [START_REF] Liero | Gradient structures and geodesic convexity for reactiondiffusion systems[END_REF][START_REF] Liero | Optimal transport in competition with reaction: the Hellinger-Kantorovich distance and geodesic curves[END_REF] and the resulting contractivity estimate WFR(ρ 1 (t), ρ 2 (t)) ≤ e t WFR(ρ 1 0 , ρ 2 0 ). This is proved rigorously in [START_REF] Chizat | A tumor growth hele-shaw problem as a gradient flow[END_REF], and therefore we retrieve convergence of the whole sequence ρ h,m → ρ in Theorem 5.1 (and not only for subsequences). Given this uniqueness, it is clearly enough to prove convergence along any discrete (sub)sequence, and this is exactly what we show below.

The strategy of proof for Theorem 5.1 is exactly as in section 3, except that we need now the estimates to be uniform in both in h → 0 and m → ∞.

Estimates and convergences

In this section, we improve the previous estimates from section 3. We start with an explicit L ∞bound: Lemma 5.2. Assume that ρ 0 1, then for all t ∈ R + ,

ρ h,m (t, •) ∞ , ρh,m (t, •) ∞ 1.
Proof. We argue by induction at the discrete level, starting from ρ 0 = ρ 0 h,m 1 by assumption.

If ρ k h,m ∞
1, Otto's maximum principle [START_REF] Otto | Double degenerate diffusion equations as steepest descent[END_REF] implies that ρ

k+1/2 h,m ∞ ρ k h,m ∞
1 in the Wasserstein step.

Assume now by contradiction that E := {ρ k+1 h,m > 1} has positive Lebesgue measure. The optimality condition (3.9) for the Fisher-Rao minimization step gives, dividing by

ρ k+1 h,m > 0 in E, ρ k+1 h,m -ρ k+1/2 h,m = h 2 ρ k+1 h,m 1 - m m -1 (ρ k+1 h,m ) m-1 Then 1 -m m-1 (ρ k+1 h,m ) m-1 1 -m m-1 < 0 in the right-hand side, hence the desired contradiction ρ k+1 h,m < ρ k+1/2 h,m 1.
Noticing that the functional 1 m-1

´ρm -´ρ corresponds to taking explicitly F 2 (z) = z m /m -1 and V 2 (x) ≡ -1 in section 3, it is easy to reproduce the computations from the proof of Lemma 3.5 and carefully track the dependence of the constants w.r.t m > 1 to obtain Lemma 5.3. There exists c > 0 such that, for all m > m 0 large enough and all h ≤ h 0 small enough,

(1 -ch)ρ k+1/2 h,m (x) ρ k+1 h,m (x) (1 + h)ρ k+1/2 h,m (x)
a.e.

(5.3)

Note that this holds regardless of any compatibility such as hm → 0. The key point is here that the lower bound c previously depended on an upper bound M on ρ k+1/2 in Lemma 3.5, but since we just obtained in Lemma 5.2 the universal upper bound ρ k+1/2 1 we end up with a lower bound which is also uniform in h, m. The proof is identical to that of Lemma 3.5 and we omit the details for simplicity.

Recalling that the Wasserstein step is mass-preserving, we obtain by immediate induction and for all 0 ≤ t ≤ T ρ h,m (t) L 1 , ρh,m (t) L 1 e T ρ 0 L 1 as well as

ρ h,m (t) -ρh,m (t) L 1 C T h. (5.4)
Testing successively ρ = ρ k h,m and ρ = ρ k+1/2 h,m in (5.2), we get

1 2h W 2 (ρ k h,m , ρ k+1/2 h,m ) + FR 2 (ρ k+1/2 h,m , ρ k+1 h,m ) F m (ρ k h,m ) -F m (ρ k+1 h,m ) + ˆΩ(ρ k+1/2 h,m -ρ k+1 h,m ).
Using Proposition 2.4 to control WFR 2 2(W 2 + FR 2 ) and the lower bound in (5.3) yields

1 4h WFR 2 (ρ k+1 h,m , ρ k h,m ) 1 2h W 2 (ρ k h,m , ρ k+1/2 h,m ) + FR 2 (ρ k+1/2 h,m , ρ k+1 h,m ) F m (ρ k h,m ) -F m (ρ k+1 h,m ) + ˆΩ(ρ k+1/2 h,m -ρ k+1 h,m ) F m (ρ k h,m ) -F m (ρ k+1 h,m ) + ch ˆΩ ρ k+1/2 h,m F m (ρ k h,m ) -F m (ρ k+1 h,m ) + che T
for all k N := T /h . Summing over k we get

1 4h N -1 k=0 WFR 2 (ρ k h,m , ρ k+1 h,m ) F m (ρ 0 ) -F m (ρ N h,m ) + C T 1 m -1 ˆΩ ρ m 0 + C T 1 m -1 ˆΩ ρ 0 + C T C T ,
where we used successively F m ≥ 0 to get rid of F m (ρ N h,m ), and ρ m 0 ≤ ρ 0 for ρ 0 ≤ 1 and m > 1. Consequently, for all fixed T > 0 and any t, s ∈ [0, T ] we obtain the classical Exploiting the explicit algebraic structure of F m (z) = 1 m-1 z m , compactness in space will be given here by Lemma 5.4.

If ρ 0 ∈ BV (Ω) then sup t∈[0,T ] ρ h,m (t, •) BV (Ω) , ρh,m (t, •) BV (Ω)
e T ρ 0 BV (Ω) .

Proof. The argument closely follows the lines of [18, prop. 5.1]. We first note from [14, thm. 1.1] that the BV -norm is nonincreasing during the Wasserstein step,

ρ k+1/2 h,m BV (Ω) ρ k h,m BV (Ω) .
Using as before the implicit function theorem, we show below that ρ k+1 h,m = R(ρ k+1/2 h,m ) for some suitable (1 + h)-Lispchitz function R. By standard Lip • BV composition [START_REF] Ambrosio | Functions of bounded variation and free discontinuity problems[END_REF] this will prove that

ρ k+1 h,m BV (Ω) (1 + h) ρ k+1/2 h,m BV (Ω)
and will conclude the proof by immediate induction. Indeed, we already know from (5.3) that ρ k+1/2 h,m and ρ k+1 h,m share the same support. In this support and from (3.9) it is easy to see that ρ = ρ k+1 h,m (x) is the unique positive solution of f (ρ, ρ

k+1/2 h,m (x)) = 0 with f (ρ, µ) = √ ρ 1 - h 2 1 - m m -1 ρ m-1 - √ µ.
For µ > 0, the implicit function theorem gives the existence of a C 1 map R such that f (ρ, µ) = 0 ⇔ ρ = R(µ), with R(0) = 0. An algebraic computation shows moreover that 0 < dR dµ = -∂µf ∂ρf |ρ=R(µ)

(1 + h) uniformly in m > 1, hence R is (1 + h)-Lipschitz as claimed and the proof is complete.

Proposition 5.5. Up to extraction of a discrete sequence h → 0, m → ∞, there holds

ρ h,m , ρh,m → ρ strongly in L 1 (Q T )
p h,m p and ph,m p weakly in all L q (Q T ) for all T > 0. If in addition mh → 0 then p = p.

Proof. The first part of the statement follows exactly as in section 3, exploiting the 1 2 -Hölder estimates (5.5) and the space compactness from Proposition 5.4 in order to apply the Rossi-Savaré theorem [START_REF] Rossi | Tightness, integral equicontinuity and compactness for evolution problems in Banach spaces[END_REF]. The fact that ρ h,m , ρh,m have the same limit comes from (5.4).

For the pressures, we simply note from ρ h,m 1 and m p in all L q (Q T ), and likewise for ph,m p. Finally, we only have to check that p = p if hm → 0. Because ρ h,m , ρh,m 1 and z → z m-1 is (m -1)-Lipschitz on [0, 1] we have for all fixed t 0 that

1 that p h,m = m m-1 ρ m-1 h,m 2ρ h,m is bounded in L 1 ∩ L ∞ (Q T ) uniformly in h,
ˆΩ |p m,h (t, •) -pm,h (t, •)| = ˆΩ m m -1 |ρ m-1 h,m (t, •) -ρm-1 h,m (t, •)| m ˆΩ |ρ h,m (t) -ρh (t)| C T hm -→ 0,
where we used (5.4) in the last inequality. Hence p = p and the proof is complete.

In order to pass to the limit in the diffusion term div(ρ∇p) we first improve the convergence of ph,m : Lemma 5.6. There exists a constant C T , independent of h and m, such that ph,m L 2 ((0,T ),H 1 (Ω)) C T for all T > 0. Consequently, up to a subsequence, ph,m converges weakly in L 2 ((0, T ), H 1 (Ω)) to p.

Proof. The proof is based on the flow interchange technique developed by Matthes, McCann and

Savaré in [START_REF] Matthes | A family of nonlinear fourth order equations of gradient flow type[END_REF]. Let η be the (smooth) solution of

∂ t η = ∆η m-1 + ε∆η, η| t=0 = ρ k+1/2 h,m .
It is well known [START_REF] Ambrosio | Gradient flows in metric spaces and in the space of probability measures[END_REF] that η is the Wasserstein gradient flow of

G(ρ) := ˆΩ ρ m-1 m -2 + ε ˆΩ ρ log(ρ).
Since G is geodesically 0-convex, η satisfies the Evolution Variational Inequality (EVI)

1 2 d + dt t=s W 2 (η(s), ρ) G(ρ) -G(η(s)),
for all s > 0 and for all ρ ∈ P ac (Ω), where d + dt f (t) := lim sup

s→0 + f (t+s)-f (t) s
. By optimality of ρ k+1/2 h,m in (5.2), we obtain that

1 2 d + dt t=s W 2 (η(s), ρ k h,m ) -h d + dt t=s F m (η(s)).
Since η is smooth due to the regularizing ε∆ term, we can legitimately integrate by parts for all s > 0

d ds F m (η(s)) = ˆΩ m m -1 η(s) m-1 (∆η(s) m-1 + ε∆η(s)) = - ˆΩ m m -1 |∇η(s) m-1 | 2 -ε ˆΩ mη(s) m-2 |∇η(s)| 2 - ˆΩ m m -1 |∇η(s) m-1 | 2 = - m -1 m ˆΩ ∇ m m -1 η(s) m-1 2 Remarking that m m-1 η(s) m-1 → m m-2 ρ k+1/2 h,m = p k+1/2 h,m
as s → 0, an easy lower semi-continuity argument gives that ˆΩ m -1

m |∇p k+1/2 h,m | 2 = ˆΩ m m -1 |∇(ρ k+1/2 h,m ) m-1 | 2 lim inf s 0 d + dt t=s F m (η(s)).
Then we have

h ˆΩ m -1 m |∇p k+1/2 h,m | 2 F m-1 (ρ k h,m ) -F m-1 (ρ k+1/2 h,m ) + ε ˆΩ ρ k h,m log(ρ k h,m ) -ˆΩ ρ k+1/2 h,m log(ρ k+1/2 h,m ) .
First arguing as in Proposition 3.7 to control

F m-1 (ρ k+1 h,m ) F m-1 (ρ k+1/2 h,m ) + C T h,
and then passing to the limit ε 0, we obtain

h ˆΩ m -1 m |∇p k+1/2 h,m | 2 F m-1 (ρ k h,m ) -F m-1 (ρ k+1 h,m ) + C T h. Summing over k gives ˆT 0 ˆΩ |∇p h,m (t, x)| 2 dxdt m m -1 (F m-1 (ρ 0 ) -F m-1 (ρ N h,m ) + C T ) 2F m-1 (ρ 0 ) + C T for all T < +∞. Due to ρ 0 1 and m 1 we can bound F m-1 (ρ 0 ) = 1 m-2 ´ρm-1 0 1 m-2 ´ρ0 ρ 0 L 1 (Ω)
and the result finally follows.

Properties of the pressure p and conclusion

We start by showing that the limits ρ, p satisfy the compatibility conditions in (5.1). Lemma 5.7. There holds 0 ρ, p 1 and p(1 -ρ) = 0 a.e. in Q T .

Proof. By Lemma 5.2 it is obvious that 0 ρ 1 and 0 p 1 are inherited from 0 ρ h,m 1 and 0 p h,m = m m-1 ρ m-1 h,m m m-1 . In order to prove that p(1 -ρ) = 0, we first observe that

p h,m (1 -ρ h,m ) → 0 a.e. in Q T .
Indeed, since ρ h,m → ρ strongly in L 1 (Q T ) we have ρ h,m (t, x) → ρ(t, x) a.e. If the limit ρ(t, x) < 1 then ρ h,m (t, x) (1-ε) for small h and large m.

Hence p h,m (t, x) = m m-1 ρ m-1 h,m m m-1 (1-ε) m-1 → 0 while 1 -ρ h,m
remains bounded, and therefore the product p h,m (1 -ρ h,m ) → 0. Now if the limit ρ(t, x) = 1 then the pressure p h,m = m m-1 ρ m-1 h,m m m-1 remains bounded, while 1 -ρ h,m (t, x) → 0 hence the product goes to zero in this case too.

A tumor growth model with nutrient

In this section we use the same approach for the following tumor growth model with nutrients, appearing e.g. in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF]

           ∂ t ρ -div(ρ∇p) = ρ ((1 -p)(c + c 1 ) -c 2 ) , ∂ t c -∆c = -ρc, 0 ρ 1, p 0 and p(1 -ρ) = 0, ρ |t=0 = ρ 0 , c |t=0 = c 0 . (6.1)
Here c 1 and c 2 are two positive constants, and the nutrient c is now diffusing in Ω in addition to begin simply consumed by the tumor ρ, according to the second equation. For technical convenience we work here on a convex bounded domain Ω ⊂ R d , endowed with natural Neumann boundary conditions for both ρ and c.

Contrarily to section 5 this is not a WFR gradient flow anymore, and we therefore introduce a semi-implicit splitting scheme. Starting from the initial datum ρ 0 h,m := ρ 0 , c 0 h,m := c 0 we construct four sequences ρ

k+1/2 h,m , ρ k h,m , c k+1/2 h,m , c k h,m , defined recursively as            ρ k+1/2 h,m ∈ argmin ρ∈M + ,|ρ|=|ρ k h,m | 1 2h W 2 (ρ, ρ k h,m ) + F m (ρ) , c k+1/2 h,m ∈ argmin c∈M + ,|c|=|c k h,m | 1 2h W 2 (c, c k h,m ) + E(ρ) , (6.2)
and

         ρ k+1 h,m ∈ argmin ρ∈M + 1 2h FR 2 (ρ, ρ k+1/2 h,m ) + E 1,m (ρ|c k+1/2 h,m ) , c k+1 h,m ∈ argmin c∈M + 1 2h FR 2 (c, c k+1/2 h,m ) + E 2 (c|ρ k+1/2 h,m ) , (6.3)
where E(ρ) := ˆΩ ρ log(ρ),

E 1,m (ρ|c) := ˆΩ (c + c 1 ) ρ m m -1 + ˆΩ(c 2 -c -c 1 )ρ, and E 2 (c|ρ) := ˆΩ ρc.
As earlier it is easy to see that these sequences are well-defined (i-e there exists a unique minimizer for each step), and the pressures are defined as before as

p k+1/2 h,m := m m -1 (ρ k+1/2 h,m ) m-1 and p k+1 h,m := m m -1 (ρ k+1 h,m ) m-1 .
We denote again by a h,m (t), ãh,m (t) the piecewise constant interpolation of any discrete quantity a k+1 h,m , a k+1/2 h,m respectively. Our main result reads: Theorem 6.1. Assume ρ 0 ∈ BV (Ω) with ρ 0 1 and c 0 ∈ L ∞ (Ω) ∩ BV (Ω). Then ρ h,m and ρh,m strongly converge to ρ in L 1 ((0, T ) × Ω) and c h,m and ch,m strongly converge to c in L 1 ((0, T ) × Ω) when h 0 and m +∞. Moreover, if mh → 0, then p h,m , ph,m converge weakly in L 2 ((0, T ), H 1 (Ω)) to a unique p, and (ρ, p, c) is a solution of (6.1).

Note that uniqueness of solutions would result in convergence of the whole sequence. Uniqueness was proved in [34, thm. 4.2] for slightly more regular weak solutions, but we did not push in this direction for the sake of simplicity. The method of proof is almost identical to section 5 so we only sketch the argument and emphasize the main differences.

We start by recalling the optimality conditions for the scheme (6. Proof. The proof of the estimates on c h,m and ch,m is obvious because one step of Wasserstein gradient flow with the Boltzmann entropy decreases the L ∞ -norm in (6.2) (see [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF][START_REF] Agueh | Existence of solutions to degenerate parabolic equations via the Monge-Kantorovich theory[END_REF]), and, because the product c k+1 h,m ρ k+1/2 h,m is nonnegative in (6.5), the L ∞ -norm is also nonincreasing during the Fischer-Rao step. The proof for ρ h,m and ρh,m is the same as in lemma 5. for some C T independent of m. Then we obtain the usual 1 2 -Hölder estimates in time with respect to the WFR distance, which in turn implies that ρ h,m , ρh,m converge to some ρ ∈ L ∞ ([0, T ], L 1 (Ω)) and c h,m , ch,m converge to some c ∈ L ∞ ([0, T ], L 1 (Ω)) pointwise in time with respect to WFR, see (3.20), Proposition 3.8,and (3.22) for details.

As before we need to improve the convergence in order to pass to the limit in the nonlinear terms. For ρ h,m and ρh,m , this follows from Proof. The proof is the same as Proposition 5.5, Lemma 5.6, and Lemma 5.7.

In order to conclude the proof of Theorem 6.1 we only need to check that ρ, p, c satisfy the weak formulation of (6.1): the strong convergence of ρ h,m , c h,m and the weak convergence of p h,m are enough to take the limit in the nonlinear terms as in section 5.2, and we omit the details.

Figure 2 :

 2 Figure 2: Mass evolution for two-species prey-predator interactions.

+ c 1 )

 1 2)-(6.3). The Euler-Lagrange equations for the tumor densities in the Wasserstein and Fisher-Rao steps are -c 2 , (6.4)where ϕ is a (backward) Kantorovich potential for W(ρk+1/2 h,m , ρ k h,m). For the nutrient, the Eulerψ a Kantorovich potential for W(ck+1/2 h,m , c k h,m). Using the optimality conditions for the Fischer-Rao steps, we obtain directly the following L ∞ bounds:Lemma 6.2. For all k 0 L ∞ (Ω) ,and at the continuous levelc h,m (t, •) L ∞ (Ω) , ch,m (t, •) L ∞ (Ω) c 0 L ∞ (Ω) ∀ t ≥ 0. Moreover, ρ h,m (t, •) ∞ , ρh,m (t, •) ∞ 1and there exists cT ≡ c T ( c 0 L ∞ ), C T ≡ C T ( c 0 L ∞ ) > 0 such that (1 -c T h)ρ k+1/2 h,m (x) ρ k+1 h,m (x) (1 + C T h)ρ k+1/2 h,m(x) a.e. in Ω. (1 -h)c k+1/2 h,m (x) c k+1 h,m (x) c k+1/2 h,m (x) a.e. in Ω. (6.6)

2 .+ c 1 )

 21 Using the fact that ρh,m (t, •) ∞ 1, we see that the term Φ(p k+1 h,m , c -c 2 in (6.4) is bounded in L ∞ uniformly in k. This allows to argue exactly as in Lemma 3.5 to retrieve the estimate (6.6) and concludes the proof.With these bounds it is easy to prove as in proposition 3.15 that F m (ρ k+1 h,m ) F m (ρ

Lemma 6 . 3 .( 1 +Proposition 6 . 4 .

 63164 For all T > 0, if ρ 0 , c 0 ∈ BV (Ω),sup t∈[0,T ] ρ h,m (t, •) BV (Ω) + c h,m (t, •) BV (Ω) e C T T ( ρ 0 BV (Ω) + c 0 BV (Ω) ) sup t∈[0,T ] ρh,m (t, •) BV (Ω) + ch,m (t, •) BV (Ω) e C T T ( ρ 0 BV (Ω) + c 0 BV (Ω) ).Proof. The argument is a generalization of Lemma 5.4, see[18, remark 5.1]. First, the BV -norm is nonincreasing during the Wasserstein step, [14, thm. 1.BV (Ω) .Arguing as in Lemma 5.4, we observe that, inside supp ρk+1/2 h,m = supp ρ k+1 h,m , the minimizer ρ = ρ k+1 h,m (x) is the unique positive solution of f (ρ, ρ (x)) = 0, with f (ρ, µ, c) = √ ρ 1 -h 2 1 -m m -1 ρ m-1 (c + c 1 ) -c 2 -√ µ.For µ > 0 the implicit function theorem gives as before aC 1 map R such that f (ρ, µ, c) = 0 ⇔ ρ = R(µ, c). An easy algebraic computation and (6.6) then gives 0< ∂ µ R(µ, c) (1 + C T h) and |∂ c R(µ, c)| C T h for some constant C T > 0 independent of h, m, k. This implies that C T h) ρ k h,m BV (Ω) + C T h c k h,m BV (Ω) .The same argument shows thatc k+1 h,m BV (Ω) (1 + C T h) c k h,m BV (Ω) + C T h ρ k h,m BV (Ω) ,and a simple induction allows to conclude. Up to extraction of a discrete sequence h → 0, m → +∞,ρ h,m , ρh,m → ρ strongly in L 1 (Q T )p h,m p and ph,m p weakly in all L q (Q T ) for all T > 0. If in addition mh → 0 then p = p ∈ L 2 ((0, T ), H 1 (Ω)) and (ρ, p) satisfies 0 ρ, p 1 and p(1 -ρ) = 0 a.e. in Q T .

  Evolution of two species with prey-predator interactions. First row: display of ρ 1 + ρ 2 . Second row: display of the prey ρ 1 . Third row: display of the predator ρ 1 .

	t = 0	t = 0.15	t = 0.35	t = 0.5		t = 0.65	t = 0.85	t = 1
	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275			IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275		IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275
	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275			IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275		IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275
	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275			IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275		IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275	IsoValue 0 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.225 0.25 0.275 0.3 0.325 0.35 0.375 0.4 0.425 0.45 0.475 0.5 0.525 0.55 0.575 0.6 0.625 0.65 0.675 0.7 0.725 0.75 0.775 0.8 0.825 0.85 0.875 0.9 0.925 0.95 0.975 1 1.025 1.05 1.075 1.1 1.125 1.15 1.175 1.2 1.225 1.25 1.275
	0 Figure 1: t 0.1 0.2 0.3 0.4 0.5 mass 0 0.02 0.04 0.06 0.08 0.1 0.12	0.6	0.7	0.8 mass of prey 0.9 mass of predator	1
								.1)

  m in any finite time interval [0, T ]. Thus up to extraction of a further sequence we have p h,m

ρ h (t) = ρ k+1 h , ρh (t) = ρ

Acknowledgements

We warmly thank G. Carlier for fruitful discussions and suggesting us the problem in section 3

Thanks to the uniform L ∞ bounds ρ h,m 1 and p h,m m m-1 2 we can apply Lebesgue's convergence theorem to deduce from this pointwise a.e. convergence that, for all fixed nonnegative ϕ ∈ C ∞ c (Q T ), there holds lim ˆQT p h,m (1 -ρ h,m )ϕ = 0.

On the other hand since ρ h,m → ρ strongly in L 1 (Q T ) hence a.e, and because 0 ρ h,m 1, we see that (1 -ρ h,m )ϕ → (1 -ρ)ϕ in all L q (Q T ). From Proposition 5.5 we also had that p h,m p in all L q (Q T ), hence by strong-weak convergence we have that

for all ϕ 0. Because p(1 -ρ) 0 we conclude that p(1 -ρ) = 0 a.e. in Q T and the proof is achieved.

We end this section with Proof of Theorem 5.1. We only sketch the argument and refer to [START_REF] Gallouët | A JKO splitting scheme for kantorovich-fischerrao gradient flows[END_REF] for the details. Fix any 0 < t 1 < t 2 and ϕ ∈ C 2 c (R d ). Exploiting the Euler-Lagrange equations (3.6)(3.9) and summing

where the remainder R(h, m) → 0 for fixed ϕ. The strong convergence ρ h,m , ρh,m → ρ and the weak convergences ∇p h,m ∇p = ∇p and p h,m p are then enough pass to the limit to get the corresponding weak formulation for all 0 < t 1 < t 2 . Moreover since the limit ρ ∈ C([0, T ]; M + WFR ) the initial datum ρ(0) = ρ 0 is taken at least in the sense of measures. This gives an admissible weak formulation of (5.1), and the proof is complete.

Numerical simulation

The constructive scheme (5.2) naturally leads to a fully discrete algorithm, simply discretizing the minimization problem in space for each W, FR step. We use again the ALG2-JKO scheme [START_REF] Benamou | An augmented lagrangian approach to wasserstein gradient flows and applications[END_REF] for the Wasserstein steps. As already mentioned the Fisher-Rao step is a mere convex pointwise minimization problem, here explicitly given by: for all x ∈ Ω,

and poses no difficulty in the practical implementation using a standard Newton method. Figure 3 depicts the evolution of the numerical solution ρ h,m for m = 100 and with a time step h = 0.005. We remark that the tumor first saturates the constraint (ρ 1) in its initial support, and then starts diffusing outwards. This is consistent with the qualitative behaviour described in [START_REF] Perthame | The Hele-Shaw asymptotics for mechanical models of tumor growth[END_REF].

IsoValue