Simone Di Marino
email: simone.dimarino@math.u-psud.fr

Augusto Gerolin
email: gerolin@mail.dm.unipi.it

Luca Nenna
email: luca.nenna@inria.fr

Optimal transportation

published or not. The documents may come

Introduction

The main objective of this survey is to present recent developments in optimal transport theory for repulsive costs with finitely many marginals. We are interested in characterizing, for a class of repulsive cost functions, the minimizers of the multi-marginal optimal transportation problem (MK) min R dN c(x 1 , . . . , x N )dγ : γ ∈ Π(R dN , µ 1 , . . . , µ N ) , where c : (R d ) N -→ R ∪ {∞} is a lower semi-continuous cost function and Π(R dN , µ 1 , . . . , µ N ) denotes the set of couplings, i.e. probability measures γ(x 1 , . . . , x N ) having µ i as marginals, for all i = 1, . . . , N . c is, typically, a sum of convex decreasing functions, as the Coulomb interaction cost [START_REF] Buttazzo | Optimal-transport formulation of electronic density-functional theory[END_REF][START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF][START_REF] Elliott | Density functionals for coulomb systems[END_REF][START_REF] Gori-Giorgi | Density-functional theory for strongly interacting electrons[END_REF][START_REF] Elliott | Density functionals for coulomb systems[END_REF][START_REF] Seidl | Strong-interaction limit of density-functional theory[END_REF]

c(x 1 , . . . , x N ) = N i=1 N j=i+1 1 |x j -x i | , 1
which is the case of most interest in DFT (see section 1). Also related to Quantum Mechanics, in the so-called weak interaction regime [START_REF] Buttazzo | Optimal-transport formulation of electronic density-functional theory[END_REF][START_REF] Friesecke | N-density representability and the optimal transport limit of the hohenberg-kohn functional[END_REF][START_REF] Paziani | Localspin-density functional for multideterminant density functional theory[END_REF][START_REF] Toulouse | Short-range exchange and correlation energy density functionals: Beyond the local-density approximation[END_REF], we will consider the repulsive harmonic potential as cost function c, cw(x 1 , . . . , x N ) = -

N i=1 N j=i+1 |x j -x i | 2 .
In the last section, we focus on the multi-marginal problem for the determinant cost c(x 1 , . . . , x N ) = det(x 1 , . . . , x N ), or c(x 1 , . . . , x N ) = | det(x 1 , . . . , x N )|, studied by Carlier and Nazaret in [START_REF] Carlier | Optimal transportation for the determinant[END_REF].

Multi-marginal optimal transport showed to be powerful in many applications as, for instance, in economics [START_REF] Carlier | Matching for teams[END_REF][START_REF] Chiappori | The roommate problem is more stable than you think[END_REF][START_REF] Chiappori | Hedonic price equilibria, stable matching, and optimal transport: equivalence, topology, and uniqueness[END_REF][START_REF] Ekeland | An optimal matching problem[END_REF][START_REF] Pass | Multi-marginal optimal transport and multi-agent matching problems: uniqueness and structure of solutions[END_REF], mathematical finance [START_REF] Beiglböck | Model-independent bounds for option prices-a mass transport approach[END_REF][START_REF] Beiglboeck | On a problem of optimal transport under marginal martingale constraints[END_REF][START_REF] Beiglboeck | Optimal transport and skorokhod embedding[END_REF][START_REF] Dolinsky | Martingale transport and robust hedging in continuous time[END_REF][START_REF] Dolinsky | Robust hedging with proportional transaction costs[END_REF][START_REF] Galichon | A stochastic control approach to no-arbitrage bounds given marginals, with an application to lookback options[END_REF], image processing [START_REF] Xia | Synthesizing and mixing stationary gaussian texture models[END_REF], tomography [START_REF] Abraham | Tomographic reconstruction from a few views: a multi-marginal optimal transport approach[END_REF] and statistics [START_REF] Carlier | Vector quantile regression[END_REF].

Although very natural in such applied problems, optimal transport problems with 3 or more marginals was only considered in literature during the 90's by Knott and Smith, Olkin and Rachev and Rüschendorf and Uckelmann [START_REF] Knott | On a generalization of cyclic monotonicity and distances among random vectors[END_REF][START_REF] Olkin | Maximum submatrix traces for positive definite matrices[END_REF][START_REF] Rüschendorf | On optimal multivariate couplings[END_REF], and later by Gangbo and Świech [START_REF] Gangbo | Optimal maps for the multidimensional mongekantorovich problem[END_REF]. In [START_REF] Gangbo | Optimal maps for the multidimensional mongekantorovich problem[END_REF], the authors proved the existence and uniqueness of Monge-type solution for the problem (MK) for the quadratic cost c(x 1 , . . . , x N ) = N i,j=1 |x j -x i | 2 . A systematic study of the structure of transport plans and general conditions in order to guarantee the existence of Monge-type solutions for the problem (MK) was made by Pass [START_REF] Pass | Structural results on optimal transportation plans[END_REF][START_REF] Pass | Uniqueness and monge solutions in the multimarginal optimal transportation problem[END_REF][START_REF] Pass | On the local structure of optimal measures in the multi-marginal optimal transportation problem[END_REF] and Kim and Pass [START_REF] Kim | Multi-marginal optimal transport on riemannian manifolds[END_REF][START_REF] Kim | A general condition for monge solutions in the multimarginal optimal transport problem[END_REF]. Several theoretical developments were made last years, as the study of partial multi-marginal transport [START_REF] Kitagawa | The multi-marginal optimal partial transport problem[END_REF], symmetric Monge-Kantorovich problem [START_REF] Colombo | Equality between monge and kantorovich multimarginal problems with coulomb cost[END_REF][START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF][START_REF] Ghoussoub | Remarks on multi-marginal symmetric mongekantorovich problems[END_REF][START_REF] Ghoussoub | Symmetric monge-kantorovich problems and polar decompositions of vector fields[END_REF][START_REF] Ghoussoub | Decoupling of degiorgi-type systems via multi-marginal optimal transport[END_REF][START_REF] Alexander | Optimal transportation of processes with infinite kantorovich distance. independence and symmetry[END_REF][START_REF] Lopes | Duality theorems in ergodic transport[END_REF][START_REF] Lopes | The dual potential, the involution kernel and transport in ergodic optimization[END_REF][START_REF] Zaev | On the monge-kantorovich problem with additional linear constraints[END_REF] and duality [START_REF] Beiglböck | A general duality theorem for the monge-kantorovich transport problem[END_REF][START_REF] De | Optimal transport with coulomb cost. approximation and duality[END_REF][START_REF] Hans | Duality theorems for marginal problems[END_REF][START_REF] Zaev | On the monge-kantorovich problem with additional linear constraints[END_REF]. Applications in Geometry and Analysis as the Wasserstein Barycenter problem and decoupling a system of PDE's can be found in [START_REF] Agueh | Barycenters in the {W}asserstein space[END_REF][START_REF] Ghoussoub | Decoupling of degiorgi-type systems via multi-marginal optimal transport[END_REF][START_REF] Kim | Wasserstein barycenters over riemannian manfolds[END_REF].

For a complete presentation of the multimarginal optimal transport theory with general costs functions, including a discussion of many applications, we refer to the survey [START_REF] Pass | Multi-marginal optimal transport: theory and applications[END_REF].

In this paper, we are rather analyzing the Multimarginal Optimal Transport Theory starting from the point of view of researchers working on Density Functional Theory.

Paola Gori-Giorgi, Michel Seidl et al were interested in computing accurately the ground state of a molecule with a large number of atoms. As we will discuss later, after doing suitable approximations, the solutions of the optimal transport (MK) for Coulomb cost (more precisely, the optimal potential), or more generally repulsive costs, can be seen as a good ansatz in order to construct solutions for Hohenberg-Kohn functional (1.12).

From a mathematical perspective the interesting case is when all marginals of the coupling γ in (MK) are absolutely continuous with respect to Lebesgue measure and are all the same. In particular, since the cost functions are "repulsive", if Monge-type solutions exists, they should follow the rule "the further, the better!", which means that we want to move the mass as much as we can. In other words, in the present case, optimal transport plans tend to be as spread as possible. But when a cost can be considered repulsive? It is difficult to find a good notion for this concept since usually it is based on the intuition that comes solving the particular problem. However a possible formal definition for repulsive cost is to require a mild condition, namely that we don't want any two particles in the same position; it should be cheaper to have them in different position. For example can write the inequality for the particles 1 and 2: c(x 1 , x 1 , x 3 , . . . , x N ) + c(x 2 , x 2 , x 3 , . . . , x N ) ≥ c(x 1 , x 2 , x 3 , . . . , x N ) + c(x 2 , x 1 , x 3 , . . . , x N ).

Then, also all other inequalities for i = j should be considered. This condition is satisfied by all the considered costs, even the ones that doesn't seem so "repulsive" at first, such as those of the form f (x 1 + . . . + xn) with f convex function.

Organization of this Survey. In section 1, we present the physical motivation in considering repulsive costs and give an introduction of Density Functional Theory from a mathematical point of view. The main goal of this section is to explain how optimal solutions of Monge-Kantorovich problem can be useful to compute the ground state energy of a quantum system of N electrons.

Section 2 introduces briefly the general theory of multimarginal optimal transport, including a discussion about the Kantorovich duality and the notion of c-cyclical monotonicity applied to this framework. In addition, we give some highlights in the interesting case of symmetric transport plans.

In sections 3, 4 and 5, respectively, we discuss the recent development for the specific cases of transportation problem with Coulomb, repulsive harmonic and determinant cost.

Finally, in section 6 we introduce an iterative method, recently used in many applications of optimal transport [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF][START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF], and we present some numerical experiments for the problems studied in previous sections.

What is new? We present in section 4 some new results we obtained for cost functions c(x 1 , . . . , x N ) = f (x 1 + • • • + x N ), f : R d → R is convex, which, in particular, include the case of repulsive harmonic costs (see corollary 4.8). In addition, we show:

(i) Existence of Fractal-like cyclical Monge solutions: In Theorem 4.6, we show the existence of couplings γ = (Id, T, T (2) , . . . , T (N -1) ) µ concentrated on the graphs of functions T, T (2) , . . . , T (N -1) , where µ a uniform probability measure on the d-dimensional cube. Moreover, T is not differentiable in any point, which differs substantially to the classical regularity theory in the classical 2 marginals case [START_REF] Luis | A localization property of viscosity solutions to the monge-ampere equation and their strict convexity[END_REF][START_REF] Luis | Some regularity properties of solutions of monge ampere equation[END_REF].

(ii) Optimal Transport plans supported on sets of Hausdorff dimension > 1: In standard optimal transportation (the 2 marginals case) the target is to understand if solutions of (MK) are well-known to be Monge-type solutions: optimal transport plans are supported by the graph of a function T . In Section 4, we give a concrete example (see example 4.14) where optimal transport plans are supported on 2-dimensional Hausdorff sets.

(iii) Numerical results: we present some interesting results for the repulsive harmonic and for the determinant cost.

1. Why Multi-marginal Transport Theory for repulsive costs?

1.1. Brief introduction to Quantum Mechanics of N -body systems.

Energies and geometries of a molecule depend sensibly on the kind of atom on a chemical environment. Semiempirical models as Lewis structure explain some aspects, but they are far from being satisfactory, mostly because they are not quantitative.

Classical behavior of atoms and molecules is described accurately, at least from a theoretical point of view, by quantum mechanics. However, in order to predict the chemical behavior of a molecule with a large number of electrons we need to deal with computational aspects and approximations are needed. For example, if we use a direct approximation of 10 grid points in R of the time-dependent Schrödinger equation in order to simulate the chemical behavior of the water molecule (H 2 O), which has 10 electrons and its position is represented in the space R 30 , then we need to compute 10 30 grid points.

Quantum Chemistry studies the ground state of individual atoms and molecules, the excited states and the transition states that occur during chemical reactions. Many systems use the so-called Born-Oppenheimer approximation and many computations involve iterative and other approximation methods.

The main goal of Quantum Chemistry is to describe accurately results for small molecular systems and to increase the size of large molecules that can be processed, which is very limited by scaling considerations.

It means that from a Quantum Model, for which the only input is the atomic number of a molecule, we want to predict the evolution of a molecular system. Classical models consider a molecule with N electrons and M nuclei. We denote by a point x i ∈ R 3 the position coordinates of the i-th electron of mass me and by s i its i-th electron spin. The charges, the masses and the positions of the α-th nucleus are represented, respectively, by

Zα ∈ R M , mα ∈ R M + and Rα ∈ (R 3 ) M
The state of the system is described by a time-dependent wave function

ψ ∈ L 2 ([0, T ]; ⊗ N i=1 L 2 (R 3 × Z 2 ) × ⊗ M α=1 L 2 (R 3 )), ψ(x 1 , s 1 , . . . , x N , s N ; R 1 , . . . , R M , t).
The group Z 2 = {↑, ↓} represents the spin of a particle.

We say that a wave function is antisymmetric or fermionic if it changes sign under a simultaneous exchange of two electrons i and j with the space coordinates x i and x j and spins s i and s j , that is,

ψ(x σ(1) , s σ(1) , . . . , x σ(N ) , s σ(N ) ) = sign(σ)ψ(x 1 , s 1 , . . . , x N , s N ) σ ∈ S N
where S N denotes the permutation set of N elements.

In physics literature, the Fermionic wave function obeys the Fermi-Dirac statistics and, in particular, the Pauli exclusion principle. On the other hand, a wave function is called symmetric or bosonic when it obeys Bose-Einstein statistics which also implies that when one swaps two bosons, the wave function of the system is unchanged.

If it is not mentioned explicitly, we will suppose that the wave function is spinless, i.e. ψ is a function depending only on time, electrons and nuclei position ψ(t, x 1 , x 2 , . . . , x N , R 1 , . . . , R M ).

The evolution of a wave function ψ boils down to a (time-dependent) Schrödinger Equation

i∂tψ = Hψ , H = Tn + Te + Vne + Vee + Vnn (1.1)
where the operators Tn, Te, Vne, Vee and Vnn are defined by

Tn = - M α=1 1 2mα
∆ Rα (Nuclei Kinetic Energy)

Te = - N i=1 1 2me ∆x i (Electron Kinetic Energy) Vne = - N i=1 v R (x i ) = - N i=1 M α=1 Zα |x i -Rα| Potential Energy of Inte- raction Nuclei-Electron Vee = N i=1 N j=1,j =i 1 |x i -x j | (Interaction Electron-Electron) Vnn = M α=1 M β=1,β =α ZαZ β |Rα -R β | (Interaction Nuclei-Nuclei)
Born-Oppenheimer Approximation: The Born-Oppenheimer approximation (named for its original inventors, Max Born and Robert Oppenheimer) is based on the fact that nuclei are several thousand times heavier than electrons. The proton, itself, is approximately 2000 times more massive than an electron. Roughly speaking, in Born-Oppenheimer approximation we suppose that nuclei behave as point classical particles. This is reasonable, since typically 2000me ≤ mα ≤ 100000me.

In other words, we suppose mα me, we fix me = 1 and we consider an ansatz of type ψ(x, R, t) = ψe(x, R)χ(R, t), which physically means that the dynamics of the electron ψe are decoupled of the dynamics of the atomic nuclei χ. Substituting that ansatz in (1.1), we can show that the electronic part ψe(x, R) solves the following eigenvalue problem, so-called Eletronic Schrödinger Equation

Heψe = λeψe, He = Te + Vne + Vee + Vnn (1.2)
where the eigenvalue λe = λe(R) depends on the position vector R = (R 1 , . . . , R M ) of the atomic nuclei. The function He = He(R) is called Eletronic Hamiltonian. From the other side, through a formal argument, the nuclei wave function χ(R, t) is a solution of a Schrödinger Equation restricted to the nuclei with potential energy λe,

i∂tχ(R) = (Tn + λe(R))χ(R).
We refer to [START_REF] Bodo | Applicazioni di meccanica quantistica: appunti per le lezioni[END_REF] for all computations and formal deduction of those formulas. Another consequence of that approximation is that the nuclear components of the wave function are spatially more localized than the electronic component of the wave function. In the classical limit, the nuclei are fully localized at single points representing classical point particles.

In Quantum Chemistry it is interesting to study the so called Geometric Optimization Problem.

The Geometric Optimization

Problem: Compute the following minimizer inf E(R 1 , . . . , R M ) + M α=1 M β=1,β =α ZαZ β |Rα -R β | : (R 1 , . . . , R M ) ∈ R 3M (1.3)
where

• The second term is the nuclei-nuclei iteration Vnn, as defined in the (1.1).

• The function E(R) = E(R 1 , . . . , R M ) corresponds to the effective potential created by the electrons and is, itself, given by a minimization problem (see Electronic minimization problem bellow).

The value of the potential E 0 = E 0 (R 1 , . . . , R M ) for given fixed positions of the nuclei is obtained by solving the electronic problem:

The Electronic minimization problem: Compute the lowest eingenvalue ("ground state energy") E 0 of the following linear operator, called Electronic Hamiltonian Definition 1.1 (Eletronic Hamiltonian). The Eletronic Hamiltonian is a linear operator H el :

L 2 anti (R 3N ) → L 2 anti (R 3N ), H el = - 2 2 N i=1 ∆x i + N i=1 v R (x i ) + N i=1 N j=i+1 f (x j -x i ),
where L 2 anti (R 3N ) denote the set of square-integrable antisymmetric functions ψ :

R 3N → R, v R : R 3 → R is a L 2 (R 3N ) function and f : R → R is a continuous function.
The first term of the Eletronic Hamiltonian is the Kinetic Energy and the second term, the function v R (x i ), depends only on the single electron. Typically, v R (x i ) = -M α=1 Zα |x i -Rα| is the interaction electron-nuclei energy of the single electron. The function f depends only on the distance of two electrons and measures the electron-electron potential interaction; typically

f is the Coulomb Interaction f (x j -x i ) = 1/|x j -x i | or the repulsive harmonic interaction f (x j -x i ) = -|x j -x i | 2 .
From both mathematical and physical viewpoint, it may be also useful to study more general convex and concave functions.

The ground state of the functional (1.1) is given by the Rayleigh-Ritz variational principle

E 0 = E(R 1 , • • • , R M ) = min{ ψ, H el ψ : ψ ∈ H},
where H = {ψ ∈ H 1 (R 3N ) : ψ is antisymmetric and ψ = 1}. Equivalently,

E 0 = min Te[ψ] + Vne[ψ] + Vee[ψ] : ψ ∈ H (1.4)
where Te[ψ] is the Kinetic energy,

Te[ψ] = 2 2 s 1 ∈Z 2 R 3 . . . s N ∈Z 2 R 3 N i=1 |∇x i ψ(x 1 , s 1 . . . , x N , s N )| 2 dx 1 . . . dx N Vne[ψ] is the electron-nuclei interaction energy, Vne[ψ] = s 1 ∈Z 2 R 3 . . . s N ∈Z 2 R 3 N i=1 v R (x i )|ψ(x 1 , s 1 . . . , x N , s N )| 2 dx 1 . . . dx N and, Vee[ψ] is the electron-electron interaction energy Vee[ψ] = s 1 ∈Z 2 R 3 . . . s N ∈Z 2 R 3 N i=1 N j=i+1 1 |x j -x i | |ψ(x 1 , s 1 . . . , x N , s N )| 2 dx 1 . . . dx N
Here, the ground state quantum N -body problem refers to the problem of finding equilibrium states for system of type (1.1). A theorem proved first by Zhislin [START_REF] Grigorii | Discussion of the spectrum of schrödinger operators for systems of many particles[END_REF], guarantees the existence of the minimum in (1.4). Some variants can be found in the literature due, for instance, to Lieb & Simon [START_REF] Elliott | The hartree-fock theory for coulomb systems[END_REF], Lions [START_REF] Lions | Solutions of hartree-fock equations for coulomb systems[END_REF] and Friesecke [START_REF] Friesecke | The multiconfiguration equations for atoms and molecules: charge quantization and existence of solutions[END_REF].

In computational chemistry, the Electronic minimization problem methods are used to compute the equilibrium configuration of molecules and solids. However, the most computationally practicable methods in this context are not numerical methods (in the sense in which this terminology is used in mathematics), but power series of analytical solutions of reduced models. 1.2. Probabilistic Interpretation and Marginals. The square norm of the wave function ψ(x 1 , s 1 , • • • , x N , s N ) can be interpreted as a N -point probability density distribution for the electrons are in the points x i with spins s i .

R 3N |ψ(x 1 , s 1 , • • • , x N , s N )| 2 dx = 1.
If ψ is a L 2 (R 3N ) function we can define the single particle density by

ρ(x 1 ) = N R 3(N -1) γ N (x 1 , s 1 , • • • , x N , s N )dx 2 , • • • dx N , (1.5)
where γ N represents the N -point position density

γ N (x 1 , • • • , x N ) = s 1 ,••• ,s N ∈Z 2 |ψ(x 1 , s 1 , • • • , x N , s N )| 2 . (1.6)
Analogously, we can define the k-density

γ k (x 1 , . . . , x k ) = N k R 3(N -k) γ N (x 1 , s 1 , • • • , x N , s N )dx k+1 . • • • dx N (1.7)
We remark that when the particle is spinless the single particle density (1.5) can be simply be written as

ρ(x 1 ) = N R 3(N -1) |ψ(x 1 , . . . , x N )| 2 dx 2 , • • • dx N .
The relevance of ρ and γ 2 to compute the ground state energy (1.4) is due to the fact that the interation nuclei-nuclei energy Vne[ψ] and electron-electron energy Vee[ψ] depends only on these.

Lemma 1.2 ([27]). If ψ, ∇ψ ∈ L 2 ((R 3 × Z 2 ) N , R), ψ antisymmetric and ψ 2 = 1, then Vne[ψ] = R 3 v(x)ρ(x)dx, Vee[ψ] = R 6 1 |x -y| γ 2 (x, y)dxdy.
. It is natural to wonder about the space of those densities arising from such a density ψ is defined, i.e, the space A of densities ρ : R 3 → R verifying (1.5) for an antisymmetric ψ such that ψ, ∇ψ ∈ L 2 ((R 3 × Z 2 ) N , R) and ψ 2 = 1. An explicit characterization exists (due to Lieb [START_REF] Elliott | Density functionals for coulomb systems[END_REF]):

A = ρ : R 3 → R : ρ ≥ 0, √ ρ ∈ H 1 (R 3 ) and R 3 ρ(x)dx = N .
(1.8)

Density Functional Theory (DFT).

The Density Functional Theory is the standard approximation to quantum mechanics in simulations of a system with more than a dozen electrons and showed to be successful in many instances but has rare drastic failures as, for instance, in predicting the behavior of Cr 2 molecules [START_REF] Cohen | Challenges for density functional theory[END_REF]. DFT theory approximates quantum mechanics via variational principles for the marginal density

ρ(x 1 ) := R d(N -1) |ψ(x 1 , . . . , x N )| 2 dx 2 . . . dx N , (1.9)
where ψ is a wave function associated to the N -body quantum problem as in equation (1.5). Roughly speaking, DFT models are semi-empirical models of the pair density

γ 2 (x 1 , x 2 ) := R 3(N -2) |ψ(x 1 , . . . , x N )| 2 dx 3 . . . dx N . (1.10)
in terms of its marginal ρ. We simply write ψ → ρ and ψ → γ 2 to denote the relation between ψ and ρ, and ψ and γ 2 . This means, respectively, that ψ has single particle density ρ and pair density γ 2 .

Concerning the ground state problem (1.4), we notice that, after the lemma 1.2, we are able to write the Electronic minimization problem (1.4) as

E 0 = inf ρ∈A F HK [ρ] + R 3 v(x)ρ(x)dx , (1.11)
with

F HK [ρ] = inf T [ψ] + Vee[γ 2 ] : ρ ∈ A, ψ → ρ . (1.12)
where F HK is the so called Hohenberg-Kohn or Levy-Lieb functional. For a mathematical sake of completeness we also consider a weak repulsive interaction given by Vee

[ρ] = K N i=1 N j=1 -|x j -x i | 2 γ N (x 1 , • • • , x N )dx 1 . . . dx N , (1.13)
in a compact set K ⊂ R 3N and +∞ otherwise. The next theorem is the main result in Density Functional Theory.

Theorem 1.3 (Hohenberg-Kohn, Levy-Lieb). Let ρ ∈ P(R 3 ) be a probability density such that ρ ∈ H 1 (R 3 ) and f : R → R be a continuous function. There exists a functional F HK : P(R 3 ) → R depending only on the single-particle density ρ such that for any potential Vnn, the exact Quantum Mechanics ground state energy (1.4) satisfies

E 0 = min F HK [ρ] + N R 3 v R (x)ρ(x)dx : ρ ∈ P(R 3 ) (1.14)
Moreover, F HK [ρ] is given itself as a minimum problem

F HK [ρ] = min ψ, - 2 2 ∆ + N i=1 N j=i+1 f (x j -x i ) ψ : ψ ∈ H, ψ → ρ
where, ψ → ρ means that ψ has single-particle density ρ and A is the set of ρ :

R 3 → R such that ρ ≥ 0, √ ρ ∈ H 1 (R 3 ) and R 3 ρ(x)dx = N .
The Hohenberg-Kohn theorem states that the functional F HK is universal in the sense that it does not depend on the molecular system under consideration. From a Physics point of view, it garantees that in a molecular system of N electrons, the single electron density ρ determines the pair density of the system (see corollary 1.4). From a mathematical perspective, the proof of this theorem is a functional analysis exercise that, for sake of completeness, we are going to present in a version we learnt from Gero Friesecke.

Proof. Firstly we remark that the non-universal part of the energy only depends on ρ:

ψ, N i=1 v(x i )ψ = R 3N N i=1 v(x i )|ψ(x 1 , . . . , x N )| 2 dx = N R 3 v(x)ρ(x)dx
Then, we rewrite problem (1.4) as a double inf problem

E 0 = inf ψ ψ, H el ψ + N v(r)ρ(r)dr = inf ρ inf ψ →ρ ψ, H el ψ + N v(r)ρ(r)dr = inf ρ F HK [ρ] + N v(r)ρ(r)dr where F HK [ρ] = inf ψ →ρ ψ, H el ψ .
The present version stated in theorem 1.3 is due by Levy and Lieb. The next corollary contains the main physical and mathematical consequence of the Hohenberg-Kohn theorem (1.3).

Corollary 1.4 (HK Theorem and Coupling problem). Let ρ ∈ P(R 3 ) be a measure and ψ ∈ H 1 (R 3N ) be an antisymmetric function with ψ 2 = 1. There exists a universal map from single-particle density ρ(x 1 ) to a pair densities γ 2 (x 1 , x 2 ) -or even k-density γ k (x 1 , . . . , x k )which gives the exact pair density of any N -electron molecular ground state ψ(x 1 , . . . , x N ) in terms of its single-particle density.

Proof. Consider ψ ∈ argmin{ ψ, H el ψ : ψ → ρ} and define by γ k the universal k-point density of that minimizer, i.e.

γ k (x 1 , • • • , x k ) = |ψ(x 1 , • • • , x N )| 2 d k+1 • • • dx N
Of course, γ k may be not unique, since the ground state may be degenerate. The existence of a map T : R 3 → R 3 which determines the density of γ 2 from ρ is highly not trivial and not feasible because still uses high dimension wave function.

1.4. "Semi-classical limit" and Optimal Transport problem. A natural approach to understand the behavior of the Hohenberg-Kohn functional F HK [ρ] (1.12) is to study separately, the contributions of the kinetic energy and the Coulomb interaction electron-electron energy. A well-know method in DFT literature, known as adiabatic connection, is to neglect the kinetic energy, in which the electron-electron interaction is rescaled with a real parameter λ while keeping the density ρ fixed [START_REF] David | The exchange-correlation energy of a metallic surface[END_REF][START_REF] Yang | Generalized adiabatic connection in density functional theory[END_REF]

F HK λ [ρ] = min ψ∈H ψ→ρ ψ, - 2 2 ∆ + N i=1 N j=i+1 λ |x j -x i | )ψ .
(1.15)

The strictly correlated electron limit (λ → ∞), up to relaxing in the space of probability measures, was first considered by two papers in physics literature: Seidl [START_REF] Seidl | Strong-interaction limit of density-functional theory[END_REF] and Seidl, Gori-Giorgi & Savin [START_REF] Seidl | Strictly correlated electrons in densityfunctional theory: A general formulation with applications to spherical densities[END_REF]. In [START_REF] Gori-Giorgi | Density-functional theory for strongly interacting electrons[END_REF], Gori-Giorgi, Seidl and Vignale interpreted the strong-interaction regime as a mass transportation problem with a Coulomb cost.

Later, an equivalent limit, so-called "semi-classical limit", was made mathematically rigorous in the two particles case by Cotar, Friesecke & Klüppelberg [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF]. They considered the Hohenberg-Kohn functional F HK = F HK as a function of both ρ and (λ = 1) and proved that -up to passage of the limit → 0 -the F HK reduces to the following functional

F [ρ] = inf R 6 1 |x -y| γ 2 (x, y)dxdy : γ 2 ∈ A 2 , γ 2 → ρ , (1.16)
where γ 2 → ρ means that γ 2 satisfies the equation (1.7), and the set A 2 of admissible pair density functions is defined by the image of A under the map ψ → γ 2 . As pointed out in [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF], instead of the corresponding single particle density case, we do not know any characterization of the space of admissible pair density function A 2 . We state in the next theorem the semi-classical limit for the 2-particles case. The general case of N particles is still an open problem. Theorem 1.5 ("Semi-classical limit" for N = 2, Cotar-Friesecke-Klüppelberg, [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF]). Let ρ : R 3 → [0, ∞] be a probability density such that √ ρ ∈ H 1 (R 3 ), N = 2 and f : R → R be the Coulomb or the repulsive harmonic interaction potentials. Then,

F HK [ρ] = min ψ∈H 1 (R 6 ),ψ→ρ ψ, - 2 2 ∆ + f (x 2 -x 1 ) ψ → →0 min γ∈Π 2 (R 6 ,ρ) R 3N f (x 2 -x 1 )dγ(x 1 , x 2 ) =: F OT [ρ], (1.17)
where the set Π 2 (R 6 , ρ) denotes the set of measures γ ∈ P(R 6 ) having ρ as marginals, i.e. (e i ) γ = ρ where, for i = 1, 2, e i : R 6 = R 3 × R 3 → R 3 are the projection maps.

The main difficulty in proving Theorem 1.5 is that for any transport plan given by a density

ψ, γ = |ψ(x 1 , . . . , x N )| 2 dx 1 . . . x N , it may be that ψ ∈ H 1 (R 3N ), ψ ∈ L 2 (R 3N ) and T [ψ] = ∞.
Moreover, smoothing the optimal γ does not work at this level, because this may change the marginals of the problem (1.17). Cotar, Friesecke & Klüppelberg developed a smoothing technique in order to deal with this problem without changing the marginals. A complete proof of the previous theorem can be found in [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF].

Let us precise somewhat the terminology. In Quantum Mechanics, the "semi-classical limit" has a precise meaning: it is an asymptotic regime for the Hamiltonian dynamics of a Quantum system defined in a Hilbert space and it is given by a Weyl-Wigner quantization (or quantization by deformation studied in the more abstract context of Poisson manifolds). In those specific cases, the limit → 0 is called "semi-classical" limit, because the first and second order terms of an asymptotic expansion of the Hamiltonian operator is given by "classical" terms, functions of the Hamiltonian function in a symplectic manifold [START_REF] Lawrence | Lectures on semiclassical analysis[END_REF][START_REF] Robert | Autour de l'approximation semi-classique[END_REF].

In DFT context, that limit seems to be up to now merely a question of re-scaling. The minimizers of F OT [ρ] are candidates of ansatz to develop approximating methods to compute the ground state energy of the Electronic Hamiltonian.

At that point, it is natural to define the DFT-Optimal Transportation ground state E DF T -OT 0 [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF],

E DF T -OT 0 [ρ] = inf T [ψ] + Vne[ψ] + E OT [ψ] : ψ ∈ H = inf T QM [ρ] + Vne[ρ] + E OT [ρ] : ρ ∈ A . (1.18)
where ρ represents the single particle density (see (1.5)), T QM and E OT are defined, respectively, by

T QM [ρ] = inf T [ψ] : ψ ∈ H, ψ → ρ and E OT [ρ] = inf R 3N N i=1 N i=1 f (x j -x i )dγ(x 1 , . . . , x N ) : γ ∈ Π((R d ) N , ρ) , (1.19)
where Π((R d ) N , ρ) denotes the set of probability measures γ

: (R d ) N → R + having ρ has marginals. The problem of minimizing E DF T -OT 0 in (1.18) is called DFT-OT problem. Notice that, in the case where f is the Coulomb potential, Vee[γ 2 ] ≥ E OT [ρ], for every probability density ρ ∈ H 1 (R 3 ).
Finally, taking the infimum in both sides of the previous equation, we have

Theorem 1.6 (Cotar-Friesecke-Klüppelberg, [27]). Consider f (|x j -x i |) = 1/|x j -x i |. For every N , and any potential v R ∈ L 3/2 + L ∞ (R 3
), the density functional with electron-electron interaction energy (1. [START_REF] Chen | Pair densities in density functional theory[END_REF]) is a rigorous lower bound of the Electronic minimization problem (1.4)

E 0 ≥ E DF T -OT 0
It turns out now that the central question of DFT-OT problem with Coulomb potentials is to characterize the minimizers γ of E OT [ρ]. Those minimizers could be used as ansatz in order to find minimizers of the Hohenberg-Kohn F HK [ρ]. In [START_REF] Seidl | Strong-interaction limit of density-functional theory[END_REF], the physicist Seidl formulated the following conjecture.

Conjecture 1.7 (Weak Seidl Conjecture [82]

). There exists a deterministic minimizer of

E OT [ρ] with Coulomb potential/cost c(x 1 , . . . , x N ) = N i=1 N j=i+1 1 |x i -x j | .
In other words, for Coulomb-type electron-electron interactions, Seidl conjectures 1 state that, at least for radially symmetric measures ρ, there exists an optimal measure γ for

E OT of form γ = (Id × T 1 × • • • × T N -1 ) ρ, with T 1 , . . . , T N -1 : R 3 → R 3
. This maps T i are called co-motion functions and, due to the symmetries of the problem, a natural group law is required for them, namely that for every i, j there exists k such that

T i • T j = T k ; this is satisfied if for example T i = T (i) (cyclical case, see Equation (2.9)).
It is natural to ask the same kind of question for other cost functions. In the following we mention some interesting results:

(i) In the 2-electrons case, the minimizers of E DF T -OT 0 are well-understood [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF] and there exists a map which correlates the density of a given minimization plan of E OT [ρ], see Theorem 3.6.

(ii) Gangbo & Świȩch (without being aware of Seidl) reinforces the conjecture 1.7 showing in [START_REF] Gangbo | Optimal maps for the multidimensional mongekantorovich problem[END_REF] that for the attractive harmonic cost c(x 1 , . . . , x N ) = N i,j=1 |x i -x j | 2 there exists a unique optimal plan, and it is a deterministic one.

(iii) Colombo, De Pascale & Di Marino [START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF] answer affirmatively (strong) Seidl conjecture in the one dimensional case (d = 1), as we will describe in section 3.

1 There is also a stronger version of the conjecture, the Strong Seidl Conjecture, where the author describes explicitly a possible optimal map in the radially symmetric case: for the precise statement see the conjecture 3.9 and the lemma 3.8 in section 3.

(iv) The authors proved that "Seidl's conjecture" is not true for the repulsive harmonic cost c(x 1 , . . . , x N ) = N i,j=1 -|x i -x j | 2 (see theorem 4.6). Moreover, for the same cost and for a particular single particle density, there exists "fractal-like" optimal transport maps, i.e. an optimal map T which are not differentiable at every point. We are going to state and prove those statements, see section 4. Definition 2.1 (Coupling). Let (X, µ), (Y, ν) be two probability spaces. A coupling γ is a measure on the product X × Y such that γ admits µ and ν as marginals, respectively, on X and Y , i.e.

DFT meets Optimal Transportation Theory

(e 1 ) γ = µ, (e 2 ) γ = ν where e 1 , e 2 are respectively the projection maps

(x, y) ∈ X × Y → x ∈ X and (x, y) ∈ X × Y → y ∈ Y .
It is easy to show that a coupling of µ and ν is equivalent to

(i) For all measurable sets A ⊂ X and B ⊂ Y , one has γ[A × Y ] = µ[A] and γ[X × B] = ν[B]. (ii) For all integrable measurable functions u, v on X, Y , X×Y (u(x) + v(y))dγ(x, y) = X u(x)dµ(x) + Y v(y)dν(y)
A trivial coupling of (X, µ) and (Y, ν) is given by the product measure µ × ν. The symbol ρ → ψ defined in eq. (1.9) means that the measure γ = |ψ(x 1 , . . . , x N )| 2 dx 1 . . . dx N is a coupling, for instance, between ρ 1 and ρ k-1 . In the classical theory of 2-marginals optimal transport [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF][START_REF] Villani | Optimal transport: old and new[END_REF] a coupling is also called a transport plan.

In the case of a 2-molecular system with Coulomb interactions, up to passage to the limit → 0, the theorem 1.5 states that there exists a special coupling γ, minimizer of (1.17), which is concentrated on the graph of a measurable function. That motivates the next definition Definition 2.2 (Deterministic Coupling). A coupling γ of two probability spaces (X, µ) and (Y, ν) is said to be deterministic if does there exists a measurable function T : X → Y such that γ = (Id, T ) µ.

In that terminology, the main question is to understand when the coupling given by the DFT-Optimal Transport (1.19) is deterministic. Moreover, Seidl's conjecture can be rephrased in the following way: Does exist a deterministic coupling in the Optimal Transportation problem (1.19) for Coulomb costs? 2.2. Multimarginal Optimal Transportation Problem. The "semi-classical limit" in the 2 particles case (thm. 1.5) and Seidl's Conjecture 1.7, motivate us to introduce a more general coupling problem.

We denote by N the number of particles (marginals), d the dimension of the space where those particles live (typically d = 3) and I = {1, . . . , N} an index set. Take x i : R d → R as the position of a single particle, then the configuration of a system of N -particles is given by a vector x = (x 1 , x 2 , . . . , x N ) or, more precisely, by a function x :

I → R d .
Let us consider the distribution of the i-th particle given by a density function

ρ i : R d → [0, 1], i ∈ I and c : R N d × R N d → [0, ∞], a continuous cost function. We set µ i = ρ i L d , ∀i ∈ I,
where L d is the Lebesgue measure on R d , and we want to characterize the optimal coupling of the Monge-Kantorovich problem

(MK) inf γ∈Π(R dN ,µ 1 ,...,µ N ) (R d ) N c(x 1 , . . . , x N )dγ(x 1 , . . . , x N ), (2.1)
where Π(R dN , µ 1 , . . . , µ N ) denotes the set of couplings γ(x 1 , . . . , x N ) having µ i as marginals, for all i ∈ I, i.e. the set Π(R dN , µ 1 , . . . , µ N ) = γ ∈ P(R dN ) : (e i ) γ = µ i , ∀i ∈ I and, for all i ∈ I, e i : R dN → R d are the projection on the i-th component and P(R dN ) denotes the set of probability measures in R dN . The elements of Π(R dN , µ 1 , . . . , µ N ) are called transport plans. In the DFT-OT problem it was mentioned that an interesting case is when all marginals are the same and equal to a measure µ on R d (it is quite natural in the DFT framework as the electrons are indistinguishable). In that case, we denote the set of transport plans simply by Π N (µ).

Notice that the existence of a minimum for (2.1) is not a big deal and it is quite standard in Optimal Transport Theory. Indeed, Π(R dN , µ 1 , . . . , µ N ) is trivially non empty, since the independent coupling µ 1 × • • • × µ N belongs to Π(R dN , µ 1 , . . . , µ N ); the set Π(R dN , µ 1 , . . . , µ N ) is convex and compact for the weak*-topology thanks to the imposed marginals; and moreover the quantity to be minimized γ → c dγ is linear with respect to µ. Hence, we can guarantee the existence of a minimum for (2.1) by imposing a very weak hypothesis on the cost function, such as lower-semicontinuity.

However, we are interested in characterizing some class of optimal transport plans or, at least, understand when that optimal coupling is deterministic. Thus, we would like to study when minimizers for (MK) correspond to minimizers the multi-marginal Monge problem

(M N ) min T =(T i ) N i=1 ∈T N (R d ) N c(x, T 2 (x), . . . , T N (x))dµ 1 (x), (2.2)
where

T N = T = (T i ) N i=1 : T i is a Borelian map such that T i µ 1 = µ i , ∀i ∈ {1, . . . , N } and T 1 = Id .
As usual, it is easy to see that the set of admissible plans in (2.1) "includes" the set of Monge transport maps in (2.2): in fact, given a transport map T = (Id, T 2 , . . . , T N ), we can consider a measures γ T on (R d ) N defined by

γ T = (Id ×T 2 × • • • × T N ) µ. Let h : (R d ) N → R be a γ T -measurable function, then (R d ) N h(x 1 , . . . , x N )dγ T (x 1 . . . , x N ) = R d h(x 1 , T 2 (x 1 ), . . . , T N (x 1 ))dµ 1 (x 1 ); and if we have γ T -measurable function f : R d → [0, ∞] ∀i ∈ I, (R d ) N f (x i )dγ T (x 1 , . . . , x N ) = R d f (x i )dµ i (x i ),
and so, γ T is a transport plan. In particular, the value of the minimization over Monge-Kantorovich couplings is less than or equal to the value of the minimization over the deterministic couplings

min(MK) ≤ inf(M N ).
Contrary to in (2.1), the existence of the optimal deterministic coupling in (2.2) is not obvious and the difficulties do not lie on the "multi-marginality" of the problem. In fact, at that point, it is not possible to apply standard methods in Calculus of Variations even in the classical 2-marginals case in order to guarantee the existence of the minimum (2.2).

Intuitively, the difference between (2.1) and (2.2) is that in M N almost every point x 1 ∈ R d is coupled with exactly one point x i ∈ R d for each i ∈ {2, . . . , N } whereas throughout the problem (2.1) we allow the splitting of mass in the "transportation" process.

The problem (2.1), when N = 2, was studied by L. Kantorovich in 1940's [START_REF] Kantorovitch | On the translocation of masses[END_REF][START_REF] Kantorovich | On a problem of monge[END_REF], without being aware of the famous article of Monge [START_REF] Gaspard Monge | Mémoire sur la théorie des déblais et des remblais[END_REF]. In 1975, L. Kantorovich won the Nobel Prize together with C. Koopmans "for their contributions to the theory of optimum allocation of resources".

An important disadvantage in using the relaxed approach (2.1) in the multi-marginal OT is that the set of optimal transportation plans could be very large and we could need to select some special classes of transportation plans.

In the 2-marginal setting, the two formulations (2.1) and (2.2) are proved to have the same value for general Polish Spaces by Pratelli in [START_REF] Pratelli | On the equality between monge's infimum and kantorovich's minimum in optimal mass transportation[END_REF].

Theorem 2.3 (Monge equals Monge-Kantorovich in the 2-marginals case, [START_REF] Pratelli | On the equality between monge's infimum and kantorovich's minimum in optimal mass transportation[END_REF]). Let (X, µ), (Y, ν) be Polish spaces, where µ ∈ P(X) and ν ∈ P(Y ) are non-atomic probability measures. If c :

X × Y → R ∪ {∞} is a continuous cost, then min X×Y c(x, y)dγ(x, y) : γ ∈ Π(X × Y, µ, ν) = = inf X c(x, T (x))dµ(x), : T µ = ν, T : X → Y Borel (2.3)
where Π(X × Y ) denotes the set of 2-marginals transport plans γ ⊂ X × Y having µ and ν as marginals.

In the particular case of the multi-marginal OT problem in (2.2) and (2.1) when all marginals are equal to ρL d , we can apply the previous theorem (2. Notice that, in general, the equivalence between (2.1) and (2.2) is not an immediate consequence of the theorem 2.3. In particular, remark that the image measure of Y = R d(N -1) is not prescribed, but only its marginals.

Finally, we conclude this section by giving an example of a particular continuous cost function where the optimal coupling γ is not deterministic, even in the 2 particles case. The nonexistence of an optimal transport maps typically happens when minimizing sequences of γ exhibit strong oscillations.

Example 2.5 (G. Carlier). Suppose µ uniformly distributed in [0, 1], ν uniformly distributed in [-1, 1] and the continuous cost c : R 2 → R given by

c(x 1 , x 2 ) = (x 2 1 -x 2
2 ) 2 . Since c ≥ 0 we have that an optimal plan for this cost is given by γ = µ ⊗ ( 1 2 δx + 1 2 δ -x ) since c dγ = 0. We want to show there doesn't exit a map T such that the cost attains zero: suppose in fact that (x 2 -T (x) 2 ) 2 dx = 0 but then one should have T (x) ∈ {x, -x} almost everywhere. Hence,

T (x) = (χ A -χ [0,1]\A )x for a measurable A ⊂ [0, 1]. But then T µ = L 1 | A + L 1 | Ã = ν, where x ∈ Ã iff -x ∈ [0, 1] \ A,
and so we obtain a contradiction since T wouldn't admissible. It is however easy to construct admissible maps Tn, linear on [i/2n, (i + 1)/2n), such that c(x 1 , Tn(x 1 ))dµ(x 1 ) tends to zero.

Remark that, we state by simplicity the example in the 2-marginal case, but the same conclusion and arguments apply in the larger interactions case under the cost

c(x 1 , . . . , x N ) = N i=1 N j=i+1 (x 2 i -x 2 j ) 2 .
Figure 1. (example 2.5 ) Left: support of the optimal coupling γ . Right: optimal coupling γ. The simulation has been performed by using the iterative method described in section 6.

Dual formulation.

We consider the formal convex dual formulation of (2.1)

(K N ) sup R d N i=1 u i (x i )dµ i (x i ) : N i=1 u i (x i ) ≤ c(x 1 , . . . , x N ) , (2.4)
where u i ∈ L 1 (R d , µ i ), called Kantorovich potentials or simply potentials, are upper semicontinuous for all i ∈ {1, . . . , N }. This is often called the dual problem of (2.2).

Notice that, when all marginals are equal to µ 1 = ρ(x)L d and the cost is symmetric, we can assume the Kantorovich potentials u i (x i ) are all the same u(x i ), and we can rewrite the constraint in (2.4) as

N i=1 u(x i ) ≤ c(x 1 , . . . , x N ).
Among all Kantorovich potentials, a particular class is of special interest.

Definition 2.6 (c-conjugate function). Let c : (R d ) N → R ∪ {∞} a Borel function. We say that the N -tuple of functions (u 1 , . . . , u N ) is a c-conjugate function for c, if u i (x i ) = inf c(x 1 , . . . , x N ) - N j=1,j =i u j (x j ) : x j ∈ X j , j = i , ∀i = 1, . . . , N.
As we will see in the next theorem 2.7, the optimal potentials u 1 , . . . , u N in (K N ) are cconjugates, and therefore semi-concave. So, they admit super-differentials ∂u i (x i ) at each point x i ∈ R d and then, at least in the compact case, we can expect that for each x i ∈ R d , there exists (x j ) j =i ∈ R d(N -1) , such that

u i (x i ) = c(x 1 , . . . , x N ) - j =i u j (x j ). Moreover, if c : (R d ) N → R is differentiable and |u i (x i )| < ∞, for some x i ∈ R d , u i is locally
Lipschitz. The well-known general result linking both (2.1) and (2.4) has been proven by Kellerer in 1984.

Theorem 2.7 (Kellerer, [START_REF] Hans | Duality theorems for marginal problems[END_REF]). Let (X 1 , µ 1 ), . . . , (X N , µ N ) be Polish spaces equipped with Borel probability measures µ 1 , . . . , µ N . Consider c : X 1 . . . X N → R a Borel cost function and assume that c = sup X c < ∞ and c = inf X c > -∞. Then, (i) There exists a solution γ of the Monge-Kantorovich problem (2.1) and a c-conjugate solution (u 1 , u 2 , ..., u N ) of the dual problem (2.4). (ii) "Duality holds",

inf X 1 ...X N cdγ : γ ∈ Π(µ 1 , . . . , µ N ) = sup N i=1 X i u i (x i )dµ i : u i ∈ D , (2.5)
where D is the set of functions u i : R d → R, i = 1, . . . , N , such that

N i=1 u i (x i ) ≤ c(x 1 , . . . , x N ), and 1 N c -(c -c) ≤ N i=1 u i (x i ) ≤ 1 N c.
(iii) For any solution γ of (2.1), any conjugate solution of (2.4) and γ a.e. (x 1 , . . . , x N ), we have

N i=1 u i (x i ) = c(x 1 , . . . , x N ).
We remark that despite its generality, Kellerer's theorem can not be applied directly to Coulomb-type costs, since Coulomb cost does not satisfies the boundedness hypothesis c = sup X c < ∞ and can not be bounded by functions in L 1 (R d , ρL d ).

A generalization of Kellerer's theorem due to M. Beiglboeck, C. Leonard and W. Schachermayer [START_REF] Beiglböck | A general duality theorem for the monge-kantorovich transport problem[END_REF], extends this duality for more general costs, but also cannot be applied in this context.

Recently, De Pascale [START_REF] De | Optimal transport with coulomb cost. approximation and duality[END_REF] extended the proof of duality which can be applied also for Coulombtype costs, see Theorem 3.1. In fact, suppose µ = ρ(x)dx a probability measure on R d which does not give mass to sets of cardinality less than or equal to d -1, then min

γ∈Π N (µ) (R d ) N N i=1 N j=i+1 1 |x i -x j | dγ = max u∈ D R d N u(x)ρ(x)dx (2.6)
where D is the set of potentials

u ∈ L 1 (R d , µ) such that N i=1 u(x i ) ≤ N i=1 N j=i+1 1 |x j -x i | , ⊗ N i=1 ρ -almost everywhere. (2.7)
Moreover, we get the conclusions (i) and (iii) of the theorem 2.7. Physically, the constraint (2.7) means that, at optimality, the allowed N d configuration space should be a minimum of the classical potential energy.

Remark 2.8. More generally, the main theorem in [START_REF] De | Optimal transport with coulomb cost. approximation and duality[END_REF] shows that the Kantorovich duality holds, for instance, for costs of form

N i=1 N j=i+1 1 |x j -x i | s , s ≥ 1.

Geometry of the Optimal Transport sets.

Our goal in this section is to present the relation between the support of a coupling γ and optimality in the Monge-Kantorovich problem (2.1). In the following, we will just summarize key results necessary to present recent developments of optimal transportation with Coulomb and repulsive harmonic-type costs.

Roughly speaking, it is well-known in optimal transport theory with 2-marginals that, for a wide class of costs c, a coupling γ is optimal if, and only if, the support of γ is concentrated in a c-cyclically monotone set [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF][START_REF] Pratelli | On the sufficiency of c-cyclical monotonicity for optimality of transport plans[END_REF][START_REF] Schachermayer | Characterization of optimal transport plans for the monge-kantorovich-problem[END_REF]: Definition 2.9 (c-cyclically monotone, 2-marginal case). Let X, Y be Polish spaces, c : X × Y → R ∪ {∞} be a cost function and Γ be subset of the product space X × Y . We say that Γ is a c-cyclically monotone set if, for any finite couples of points

{(x i , y i )} M i=1 ⊂ Γ, M i=1 c(x i , y i ) ≤ M i=1 c(x i , y σ(i) ).
In [START_REF] Pass | Structural results on optimal transportation plans[END_REF], Pass suggested a possible extension of the concept of c-cyclically monotone sets for the multi-marginal case. For our proposes, it will be enough to consider this notion in R d . Definition 2.10 (c-cyclically monotone set). Let c : R dN → R be a cost function. A subset Γ ⊂ R dN is said to be c-monotone with respect to a partition (p, p c ), p ⊂ I, if, for all x = (x 1 , . . . , x N ), y = (y 1 , . . . , y N ) ∈ Γ c(x) + c(y) ≤ c(X(x, y, p)) + c(Y (x, y, p)), where X(x, y, p), Y (x, y, p) ∈ R N are functions obtained from x and y by exchanging their coordinates on the complement of p, namely

X i (x, y, p) = x i , if i ∈ p y i , if i ∈ p c and Y i (x, y, p) = y i , if i ∈ p x i , if i ∈ p c ∀i ∈ I
We say that Γ ⊂ R dN is c-cyclically monotone, if it is c-monotone with respect to any partition (p, p c ), p ⊂ I.

The notion of c-cyclical monotonicity in the multi marginal setting was studied by Ghoussoub, Moameni, Maurey, Pass [START_REF] Ghoussoub | Remarks on multi-marginal symmetric mongekantorovich problems[END_REF][START_REF] Ghoussoub | Symmetric monge-kantorovich problems and polar decompositions of vector fields[END_REF][START_REF] Pass | Structural results on optimal transportation plans[END_REF][START_REF] Pass | On the local structure of optimal measures in the multi-marginal optimal transportation problem[END_REF] and successfully applied in the study of decoupling PDE systems [START_REF] Ghoussoub | Decoupling of degiorgi-type systems via multi-marginal optimal transport[END_REF]. For Coulomb costs, it was useful in order to characterize the optimal transport maps in the one-dimensional case [START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF].

More recently, Beiglboeck and Griessler [START_REF] Beiglböck | An optimality principle with applications in optimal transport[END_REF], presented a new notion called finistic optimality inspired by the martingale optimal transport, which is equivalent, in the 2-marginals case, to the definition 2.9.

The next proposition gives a necessary condition on the support for an optimal transport plan γ. For the proof, which is based on the 2-marginal result [START_REF] Smith | On hoeffding-fréchet bounds and cyclic monotone relations[END_REF], we refer the reader to Lemma 3.1.2 of [START_REF] Pass | Structural results on optimal transportation plans[END_REF].

Proposition 2.11 (Support of transport plan are c-cyclically monotonic). Let c : (R d ) N → [0, ∞] be a continuous cost and µ 1 , . . . , µ N absolutely continuous probability measure on R d . Suppose that γ ∈ Π((R d ) N , µ 1 , . . . , µ N ) is an optimal transport plan for Multimarginal Monge-Kantorovich Problem (2.1) and assume (MK) < ∞. Then γ is concentrated in a c-cyclically monotone set.

Symmetric Case.

We are going to show in this section that we can reduce our study of the Monge-Kantorovich (2.1) and the Monge problems (2.2), respectively, to a suitable class of transport maps (called cyclic maps) and transport plans (symmetric plans). Having in mind the DFT-OT problem, it is natural to suppose that the cost functions does not change if there are some symmetries among the configuration system of particles or, equivalently, we can suppose the couplings γ ∈ Π N (µ) of those particles are cyclically symmetric. Definition 2.12 (Cyclically symmetric probability measures). A probability measure γ ∈ P(R dN ) is cyclically symmetric if

R dN f (x 1 , . . . , x N )dγ = R dN f (σ(x 1 , . . . , x N ))dγ, ∀f ∈ C(R N d )
where σ is the cyclical permutation σ(x 1 , x 2 , . . . , x N ) = (x 2 , x 3 , . . . , x N , x 1 ).

We will denote by Π sym N (µ 1 ), the space of all symmetric probability measures on R dN having the same first marginal µ 1 .

Notice that if (µ 1 , . . . , µ N ) are the marginals of γ ∈ Π sym (R dN , µ 1 ), then all marginals µ 1 , . . . , µ N are necessarily the same:

∀ f ∈ C(R d ), ∀i ∈ I, R d f (x i )dµ i (x i ) = R dN f (x 1 , . . . , x N )dγ(x 1 , . . . , x N ) = R dN f (σ 1-i (x 1 , . . . , x N ))dγ = R d f (x 1 )dµ i (x 1 ), where f (x 1 , . . . , x N ) = f (x i ). Briefly, (e i ) γ = µ 1 , ∀i ∈ I,
where e i is the canonical projection.

As previously, we can define the multi marginal symmetric Monge-Kantorovich Problem as follows:

(MKsym) inf

γ∈Π sym N (µ) R dN c(x 1 , . . . , x N )dγ(x 1 , . . . , x N ), (2.8)
and the associated multi marginal cyclically Monge Problem

(M cycl ) inf T µ=µ, T (N ) =I R d c(x, T (1) (x), T (2) (x), . . . , T (N -1) (x)), (2.9)
Notice that, (M cycl ) is related to (M) by choosing T µ = µ, T i+1 (x) = T (i) (x) for every i ∈ {1, . . . , N -1}, T (N ) = Id.

Proposition 2.13. Suppose µ ∈ P(R d ) is an absolutely continuous probability measure with respect to the Lebesgue measure and c

: (R d ) N → R is a continuous symmetric cost. Then, inf(MK) = inf(MKsym).
Proof. The fact that the infimum of (2.1) is less than or equal to (2.8) is obvious. Now, we need to show that for every transport plan γ ∈ Π N (µ) we can associate a symmetric plan in Π sym N (µ). Indeed, given γ ∈ Π N (µ), we define

γ = 1 N ! σ∈S N σ γ,
and finally we notice that γ has the same cost as γ.

We remark that prove inf(M cycl ) ≥ inf(M N ) is not obvious . However, we can avoid this problem in order to prove the equivalence between the Multi-marginal cyclic problem and Multimarginal OT problem, by noticing that inf(Msym ≤ inf(M cycl ). The last part was proved by M. Colombo and S. Di Marino [START_REF] Colombo | Equality between monge and kantorovich multimarginal problems with coulomb cost[END_REF]. We give the precise statement of this theorem. The last result of this section is due Pass [START_REF] Pass | Remarks on the semi-classical hohenberg-kohn functional[END_REF], which states that optimal transport plans for symmetric costs can be supported on sets of Hausdorff dimension equal or bigger than 2N -2.

Theorem 2.16 (Pass, [START_REF] Pass | Remarks on the semi-classical hohenberg-kohn functional[END_REF]). Let c : (R d ) N → R be a symmetric cost and {µ i } N i=1 be radially symmetric and absolutely continuous with respect to Lebesgue measure in R d . Suppose that for every radii (r 1 , ...r N ) the minimizers (x 1 , . . . , x N ) ∈ argmin |y i |=r i c(y 1 , y 2 , . . . y N ) are not all co-linear.

Then, there exists solutions γ for the Monge-Kantorovich symmetric problem (2.8) whose support is at least (2N -2)-dimensional.

Multi-marginal OT with Coulomb Cost

This section is devoted to the summary of the results present in literature on the multi marginal optimal transport with Coulomb cost. We already highlighted the fact that the interesting case is when all the marginals are equal, since this is the physical case. We point out that, as always, the focus will be on characterizing the optimal plans as well as looking for cyclical transport maps: in this direction an important conjecture has been made by Seidl in its seminal paper [START_REF] Seidl | Strong-interaction limit of density-functional theory[END_REF], where he gave an explicit construction of a map for measures with radial symmetry in R d . In the last years the conjecture was proven to be right: first in the 1D case when µ is the uniform measure on an interval, and then generalized to any diffuse measure on the real line. However quite recently various authors disproved the conjecture in the case R d for d ≥ 2 and for N = 3 marginals, yet the Seidl map is still optimal for some class of measures.

We begin by analyzing the general problem, and then we will proceed to look at the radial case.

3.1. General theory: duality, equivalent formulations and many particles limit. We begin by stating the duality theorem in this case proved by De Pascale [START_REF] De | Optimal transport with coulomb cost. approximation and duality[END_REF], by means of Γconvergence of finite dimensional problems for which the duality is classical. Theorem 3.1. Let µ ∈ P(R d ) be a measure that is not concentrated on a set of cardinality smaller or equal than N -1. Then the duality (2.6) holds and the dual problem has a bounded maximizer.

The boundedness of the maximizer u can let us prove also some regularity properties. The first one is the fact that any optimal plan is supported away from the diagonals, while the second one proves second order regularity of u. Lemma 3.2 (Corollary 3.13 in [START_REF] De | Optimal transport with coulomb cost. approximation and duality[END_REF]). Let µ as in Theorem 3.1; then there exists α > 0 such that for every optimal plan γ we have supp(γ) ⊆ {(x 1 , . . . , x N ) :

|x i -x j | ≥ α, ∀ i = j}.
Lemma 3.3. Let µ as in Theorem 3.1 and u be a bounded maximizer in the dual problem in (2.6); then u has a µ-representative that is Lipschitz and semi concave.

Proof. As in the classical case we first prove a structural property of the maximizer u, namely that the constraint (2.7) is saturated, that is (u, . . . , u) is a c-conjugate N -tuple. Let us consider w(x 1 ) = ess-inf

(x 2 ,...,x N )∈R d(N -1)    1≤i<j≤N 1 |x i -x j | - N i=2 u(x i )    ,
where the ess-inf is made with respect to the measure µ ⊗N -1 . It is obvious that w(x 1 ) ≥ u(x 1 ) for µ-a.e x 1 , thanks to (2.7). Suppose that w > u in a set of positive measure; but then (w, u, . . . , u) would be a better competitor in the (not symmetric) dual problem contradicting the fact that u is a maximizer (we use the fact that the symmetric dual problem and the not symmetric one have the same values). Now we have that w is defined everywhere and so we can talk about its regularity.

Let d > 0 be a number such that there exists points p 1 , . . . p N in the support of µ such that |p i -p j | ≥ d; this number exists as long as supp(µ) ≥ N . Let us fix ε > 0 such that

1 2ε > 2N (N -1) d + 2N u ∞ and ε ≤ d/8
; in this way it is true that if we define

wx 0 (x 1 ) = ess-inf |x i -x 0 |≥ε, i=2,...,N    1≤i<j≤N 1 |x i -x j | - N i=2 u(x i )    ,
we have wx 0 = w on B(x 0 , ε). In fact in the definition of w we can choose x i among the p i (or very close to them, as they belong to the support) such that |x i -

x j | > d/2 -2ε ≥ d/4
for every i = j, and so we have that w ≤ 2N (N -1) d

+ N u ∞ but then it is clear that for any (x 2 , . . . , x N ) such that |x i -x 0 | ≤ ε we will have

1≤i<j≤N 1 |x i -x j | - N i=2 u(x i ) ≥ 1 2ε > w(x 1 ) ∀x 1 ∈ B(x 0 , ε);
this proves that in fact wx 0 = w on the set B(x 0 , ε) and so in particular also in B(x 0 , ε/2). But in this set wx 0 is Lipschitz and semi concave since it is an infimum of uniformly C ∞ functions on B(x 0 , ε/2). Moreover the bounds on the first and second derivatives don't depend on x 0 but only on ε, that is fixed a priori, and so by a covering argument we obtain the thesis.

Another interesting reformulation of the Coulomb-like problem, or more generally when we have only interaction between two particles, can be found in [START_REF] Friesecke | N-density representability and the optimal transport limit of the hohenberg-kohn functional[END_REF] where, seeking a dimensional reduction of the problem, the authors use the fact that if γ ∈ Π N (µ) is a symmetric plan then

R N d 1≤i<j≤N 1 |x i -x j | dγ = N 2 R 2d 1 |x -y| d(e 1 , e 2 ) γ,
that is, since the electron are indistinguishable, it is sufficient to look at the potential energy of a couple of electrons and then multiply it by the number of couples of electron. It is clear that γ 2 = (e 1 , e 2 ) γ is a 2-plan whose marginal are still µ; we will say that η ∈ Π 2 (µ) is N -representable whenever it exists γ ∈ Π N (µ) such that η = (e 1 , e 2 ) γ. The equivalent formulation they give is

min γ∈Π N (µ) R N d 1≤i<j≤N 1 |x i -x j | dγ = min η∈Π 2 (µ) η is N -representable N 2 R 2d 1 |x -y| dη.
Unfortunately the conditions for being N -representable are not explicit and they are very difficult to find; this is correlated to the N -representability for matrices, but here, since we are in the semiclassical limit, we deal with densities instead. In [START_REF] Friesecke | N-density representability and the optimal transport limit of the hohenberg-kohn functional[END_REF] the authors propose, as a method for reducing the dimension of the problem, to substitute the condition of being N -representable with that of being k-representable with k ≤ N ; the resulting problem will give a lower bound to the SCE functional .

We present also a general theorem regarding the many particles limit, that embodies the wellknown fact that when N → ∞ then the solution is the mean field one Theorem 3.4 ( [START_REF] Cotar | Infinite-body optimal transport with coulomb cost[END_REF]). Let µ ∈ P(R d ). Then we have that lim

N →∞ 1 N 2 min γ∈Π N (µ) R N d 1≤i<j≤N 1 |x i -x j | dγ = R 2d 1 |x -y| dµ ⊗ µ.
In terms of DFT we are saying that

E OT (µ) = N 2 E M F (µ) + o(N 2 )
, where E M F is the normalized mean field energy E M F (µ) = 1 |x-y| dµ(x) dµ(y). Remark 3.5. The statement about the Coulomb cost in the physical case is quite classical. In fact, for a measure ρ ∈ L 4/3 (R 3 ) the Lieb-Oxford bound holds [START_REF] Elliott | A lower bound for coulomb energies[END_REF][START_REF] Elliott | Improved lower bound on the indirect coulomb energy[END_REF]:

N 2 E M F (ρ) ≥ E OT (ρ) ≥ N 2 E M F (ρ) -CN 4/3 R 3 ρ 4/3 (x) dx,
and so the conclusion is immediate. However, in [START_REF] Cotar | Infinite-body optimal transport with coulomb cost[END_REF] the proof is completely different and relies on the fact that a measure γ ∈ Π 2 (µ) that is N -representable for every N must be in the convex envelope of measures of the type ν ⊗ ν, and then on a direct computation using the Fourier transform. In particular, aside from the Coulomb cost in dimension 3 they prove the theorem for more general costs with binary interaction of the form 1≤i<j≤N l(x i -x j ), where l ∈ C b (R d ) ∩ L 1 (R d ) satisfies l(z) = l(-z) and also l ≥ 0.

3.2.

The Monge problem: deterministic examples and counterexamples. In this section we will illustrate the cases in which we know that there exists a deterministic solution; in those cases moreover there is also a result of uniqueness. We remark that these are the two extreme case, namely the case N = 2 (and any d) and the case d = 1 (and any N ). We will just sketch the proofs in order to make clear the method used here, and why it is difficult to generalize it. We begin with the 2-particles case, in every dimension d: this result was proved in [START_REF] Cotar | Density functional theory and optimal transportation with coulomb cost[END_REF] by means of standard optimal transport techniques. Theorem 3.6. Let µ ∈ P(R d ) be a probability measure that is absolutely continuous with respect to the Lebesgue measure. Then there exists a unique optimal plan γ O ∈ Π 2 (µ) for the problem

min γ∈Π 2 (µ) R 2d 1 |x 1 -x 2 | dγ.
Moreover this plan is induced by an optimal map T , that is, γ = (Id, T ) µ, and

T (x) = x+ ∇φ |∇φ| 3/2
µ-almost everywhere, where φ is a Lipschitz maximizer for the dual problem.

Proof. Let us consider γ a minimizer for the primal problem and φ a bounded and Lipschitz maximizer of the dual problem (it exists thanks to Lemma 3.3). Then we know that

F (x 1 , x 2 ) = φ(x 1 ) + φ(x 2 ) - 1 |x 1 -x 2 | ≤ 0 for µ ⊗ µ-almost every (x 1 , x 2 ),
and we know also that F = 0 γ-almost everywhere. But then F has a maximum on the support of γ and so ∇F = 0 in this set; in particular we have that ∇φ(x 1 ) = -

x 1 -x 2 |x 1 -x 2 | 3 on the support of γ. Finding x 2 we have that x 2 = x 1 -∇φ |∇φ| 3/2 (x 1 ) = T (x 1
) on the support of γ and this implies γ = (Id, T ) µ as we wanted to show.

The first positive N -marginal result for the Coulomb cost is instead shown in [START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF] where, in dimension d = 1, the authors can prove that for non-atomic measure an optimal plan is always "induced" by a cyclical optimal map T . Theorem 3.7. Let µ ∈ P(R) be a diffuse probability measure. Then there exists a unique optimal symmetric plan γ O ∈ Π sym 2 (µ) that solves

min γ∈Π sym N (µ) R N 1≤i<j≤N 1 |x j -x i | dγ.
Moreover this plan is induced by an optimal cyclical map T , that is,

γ O = 1 N !
σ∈S N σ γ T , where γ T = (Id, T, T (2) , . . . , T (N -1) ) µ. An explicit optimal cyclical map is

T (x) = F -1 µ (Fµ(x) + 1/N ) if Fµ(x) ≤ (N -1)/N F -1 µ (Fµ(x) + 1 -1/N ) otherwise.
Here Fµ(x) = µ(-∞, x] is the distribution function of µ, and F -1 µ is its lower semicontinuous left inverse.

Proof. We begin by observing that if γ is a symmetric optimal plan then a stronger statement than Proposition 2.11 holds, namely we have that for every x, y ∈ supp(γ):

c(x) + c(y) = min{c(X(x, σ(y), p)) + c(Y (x, σ(y), p)) : ∀p ⊂ {1, . . . , N }, ∀σ ∈ S N },
where c is the Coulomb cost on the N -tuple of points x = (x 1 , . . . , x N ) and σ acts on the indices. The key point here is that one can show (Proposition 2.4 [START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF]) that this property holds if and only if a more geometric condition holds true for x and y: they are well-ordered. This property amounts to the fact that, up to permute the coordinates of both points, we have

x 1 ≤ y 1 ≤ x 2 ≤ • • • ≤ x N ≤ y N ,
or the other way around. Once this property is proven the rest of proof is rather straightforward, since we proved that the support of (e 1 , e 2 ) γ| D , where D = {x i ≤ x j and y i ≤ y j , ∀ i < j}, is monotone in the usual sense. We remark that the key property does not easily generalize to d > 1 since it (heavily) uses the ordered structure of R.

Let us notice that in this case we don't have uniqueness of optimal plan neither of cyclical optimal maps, as pointed out in Remark 1.2 in [START_REF] Di Marino | Multimarginal optimal transport maps for 1-dimensional repulsive costs[END_REF].

For the general case we thus have to assume d > 1 and N > 2: from now on we will look at the case in which the measure µ has radial symmetry. We can easily see that in this case we can reduce our problem to a multi marginal problem in R. where c 1 is the reduced Coulomb cost

c 1 (r 1 , . . . , r N ) = min    1≤i<j≤N 1 |x j -x i | : |x i | = r i ∀i = 1, . . . , N    . Moreover R dN 1≤i<j≤N 1 |x i -x j | dγ = R d c 1 (r 1 , . . . r N ) dγr.
Given the one dimensional character of this problem, one could expect that the solution is similar to the one depicted in Theorem 3.7. In fact in [START_REF] Seidl | Strictly correlated electrons in densityfunctional theory: A general formulation with applications to spherical densities[END_REF] the authors conjecture a similar structure: Conjecture 3.9 (Strong Seidl Conjencture). Let µ ∈ P(R d ) be an absolutely continuous measure with respect to the Lebesgue measure, with radial symmetry, and let µr

= | • | µ. Let 0 = r 0 < r 1 < . . . < r N -1 < r N = ∞ such that the intervals A i = [r i , r i+1
) have all the same radial measure µr(A i ) = 1/N . Then let F (r) = µr(0, r] be the cumulative radial function and let S : [0, ∞) → [0, ∞) be defined piecewise such that the interval A i is sent in the interval A i+1 in an anti monotone way:

S(r) = F -1 (2i/N -F (r)) if r i-1 ≤ r < r i and i < N F (S(r)) = F -1 (F (r) + 1/N -1) if N is even F -1 (1 -F (r)) if N is odd, if r N -1 ≤ r < r N .
Then S is an optimal cyclical map for the problem (3.1).

However in the recent papers [START_REF] Colombo | Counterexamples to multimarginal optimal transport maps with coulomb cost and radial measures[END_REF][START_REF] Di Marino | The strictly-correlated electron functional for spherically symmetric systems revisited[END_REF] this conjecture is proven to be wrong, looking at radial measures µ concentrated in some thin annulus {1 -δ ≤ |x| ≤ 1 + δ}; in particular in [START_REF] Di Marino | The strictly-correlated electron functional for spherically symmetric systems revisited[END_REF] the 1D problem is proven to be equivalent, when δ → 0, to the repulsive harmonic one, for which high non-uniqueness holds and the Seidl map is not optimal (see the following section). However in [START_REF] Colombo | Counterexamples to multimarginal optimal transport maps with coulomb cost and radial measures[END_REF] also a positive example is found, namely a class of measures for which the conjecture holds.

These results show that the solution to the multi marginal problem with Coulomb cost is far from being understood: in particular there is no clear condition on the marginal µ for the strong Seidl conjecture to hold.

Multi-marginal OT with repulsive Harmonic Cost

This section is devoted to the study of the repulsive harmonic cost. In DFT-OT problem (see section 1), we replace the electron-electron Coulomb interaction by a weak force which decreases with the square distance of the particles. The term weak comes from the fact that, in the repulsive harmonic cost the, interaction cost function has value zero when the particles overlap, instead of ∞ in the Coulomb case.

More precisely, we are interested in characterizing the minimizers of the following problem

(MK weak ) min γ∈Π((R d ) N ,µ 1 ,...,µ N ) (R d ) N N i,j=1 -|x j -x i | 2 dγ(x 1 , . . . , x N ), (4.1)
where µ 1 , . . . , µ N are absolutely continuous probability measures in R d . From a mathematical viewpoint this cost has some advantages compared to the Coulomb one, since here we can do explicit examples. We can interpret the solutions of (4.1) as an ansatz for DFT problem (Hohenberg-Kohn functional (1.12)) for particles interacting under the repulsive harmonic potential. From a technical aspect, this could be an interesting toy model to approach the case of Coulomb cost. From applications, some minimizers of (4.1) seem to have no particular relevance in physics because, as we will see in the examples 4.10, 4.11 and 4.12 below, certain optimal γ of (4.1) allow particles to overlap.

We will see that this problem has a very rich structure, that is very different from the classical 2-marginal case. First of all we notice that minimizers of this problem are also minimizers of the problem with the cost c(x 1 , . . . , xn) = |x 1 + . . . + xn| 2 ; in fact we have that

R dN N i,j=1 -|x i -x j | 2 dγ = 2 R dN c dγ -(N + 1) N i=1 R d |x| 2 dµ i ,
but this last additive term depends only on the marginals and not on the specific plan γ. The cost c is very particular since it has a wide class of "trivial" optimal plans, that is the ones that are concentrated on x 1 + . . . + xn = 0; however the structure is very rich, see Lemma 4.3.

Concerning the existence of minimizers of (4.1), we will assume that the measures µ i have finite second moments; then existence follows immediately from the equality

argmin γ - N i=1 N j=i+1 |x j -x i | 2 dγ = argmin γ |x 1 + • • • + x N | 2 dγ.
In particular, Corollary 2.15 holds in this case too. Notice that the fact that the repulsive harmonic cost is smooth and has linear gradient does not make the Multi-marginal Optimal Transportation problem easier compared to the Coulomb cost. In fact, in this case the problem is that if we write down the optimality conditions for the potentials, in the case described in Lemma 4.3, we simply find the condition x 1 + . . . + x N = const, since in this case the potentials are linear functions.

However, we can enjoy the symmetries of this problems and build easy Monge solutions for some particular cases of (4.1), see examples 4.9, 4.10, 4.11, 4.12 and 4.14 below.

Before stating the main result in the multi-marginal setting, we start analyzing the problem (4.1) in the 2-marginals case, where everything seems to work fine, just as in the square distance case.

Proposition 4.1. Let µ, ν ∈ P(R d ) and cw(x, y) = -|x -y| 2 be the opposite of the squaredistance in R d . Suppose that µ is an absolutely continuous with respect to the Lebesgue measure in R d and ν has no atoms. Then, there exists a unique optimal transport map T : R d → R d for the problem

min cw(x, y)dγ : γ ∈ Π 2 (R 2d , µ, ν) = inf cw(x, T (x))dµ : T µ = ν
Moreover, T = ∇φ, where φ : R d → R is a concave function and there exists a unique optimal transport plan γ, that is γ = (I d × T ) µ .

Proof. This result is a easy consequence of the Brenier's theorem. Indeed, it is enough to verify that, taking C = 2 |x| 2 dµ + 2 |x| 2 dν, we have

C + inf T µ=ν -|x -T (x)| 2 dµ(x) = inf G µ=ν |x -G(x)| 2 dµ(x)
where G = -T and ν = (-Id) ν. Then by Brenier's theorem there exists an unique optimal map G which can be written as G(x) = ∇ψ(x), and ψ : R d → R is a convex map. In other words, T is a gradient of a concave function.

Notice that if we suppose µ = f (x)dx and ν = g(x)dx are probability measures with densities f, g concentrated, respectively, in convex sets Ω f , Ωg and assume there exits a constant λ > 0 such that λ ≤ f, g ≤ 1/λ. Then, T is a C 1,α function inside Ω f [START_REF] Luis | A localization property of viscosity solutions to the monge-ampere equation and their strict convexity[END_REF][START_REF] Luis | Some regularity properties of solutions of monge ampere equation[END_REF].

Example 4.2 (2 marginals case, uniform measure in the d-dimensional cube). Suppose that µ = ν = L [0,1] d . In this case, we can verify easily that the optimal map T : [0, 1] d → [0, 1] d is given by the anti-monotone map T (x) = (1, . . . , 1) -x.

Surprisingly, the next theorem says that we can not always expect Caffarelli's regularity for the optimal transport maps for the repulsive harmonic cost with finitely many marginals N > 2 even when the support of the measures has convex interior (see corollary 4.8 below). Before stating the theorem we will prove a lemma that characterizes many optimal plans: Lemma 4.3. Let {µ i } N i=1 be probability measures on R d and h : R d → R be a strictly convex function and suppose c :

([0, 1] d ) N → R a cost function of the form c(x 1 , . . . , x N ) = h(x 1 + • • • + x N ).
Then, if there exists a plan γ ∈ Π(µ 1 , . . . , µ N ) concentrated on some hyperplane of the form x 1 + . . . x N = k, this plan is optimal for the multi marginal problem with cost c and supp(γ) ⊂ {x 1 + . . . + x N = k} is a necessary and sufficient condition for γ to be optimal. In this case we will say that γ is a flat optimal plan and {µ i } N i=1 is a flat N -tuple of measures.

Proof. First of all we show that k is fixed by the marginals µ i , and it is in fact the sum of the barycenters of these measures. Let c i = x dµ i ; then let us suppose that there exists γ ∈ Γ(µ 1 , . . . , µ N ) that is concentrated on {x 1 + . . . + x N = k}. Then, using that the i-th marginal of γ is µ i , we can compute

k = (x 1 + . . . + x N ) dγ = N i=1 x i dγ = N i=1 x dµ i = N i=1 c i .
In particular we notice also that for every admissible plan γ we have (x 1 + . . . + x N ) dγ = k and so by Jensen inequality we have

h(x 1 + . . . + x N ) dγ ≥ h (x 1 + . . . + x N ) dγ = h(k) = h(x 1 + . . . + x N ) dγ.
This proves that γ is an optimal plan. Thanks to the strict convexity of h, this shows also that if γ is optimal then γ-a.e. we should have x 1 + . . .

+ x N = k.
This reveals a very large class of minimizers in some cases, as we will see later. However not every N -tuple of measures is flat: it is clear that we can have marginals such that there is no plan with such property: Remark 4.4. Let N = 3 and µ i = µ for every i = 1, 2, 3 with µ = ν 1 + ν 2 with ν 1 = (-Id) ν 2 and ν 1 concentrated on [START_REF] Agueh | Barycenters in the {W}asserstein space[END_REF][START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF]. Now it is clear that the barycenter of µ is 0 but for every 3 points x 1 , x 2 , x 3 in the support of µ we cannot have x 1 + x 2 + x 3 = 0: two of them have the same sign, let's say x 1 and x 2 , but then we have

|x 1 + x 2 | ≥ 4 > 3 ≥ |x 3 |, which contradicts x 1 + x 2 = -x 3 .
So there is no an admissible plan concentrated on H 0 = {x 1 + x 2 + x 3 = 0}; in fact we showed that for every γ ∈ Π(µ) we have supp(γ) ∩ H 0 = ∅. Remark 4.5. In the case N = 2 the condition for which there exists an admissible plan γ concentrated on an hyperplane of the form x 1 + x 2 = k implies that the two measures µ and ν = (-Id) ν are equal up to translation. This condition is in fact very restrictive. However when we think at the DFT problem and in particular we assume that µ = ν then the said condition amounts to have that µ is centrally symmetric about its barycenter and in the context of density of 2 electrons around a nucleus this seems a fairly natural condition.

Theorem 4.6. Let µ i = µ = L d [0,1] d , ∀i = 1, . . . , N be the uniform measure on the d- dimensional cube [0, 1] d ⊂ R d , h : R d → R a convex function and suppose c : ([0, 1] d ) N → R a cost function such that c(x 1 , . . . , x N ) = h(x 1 + • • • + x N ).
Then, there exits a transport map

T : [0, 1] d → [0, 1] d such that T N (x) = x and min γ∈Π N (µ) cdγ = min T µ=µ, T (N ) =I c(x, T (x), . . . , T (N -1) (x))dµ
Moreover, T is not differentiable at any point and it is a fractal map, meaning that it is the unique fixed point of a "self-similar" linear transformation F acting on

L ∞ ([0, 1] d ; [0, 1] d ).
Proof. We express every z ∈ [0, 1] by its base

-N system, z = ∞ k=1 a k N k with a k ∈ {0, 1 . . . , N -1}. Consider the map given by S(z) = ∞ k=1 S(a k )
N k , where S is the permutation of N symbols such that S(i) = i + 1, ∀{i = 0, 1, . . . , N -2} and S(N -1) = 0. A straightforward computation shows that

z + N i=1 S i (z) = N 2 (4.2)
Let T : R d → R d be a map defined by T (x) = T (z 1 , . . . , z d ) = (S(z 1 ), . . . , S(z d )) and denote by T (j) (x) = (S (j) (z 1 ), . . . , S (j) (z d )), j = 1, . . . , N -1. We will first show that S is a measurepreserving map. In fact, we can show that there exist functions S k : [0, 1] → [0, 1] defined recursively by S 0 (x) = x, and S k+1 (x) = T k+1 (S k (x))

where T k acts only of the k-th digit:

T k (x) = x -(N -1) • N -k if x ∈ C k and T k (x) = x + N -k if x ∈ [0, 1] \ C k . The sets C k are defined by C k = N k-1 j=1 j N k-1 - 1 N k , j N k-1 .
Moreover, it is easy to see that

(T k ) L [0,1] = L [0,1] ∀k ∈ N,
and S k → S uniformly.

Hence

f (x)S L d = f (x)dx, ∀f ∈ C 0 ([0, 1]).
Now, it remains to show that T is optimal. But this is true thanks to the fact that (4.2) implies that the plan induced by T satisfies the hypothesis of Lemma 4.3.

It is clear that we can reduce to prove the non-differentiability and the fractal properties only for the map S. The non-differentiability comes from the fact that for each z ∈ [0, 1], we have that the base-

N representation z = ∞ k=1 a k
N k is such that a k can't be definitively N -1. Now it is sufficient to choose those k j such that a k j = N -1 and consider the numbers z j = z + (N -1 -a k j )/N k j and z j = z ± 1/N k j (depending on whether a k j = 0 or not) which have the same digits as z apart from the j k -th digit. Then it is straightforward to see that S(z j ) -S(z) = -(a k j + 1)/N k j while S(z j ) -S(z) = z j -z; in particular we have

S(z j )-S(z) z j -z = 1 while S(z j )-S(z) z j -z
≤ -1 N and so, letting j → ∞, we get that S is not differentiable at the point z.

As for the fractal property: we consider the transformation

F (g)(x) = 1 N g(N x -i) + i+1 N for i N ≤ x < i+1 N , with i = 0, . . . , N -2 1 N g(N x -N -1) for N -1 N ≤ x < 1.
In order to see what the construction is doing we imagine to divide [0, 1]2 in a grid N × N of squares and then we are putting scaled copies (by a factor N ) of the original function in the above diagonal squares and in the rightmost bottom one. Now it is easy to see that S is the unique 2 fixed point of F , and this gives the property of self-similarity. Remark 4.7. Another proof of Theorem 4.6 can be done noticing that the transformation F (the one defined in the proof) acts as a 1/N -contraction in the space of measure preserving L ∞bijections from [0, 1) to [0, 1).

Corollary 4.8. Let µ = L d [0,1] d be the uniform measure on the d-dimensional cube [0, 1] d in R d and suppose c : ([0, 1] d ) N → R the N -dimensional repulsive harmonic cost c(x 1 , . . . , x N ) = - N i=1 N j=i+1 |x j -x i | 2 , (x 1 , . . . , x N ) ∈ ([0, 1] d ) N
Then, there exits an optimal cyclical transport map T :

[0, 1] d → [0, 1] d : in particular min γ∈Π N (µ) cdγ = min T µ=µ, T (N ) =I c(x, T (x), T (2) (x), . . . , T (N -1) (x))dµ.
Moreover, T is not differentiable at any point.

Proof. As we already observed, the problem with the cost c is equivalent to the problem with the cost |x 1 + . . . + x N | 2 and the result follows from the Theorem 4.6.

We remark the construction of T in the proof of the Theorem 4.6 also works if N = 2 and T is exactly the optimal transport map described in the example (4.2).

The unexpected aspect of Corollary 4.8 is the existence -for N ≥ 3 -of an optimal transport map which is not differentiable almost everywhere. It turns out that this optimal map T could be not unique if d > 1 and N ≥ 3 and, in that case, we can construct explicitly a regular optimal map. This kind of richness appears already in the case d = 1 with the Coulomb cost; however in that case we have uniqueness if we restrict ourselves to the symmetric optimal plan. Here this is not the case as it is easy to see that modifying the action S on the digits (the important thing is that when we see S as a permutation, it is a cycle), we obtain another map, and the symmetrized plan is not equal to the one generated by the map described before.

In the following, we are going to present some concrete examples where we have other explicit solutions for such kind of optimal maps for (MK weak ) in (4.1). For these examples the goal is to show that there can be smooth optimal maps but they don't necessarily satisfy the "group rule", that is, we can find maps T 2 , . . . T N such that (Id, T 1 , T 2 , . . . , T N -1 ) µ 1 is an optimal plan but T 1 • T 1 = T i for any index i ∈ {1, . . . , N -1}.

Example 4.9 (3 particles, asymmetric). Consider the case when three particles are distributed in R 3 as Gaussians,

µ 1 = µ 2 = 1 (2π) 3/2 exp(-1 2 (x 2 1 + x 2 2 + x 2 3 )) and µ 3 = 1 2π 3/2 exp(-(x 2 1 + x 2 2 + x 2 3
)). In this case, we can verify that the couple (T, S) of maps T, S : R d → R d , T (x) = x a.e. and S(x) = -2x a.e. is admissible and it is optimal since x + T (x) + S(x) = 0.

Example 4.10 (2N particles on S 1 , ). Suppose µ 1 , . . . , µ 2N uniform probability measures on the circle S 1 . The rotation map R θ : R 2 → R 2 with angle θ = π/N is an optimal transport map. Also, the maps R kθ , k = 2, . . . , N , are optimal transport maps for the repulsive harmonic cost.

Example 4.11 (2N particles, breathing map). Suppose µ 1 , . . . , µ 2N uniform probability measures on S 2 ⊂ R 3 . Consider the vector v, A : R 3 → R 3 the antipodal map A(x) = -x and R v θ : R 3 → R 3 the rotation of angle θ = π/N and direction v. Then, T = A • R v θ is a cyclic optimal transport map.

The optimal solution γ = (x, T, T (2) , . . . , T (N -1) ) µ 1 is called "breathing" solution [START_REF] Seidl | Strictly correlated electrons in densityfunctional theory: A general formulation with applications to spherical densities[END_REF]: the coupling γ represent the configuration where the 2N electrons are always at the same distance from the center, opposite to each other in the equilibrium configuration. Notice that, the map G = T 2 is also an cyclical optimal transportation map. Moreover, T is a C ∞ function and it is a gradient of the convex-concave function φ(x 1 , x 2 , . . . ,

x 2N -1 , x 2N ) = 1 2 (x 2 1 + x 2 3 + • • • + x 2 2N -1 ) - 1 2 (x 2 2 + x 2 4 + • • • + x 2 2N ).
Example 4.12 (2N particles in R d , symmetric ρ). In this case, on can consider the maps T (x) = x a.e. and S(x) = -x a.e. which are such that c(T (x), S(x), . . . , T (x), S(x)) = 0.

We notice that in this case we have S (2) (x) = x and so in particular this solution has the cyclic structure (T (x), S(x), . . . , T (x), S(x)) = (x, S(x), S (2) (x), . . . , S (2N -1) ).

Example 4.13 (Optimal Maps which doesn't satisfy a group law). Let us consider the case d = 1 and N = 3 where the measures are

µ 1 = µ 2 = µ 3 = 1 2 L| [-1,1]
. Then we define the maps

T (x) = x + 1 if x ≤ 0 x -1 if x > 0, S(x) = -1 -2x if x ≤ 0 1 -2x if x > 0.
We have that T µ = µ and S µ = µ, and moreover x + T (x) + S(x) = 0; in particular (Id, S, T ) µ is a flat optimal plan and so the thesis. , where H is defined as H = {x + y + z = 0} ∩ {|x| ≤ 1, |y| ≤ 1, |z| ≤ 1} and f is a function to be chosen later. This is a symmetric flat optimal plan; now we compute the marginals. Since it is clear that γ is invariant under x → -x, it is sufficient to consider the marginal on the set x > 0. But then we can make the computation

+ √ 3 1 0 -x -1 φ(x)f (|y|) dy dx = 1 0 √ 3φ(x) xf (x) + 2 1 x f (t) dt dx.
In particular the choice f (x) = √ 3 6 x gives the marginals equal to 1 2 L| [-1 ,1] . We notice also that any even density ρ = h(|x|) for some decreasing function h : [0, 1] → [0, ∞) can be represented in this way: in fact it is sufficient to choose f

(x) = 1 √ 3 h(x) x -2x 1 x h(t) t 3 dt .
Example 4.15 (A Counterexample on the uniqueness N > 3). As mentioned in theorem 2.16 and in [START_REF] Friesecke | N-density representability and the optimal transport limit of the hohenberg-kohn functional[END_REF], it was already understood by Pass that in these high dimensional cases, the solution of repulsive costs may also be non unique, as opposite to the two marginals case. On a higher dimensional surface there can be enough wiggle room to construct more than one measure with common marginals, as shown in the examples 4.9-4.14. In most of the cases, the non-uniqueness seems to be given by the symmetries of the problem, but in Example 4.14 and Corollary 4.8 this is not the case, as we exploit the fact that the dimension of the set c(x+y+z)-φ(x)-φ(y)-φ(z) = 0, is greater than the minimal one.

Finally, the last proposition of this section states that when d = 1 and odd N we have no hope in general to find piecewise regular cyclic optimal transport maps. is not attained by a map T which is differentiable almost everywhere.

Proof. First of all we notice that if T is differentiable almost everywhere then also T (2) has the same property, thanks to the fact that T µ = µ. In particular, since the Lusin property holds true for T (i) for every i = 1, . . . , N and T µ = µ, the change of variable formula holds and in particular we have that (T (i) ) (x) = ±1 for almost every x (notice also that T is bijective almost everywhere since T (N ) (x) = x). Since µ is N -flat we have that the condition on T in order to be an optimal cyclical map is x + T (x) + T (2) (x) + . . . + T (N -1) (x) = N/2; now we can differentiate this identity and so we will get 1 + T (x) + (T (2) ) (x) + . . . + (T (N -1) ) (x) = 0 for a.e. x.

But this is absurd since on the left hand side we have an odd number of ±1 and their sum will be an odd number.

In conclusion, also in the case of the repulsive harmonic cost, the picture is far from being clear: an interesting structure appears when {µ i } N i=1 is a flat N -tuple of measures but we still can't characterize this property. Moreover, in the flat case in which µ i = µ, for example when

µ = 1 2 L 1 | [-1,1]
, we have both diffuse optimal plan and a cyclical optimal map. An interesting open problem is whether for any N -flat measure µ, say absolutely continuous with respect to the Lebesgue measure, we have a cyclical optimal map and a diffuse plan.

Multi-marginal OT for the Determinant

We are going to give a short overview of the main results in [START_REF] Carlier | Optimal transportation for the determinant[END_REF], where Carlier and Nazaret consider the following optimal transport problems for the determinant:

(MK Det ) sup γ∈Π((R d ) d ,
d i=1 R d |x i | p i p i dµ(x i ) < +∞.
Notice that for this particular cost the problem (2.1) makes sense only when N = d. In the following, we will focus on problem (5.1) and exhibit explicit minimizers γ in the radial case. Clearly, the difference between (5.1) and (5.2) is that the second one admits positively and negatively oriented basis of vectors, while the first one "chooses" only the positive ones. Moreover, if we assume that among marginals µ 1 , . . . , µ d there exist two symmetric probability measures µ i , µ j , i = j, i.e. µ i = (-Id) µ i and µ j = (-Id) µ j , then any solution γ of (5.1) satisfies det(x 1 , . . . , x d ) ≥ 0 γ-almost everywhere and so solves also (5.2) (Proposition 6, [START_REF] Carlier | Optimal transportation for the determinant[END_REF]).

Similarly to the Gangbo-Świȩch cost [START_REF] Gangbo | Optimal maps for the multidimensional mongekantorovich problem[END_REF], the Monge-Kantorovich problem for the determinant (5.1) can be seen as a natural extension of classical optimal transport problem with two marginals and so, it is equivalent to the 2-marginals repulsive harmonic cost (4.1).

Indeed, we can write in the two marginals case, det(x 1 , x 2 ) = x 1 , Rx 2 , where R : R 2 → R 2 is the rotation of angle -π/2. Hence, since µ 1 and µ 2 have finite second moments, up to a change of variable x2 = Rx 2 , the problem (MK Det ) in (5.1) is equivalent to the classical Brenier's optimal transportation problem:

argmax γ∈Π(µ 1 ,µ 2 ) R 2 det(x 1 , x 2 )dγ(x 1 , x 2 ) = argmax γ∈Π(µ 1 ,µ 2 ) R 2 x 1 , x2 dγ(x 1 , x 2 ) = argmax γ∈Π(µ 1 , μ2 ) R 2 x 1 , x 2 dγ = argmin γ∈Π(µ 1 , μ2 ) R 2 |x 1 -x 2 | 2 2 dγ(x 1 , x 2 ) -C where C = 1/2( |x 1 | 2 dµ 1 + |x 2 | 2 dµ 2 ) and μ2 = R µ 2 .
In the sequel, we are going to construct maximizers for (5.1), thanks to some properties of the Kantorovich potentials of the dual problem associated to (5.1) (see theorem 5.1 bellow),

(K Det N ) inf R d d i=1 u i (x i )dµ i (x i ) : det(x 1 , . . . , x d ) ≤ N i=1 u i (x i ) . (5.3)
In [START_REF] Carlier | Optimal transportation for the determinant[END_REF], the authors provide a useful characterization of optimal transport plans through the potentials u i , given by theorem 5.1 . In addition, by means of a standard convexification trick we obtain regularity results on the Kantorovich potentials.

Theorem 5.1. A coupling γ ∈ Π((R d ) d , µ 1 , . . . , µ d ) is optimal in (5.1) if and only if there exists lower semi-continuous convex functions u i : R d → R ∪ {∞} such that for all i ∈ {1, . . . , d},

d i=1 u i (x i ) ≤ d i=1 u * i ((-1) i+1 i =j x j ), on (R d ) d ; d i=1 u i (x i ) = d i=1 u * i ((-1) i+1 i =j
x j ), γ -almost everywhere;

(-1) i+1 i =j

x j ∈ ∂u i (x i ), γ -almost everywhere.

where d i=1 x j denotes the wedge product and, for every i, u * i is the convex dual of the Kantorovich potential u i . Now, the main idea is to use the geometrical constraints on the Kantorovich potentials u i (5.3), given by the theorem (5.1), in order to construct an explicit solution.

We illustrate the theorem (5.1) and explain how to construct a particular optimal γ for (5.1) by an example in the three marginals case. Let µ i = ρ i L 3 , i = 1, 2, 3 radially symmetric probability measures on the 3-dimensional ball B.

In this particular situation, the optimizers of (5.1) and (5.2) have a natural geometric interpretation: what is the best way to place three random vectors x, y, z, distributed by probability measures µ 1 , µ 2 , µ 3 on the sphere, such that the simplex generated by those three vectors (x, y, z) has maximum average volume? Suppose γ ∈ Π(B, µ 1 , µ 2 , µ 3 ) optimal in (5.1) when d = 3. From optimality of γ, we have

u 1 (x) + u 2 (y) + u 3 (z) = det(x, y, z), γ -almost everywhere.
Applying the theorem (5.1), we get

u 2 (y) + u 3 (z) = u * 1 (y ∧ z) u 1 (x) + u 3 (z) = u * 2 (-x ∧ z) u 1 (x) + u 2 (y) = u * 3 (x ∧ y)
, γ -almost everywhere, and,

∇u 1 (x) = y ∧ z ∇u 2 (y) = -x ∧ z ∇u 3 (z) = x ∧ y
, γ -almost everywhere.

(5.4)

It follows from (5.4), given a vector x on in the ball, the conditional probability of y given x is supported in a "meridian" M (x)

M (x) = {y ∈ S 2 : ∇u 1 (x), y = 0},
where S 2 is the 2-sphere. Finally, assuming that x, ∇u 1 (x) = 0, the conditional probability of z given the pair (x, y) is simply given by a delta function on z z = ∇u 1 (x) ∧ ∇u 2 (y)

x, ∇u 1 (x) ;

In particular, we have

x, ∇u 1 (x) + y, ∇u 2 (y) + z, ∇u 3 (z) = det(x, y, z) = det(∇u 1 (x), ∇u 2 (y), ∇u 3 (z)).

Example 5.2 (An explicit solution in the ball B ⊂ R 3 ). Suppose µ i = L 3 B , i = 1, 2, 3, the 3-dimensional Lebesgue measure in the ball B ⊂ R 3 . The following coupling γ *

B 3 f dγ * = 1 L 3 (B) B M (x) f (x, |x|y, x ∧ y) dH 1 (y) 2π dx, ∀f ∈ C(B 3 , R). (5.5)
is an optimizer for (5.1) with d = 3. Indeed, from we can show explicit potentials u

* 1 (x) = u * 2 (x) = u * 3 (x) = |x| 3 /3; clearly we have det(x, y, z) ≤ |x||y||z| ≤ |x| 3 /3 + |y| 3 /3 + |z| 3 /3, ∀ (x, y, z) ∈ B,
with equality when |x| = |y| = |z| and x, y, z are orthogonal. Since γ * is concentrated on this kind of triples of vector we have the optimality. Finally, by a suitable change of variable, it is easy to see that γ * ∈ Π 3 (L 3 B ). Some comments on the radially symmetric d-marginals case: In [START_REF] Carlier | Optimal transportation for the determinant[END_REF], for d radially symmetric probability measures, the authors exhibit explicit optimal couplings γ * . In their proof, two aspects were crucial: the first one is remark that if {µ i } d i=1 are radially symmetric measures in R d , then the optimal Kantorovich potentials u i (x i ) = u i (|x i |) are also radially symmetric; in particular the system (5.4) for general d implies that the support of γ * is contained in the set of an orthogonal basis. The second observation is to notice that in the support of γ * we have

H i (|x 1 |) = |x i |, where H i is the unique monotone increasing map such that (H * i ) µ 1 = µ i , where H * i (x) = x |x| H i (|x|)
. This is done analyzing the correspondent radial problem (with cost c(r 1 , . . . , r d ) = r 1 • • • r d ), using the optimality condition φ i (r i ) = ∂ i c and the fact that in this case r i φ i (r i ) = r 1 φ 1 (r 1 ) = c ≥ 0; then, exploiting the convexity of φ i we get r i = H i (r 1 ) for some increasing function H i , that is uniquely determined.

Existence of Monge-type solutions: In the 3-marginals case in the unit ball, by construction of the the coupling γ * in (5.5) or, more generally, the optimal coupling in the d-marginal case (see theorem 4 in [START_REF] Carlier | Optimal transportation for the determinant[END_REF]), we can see that their support are not concentrated in the graph of cyclic maps T, T 2 , . . . , T d-1 or simply on the graph of maps T 1 , . . . , T d-1 as we could expect from corollary (2.15). In other words, γ * in (5.5) is not Monge-type solution.

The existence of Monge type solutions for the determinant cost is still an open problem for odd number of marginals. From the geometric conditions we discussed above, in the case in which µ i = µ a radial measure, if Monge solutions exists then, for every x ∈ R d , (x, T 1 (x), . . . , T d-1 (x)) should be an orthogonal basis, and

|T i (x)| = |x| , i = 1, . . . , d -1.
For the interesting even dimensional case, we can observe a similar phenomena remarked in the repulsive harmonic costs, concerning the existence of trivial even dimensional solutions for the Monge problem in (5.1). We expect the existence of non-regular optimal transport map also to this case.

Example 5.3 (Carlier & Nazaret, [START_REF] Carlier | Optimal transportation for the determinant[END_REF]). The even dimensional phenomena: As in the repulsive harmonic cost, it is easy to construct Monge minimizers for the determinant cost for even number of marginals ≥ 4. For instance, suppose c(

x 1 , x 2 , x 3 , x 4 ) = det(x 1 , x 2 , x 3 , x 4 ), define transport maps T 1 , T 2 , T 3 : B → R 4 by, for x = (x 1 , x 2 , x 3 , x 4 ) ∈ B T 1 (x) =     -x 2 x 1 -x 4 x 3     , T 2 (x) =     -x 3 x 4 x 1 -x 2     , T 3 (x) =     -x 4 -x 3 x 2 x 1     .
We can see that γ T = (Id, T 1 , T 2 , T 3 ) γ is a Monge-type optimal transport plan for (5.1) and (5.2).

Numerics

Numerics for the multi-marginal problems have so far not been extensively developed. Discretizing the multi-marginal problem leads to a linear program where the number of constraints grows exponentially in the number of marginals. In [START_REF] Carlier | Numerical methods for matching for teams and wasserstein barycenters[END_REF] Carlier, Oberman and Oudet studied the matching of teams problem and they were able to reformulate the problem as a linear program whose number of constraints grows only linearly in the number of marginals. More rencently a numerical method based on an entropic regularization has been developed and it has been applied to various optimal transport problem in [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF][START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. Let us mention that the Coulomb cost has been treated numerically in various paper as in [START_REF] Malet | Strong Correlation in Kohn-Sham Density Functional Theory[END_REF][START_REF] Gori-Giorgi | Density-functional theory for strongly interacting electrons[END_REF], where the solutions are based on the analytical form of transport maps given in [START_REF] Gori-Giorgi | Density-functional theory for strongly interacting electrons[END_REF], in [START_REF] Mendl | Kantorovich dual solution for strictly correlated electrons in atoms and molecules[END_REF], where parameterized functional forms of the Kantorovich potential are used, in [START_REF] Chen | Numerical methods for a kohn-sham density functional model based on optimal transport[END_REF][START_REF] Chen | Pair densities in density functional theory[END_REF], where a linear programming approach has been developed to solve the 1-dimensional problem, and, as already mentioned, in [START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF]. In this section we focus on the regularized method proposed in [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF][START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] and we, finally, present some numerical experiments for the costs studied above.

γ∈Π(R dN ,µ 1 ,...,µ N ) R dN c(x 1 , . . . , x N )γ(x 1 , . . . , x N )dx 1 , . . . , dx N , (6.1)
where N is the number of marginals µ i which are probability distributions over R d , c(x 1 , ..., x N ) is the cost function, γ the coupling, is the probability distribution over R dN Remark 6.1. From now on the marginals µ i and the coupling γ are densities and when the optimal coupling γ is induced by maps T i , we write it as γ = µ

1 (x 1 )δ(x 2 -T 2 (x 1 )) • • • δ(x N -T N (x 1 )).
In order to discretize (6.1), we use a discretisation with M d points of the support of the kth marginal as

{x j k } j k =1,••• ,M d . If the densities µ k are approximated by µ k ≈ j k µ j k δx j k , we get min γ∈ j 1 ,•••j N c j 1 ,••• ,j N γ j 1 ,••• ,j N , (6.2)
where Π k is the discretization of Π, c j 1 ,••• ,j N = c(x j 1 , • • • , x j N ) and the coupling support for each coordinate is restricted to the points {x

j k } j k =1,••• ,M d thus becoming a (M d ) N matrix again denoted γ with elements γ j 1 ,••• ,j N . The marginal constraints C i (such that Π k = N k=1 C k ) becomes C k := γ ∈ R (M d ) N + : j 1 ,...,j k-1 ,j k+1 ,...,j N γ j 1 ,...,j N = µ j k , ∀j k = 1, • • • , M d . (6.3)
As in the continous framework the problem (6.1) admits a dual formulation max

u j k N k=1 M j=1 u j k µ j k s.t. N k=1 u j k ≤ c j 1 •••j N ∀ j k = 1, • • • , M d , (6.4)
where u j k = u k (x j k ) is the kth Kantorovich potential. One can notice that the primal (6.1) has (M d ) N unknown and M d × N linear constraints and the dual problem (6.4) has M d × N unknown, but (M d ) N constraints. This actually makes the problems computationally unsolvable with standard linear programming methods even for small cases.

Remark 6.2. We underline that in many applications we have presented (as in DFT) the marginals µ j k are equal (µ j k = µ j , ∀ k ∈ {1, • • • , N }). Thus the dual problem can be re-written in a more convenient way (but still computationally unfeasible for many marginals) max

u j M j=1 N u j µ j s.t. N k=1 u j k ≤ c j 1 •••j N ∀j k = 1, • • • , M d , (6.5)
where u j = u j k = u(x j k ). Now the dual problem has M d unknown, but (M d ) N linear constraints.

A different approach consists in computing the problem (6.1) regularized by the entropy of the joint coupling. This regularization dates to E. Schrödinger [START_REF] Schrödinger | Über die umkehrung der naturgesetze[END_REF] and, as mentioned above, it has been recently introduced in many applications involving optimal transport [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF][START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF][START_REF] Galichon | Matching with trade-offs: Revealed preferences over competing characteristics[END_REF]. Thus, we consider the following discrete regularized problem where E(γ) is defined as follows

E(γ) = j 1 ,•••j N γ j 1 ,••• ,j N log(γ j 1 ,••• ,j N ) if γ ≥ 0 +∞ otherwise, (6.7)
and C is the intersection of the set associated to the marginal constraints (we remark that the entropy is a penalization of the non-negative contraint on γ). After elementary computations, we can re-write the problem as min As explained in section 2, when the transport plan γ is concentrated on the graph of a transport map which solves the Monge problem, after discretisation of the densities, this property is lost along but we still expect the matrix γ to be sparse. The entropic regularization spreads the support and this helps to stabilize the computation as it defines a strongly convex program with a unique solution γ . Moreover the solution γ can be obtained through elementary operations. The regularized solutions γ then converge to γ (see figure 2), the solution of (6.1) with minimal entropy, as → 0 (see [START_REF] Cominetti | Asymptotic analysis of the exponential penalty trajectory in linear programming[END_REF] for a detailed asymptotic analysis and the proof of exponential convergence). The aim of the IPFP is to find the KL projection of γ on the set

C = {γ ij ∈ R M d × R M d : j γ ij = µ i } ∩ {γ ij ∈ R M d × R M d : i γ ij = ν j }.
By writing down the lagrangian associated to (6.10) and computing the optimality condition, we find that γ ij can be written as γ ij = a i b j γij with a i = e u i / , b j = e v j / , (6.11) where u i and v j are the regularized Kantorovich potential. Then, a i and b j can be uniquely determined by the marginal constraint Remark 6.3. In [START_REF] Ruschendorf | Convergence of the iterative proportional fitting procedure[END_REF] Rüschendorf proves that the iterative method (6.13) converges to the KLprojection of γ on C. Remark 6.4. Rüschendorf and Thomsen (see [START_REF] Rüschendorf | Closedness of sum spaces andthe generalized schrödinger problem[END_REF]) proved, in the continous measure framework, that a unique KL-projection exists and takes the form γ(x, y) = a(x) ⊗ b(y)γ(x, y) (where a(x) and b(y) are non-negative functions).

a i = µ i j b j γij , b j = ν j i a i γij . ( 6
The extension to the multi-marginal framework is straightforward but cumbersone to write. It leads to a problem set on N M d -dimensional vectors a j,i (•) , j = 1, • • • , N, i (•) = 1, • • • , M d . Each update takes the form a n+1 j,i j = ρ i j i 1 ,i 2 ,...i j-1 ,i j+1 ,...,i N γi 1 ,...,i N a n+1 1,i 1 a n+1 2,i 2 ...a n+1 j-1,i j-1 a n j+1,i j+1 ...a n N,i N . optimal coupling γ12 (x, y). Bottom-Left: support of the optimal coupling γ13 (x, z). Bottom-Right: optimal coupling γ13 (x, z).

is simple to show that the second order Γ-limit of the regularized problem as → 0 is the relative entropy of γ with respect to H 2 | {x+y+z=0} ; we expect that the Γ-limit is the same also when we add the marginal constraint.

In figure 7 we present the projection γ12 (x, y) = (e 1 , e 2 )γ(x, y, z) of the computed coupling (it is enough to visualize only this projection because of the symmetries of the problem) and, as expected, we observe a diffuse plan. . The simulation has been performed on a discretization of [-1, 1] with M d = 1000 gridpoints, = 0.0005. 6.4. Numerical experiments: Determinant cost. Numerical simulations for the multi-marginal problem with the determinant cost present an obvious computational difficulty: in order to compute the solution just for the 3-marginals case, we have to introduce a discretization of R 3 (as the cost function now is c(x 1 , • • • , x N ) = det(x 1 , • • • , x N ) with x i ∈ R N ). However if we take radially symmetric densities µ i as marginals (see section 5 for a more complete description of the problem in the measure framework), Carlier and Nazaret show that (MK) can be reduced to a 1-dimensional problem: the only unknowns are the relations between the norms r i = x i of each vector. These relations can be obtained by solving the following problem min where λ i (r) = σ(r)µ i (r) (e.g. N = 3 then r = (r, θ, ϕ) and σ(r) = 4πr 2 ). Moreover, for this problem we know that there exists a unique (deterministic) optimal coupling (see Proposition 5 in [START_REF] Carlier | Optimal transportation for the determinant[END_REF]). The discretized cost now reads c j 1 •••j N = N k=1 r j k . Let us consider the 3-marginals case and all densities λ i equal to the uniform density on the ball B 0.5 = {x ∈ R 3 | x ≤ 0.5}. As proved in [START_REF] Carlier | Optimal transportation for the determinant[END_REF], the optimal coupling is actually supported by the graph of maps which rearrange measure (see figure 8 ). The simulation in figure 8 has been performed on a discretization of [0, 0.5] with M d = 100 and = 0.01. In the same way we can take all marginals λ i (r) = 4πr 2 e -4r We, finally, present a simulation (see figure 11) with three different marginals λ i (see figure 10). As one can observe we obtain that the projection of the optimal coupling are concentrated on the graph of a map which rearranges the densities in a monotone way. The simulation has been performed on a discretization of [0, 4] with M d = 300 and = 0.07. ). Top-Center: optimal coupling γ13 (r 1 , r 3 ). Top-Right: optimal coupling γ23 (r 2 , r 3 ). Bottom: supports of those couplings.

Conclusion

In this survey we presented a general famework for the multi-marginal optimal transportation with repulsive cost. We also gave some new results which highlight the importance of the multimarginal approach. As mentioned there are still a lot of open questions about the solutions for the problems presented and we hope to treat them in future works.

2. 1 .

 1 Couplings and Multi-marginal Optimal Transportation problem.

  3) on the Polish spaces X = R d , Y = (R d ) N -1 with measures µ = ρL d and ν = (e 2 , . . . , en) γ, for every γ ∈ Π(R dN , ρL d ), and obtain the following corollary Corollary 2.4. Let µ 1 = • • • = µ N = ρL d be probability measures on R d with density ρ and c : (R d ) N → R be a continuous cost function. Then, min(MK) = inf(M N ).

Theorem 2 .

 2 14 (M. Colombo and S. Di Marino,[START_REF] Colombo | Equality between monge and kantorovich multimarginal problems with coulomb cost[END_REF]). Let c : (R d ) N → R be a continuous function and bounded from below. If µ has no atoms, then(MKsym) = (M cycl ).Finally, we have the equality between the Multi-marginal cyclic problem and Multi-marginal OT problem. Corollary 2.15. If µ has no atoms and c : (R d ) N → R is a continuous function and bounded from below, then inf(MKsym) = inf(MK) = inf(M N ) = inf(M cycl ).

Lemma 3 . 8 (

 38 Radial case). Let µ ∈ P(R d ) be a radially symmetric measure. In particular µ is determined by µr = | • | µ. Then for every optimal plan γ ∈ Π N (µ), if we consider γr = R γ, where R : (x 1 , . . . , x N ) → (|x 1 |, . . . , |x N |), we have that γr is an optimal plan for the 1D multi marginal problem min η∈Π N (µr ) R N c 1 (r 1 , . . . , r N ) dη. (3.1)

Example 4 . 14 (

 414 Optimal Transport diffuse plan). Let us consider the same problem as in Example 4.13. Now we consider a general symmetric plan γ = 1 2 H 2 | H f (max{|x|, |y|, |z|})

  x≥0 φ(x) dγ(x, y, z) = |x|,|y|,|x+y|≤1, x≥0 √ 3φ(x)f (max{|x|, |y|, |x + y|}) dxdy )f (x) dy dx

Proposition 4 . 16 .

 416 Let µ = L| [0,1] and N ≥ 3 be an odd number. Then, the infimum inf c(x, T (x), T (2) (x), . . . , T (N -1) )dµ : T µ = µ, T N = I .

6. 1 .

 1 The regularized problem and the Iterative Proportional Fitting Procedure. Let us consider the problem (MK) inf

min γ∈C j 1 ,

 1 •••j N c j 1 ,••• ,j N γ j 1 ,••• ,j N + E(γ) (6.6)

- c j 1 ,

 1 γ|γ) = i 1 ,...,i N γ i 1 ,...,i N log( γ i 1 ,...,i N γi 1 ,...,iN) is the Kullback-Leibler distance and γi 1 ,...,i N = e

Figure 2 .

 2 Figure 2. Support of the coupling γ for the Coulomb cost and µ 1 = µ 2 = (1 + cos(πx/2))/2. The simulation has been performed on a discretization of [-2, 2] with M d = 1000. In order to introduce the Iterative Proportional Fitting Procedure (IPFP), we consider the two marginals problem min γ∈C i,j c ij γ ij . (6.10)

Figure 3 .

 3 Figure 3. On the Left: support of the optimal coupling γ12 (x, y), on the Right: graph of the optimal coupling γ12 (x, y) , in the cases N = 3 (1st Row), N = 4 (2nd Row), N = 5 (3rd Row). The dashed lines delimit the supports of μi , with i = 1, • • • , N .

Figure 4 .

 4 Figure 4. (logarithmic cost) Left: support of the optimal coupling γ12 (x, y) for N = 3. Right: support of the optimal coupling γ12 (x, y) for N = 3. The dashed lines delimit the intervals where μi , with i = 1, • • • , 3, are defined.

Figure 5 .

 5 Figure 5. Left: support of the optimal coupling. Right: optimal coupling.

Figure 6 .

 6 Figure 6. Top-Left: support of the optimal coupling γ12 (x, y). Top-Right: optimal coupling γ12 (x, y). Bottom-Left: support of the optimal coupling γ13 (x, z). Bottom-Right: optimal coupling γ13 (x, z).

γ∈Π(λ 1 Figure 7 .

 17 Figure 7. Left: support of the optimal coupling γ12 (x, y). Right: optimal coupling γ12 (x, y).

Figure 8 .

 8 Figure 8. (Uniform density) Left: support of the optimal coupling γ12 (r 1 , r 2 ).Right: optimal coupling γ12 (r 1 , r 2 ).

Figure 9 .

 9 Figure 9. (Exponential density)Left: support of the optimal coupling γ12 (r 1 , r 2 ). Right: optimal coupling γ12 (r 1 , r 2 ).

λ 1 λ 2 λ 3 Figure 10 .

 310 Figure 10. Densities λ 1 , λ 2 and λ 3 .

Figure 11 .

 11 Figure 11. Top-Left: optimal coupling γ12 (r 1 , r 2 ). Top-Center: optimal coupling γ13 (r 1 , r 3 ). Top-Right: optimal coupling γ23 (r 2 , r 3 ). Bottom: supports of those couplings.

  µ 1 ,...,µ d ) (R d ) d det(x 1 , . . . , x d )dγ(x 1 , . . . , x d ) , . . . , µ d are the uniform probability measures in R d . In addition, in order to guarantee existence of a solution, we assume that there exist p 1 , . . . , p d ∈ [1, ∞[ such that

	and
	(MK |Det| )	sup γ∈Π((R d ) d ,µ 1 ,...,µ d ) (R d ) d	| det(x 1 , . . . , x d )|dγ(x 1 , . . . , x d ),	(5.2)
	where µ 1 d i=1	1 p i	= 1, and
					(5.1)

We notice that F is a contraction:F (g) -F (g ) ∞ ≤ 1 N g -g ∞

Acknowledgements

Acknowledgements

We would like to thank Paola Gori-Giorgi for many helpful and stimulating discussions as well as Jean-David Benamou and Guillaume Carlier for a critical reading of the manuscript.

Simone Di Marino and Luca Nenna gratefully acknowledge the support of the ANR, through the project ISOTACE (ANR-12-MONU-0013). Luca Nenna also ackwoledges INRIA through the "action exploratoire" MOKAPLAN.

Augusto Gerolin is member of the Gruppo Nazionale per l'Analisi Matematica, la Probabilità e le loro Applicazioni (GNAMPA) of the Instituto Nazionale di Alta Matematica (INdAM). He also acknowledges and thanks CNPq for the financial support of his Ph.D, when this work was done.

Example 6.5 (IPFP and 3 marginals). In order to clarify the extension of the IPFP to the multimarginal case , we consider 3 marginals and we write down the updates of the algorithm

In the following sections we present some numerical results obtained by using the IPFP.

6.2.

Numerical experiments: Coulomb cost. We present now some results for the multimarginal problem with Coulomb cost in the real line. We consider the case where the N marginals are equal to a density µ (we work in a DFT framework so the marginals rapresent the electrons which are indistinguishable). We recall that if we split µ into N μi with equal mass ( μi (x)dx = 1 N µ(x)dx), then we expect an optimal plan induced by a cyclical map such that

The simulations in figure 3 are all performed on a discretization of [-5, 5] with

02. If we focus on the support of the coupling γ12 (x, y) = (e 1 , e 2 )γ(x, y, z) we can notice that the numerical solution correctly reproduces the prescribed behavior: the transport plan is induced by a cyclical optimal map (see section 3). We refer the reader to [START_REF] Benamou | A numerical method to solve optimal transport problems with coulomb cost[END_REF] and [START_REF] Di Marino | The strictly-correlated electron functional for spherically symmetric systems revisited[END_REF] for examples in higher dimension. Remark 6.6. Theorem 3.7 actually works also for other costs functions as c(x

This means that if we use the density µ(x) = N 10 (1 + cos( π 5 x)), we expect to obtain the same solution as for the Coulomb cost. Thus if we now consider the 3-marginals case and the discretized cost c j 1 ,j 2 ,j 3 = -log(|x

we can notice (see figure 4) that we recover the same optimal coupling as the corresponding case for Coulomb (first row of figure 3). The simulation is performed on a discretizetion [-5, 5] with

, so for the discrete problem (6.2) we take

Let us first consider the two marginal case and the uniform density on [0, 1] (see example 4.2), and, as expected, we find a deterministic coupling given by γ(x, y) = µ(x)δ(y -T (x)) where T (x) = 1 -x, see figure 5. The simulation in figure 5 has been performed on a discretization of [0, 1] with M d = 1000 gridpoints and = 0.005. The multi-marginals case is more delicate to treat: the original OT problem does not admit a unique solution whereas the regularized problem does. As we have explained in section 6.1 the regularized problem is a strongly convex problem which admits a unique solution and the resulting coupling is the one with the minimal entropy. However we are able to make the IPFP algorithm converge to a selected coupling among the optimal ones for the original problem. Let us focus on the example 4.13. In this case we have all marginals equal to µ = 1 2 L| [-1,1] and we can find a deterministic coupling given by γ(x, y, z) = µ(x)δ(y -T (x))δ(z -S(x)) for the maps T (x) and S(x)). In order to select this coupling, the idea is to modify the cost function by adding a penalization term p(x, y) = τ |x-y| which makes µ 1 and µ 2 be as far as possible (if we consider µ 1 and µ 2 as particles). Thus, the discretized cost now reads as

.14)

Then we expect, as shown in figure 6, that the projection γ12 (x, y) = (e 1 , e 2 )γ(x, y, z) and γ13 (x, z) = (e 1 , e 3 )γ(x, y, z) of the computed coupling are induced by map T and S respectively. The simulation has been performed on a discretization of [-1, 1] with M d = 1000 gridpoints, = 0.0005 and η = 0.1. We, finally, take the same marginals as in the previous example, but we do not add the penalization. In this case we expect a diffuse plan since, if we forget about the marginal constraint, it