Nicolas Fournier
email: nicolas.fournier@sorbonne-universite.fr

Benoît Perthame
email: benoit.perthame@sorbonne-universite.fr.

Transport costs for PDEs: the coupling method

Keywords: 2010 Mathematics Subject Classification: 35A05, 35K55, 60J99, 28A33 Transport cost, Monge-Kantorovich distance, Coupling, Fokker-Planck equation, Fractional Laplacian, Homogeneous Boltzmann equation, Integro-differential equations, Porous media equation

We informally review a few PDEs for which some transport cost between pairs of solutions, possibly with some judicious cost function, decays: heat equation, Fokker-Planck equation, heat equation with varying coefficients, fractional heat equation with varying coefficients, homogeneous Boltzmann equation for Maxwell molecules, and some nonlinear integro-differential equations arising in neurosciences. We always use the same method, that consists in building a coupling between two solutions. This means that we double the variables and solve, globally in time, a well-chosen PDE posed on the Euclidian square of the physical space. Finally, although the above method fails, we recall a simple idea to treat the case of the porous media equation. We also introduce another method based on the dual Monge-Kantorovich problem.

Introduction

It is usual to study the well-posedness, stability and large-time behavior of stochastic processes (e.g. solutions to Stochastic Differential Equations) by using coupling methods: we consider two such processes, with different initial conditions, driven by suitably correlated randomness, and we measure the -transport cost T between their distributions.

We work in R d and we always assume that the cost function : R d × R d → R satisfies (x, x) = 0 and (x, y) = (y, x) > 0 for x = y. We recall that for two probability densities u 1 , u 2 on R d ,

   T (u 1 , u 2 ) = inf v∈K(u 1 ,u 2 )
(x, y)v(x, y)dx dy, K(u 1 , u 2 ) = {v : R d × R d → R + such that v(x, y)dy = u 1 (x), v(x, y)dx = u 2 (y)}.

(

) 1
When for some p ≥ 1,

p (x, y) = |x -y| p p ,
and we put T p = T p . The distances d p = T 1/p p are also called Wasserstein distances, and one refers to the Monge-Kantorovich distance when p = 1.

The probabilistic coupling method consists in finding, for two solutions u 1 and u 2 of a given Partial Differential Equation (PDE in short), two coupled stochastic processes X 1 and X 2 , with time-marginals u 1 and u 2 , so that X 1 and X 2 remain as close as possible. One then controls T (u 1 (t), u 2 (t)) by E[ (X 1 (t), X 2 (t))]. The goal of the present survey paper is to describe, in an informal way, this method, using only arguments based PDEs. This is what we call the coupling method: for each problem, we introduce a PDE, with doubled variables, describing the time-evolution of density v(x, y, t) of the law at time t ≥ 0 of two underlying coupled processes, in such a way that its two marginals solve the original equation, and such that v(x, y, t) is as concentrated as possible near the diagonal x = y. This is possible only when the underlying pair of stochastic coupled processes is Markov. We then study the quantity (x, y)v(x, y, t), which controls the transport cost T between two solutions. This can be seen as a translation, in terms of PDEs, of the probabilistic coupling method. For example, our proof of Tanaka's theorem [START_REF] Tanaka | Probabilistic treatment of the Boltzmann equation of Maxwellian molecules[END_REF] for the Boltzmann equation really relies on the same main arguments as Tanaka (who is using some Poisson-driven stochastic differential equations), but we avoid introducing any stochastic process.

The difficulties and novelties rely on the choice of the cost function and on the choice of coupling between two solutions by solving a well-chosen PDE posed on the Euclidian square of the physical space, R 2d in general. Each time, we try to emphasize the main technical difficulties that would allow one to justify the computations.

Of course, any way to produce some non-expansion estimates along solutions of some PDEs for some transport cost relies on some coupling. However, for example, the deterministic methods in Carrillo [START_REF] Carrillo | EDPs de difusión y transporte óptimo de masa[END_REF] for the heat equation or Villani [START_REF]Topics in optimal transportation[END_REF]Section 7.5.6], see also [START_REF] Bolley | Tanaka theorem for inelastic Maxwell models[END_REF], for Tanaka's theorem, are really different in spirit. In [START_REF] Carrillo | EDPs de difusión y transporte óptimo de masa[END_REF], the main tool is that solutions to the heat equation can be represented by a convolution formula. In [START_REF]Topics in optimal transportation[END_REF]Section 7.5.6] and [START_REF] Bolley | Tanaka theorem for inelastic Maxwell models[END_REF], everything relies on the contractive property of the gain operator and on the Duhamel principle.

For example, considering the Brownian motion leads to the heat equation. We first give a simple proof that the heat equation is non-expansive (weak contraction) for any smooth cost function of the form (x, y) = r(|x-y|). This is standard but the PDE literature seems to ignore this simple approach. The method can be extended to various cases. The Fokker-Planck equation is the simplest extension. The case of the heat equation with variable coefficients, of the form ∂ t u -∆(a(x)u) = 0, is more involved: in one dimension, the Monge-Kantorovich distance T 1 plays a central role and is always non-expansive (under technical conditions); we illustrate the general structure in higher dimension and show that if the cost function satisfies some elliptic PDE, which does seem to enter a class with generic existence results, then T is non-expansive along solutions. The method also applies to some jump processes: fractional heat equation with variable coefficients in dimension one, scattering equations, kinetic scattering equations, Boltzmann equation for Maxwell molecules.

For the porous media equation, the situation is more intricate and the above method does not seem to apply. However, we recall from [START_REF] Bolley | Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction[END_REF] another, somehow related and rather simple, path to treat this equation.

Concerning piecewise deterministic jump processes and (inhomogeneous) kinetic scattering equations, we present a new result showing that some transport costs are non-expansive.

Finally, concerning jump processes, in particular those related to the discretized heat equation, we present another approach, based on the dual formulation of the transport costs.

Maybe the first example of use of the coupling method can be found in Dobrushin [START_REF] Dobrušin | Vlasov equations[END_REF], where the Vlasov equation is derived as mean-field limit of a deterministic system of interacting particles, making use of some transport cost. No PDE is written for the coupling in [START_REF] Dobrušin | Vlasov equations[END_REF], because eveything may be written in terms of characteristics. See [START_REF] Golse | On the mean field and classical limits of quantum mechanics[END_REF]Section 3] for a PDE analogue of Dobrushin's argument. In the same spirit, the Euler equation is derived from a deterministic system of interacting vortices in Marchioro-Pulvirenti [START_REF] Marchioro | Mathematical theory of incompressible nonviscous fluids[END_REF]Section 5.3], using also a coupling argument. See [START_REF] Hauray | Wasserstein distances for vortices approximation of Euler-type equations[END_REF] for a result with the strong transport distance d ∞ but for a problem slightly less singular than the Euler equation.

Recently, the topic of transport costs has developed quickly for PDEs and integro-differential equations (IDEs) after new understanding of optimal transportation and of the Brenier-Kantorovich map by [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF][START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]. There are several approaches to use the transport costs in PDEs. A geometrical approach based on gradient flow structures has been introduced in [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF] and extended in [START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF][START_REF]Uniform convergence to equilibrium for granular media[END_REF], in particular for the porous media equation (with T 2 ), for interacting particle systems and for granular flows. Also, many results on PDEs have been derived from the discretization algorithm named JKO after [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF]. See the book [START_REF]Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF] for a complete presentation of these results. See also [START_REF] Vázquez | The porous medium equation. New contractivity results[END_REF] who showed that contractivity for some other transport cost may fail for the porous media equation. Let us also mention that the special structure associated with dimension 1 has been used to prove strict contraction for the porous media equation [START_REF] Carrillo | Strict contractivity of the 2-Wasserstein distance for the porous medium equation by mass-centering[END_REF] for the cost T 2 , and to treat other equations as scalar conservation laws [START_REF] Bolley | Contractive metrics for scalar conservation laws[END_REF], the Keller-Segel system [START_REF] Calvez | Refined asymptotics for the subcritical Keller-Segel system and related functional inequalities[END_REF] or the granular media equation [START_REF] Li | Long-time asymptotics of kinetic models of granular flows[END_REF]. See also the general reference [START_REF] Carrillo | Wasserstein metric and large-time asymptotics of nonlinear diffusion equations[END_REF] for nonlinear diffusion equations and [START_REF] Blanchet | Convergence of the mass-transport steepest descent scheme for the subcritical Patlak-Keller-Segel model[END_REF] for the multi-dimensional Keller-Segel system. Methods based on optimal transportation have also been recently used to treat singular congestion (incompressible) equations arising in crowd modeling, see for instance [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations [corrected version of mr2584740[END_REF][START_REF] Maury | Congestion-driven dendritic growth[END_REF][START_REF] Craig | Congested aggregation via Newtonian interaction[END_REF].

Most of the recent papers using some transport cost for PDEs have been using the gradient flow structure which is closely related to a variational formulation of the fluxes. Here, with several examples of conservative equations, which do not necessarily have a gradient flow structure, we control the transport cost using the coupling method. We often borrow our examples from the stochastic processes which represent the PDEs thanks to their Kolmogorov equation. The cases of variable coefficients are particularly interesting because they often require some special choice of the cost function.

We organize our examples as follows. We begin with three simple examples: heat equation, Fokker-Planck equation, and a class of nonlinear transport equations. We show directly that the transport costs are non-expansive along these equations. Then we turn, in Section 2, to the heat equation with variable coefficients. In Section 3 we consider some IDEs: scattering equations, including kinetic scattering and inhomogeneous fractional heat equation. Zero-th order terms, describing absorption and re-emission, as they appear in models of neural networks, can also be treated by adapting the method; this is explained in Section 4. The famous Tanaka theorem for the homogeneous Boltzmann equation can be included in our framework and this is done in Section 5. We treat the porous media equation in Section 6. Finally, we exemplify in Section 7 how the same results can be proved using the dual formulation of the transport costs.

Heat, Fokker-Planck and transport equations

In order to explain the coupling method in a very simple, but still relevant, framework, we begin with the heat equation. Then we turn to drift and transport terms.

Heat equation

Here is the well-known result, see e.g. [START_REF]Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF], we want to quickly recall.

Theorem 1 Consider any increasing function r : [0, ∞) → [0, ∞) such that the cost function :

R d × R d → R + defined by (x, y) = r(|x -y|) is of class C 2 . Consider two probability densities u 0 1 , u 0 2
on R d , and the corresponding solutions u 1 , u 2 to the heat equation

∂ t u -∆u = 0, x ∈ R d , t ≥ 0. ( 2
)
For any t ≥ 0, one has

T (u 1 (t), u 2 (t)) ≤ T (u 0 1 , u 0 2 ).
Proof. We consider an initial density v 0 : R d × R d → R + with marginals u 0 1 and u 0 2 , that is such that v 0 ∈ K(u 0 1 , u 0 2 ). We next consider the solution v(x, y, t) to the degenerate heat equation ∂v ∂t

-∆ x v -∆ y v -2∇ x • ∇ y v = 0, x, y ∈ R d , t ≥ 0 (3)
starting from v 0 . Observe that (3) may be rewritten as

∂v ∂t -(∇ x + ∇ y ) • (∇ x + ∇ y )v = 0, x, y ∈ R d , t ≥ 0. (4)
Clearly, it holds that v(x, y, t) ≥ 0, because of the non-negativity of the operator in (3), which can be written in the variables (x + y, x -y) as -∆ x+y .

We then define the marginals

v 1 (x, t) = v(x, y, t)dy, v 2 (y, t) = v(x, y, t)dx
and show that v 1 = u 1 and v 2 = u 2 : for instance, integrating (3) with respect to y, one finds

∂v 1 (x,t) ∂t -∆ x v 1 (x, t) = 0, x ∈ R d , t ≥ 0, v 1 (x) = u 0 1 (x), x ∈ R d
and uniqueness of the solution of the heat equation gives us v 1 = u 1 .

Recalling (1), we conclude that

T (u 1 (t), u 2 (t)) ≤ (x, y)v(x, y, t)dxdy = r(|x -y|)v(x, y, t)dxdy.
Finally, we may also compute, using (4) and integrating by parts,

d dt r(|x -y|)v(x, y, t)dxdy = v(x, y, t) (∇ x + ∇ y ) • (∇ x + ∇ y )[r(|x -y|)] dxdy = 0.
Therefore, for any initial data v 0 ∈ K(u 0 1 , u 0 2 ),

T (u 1 (t), u 2 (t)) ≤ r(|x -y|)v 0 (x, y)dxdy = (x, y)v 0 (x, y)dxdy
and minimizing among such v 0 completes the proof.

The only technical question is to justify the integration by parts, which is immediate if we assume enough moments initially (otherwise there is nothing to prove), at least when we restrict ourselves to power cost functions p (x, y) = |x -y| p /p with p ≥ 2. Notice that the well-posedness for (3) follows from the observation that we actually deal with -∆ x+y . It is also possible, under some conditions, to treat the case of some non smooth cost functions, e.g. p for some p ∈ [1, 2): this issue is discussed in Section 2.

Remark that the solution to (3) can be represented as v(x, y, t

) = E[v 0 (x + √ 2B t , y + √ 2B t )]
, for (B t ) t≥0 a Brownian motion. Note also that in (3), diffusion occurs only in the x + y variable, while x -y remains invariant under the dynamics.

Another way to understand [START_REF] Blanchet | Convergence of the mass-transport steepest descent scheme for the subcritical Patlak-Keller-Segel model[END_REF]

is to see that v(t) is the law of the couple (X 0 + √ 2B t , Y 0 + √ 2B t )
, where (X 0 , Y 0 ) is v 0 -distributed and where we use the same Brownian motion (B t ) t≥0 for both coordinates. If using e.g. independent Brownian motions (B t ) t≥0 and (C t ) t≥0 , i.e. if considering the coupling

(X 0 + √ 2B t , Y 0 + √ 2C t )
, one would find v with the good marginals but solving ∂ t v -(∆ x + ∆ y )v = 0, leading for example to the far from optimal estimate T 2 (u 1 (t), u 2 (t)) ≤ T 2 (u 0 1 , u 0 2 ) + 2dt.

Fokker-Planck equation

The coupling method can be extended to the Fokker-Planck equation, see [START_REF] Bolley | Convergence to equilibrium in Wasserstein distance for Fokker-Planck equations[END_REF] for some more elaborate consequences. The result can be stated as follows

Theorem 2 Consider some function V : R d × R + → R d such that, for some α ∈ R, (V (x, t) -V (y, t)) • (x -y) ≤ α|x -y| 2 , x, y ∈ R d , t ≥ 0. ( 5
)
Consider two probability densities u 0 1 , u 0 2 on R d and the corresponding solutions u 1 , u 2 to the Fokker-Planck equation

∂ t u -∆u + div(V (x, t)u) = 0, x ∈ R d , t ≥ 0. ( 6
)
For any t ≥ 0, any p ≥ 1, one has

T p (u 1 (t), u 2 (t)) ≤ T p (u 0 1 , u 0 2 ) exp(αpt).
This inequality is well-known, see for example [START_REF]Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF] §9.1.5 and the references therein. One can also find relations to several deep and recent functional analysis tools. This goes far beyond our present purpose.

Proof. This is the same proof as for the heat equation, with longer expressions. We consider any v 0 , with marginals u 0 1 and u 0 2 , and the solution v to the equation

∂ t v -∆ x v -∆ y v -2∇ x • ∇ y v + div x (V (x, t)v) + div y (V (y, t)v) = 0, x, y ∈ R d , t ≥ 0 (7)
starting from v 0 . One easily checks that v(x, y, t) ≥ 0 and, integrating [START_REF] Bolley | Convergence to equilibrium in Wasserstein distance for Fokker-Planck equations[END_REF] with respect to y, that v 1 (x, t) := v(x, y, t)dy solves ( 6) and starts from u 0 1 , whence v 1 = u 1 . The second marginal is treated similarly, and we conclude that

T p (u 1 (t), u 2 (t)) ≤ p -1
|x -y| p v(x, y, t)dxdy. Finally, using the same computation as for the heat equation, with some additional terms, we see that

d dt |x -y| p p v(x, y, t)dxdy = 0 + v(x, y, t)|x -y| p-2 (x -y) • (V (x, t) -V (y, t))dxdy ≤ α |x -y| p v(x, y, t)dxdy
by assumption [START_REF] Bolley | Tanaka theorem for inelastic Maxwell models[END_REF]. The result follows using the Gronwall lemma

T p (u 1 (t), u 2 (t)) ≤ |x -y| p p v(x, y, t)dxdy ≤ |x -y| p p v 0 (x, y)dxdy e αpt
and minimizing in v 0 .

Observe that one can treat in a similar way the case with variable diffusion coefficients, of the form

∂ t u -d i,j=1 ∂ ij (a ij (x, t)u) + div(V (x, t)u) = 0
, under suitable conditions on the nonnegative symmetric matrix a and on V . For example, we will find that

T 2 (u 1 (t), u 2 (t)) ≤ T 2 (u 0 1 , u 0 2 ) exp(αt) as soon as Tr [σ(x, s) -σ(y, s)][σ(x, s) -σ(y, s)] * + (x -y) • (V (x, t) -V (y, t)) ≤ α|x -y| 2 , where σ(x, t) is, for each (x, t) ∈ R d × [0, ∞), a matrix such that σ(x, t)[σ(x, t)] * = a(x, t).
To our knowledge, a divergence form equation ∂ t u -d i,j=1 ∂ i (a ij (x, t)∂ j u) = 0 does not enjoy particular properties from this point of view.

A nonlinear transport equation

We next consider a fully deterministic problem which arises in several types of modeling, such as polymers, cell division, neuron networks, etc :

∂ t u + div[V (x, I(t))u] = 0, x ∈ R d , t ≥ 0, (8)
where the nonlinearity stems from the quantity I(t) defined, with a given weight ψ : R d → R,

I(t) = R d ψ(x)u(x, t)dx. (9)
We again complement this equation with an initial condition u 0 ≥ 0 with mass u 0 = 1.

Theorem 3 Assume that V : R d × R → R d and ψ : R d → R are of class C 1
, and that for some α > 0

(x -y) • V (x, I) -V (y, I) ≤ -α|x -y| 2 , ∀x, y ∈ R d , I ≥ 0. ( 10
)
Setting x(t) = xu(x, t)dx, we have

T 2 (u(t), δ x(t) ) = |x -x(t) | 2 2 u(x, t)dx ≤ e -2αt |x -x(0) | 2 2 u 0 (x, t)dx = e -2αt T 2 (u(0), δ x(0) ).
Assume additionally that

β = D I V ∞ Dψ ∞ < α
and fix any initial point X 0 ∈ R d . Consider the solution X to X (t) = -V (X(t), ψ(X(t)) starting from X 0 . For all t ≥ 0, one has

T 2 (u(t), δ X(t) ) = |x -X(t)| 2 2 u(x, t)dx ≤ e 2(β-α)t |x -X 0 | 2 2 u 0 (x)dx = e 2(β-α)t T 2 (u 0 , δ X 0 ).
It holds that (δ X(t) ) t≥0 solves (8) in a weak sense. A more general result, involving any pair of solutions, can be found in [START_REF]Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF]. One could prove Theorem 3 without using a PDE for the coupling, using only the characteristics, in the spirit of Dobrushin [START_REF] Dobrušin | Vlasov equations[END_REF].

Proof. We consider two solutions u 1 and u 2 to (8), and denote by I 1 (t) and I 2 (t) the corresponding functions, see [START_REF] Brasco | Congested traffic dynamics, weak flows and very degenerate elliptic equations [corrected version of mr2584740[END_REF]. As we are interested in the case where one of the two solutions is a Dirac mass (for each t ≥ 0, u 2 (t) = δ X(t) ), we can consider the trivial coupling v(x, y, t) = u 1 (x, t)u 2 (y, t), which of course has the correct marginals, and satisfies

∂ t v + div x [V (x, I 1 (t))v] + div y [V (y, I 2 (t))v] = 0.
Therefore, we may compute

d dt |x -y| 2 2 u 1 (x, t)u 2 (y, t)dxdy = (x -y) • V (x, I 1 (t)) -V (y, I 2 (t)) u 1 (x, t)u 2 (y, t)dxdy = (x -y) • V (x, I 1 (t)) -V (y, I 1 (t)) u 1 (x, t)u 2 (y, t)dxdy + (x -y) • V (y, I 1 (t))) -V (y, I 2 (t)) u 1 (x, t)u 2 (y, t)dxdy ≤ -α |x -y| 2 u 1 (x, t)u 2 (y, t)dxdy + D I V ∞ |I 1 (t) -I 2 (t)| |x -y| 2 u 1 (x, t)u 2 (y, t)dxdy 1/2 .
We first apply this in the case of single solution u := u 1 = u 2 , whence I 1 = I 2 , and we directly conclude by Gronwall's lemma that

|x -y| 2 2 u(x, t)u(y, t)dxdy ≤ e -2αt |x -y| 2 2 u 0 (x, t)u 0 (y, t)dxdy.
This classically rewrites as

|x -x(t) | 2 u(x, t)dx ≤ e -2αt |x -x(0) | 2 u 0 (x, t)dx.
as desired. Next, when considering two solutions, we notice that

I 1 (t) -I 2 (t) = [ψ(x) -ψ(y)]u 1 (x, t)u 2 (y, t)dxdy, whence I 1 (t) -I 2 (t) ≤ Dψ ∞ |x -y|u 1 (x, t)u 2 (y, t)dxdy ≤ Dψ ∞ |x -y| 2 u 1 (x, t)u 2 (y, t)dxdy 1/2 .
Therefore,

d dt |x -y| 2 2 u 1 (x, t)u 2 (y, t)dxdy ≤ (β -α) |x -y| 2 u 1 (x, t)u 2 (y, t)dxdy.
Applying this to the case where u 2 (t) = δ X(t) concludes the proof.

For consistency with the other presentations in this section, we have written this result for an L 1 density v with a finite second moment, but the extension to a probability measure is immediate.

Agregation type equation

Among the nonlinear transport equations let us mention the agregation equation

∂ t u -div[u (K * u)] = 0, x ∈ R d , t ≥ 0, (11)
where the convolution term is defined by

K * u(x, t) = R d K(x -y)u(y, t)dy.
It is well-known, under various assumptions on K and for several variants of the equation, that it is contractive for transport costs and we refer to the papers [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF][START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF][START_REF] Carrillo | Wasserstein metric and large-time asymptotics of nonlinear diffusion equations[END_REF][START_REF] Li | Long-time asymptotics of kinetic models of granular flows[END_REF][START_REF] Golse | On the mean field and classical limits of quantum mechanics[END_REF].

Here we assume that the smooth kernel satisfies, for some α ≥ 0,

K ∈ C 1 (R d ; R d ), (x -y) • K(x) -K(y) ≥ -α|x -y| 2 K(x) = -K(-x). (12)
We also assume that the L 1 (R d )-initial conditions u 0 i ≥ 0 with unit mass u 0 i = 1, for i = 1, 2 satisfy

R d |x| 2 u 0 i dx < ∞, R d xu 0 i (x)dx = 0. ( 13
)
These imply that R d xu i (x, t)dx = 0 for all t > 0.

Theorem 4 Assume [START_REF] Carrillo | EDPs de difusión y transporte óptimo de masa[END_REF]. Consider two initial data u 0 1 and u 0 2 satisfying (13) and u 1 , u 2 the corresponding solutions to [START_REF] Calvez | Refined asymptotics for the subcritical Keller-Segel system and related functional inequalities[END_REF]. One has for all t ≥ 0,

T 2 (u 1 (t), u 2 (t)) ≤ e 2αt T 2 (u 0 1 , u 0 2 ).
Notice that, combining the coupling below with that in Section 1.1, it is immediate to also treat the case with diffusion

∂ t u -∆u -div[u K * u] = 0, Proof. Again the simple coupling ∂ t v(x, y, t) -div x [v K * u 1 ] -div y [v K * u 2 ] = 0, x, y ∈ R d , t ≥ 0 ( 14
)
is enough. Indeed, we can verify immediately the properties that it preserves nonnegativity and the marginales satisfy the correct equations. Finally, one can compute

d dt |x -y| 2 2 v(x, y, t)dxdy = - (x -y) • K(x -x )u 1 (x )dx -K(y -y )u 2 (y )dy v(x, y, t)dxdy = - (x -y) • K(x -x ) -K(y -y ) v(x , y , t)v(x, y, t)dx dy dxdy = - 1 2 (x -x -y + y ) • K(x -x ) -K(y -y ) v(x , y , t)v(x, y, t)dx dy dxdy ≤ α 2 |x -x -y + y | 2 v(x , y , t)v(x, y, t)dx dy dxdy = α |x -y| 2 v(x, y, t)dxdy.
.

We have used a symmetry argument and that (x -y)v(x, y, t)dxdy = 0. We complete the proof as in Section 1.1.

Heat equation with variable coefficients

We consider the heat equation with variable coefficient. This is much more intricate than the previous examples. In 1 dimension, we use the T 1 distance and recover a result implicitly included in [START_REF] Gall | Applications du temps local aux équations différentielles stochastiques unidimensionnelles[END_REF].

In higher dimension, we indicate a general way to construct cost functions. This leads to a poorly explored degenerate elliptic PDE, see however [START_REF] Pierre | Global existence in reaction-diffusion systems with control of mass: a survey[END_REF] and the references therein.

One-dimensional case

We consider some a : R → R + and the following heat equation.

∂u ∂t - ∂ 2 ∂x 2 [a(x)u] = 0, x ∈ R, t ≥ 0. ( 15
)
Theorem 5 Assume that d = 1 and that a = σ 2 for some σ ∈ C 1/2 (R). Consider two probability densities u 0 1 , u 0 2 on R and the corresponding solutions u 1 , u 2 to [START_REF] Carrillo | Wasserstein metric and large-time asymptotics of nonlinear diffusion equations[END_REF]. For all t ≥ 0, one has

T 1 (u 1 (t), u 2 (t)) ≤ T 1 (u 0 1 , u 0 2 ).
Proof. We give a proof for σ ∈ C α (R), with α > 1/2, the remark below explains how to treat α = 1/2. We consider any probability density v 0 (x, y) with marginals u 0 1 and u 0 2 and consider the coupling equation

∂ t v -∂ xx (σ 2 (x)v) -∂ yy (σ 2 (y)v) -2∂ x ∂ y [σ(x)σ(y) v] = 0, x, y ∈ R, t ≥ 0 ( 16
)
starting from v 0 . This equation preserves non-negativity. A simple way to see this is the following computation: mutliplying ( 16) by -v -, integrating on R 2 and using some integrations by parts, one can check that 1 2

d dt v 2 -(x, y, t)dx dy = - |σ(x)∂ x v -(x, y, t)| 2 -2σ(x)σ(y)∂ x v -(x, y, t)∂ y v -(x, y, t) + |σ(y)∂ y v -(x, y, t)| 2 dx dy + 1 2 v 2 -(x, y, t)[∂ xx (σ 2 (x)) + ∂ yy (σ 2 (y)) -2∂ x σ(x)∂ y σ(y)]dx dy ≤ 1 2 v 2 -(x, y, t)[∂ xx (σ 2 (x)) + ∂ yy (σ 2 (y)) -2∂ x σ(x)∂ y σ(y)]dx dy. Since v 2 -(
x, y, 0)dx dy = 0, the result follows from the Gronwall lemma if σ is smooth. Otherwise, one can work by approximation.

Integrating [START_REF]Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF] with respect to y, we see that v 1 (x, t) := v(x, y, t)dy solves [START_REF] Carrillo | Wasserstein metric and large-time asymptotics of nonlinear diffusion equations[END_REF] and starts from u 0 1 , whence v 1 = u 1 . The second marginal is treated similarly, and we conclude that T 1 (u 1 (t), u 2 (t)) ≤ |x -y|v(x, y, t)dxdy. Because of its singularity, we need to regularize the absolute value as a W 2,∞ function and define

ω ε (r) = r 2 2ε for r ≤ ε, r -ε 2 for r ≥ ε.
Using the Hölder constant C σ of σ(•), we see that

d dt ω ε (|x -y|)v(x, y, t)dxdy = v(x, y, t) ω ε (|x -y|) [σ(x) -σ(y)] 2 dxdy ≤ C 2 σ v(x, y, t) 1 I {|x-y|≤ε} ε |x -y| 2α dxdy ≤ C 2 σ ε 2α-1 ,
because v(t) is a probability measure. Since now 2α -1 > 0, we may let ε → 0 and we find that

T 1 (u 1 (t), u 2 (t)) ≤ |x -y|v(x, y, t)dxdy ≤ |x -y|v 0 (x, y)dxdy.
We conclude, as usual, by minimizing in v 0 .

Remark 6

The condition σ ∈ C 1/2 (R) is enough. To treat this exponent, a better construction of the regularization is required, using the so-called Yamada function:

ω ε (r) = 0 for r ≤ ε 3/2 , ω ε (r) = 2 r| ln(ε)| for ε 3/2 ≤ r ≤ ε, ω ε (r) = 1 for r ≥ ε.
There are other technical issues here. For example, the well-posedness of ( 16), which is necessary to identify the marginals of the solution v to the coupling equation, is not so easy. A possible direction is to use results established in [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF], in the spirit of [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF].

A general construction of the weight

In order to unravel the algebraic structure behind the choice of the weight , we now consider the general case of dimension d. We assume that a : R d → M d×d (R) is everywhere symmetric and nonnegative, of the form

a ij (x) = K k=1 σ ik (x)σ jk (x), (17)
for some σ : R d → M d×K (R), and we consider the heat equation

∂u ∂t - d i,j=1 ∂ 2 ∂x i ∂x j [a ij (x)u] = 0, x ∈ R d , t ≥ 0, (18)
completed with an initial probability density u 0 on R d .

Proposition 7 Assume that σ is regular enough and consider two probability densities u 0 1 , u 0 2 on R d and the corresponding solutions u 1 , u 2 to [START_REF] Daneri | Eulerian calculus for the displacement convexity in the Wasserstein distance[END_REF]. For all t ≥ 0, one has

T (u 1 (t), u 2 (t)) ≤ T (u 0 1 , u 0 2 ),
for any smooth cost :

R d → R + satisfying d i,j=1 a ij (x) ∂ 2 (x, y) ∂x i ∂x j + d i,j=1 a ij (y) ∂ 2 (x, y) ∂y i ∂y j + 2 d i,j=1 K k=1 σ ik (x)σ jk (y) ∂ 2 (x, y) ∂x i ∂y j ≤ 0, x, y ∈ R d . ( 19
)
When a is constant, we recover that any C 2 cost function of the form (x, y) = r(|x -y|) works. In dimension 1, (x, y) = |x -y| is indeed a (weak) solution to [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF]. We do not know of a theory to solve [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF], in dimension d ≥ 2, for a general coefficient a, so that we do not know if this result is useful. Notice that equation [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] should be completed by the boundary value (x, x) = 0 with some growth condition to mimic |x -y| p .

Proof. We consider any probability density v 0 (x, y) with marginals u 0 1 and u 0 2 and consider the coupling equation

∂v ∂t - d i,j=1 ∂ 2 ∂x i ∂x j [a ij (x)v] - d i,j=1 ∂ 2 ∂y i ∂y j [a ij (y)v] -2 d i,j=1 K k=1 ∂ 2 ∂x i ∂y j [σ ik (x)σ jk (y)v] = 0
starting from v 0 . We show as usual that v(x, y, t)dy = u 1 (x, t) and that v(x, y, t)dx = u 2 (y, t). Moreover, we have v(x, y, t) ≥ 0: we multiply the coupling equation by -v -and integrate, finding which can also be written

1 2 d dt v 2 -(x, y, t)dxdy = - K k=1 I k -J, with I k = d i,j=1 σ ik (x) ∂v -(x,
I k = d i=1 σ ik (x) ∂v -(x, y, t) ∂x i + d i=1 σ ik (y) ∂v -(x, y, t) ∂y i 2 dxdy ≥ 0.
The other term is

J = d i,j=1 ∂v -(x, y, t) ∂x i ∂a ij (x) ∂x j v -(x, y, t) + ∂v -(x, y, t) ∂y i ∂a ij (y) ∂y j v -(x, y, t)+ 2 ∂v -(x, y, t) ∂x i v -(x, y, t) ∂ ∂y j K k=1 σ ik (x)σ jk (y) dxdy,
which can also be written after integration by parts

J = - 1 2 (v -(x, y, t)) 2 d i,j=1 ∂ 2 a ij (x) ∂x i ∂x j + ∂ 2 a ij (y) ∂y i ∂y j + 2 ∂ 2 ∂x i ∂y j K k=1 σ ik (x)σ jk (y) dxdy.
Assuming that the entries σ ik are bounded with two bounded derivatives, we conclude by Gronwall's lemma that v -≡ 0, since we initially have v 2 -(x, y, 0)dx dy = 0. Recalling (1), we conclude that T (u We leave open the question to formalize this approach rigorously, in particular for degenerate coefficients σ, and to build other examples where one can prove the existence of a weight .

1 (t), u 2 (t)) ≤ (x, y)v(x,

Scattering and integral kernels

We now turn to equations that describe the probability law of various jump processes. These are well-known results except the case of kinetic scattering in Subsection 3.2 which seems to be new.

Simple scattering

For x ∈ R d , we parameterize the pre-jump location X = Φ(x, h) by h ∈ R d , distributed according to a bounded measure µ. We assume that for all fixed h ∈ R d ,

x → X = Φ(x, h) is invertible on R d and D x Φ(x, h) is an invertible matrix, (20)
and we use the notation X → x = Φ -1 (X, h) for the inverse in x (with h fixed).

We consider the scattering problem

∂ t u(x, t) = u(Φ(x, h), t) det(D x Φ(x, h)) -u(x, t) dµ(h), (21)
with initial condition u 0 , a probability density on R d . Actually, this equation is to be understood in the weak sense: integrating the right hand side against a test function ϕ(x), we see that

ϕ(x) u(Φ(x, h), t) det(D x Φ(x, h)) -u(x, t) dµ(h) = u(X, t)[ϕ(Φ -1 (X, h)) -ϕ(X)]dXdµ(h),
which shows that the determinant det(D x Φ(x, h)) is only used informally. We briefly prove the following result, which is classical, see for instance [START_REF] Bal | Generalized stability estimates in inverse transport theory[END_REF].

Theorem 8 Assume (20), fix p ∈ [1, ∞) and suppose there is δ ∈ R such that for all X, Y ∈ R d , |Φ -1 (X, h) -Φ -1 (Y, h)| p dµ(h) ≤ KL|X -Y | p , where K = µ(R d ). ( 22
)
Consider two probability densities u 0 1 , u 0 2 on R d and the corresponding solutions u 1 , u 2 to [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF]. For all t ≥ 0, one has

T p (u 1 (t), u 2 (t)) ≤ e K(L-1)t T p (u 0 1 , u 0 2 ),
The homogeneous scattering corresponds to Φ(x, h) = x + h and obviously fulfills the above assumptions.

Proof. For a probability density v 0 on R d ×R d with marginals u 0 1 and u 0 2 , we consider the solution v of the coupled equation built in such a way that the jumps parameter h is common to the two variables. Namely, we choose

∂ t v(x, y, t) = [v(Φ(x, h), Φ(y, h), t) det(D x Φ(x, h)) det(D x Φ(y, h)) -v(x, y, t)]dµ(h), (23)
starting from v 0 . We clearly have v ≥ 0, and integrating in y and using the change of variable y → Φ(y, h), we find that v 1 (x, t) = v(x, y, t)dy satisfies [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF]. Since it starts from u 0 1 , we conclude that v 1 = u 1 . The second marginal is treated similarly, and we conclude as usual that T p (u 1 (t), u 2 (t)) ≤ p -1 |x -y| p v(x, y, t)dxdy. Next, we compute, using [START_REF] Fournier | On a toy model of interacting neurons[END_REF]:

d dt |x -y| p v(x, y, t)dxdy + K |x -y| p v(x, y, t)dxdy = |x -y| p v(Φ(x, h), Φ(y, h), t) det(D x Φ(x, h)) det(D y Φ(y, h))dµ(h)dxdy = |Φ -1 (X, h) -Φ -1 (Y, h)| p v(X, Y, t)dµ(h)dXdY.
We used the changes of variables X = Φ(x, h) and Y = Φ(y, h) (with h fixed). Recalling [START_REF] Fournier | On pathwise uniqueness for stochastic differential equations driven by stable Lévy processes[END_REF], we conclude that

d dt |x -y| p v(x, y, t)dxdy ≤ K(L -1) |x -y| p v(x, y, t)dxdy.
Using the Gronwall lemma, we thus find that

T p (u 1 (t), u 2 (t)) ≤ p -1
|x -y| p v(x, y, t)dxdy ≤ p -1 e K(L-1)t |x -y| p v 0 (x, y)dxdy and we conclude as usual, minimizing in v 0 .

The most general scattering equation reads

∂ t u(x, t) = π(x, x * )u(x * ) -π(x * , x)u(x, t) dx * , (24)
and equation ( 21) corresponds to the homogeneous cases when π(x * , x)dx * = 1, and the above method can easily be adapted. For the inhomogeneous case, see Section 4.

Kinetic scattering

We next consider some kinetic scattering models, that means we work in the phase space. We consider some finite measure µ on R d , some application

V : R d → R d such that, for all h ∈ R d , v → V = Φ(v, h) is invertible and D v Φ(v, h) is an invertible matrix, ( 25
)
and the kinetic scattering equation

∂ t f (x, v, t) + v.∇ x f = [f (x, Φ(v, h), t) det(D v Φ(v, h)) -f (x, v, t)]dµ(h) (26)
completed with an initial data f 0 (x, v) ≥ 0 with f 0 dxdv = 1.

Theorem 9 Assume (25). Set K = µ(R d ) and suppose that for some L ∈ R + , for all v, w ∈ R d ,

|Φ -1 (V, h) -Φ -1 (W, h)|dµ(h) ≤ KL|V -W |. ( 27
)
Suppose that K ≥ KL + 1. Consider two probability densities f 0 1 , f 0 2 on R d × R d and the corresponding solutions f 1 , f 2 to [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF]. It holds that for all t ≥ 0, (here T 1 is associated to the cost function

((x, v), (y, w)) = |x -y| + |v -w|) T 1 (f 1 (t), f 2 (t)) ≤ T 1 (f 0 1 , f 0 2 ).
Proof. As usual, we consider any probability density F 0 ((x, v), (y, w)) on (R d × R d ) 2 with marginals f 0 1 and f 0 2 , and we consider F ((x, v), (y, w), t) starting from F 0 and solving

∂ t F + v.∇ x F + w.∇ y F = F ((x, Φ(v, h)), (y, Φ(w, h)), t) det(D v Φ(v, h)) det(D v Φ(w, h)) -F ((x, v), (y, w), t) dµ(h).
This function is clearly nonnegative and has the correct marginals. For example, with F 1 (x, v, t) = F (x, y, v, w, t)dydw, we see that

∂ t F 1 + v.∇ x F 1 = [F 1 (x, Φ(v, h), t) det(D v Φ(v, h)) -F 1 (x, v, t)]dµ(h) because F ((x, Φ(v, h)), (y, Φ(w, h)), t) det(D v Φ(w, h))dydw = F 1 (x, Φ(v, h), t): use the substitution V = Φ(w, h) (with h fixed). Since F 1 (0) = f 1 (0), we conclude that F 1 (t) = f 1 (t). Hence we conclude that T 1 (f 1 (t), f 2 (t)) ≤ (|x -y| + |v -w|)F ((x, v), (y, w), t)dxdydvdw.
Next, using the equation for F , we find with V = Φ(v, h) and W = Φ(w, h),

d dt (|x -y| + |v -w|)F (x, y, v, w, t)dxdydvdw = x -y |x -y| • (v -w)F (x, y, v, w, t)dxdydvdw -K (|x -y| + |v -w|)F (x, y, v, w, t)dxdydvdw + |x -y| + |Φ -1 (V, h) -Φ -1 (W, h)| F ((x, V ), (y, W ), t)dxdydV dW dµ(h) ≤ (1 -K + KL) |v -w|F (x, y, v, w, t)dxdydvdw.
by [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF]. Since now K ≥ 1 + KL by assumption, we deduce that

T 1 (f 1 (t), f 2 (t)) ≤ (|x -y| + |v -w|)F 0 ((x, v), (y, w))dxdydvdw
and complete the proof as usual, minimizing in F 0 .

Remark 10 Fix a > 0. If using the transport cost with weight = a|x -y| + |v -w|, the condition K ≥ 1 + KL is replaced by the condition K > a + KL.

Fractional heat equation with variable coefficients

Informally, the fractional Laplacian is a variant of the integral equation treated in Subsection 3.1. However there is a particular interest when the coefficients depend on space, an example we borrow from [START_REF] Komatsu | On the pathwise uniqueness of solutions of one-dimensional stochastic differential equations of jump type[END_REF][START_REF] Fournier | On pathwise uniqueness for stochastic differential equations driven by stable Lévy processes[END_REF]. Consider the parabolic equation with derivatives of order α ∈ (0, 2)

     ∂ t u(x, t) = L α [u], x ∈ R, t ≥ 0, L * α [ϕ](x) := [ϕ(x + σ(x)h) -ϕ(x) -hσ(x)ϕ (x)] dh |h| 1+α . (28)
Theorem 11 Assume that α ∈ (1, 2) and that σ ∈ C 1/α and consider two initial probability densities u 0 1 and u 0 2 on R and the corresponding solutions u 1 and u 2 to (28). For all t ≥ 0,

T α-1 (u 1 (t), u 2 (t)) ≤ T α-1 (u 0 1 , u 0 2 ).
Proof. We consider an initial probability density v 0 on R 2 with marginals u 0 1 and u 0 2 and the solution v to the problem (written in weak form): for all smooth ϕ : R 2 → R,

d dt ϕ(x, y)v(x, y, t)dxdydt (29)
= v(x, y, t) ϕ(x + σ(x)h, y + σ(y)h) -ϕ(x, y) -h[σ(x)∂ x ϕ(x, y) + σ(y)∂ y ϕ y (x, y)] dh |h| 1+α dxdy
starting from v 0 . The solution is clearly nonnegative and one checks as usual that for each t ≥ 0, the marginals of v(t) are u 1 (t) and u 2 (t): for example, we apply the above formula with ϕ depending only on x and deduce that v 1 (x, t) = v(x, y, t)dy solves the weak form of (28), whence

v 1 = u 1 since v 1 (0) = u 1 (0). Consequently, we have T α-1 (u 1 (t), u 2 (t)) ≤ (α -1) -1
|x -y| α-1 v(x, y, t)dxdy, and, using the same arguments as usual, it suffices to show that

d dt |x -y| α-1 v(x, y, t)dxdydt ≤ 0.
This follows from the fact that for all x, y ∈ R, setting u = σ(x)-σ(y)

x-y

, R |x + σ(x)h -y -σ(y)h| α-1 -|x -y| α-1 -(α -1)h[σ(x) -σ(y)]|x -y| α-3 (x -y) dh |h| 1+α = |x -y| α-1 R |1 + hu| α-1 -1 -(α -1)hu dh |h| 1+α = |x -y| α-1 |u| α R |1 + h| α-1 -1 -(α -1)h dh |h| 1+α = 0.
The proof of this last equality can be found in [22, Lemma 9-(ii)], case a + = a -and

β = α -1. Observe that |x -y| α-1 |u| α = |σ(x) -σ(y)| α |x -y| ≤ C σ
so that (29) makes sense with ϕ(x, y) = |x -y| α , thanks to our regularity assumption on σ.

Here again, as in Section 2.1, the main technical difficulty is to prove the well-posedness of [START_REF] Komatsu | On the pathwise uniqueness of solutions of one-dimensional stochastic differential equations of jump type[END_REF], in particular when σ may degenerate. This is useful to check that the solution v to the coupled equation has the correct marginals.

Inhomogeneous integral equations

Our next purpose is to give an example on the way to take into account x-dependency in IPDE models, for instance when considering a measure µ(x, h) in the scattering equation [START_REF] Figalli | Existence and uniqueness of martingale solutions for SDEs with rough or degenerate coefficients[END_REF]. We exemplify this issue with a simple equation we borrow from [START_REF] Fournier | On a toy model of interacting neurons[END_REF]. Consider an interval I of R, a rate function d ≥ 0 defined on I and some probability density b on I. We consider the conservative equation

∂ t u(x, t) + d(x)u = b(x)A(t), A(t) = I d(x)u(x, t)dx (30)
starting from an initial probability density u 0 on I. We notice at once that this equation makes sense for probability measures u(dx, t) (for each t ≥ 0, u(dx, t) is a probability measure on I) in the following weak sense: for all smooth ϕ : I → R,

d dt ϕ(x)u(dx, t) = [ϕ(z) -ϕ(x)]b(z)d(x)u(dx, t)dz. (31)
Theorem 12 Consider two probability densities u 0 1 , u 0 2 on R d and the corresponding solutions u 1 , u 2 to [START_REF] Marchioro | Mathematical theory of incompressible nonviscous fluids[END_REF]. Under one of the two conditions (a) or (b) below, for all t ≥ 0,

T (u 1 (t), u 2 (t)) ≤ T (u 0 1 , u 0 2 ).

Homogeneous Boltzmann equation

In his seminal paper [START_REF] Tanaka | Probabilistic treatment of the Boltzmann equation of Maxwellian molecules[END_REF], Tanaka observed, using a probabilistic approach based on nonlinear Poissondriven stochastic differential equations, that the homogeneous Boltzmann equation for Maxwell molecules is non-expansive for the 2-transport cost. A deterministic proof was provided by Villani in [41, Section 7.5.6] and was extended to inelastic collisions in [START_REF] Bolley | Tanaka theorem for inelastic Maxwell models[END_REF]. A survey of results concerning homogeneous kinetic equations can be found in [START_REF]Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF]. The approach in [START_REF]Topics in optimal transportation[END_REF][START_REF] Bolley | Tanaka theorem for inelastic Maxwell models[END_REF][START_REF]Contractive probability metrics and asymptotic behavior of dissipative kinetic equations[END_REF] is to first derive a contractivity result for the gain operator by coupling, and then to conclude via Duhamel's principle.

Here we follow Tanaka's original approach, or rather we show how to write down his main arguments without introducing stochastic processes.

Let us also mention that [START_REF] Rousset | A N -uniform quantitative Tanaka's theorem for the conservative Kac's N -particle system with Maxwell molecules[END_REF] managed to study the corresponding dissipation in order to quantify the convergence to equilibrium of the solutions and, even more interesting, to prove the convergence to equilibrium of the Kac particle system, with a rate of convergence not depending on the number of particles.

The homogeneous Boltzmann writes

           ∂ t f (v, t) = Q(f ) := R 3 S 2 [f (v , t)f (v * , t) -f (v, t)f (v * , t)]B(θ)dv * dσ, v = 1 2 (v + v * ) + 1 2 |v -v * |σ, v * = 1 2 (v + v * ) -1 2 |v -v * |σ, cos(θ) = v-v * |v-v * | • v -v * |v-v * | . (34)
The collision kernel B : (0, π) → R + is assumed to satisfy π 0 B(θ)dθ = 1. As is well-known, this equation writes, in weak form, for all mapping ϕ : R 3 → R,

d dt R 3 φ(v)f (v, t)dv = R 3 ×R 3 S 2 φ(v ) + φ(v * ) -φ(v) -φ(v * ) B(θ)f (v, t)f (v * , t)dvdv * dσ. ( 35
)
Theorem 13 Consider two initial probability densities f 0 1 , f 0 2 on R 3 with a finite moment of order 2 and the corresponding solutions f 1 , f 2 to (34). Then, for all t ≥ 0, one has

T 2 (f 1 (t), f 2 (t)) ≤ T 2 (f 0 1 , f 0 2 ).
Proof. We fix a probability density F 0 on (R 3 ) 2 with marginals f 0 1 and f 0 2 and build a coupled equation with the same principle as for scattering, that is the jump parameters are taken in common to the two variables, in such a way that the post-collisional velocities are as close as possible. We consider the solution F (v, w, t), starting from F 0 , to the following coupling equation written in weak form: for all mapping Ψ : R

3 × R 3 → R, d dt (R 3 ) 2 Ψ(v, w)F (v, w, t)dvdw = (R 3 ×R 3 ) 2 π 0 2π 0 Ψ(v , w ) + Ψ(v * , w * ) -Ψ(v * , w * ) -Ψ(v, w) B(θ)F (v, w, t)F (v * , w * , t)dvdwdv * dw * dθdϕ, (36)
where, for v, w, v * , w * ∈ R 3 , θ ∈ (0, π) and ϕ ∈ (0, 2π), we have set

σ = cos(θ) v -v * |v -v * | +sin (
v = 1 2 (v + v * ) + 1 2 |v -v * |σ, v * = 1 2 (v + v * ) - 1 2 |v -v * |σ, w = 1 2 (w + w * ) + 1 2 |w -w * |ω, w * = 1 2 (w + w * ) - 1 2 |w -w * |ω. It is not difficult to show that (36) is well-posed in L ∞ loc ([0, ∞), L 1 (R 3 × R 3
)) and that the solution remains nonnegative for all times. Indeed, applying [START_REF] Pierre | Global existence in reaction-diffusion systems with control of mass: a survey[END_REF] with the test function Ψ(v, w) = -1 I {F (v,w)<0} , using that

-1 I {F (v ,w )<0} -1 I {F (v * ,w * )<0} + 1 I {F (v,w)<0} + 1 I {F (v * ,w * )<0} F (v, w)F (v * , w * ) ≤4 F -(v, w)|F (v * , w * )| + |F (v, w)|F -(v * , w * ) and that π 0 B(θ)dθ = 1, we conclude that d dt (R 3 ) 2 F -(v, w, t)dvdw ≤ 8 (R 3 ) 2 F -(v, w, t)dvdw × (R 3 ) 2 |F (v, w, t)|dvdw,
whence the conclusion by the Gronwall lemma.

Also, it holds that

R 3 F (v, w, t)dw = f 1 (v, t) and R 3 F (v, w, t)dv = f 2 (w, t).
For example concerning f 1 , we apply the weak coupling equation to some Ψ depending only on v and we show that R 3 F (v, w, t)dw solves [START_REF] Otto | Eulerian calculus for the contraction in the Wasserstein distance[END_REF]. This follows from the fact that, when fixing (v, w) and (v * , w * ), the expression between brackets in (36) only depends on σ, so that for any function H : S 2 → R, we may write

S 2 H(σ)B(θ)dσ = π 0 2π 0 H cos(θ) v -v * |v -v * | + sin(θ)[I cos(ϕ) + I 1 sin(ϕ)] B(θ)dϕdθ.
We conclude that

(R 3 ×R 3 ) 2 π 0 2π 0 Ψ(v ) + Ψ(v * ) B(θ)F (v, w, t)F (v * , w * , t)dvdwdv * dw * dθdϕ = (R 3 ×R 3 ) 2 S 2 Ψ(v ) + Ψ(v * ) B(θ)F (v, w, t)F (v * , w * , t)dvdwdv * dw * dσ = (R 3 ×R 3 ) 2 S 2 Ψ(v ) + Ψ(v * ) B(θ)f 1 (v, t)f 1 (v * , t)dvdv * dσ.
Consequently, it holds that T 2 (f 1 (t), f 2 (t)) ≤ |v -w| 2 F (v, w, t)dvdw =: h(t), and it suffices, as usual, to show that h (t) ≤ 0. For this, it suffices to verify that for all fixed v, w, v * , w * ∈ R 3 , all θ ∈ (0, π),

∆ = 2π 0 [|v -w | 2 + |v * -w * | 2 -|v -w| 2 -|v * -w * | 2 ]dϕ ≤ 0. A simple computation, using that 2π 0 σdϕ = 2π 0 ωdϕ = 0, shows that ∆ = 2π 0 [(v -w) • (v * -w * ) -|v -v * ||w -w * |σ • ω]dϕ = 2π 0 [(1 -cos 2 θ -sin 2 θ sin 2 ϕ)(v -w) • (v * -w * ) -sin 2 θ cos 2 ϕ|v -v * ||w -w * |]dϕ. We used that |v -v * ||w -w * |I 1 • I 2 = (v -w) • (v * -w * ). All in all, we arrive at ∆ = [(v -w) • (v * -w * ) -|v -v * ||w -w * |] sin 2 θ 2π 0 cos 2 ϕdϕ,
and the proof is complete.

Porous media equation

We now consider the generalized porous media equation written, with A : R + → R of class C 2 , as

∂ t u -div(u∇[A (u)]) = 0, x ∈ R d , t ≥ 0. ( 37
)
Using a gradient flow approach taking advantage of displacement convexity, as introduced in [START_REF] Mccann | A convexity principle for interacting gases[END_REF], it was discovered by Otto [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF], see also [START_REF] Carrillo | Contractions in the 2-Wasserstein length space and thermalization of granular media[END_REF][START_REF] Otto | Eulerian calculus for the contraction in the Wasserstein distance[END_REF][START_REF] Daneri | Eulerian calculus for the displacement convexity in the Wasserstein distance[END_REF], that this equation is non-expansive for T 2 , under a few conditions on A, including convexity. The method was improved by Bolley and Carrillo [START_REF] Bolley | Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction[END_REF] who make the contraction property an equivalence to displacement convexity for the energy functional. The coupling method, as defined in the introduction and used in the whole present paper, does not seem to apply directly. However, using Brenier's map, the argument in [START_REF] Bolley | Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction[END_REF] can somehow be presented in relation with a coupling. We present this argument, staying at an informal level.

Theorem 14 ([34, 6]) Consider some C 2 function A : R + → R such that B(r) = r 0 wA (w)dw ≥ 0 for all r ≥ 0 and such that r → r 1/d-1 B(r) is non-decreasing. Consider two probability densities u 0 1 , u 0 2 on R d and the corresponding solutions u 1 , u 2 to (37). Then for all t ≥ 0, one has

T 2 (u 1 (t), u 2 (t)) ≤ T 2 (u 0 1 , u 0 2 ).
This applies to the porous media equation, i.e. with A(u) = m -1 u m , as soon as m ≥ 1. The justification of the computation requires at least that R d B(u 1 (x, t)) < ∞. See [START_REF] Otto | The geometry of dissipative evolution equations: the porous medium equation[END_REF][START_REF] Bolley | Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction[END_REF] for the rigorous proof, which uses that general weak solutions of the porous media equation can be approximated by smooth and positive solutions of suitably regularized nonlinear diffusion equations.

Proof. We consider Brenier's map [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF] for u 0 1 and u 0 2 , i.e. a convex function Φ :

R d → R such that T 2 (u 0 1 , u 0 2 ) = 1 2 R d |x -∇Φ(x)| 2 u 0 1 (x)dx and ∇Φ#u 0 1 = u 0 2 .
We next consider the probability measure v(dx, dy, t) (for each t ≥ 0, v(t) ∈ P((R d ) 2 )) solving the coupling equation ∂v ∂t = div x v∇ x A (u 1 (x, t) + div y v∇ y A (u 2 (y, t)

and starting from v 0 ∈ K(u 0 1 , u 0 2 ) defined by the formula v 0 (A) = R d 1 {(x,∇Φ(x))∈A} u 0 1 (x)dx for all Borel set A ⊂ R d , in short notation v 0 (x, y) = u 0 1 (x)δ ∇Φ(x) (y). Because of this specific initial data, the equation on v will not provide a global coupling as in the equations treated so far. Again, we only use the weak form, for all smooth ϕ : R d × R d → R, d dt ϕ(x, y)v(dx, dy, t) = -∇ x ϕ(x, y)•∇ x [A (u 1 (x, t))]+∇ y ϕ(x, y)•∇ y [A (u 2 (y, t))] v(dx, dy, t).

One easily verifies, as usual, that for each t ≥ 0, v(t) has u 1 (t) and u 2 (t) for marginals: for example, applying the weak equation of v to some ϕ depending only on x shows that y∈R d v(dx, dy, t) is a Since r → r 1/d-1 B(r) is non-decreasing by assumption, we conclude that I (0) ≤ 0.

The above considerations hold true at any time, and not only at t = 0. In other words, for all t ≥ 0, we can find a function I t : [t, ∞) → R such that I t (t) = T 2 (u 1 (t), u 2 (t)), I t (t) ≤ 0 and T 2 (u 1 (s), u 2 (s)) ≤ I t (s) for all 0 ≤ t ≤ s.One immediately concludes that for all t ≥ 0, lim sup h↓0 T 2 (u 1 (t + h), u 2 (t + h)) -T 2 (u 1 (t), u 2 (t)) h ≤ I t (t) ≤ 0, so that t → T 2 (u 1 (t), u 2 (t)) is non-increasing.

An approach by duality

In order to complete the presentation, we quickly mention another possible and original approach, based on duality. We consider the simplest model, i.e. the heat equation in dimension 1, but all the models treated in the present paper, except the porous media equation, may be treated similarly, with more complicated discretization procedures and more involved computations. See also Villani [40, pages 41-43] for a similar approach concerning the Vlasov equation, in the case of the Monge-Kantorovich distance T 1 , of which the dual expression is particularly simple.

Proof of Theorem 1 when d = 1 for ρ(x, y) = |x -y| p with p ≥ 1. We consider two solutions u 1 , u 2 to (2), starting from probability measures u 0 1 , u 0 2 with finite p-moment. For h > 0, we consider the solutions u 1,h , u 2,h , starting from u 0 1 , u 0 2 , to the discrete heat equation It is standard that u 1,h → u 1 and u 2,h → u 2 as h → 0 (in the weak topology of measures for instance). We will verify that for each h > 0, it holds that T p (u 1,h (t), u 2,h (t)) ≤ T p (u 0 1 , u 0 2 ), for any p ≥ 1, and this will complete the proof.

∂
We fix p ≥ 1 and introduce the set Q p of pairs (ϕ, ψ) of functions from R d to R such that for all x, y ∈ R d , ϕ(x) + ψ(y) ≤ |x -y| p . For any pair of probability densities f, g on R d , the transport cost can also be expressed by duality, see [START_REF]Grundlehren der Mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences[END_REF], as T p (f, g) = 1 p sup (ϕ,ψ)∈Qp ϕ(x)f (x)dx + ψ(y)g(y)dy .

For (ϕ, ψ) ∈ Q p , we set ∆ ϕ,ψ (t) = ϕ(x)u Taking the supremum over all pairs (ϕ, ψ) in Q p and dividing by p, we conclude that e 2h -2 t T p (u 1,h (t), u 2,h (t)) ≤ T p (u 0 1 , u 0 2 ) + 2h -2 t 0 e 2h -2 s T p (u 1,h (s), u 2,h (s))ds.

By the Gronwall lemma, we conclude that e 2h -2 t T p (u 1,h (t), u 2,h (t)) ≤ e 2h -2 t T p (u 0 1 , u 0 2 ) as desired. Unfortunately, we are not able to use a similar procedure directly on the (non discretized) heat equation.

  (a) I = R + , d(0) = 0, d is increasing, b = δ 0 , and (x, y) = |d p (x) -d p (y)| for some p ≥ 1. (b) I = R + , d(x) = αx p + β for some α, β ≥ 0 and p ≥ 1, with (x, y) = |x p -y p |, under the condition that β ≥ α ∞ 0 z p b(z)dz.

  θ)[I cos(ϕ)+I 1 sin(ϕ)] and ω = cos(θ) w -w * |w -w * | +sin(θ)[I cos(ϕ)+I 2 sin(ϕ)], where I = (v-v * )∧(w-w * ) |(v-v * )∧(w-w * )| , and I 1 , I 2 are chosen so that ( v-v * |v-v * | , I, I 1 ) and ( w-w * |w-w * | , I, I 2 ) are two direct orthonormal bases, and where

  1,h (x, t)dx+ ψ(y)u 2,h (y, t)dy. Using that (ϕ(•+h), ψ(•+h)) and (ϕ(• -h), ψ(• -h)) both belong to Q p , we find d dt ∆ ϕ,ψ (t) ≤ -2h -2 ∆ ϕ,ψ (t) + 2ph -2 T p (u 1,h (t), u 2,h (t)).

	This implies that

e 2h -2 t ∆ ϕ,ψ (t) ≤ ∆ ϕ,ψ (0) + 2ph -2 t 0 e 2h -2 s T p (u 1,h (s), u 2,h (s))ds ≤ pT p (u 0 1 , u 0 2 ) + 2ph -2 t 0 e 2h -2 s T p (u 1,h

(

s), u 2,h (s))ds.

B.P. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 740623).

Other assumptions on I, b, d are possible: it suffices that , b and d satisfy the dual inequality [START_REF] Mccann | A convexity principle for interacting gases[END_REF] below, which corresponds to [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF] for the heat equation with variable coefficients.

Proof. We consider some probability density v 0 on I 2 with marginals u 0 1 and u 0 2 and define the probability measure v(dx, dy, t) as solving, for all smooth ϕ :

It holds true that v(t) is a probability measure on I 2 for each t ≥ 0, and that its marginals are u 1 (t) and u 2 (t). For example, applying the coupling equation with ϕ depending only on x and using that min(d(x), d(y))

one verifies that y∈I v(dx, dy, t) solves [START_REF] Marchioro | Mathematical theory of incompressible nonviscous fluids[END_REF], whence y∈I v(dx, dy, t) = u 1 (dx, t) by uniqueness. Hence for any cost function :

. Furthermore, we easily compute, using that (z, z) = 0 for all z ∈ I, that b is a probability density, and that min(r, s)

Therefore, using the same arguments as usual, the result will follows from the fact that for all x, y ∈ I,

(a) Assuming that I = R + , that d(0) = 0, that d is increasing, that ρ(x, y) = |d(x) -d(y)| p for some p ≥ 1 and that b = δ 0 , we check that, when e.g. x ≥ y ≥ 0,

which holds true since indeed, for any s ≥ t ≥ 0, (s p -t p )s ≥ t p (s -t) because p ≥ 1.

(b) Assume next that I = R + , d(x) = αx p + β for some α, β ≥ 0 and p ≥ 1 and that (x, y) = |x p -y p |. We have to verify that, for all x ≥ y ≥ 0,

Observe that the strong equation corresponding to the weak form (32) is nothing but

(weak) solution to [START_REF] Rousset | A N -uniform quantitative Tanaka's theorem for the conservative Kac's N -particle system with Maxwell molecules[END_REF] and since it starts from u 0 1 , we conclude by uniqueness. As a conclusion, for all t ≥ 0, T 2 (u 1 (t), u 2 (t)) ≤ I(t), where I(t) = 1 2 |x -y| 2 v(dx, dy, t). Next, following [START_REF] Bolley | Nonlinear diffusion: geodesic convexity is equivalent to Wasserstein contraction[END_REF], we observe that

where

In particular, by definition of v 0 ,

where we recall that B(r) = r 0 wA (w)dw. Integrating by parts, we thus find

This uses that for any convex function Φ :

, we have, for any ϕ : R d → R, ϕ(x)u 0 1 (x)dx = ϕ((∇Φ) -1 (y))u 0 2 (y)dy = ϕ(x)u 2 (∇Φ(x)) det D 2 Φ(x)dx, so that det D 2 Φ(x) = u 0 1 (x)/u 2 (∇Φ(x)) (see [START_REF] Gutiérrez | Progress in Nonlinear Differential Equations and their Applications[END_REF] for an account on this Monge-Ampère equation). All in all, we have checked that

Proceeding similarly, we see that

Performing the substitution x = (∇Φ) -1 (y), we end with

dx.

We thus find, with the notation y = ∇Φ(x)