Luigi Ambrosio

AND Michael Goldman
email: goldman@u-paris.fr

Dario Trevisan
email: dario.trevisan@unipi.it

On the quadratic random matching problem in two-dimensional domains

Keywords: 2010 Mathematics Subject Classification. 60D05, 90C05, 39B62, 60F25, 35J05 Matching problem, optimal transport, geometric probability. L

We investigate the average minimum cost of a bipartite matching, with respect to the squared Euclidean distance, between two samples of n i.i.d. random points on a bounded Lipschitz domain in the Euclidean plane, whose common law is absolutely continuous with strictly positive Hölder continuous density. We confirm in particular the validity of a conjecture by D. Benedetto and E. Caglioti stating that the asymptotic cost as n grows is given by the logarithm of n multiplied by an explicit constant times the volume of the domain. Our proof relies on a reduction to the optimal transport problem between the associated empirical measures and a Whitney-type decomposition of the domain, together with suitable upper and lower bounds for local and global contributions, both ultimately based on PDE tools. We further show how to extend our results to more general settings, including Riemannian manifolds, and also give an application to the asymptotic cost of the random quadratic bipartite travelling salesperson problem.

Introduction

The minimum weight perfect matching problem on bipartite graphs, also called assignment problem, is a combinatorial optimization problem which has been historically subject of intense research by several communities, well beyond operation research and algorithm theory, including combinatorics and graph theory [START_REF] Lovász | Matching theory[END_REF], probability and statistics [START_REF] Talagrand | Upper and lower bounds for stochastic processes: modern methods and classical problems[END_REF] and even theoretical physics [START_REF] Mézard | Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications[END_REF][START_REF] Mezard | Information, physics, and computation[END_REF]. Applications related to planning and allocation of resources are classical, but recently they have seen an increased interest, e.g. in the online version of the problem, related to Internet advertising [START_REF] Mehta | Online matching and ad allocation[END_REF]. The assignment problem and its common linear programming relaxation, the optimal transport problem, provide also useful tools in machine learning and data science [START_REF] Arjovsky | Wasserstein generative adversarial networks[END_REF][START_REF] Peyré | Computational optimal transport: With applications to data science[END_REF], mostly because of its efficiency and versatility at discriminating between empirical distributions.

In this article, we focus on the Euclidean formulation of the problem, where we are given two finite families of points (x i ) n i=1 , (y j ) m j=1 ⊆ R d , with m ≥ n, and the matching cost is defined, for a given parameter p > 0, min σ∈Sn,m n i=1

|x i -y σ(i) | p , (1.1)
where S n,m denotes the set of injective maps σ : {1, . . . , n} → {1, . . . , m}. This corresponds to a minimum weight perfect matching problem on the complete bipartite graph and edge weights w ij = |x i -y j | p . In particular, when n = m, S n,n = S n is the symmetric group of permutations over n elements.

A linear programming approximation of (1.1) is given by the Monge-Kantorovich optimal transport problem, where σ is replaced by a transport plan or coupling, i.e., a matrix (π ij ) j=1,...,m i=1,...,n ∈ [0, 1] n×m with i π ij = 1/m, j π ij = 1/n. Minimization among all couplings yields a particular instance of the Wasserstein cost of order p > 0 between the empirical measures 1 (1.

2)

The Birkhoff-von Neumann theorem [START_REF] Lovász | Matching theory[END_REF] provides an exact correspondence between (1.1) and (1.2) when n = m, i.e., both the costs and optimizers are the same, an optimal plan is as a permutation matrix (both up to a factor 1/n). As a consequence, one can exploit the rich analytic structure of optimal transport and use tools from convex analysis and partial differential equations. We further assume that the points x i = X i , y j = Y j are samples of independent identically distributed (i.i.d.) random variables with a common law. This is the paradigm of many problems in geometric probability, in particular related to the theory of Euclidean additive functionals [START_REF] Steele | Probability theory and combinatorial optimization[END_REF][START_REF] Yukich | Probability theory of classical Euclidean optimization problems[END_REF]. The random Euclidean bipartite matching problem is in fact known to be quite challenging: general results and techniques that give precise results in similar problems (e.g. the non-bipartite matching problem) fail here, most notably when the dimension of the space is small.

The "critical" case turns out to be d = 2, which we also consider in this work. In this case, Ajtai, Komlós and Tusnády [START_REF] Ajtai | On optimal matchings[END_REF] first showed that for i.i.d. uniform samples on the square (0, 1) 2 the asymptotics in the case n = m → ∞ reads, for 1 p ≥ 1,

E W p p 1 n n i=1 δ X i , 1 n n i=1 δ Y i ∼ log n n p/2
.

(1.3)

Talagrand [START_REF] Talagrand | The Ajtai-Komlós-Tusnády matching theorem for general measures[END_REF] subsequently investigated the case of general laws, possibly not even absolutely continuous, providing a universal upper bound for p = 1, with the same rate. The parameter p > 0 plays a relevant role, since the solution to both problems (1.1) and (1.2) in general depends on p, and for d = 2, p < 1 the correct rate in (1.3) turns out to be n -p/2 , i.e. no logarithmic corrections appear [START_REF] Barthe | Combinatorial optimization over two random point sets[END_REF]. The classical literature focused mostly on the case p = 1 -but also on the case p = ∞, see [START_REF] Leighton | Tight bounds for minimax grid matching with applications to the average case analysis of algorithms[END_REF][START_REF] Trillos | On the rate of convergence of empirical measures in ∞transportation distance[END_REF].

A recent breakthrough in the quadratic case p = 2 was obtained by the statistical physics community, starting from the seminal work [START_REF] Caracciolo | Scaling hypothesis for the Euclidean bipartite matching problem[END_REF] and further developed in [START_REF] Caracciolo | Scaling hypothesis for the euclidean bipartite matching problem. ii. correlation functions[END_REF][START_REF] Sicuro | The Euclidean Matching Problem[END_REF]. By formally linearising the Monge-Ampère equation around the constant density to obtain the Poisson equation, they argued in particular that the optimal transport cost should be well-approximated by a negative Sobolev norm of the difference of the empirical measures. After a renormalization procedure to cut-off divergences, this led for uniform points on (0, 1) 2 to the conjecture

lim n→∞ n log n E W 2 2 1 n n i=1 δ X i , 1 n n i=1 δ Y i = 1 2π . (1.4)
A rigorous mathematical proof of (1.4) was later obtained in [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]. Further simplifications and improvements of this method have been done in [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF], leading to quantitative bounds for optimizers [START_REF] Ambrosio | On the optimal map in the 2-dimensional random matching problem[END_REF]. See also [START_REF] Goldman | A fluctuation result for the displacement in the optimal matching problem[END_REF] for a justification of the linearisation ansatz of [START_REF] Caracciolo | Scaling hypothesis for the Euclidean bipartite matching problem[END_REF] down to mesoscopic scales based on a large-scale regularity theory for the Monge-Ampère equation [START_REF] Goldman | Quantitative linearization results for the mongeampère equation[END_REF]. Let us point out that for p ≥ 1, with the exception of p = 2, even in the case of uniform points on (0, 1) 2 it is still an open problem [START_REF] Talagrand | Upper and lower bounds for stochastic processes: modern methods and classical problems[END_REF]Section 4.3.3] to determine if the limit in (1.3) exists, i.e. prove the existence of

lim n→∞ n log n p/2 E W p p 1 n n i=1 δ X i , 1 n n i=1 δ Y i .
A natural question is how (1.4) should be modified when the uniform law is replaced by a different density. As a first step in this direction, it is proven in [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF][START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF] that (1.4) holds on general closed compact Riemannian manifolds (M, g), when the cost is the square of the Riemannian distance and the law of the samples is the normalized Riemannian volume measure. This however does not even cover the case of a uniform measure on a general convex set M = Ω ⊆ R 2 , because of the presence of the boundary. For non-convex sets, an additional problem is that the Riemannian distance is different from the Euclidean one. On the other side, it was recently predicted in [9, Conjecture 2] that for every bounded connected open set Ω ⊂ R 2 with smooth boundary and every smooth, uniformly positive and bounded density with respect to the Lebesgue measure, the limit (1.4) holds true with the right-hand side multiplied by |Ω|. A rigorous proof of the upper bound lim sup

n→∞ n log n E W 2 2 1 n n i=1 δ X i , 1 n n i=1 δ Y i ≤ |Ω| 2π , (1.5)
has been obtained in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF]Theorem 1] under the hypothesis that Ω = (0, 1) 2 and that the common law of the points is absolutely continuous with a uniformly strictly positive and bounded Lipschitz density. Further conjectures on higher order terms in the asymptotic expansion of (1.4) can be found in [START_REF] Benedetto | Random assignment problems on 2d manifolds[END_REF], relying on a clever comparison between the expansions in any two smooth compact manifolds. Unfortunately, a rigorous mathematical justification of these predictions seems currently out of reach.

Main results. Our main results fully settle the validity of [9, Conjecture 2], allowing in fact for densities that are not necessarily smooth, but only Hölder continuous and strictly positive on a bounded connected domain with Lipschitz boundary. We state and prove separately first the case of the quadratic Wasserstein distance between an empirical measure and the common law of the sampled points, and then that of optimal transport between two independent samples. In the former setting we have the following Theorem 1.1. Let Ω ⊆ R 2 be a bounded connected domain, with Lipschitz boundary and ρ be a Hölder continuous probability density on Ω uniformly strictly positive and bounded from above. Given i.i.d. random variables (X i ) ∞ i=1 with common distribution ρ, we have

lim n→∞ n log n E W 2 2 1 n n i=1 δ X i , ρ = |Ω| 4π .
In the latter case, we also allow for possibly a different number of points n, m among two (jointly) i.i.d. families of random variables, extending [2, Remark 7.1].

Theorem 1.2. Let Ω ⊆ R 2 be a bounded connected domain with Lipschitz boundary and ρ be a Hölder continuous probability density on Ω uniformly strictly positive and bounded from above. Given i.i.d. random variables

(X i , Y i ) ∞ i=1 with common distribution ρ, for every q ∈ [1, ∞], lim n,m→∞ m/n→q n log n E   W 2 2   1 n n i=1 δ X i , 1 m m j=1 δ Y j     = |Ω| 4π 1 + 1 q .
Building on these results for Euclidean planar domains, it is possible using similar subadditivity arguments to generalize them to some metric measure spaces, including e.g., smooth connected domains of two dimensional Riemannian manifolds with densities, see Section 6 for more precise statements.

Notice that when q → ∞ we obtain in Theorem 1.2 the same limit as in Theorem 1.1. When n = m, Theorem 1.2 gives the asymptotic value of the minimum bipartite matching. In fact, using simple upper and lower bounds we can obtain a related result for n, m → ∞ with m -n that does not grow too fast.

Corollary 1.3. Let Ω ⊆ R 2 be a bounded connected domain, with Lipschitz boundary and ρ be a Hölder continuous probability density on Ω uniformly strictly positive and bounded from above. Given i.i.d. random variables

(X i , Y i ) ∞ i=1 with common distribution ρ, for any sequence m = m(n) ≥ n with lim n→∞ (m -n)/ log n = 0, it holds lim n→∞ 1 log n E min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 = |Ω| 2π .
We further deduce a result for the average cost of the random quadratic Euclidean bipartite travelling salesperson problem. The Euclidean travelling salesperson problem searches for the shortest cycle that visits a given set of points in R d , and the study of its random version is classical [START_REF] Beardwood | The shortest path through many points[END_REF][START_REF] Redmond | Asymptotics for Euclidean functionals with power-weighted edges[END_REF][START_REF] Yukich | Asymptotics for weighted minimal spanning trees on random points[END_REF][START_REF] Yukich | Asymptotics for the Euclidean TSP with power weighted edges[END_REF]. Its random bipartite variant requires that the cycle must alternatively connect points from two given sets of points. By extending the argument from [START_REF] Capelli | Exact value for the average optimal cost of the bipartite traveling salesman and two-factor problems in two dimensions[END_REF] from the case of the square to a general domain, we deduce that the cost of the random quadratic Euclidean bipartite travelling salesperson problem in two dimensions is asymptotically twice that of the bipartite matching problem.

Corollary 1.4. Let Ω ⊆ R 2 be a bounded connected domain, with Lipschitz boundary and ρ be a Hölder continuous probability density on Ω uniformly strictly positive and bounded from above. Given i.i.d. random variables

(X i , Y i ) ∞ i=1 with common distribution ρ, then lim n→∞ 1 log n E min σ,τ ∈Sn n i=1 |X σ(i) -Y τ (i) | 2 + |Y τ (i) -X σ(i+1) | 2 = |Ω| π ,
where we use the convention σ(n + 1) = σ(1), for σ ∈ S n .

Finally, let us comment that our main focus on limit results for the expected costs is justified by general concentration arguments that can improve to convergence in probability of the sequence of renormalized cost -see [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]Remark 4.7], where the only assumption is the validity of a L 2 -Poincaré-Wirtinger inequality, satisfied on such regular domains.

Idea of the proof. We follow the sub-additivity method introduced in this context in [START_REF] Barthe | Combinatorial optimization over two random point sets[END_REF][START_REF] Dereich | Constructive quantization: approximation by empirical measures[END_REF] and later improved in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF][START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF]. This splits the energy in a local and a global part. A central point is to prove that the global part, which encodes the defect in sub-additivity, has asymptotically only a vanishing contribution. While this method has already been implemented in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF] to obtain the upper bound (1.5) one of our main achievements is to prove that it may indeed be used also to get the corresponding lower bound. To this aim, we rely on a "boundary" variant of the optimal transport which has super-additivity properties. Similar functionals have been widely used in the theory of Euclidean additive functionals [START_REF] Barthe | Combinatorial optimization over two random point sets[END_REF][START_REF] Yukich | Probability theory of classical Euclidean optimization problems[END_REF]. The main point is to prove that when Ω = (0, 1) 2 and the law is the uniform one, the asymptotic costs for the usual optimal transport and the "boundary" versions are equal (see Proposition 3.1 and Proposition 3.2). Our proof relies on the PDE approach from [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF][START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]. To the best of our knowledge, this is the first example where the so-called Dirichlet-Neumann bracketing method is used in the context of optimal transport (see [START_REF] Armstrong | Local laws and rigidity for Coulomb gases at any temperature[END_REF] for a recent application of these ideas for Coulomb gases).

On a more technical side, we introduce two main ideas. First, in order to deal with domains which are not (0, 1) 2 or a finite union of cube, we consider a Whitney partition of our domain. While this does not affect too much the treatment of the local term, it requires a finer estimate of the global one. Indeed, as in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF][START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF] (see also [START_REF] Peyre | Comparison between W2 distance and Ḣ-1 norm, and Localization of Wasserstein distance[END_REF][START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF][START_REF] Goldman | Quantitative linearization results for the mongeampère equation[END_REF][START_REF] Ledoux | On optimal matching of Gaussian samples[END_REF]) we first rely on the estimate

W 2 2 (f, 1) f -1 2 H -1
to reduce ourselves to an estimate in H -1 . In [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF][START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF] (and also [START_REF] Goldman | Quantitative linearization results for the mongeampère equation[END_REF]) the right-hand side is then estimated thanks to Poincaré inequality by an L 2 norm. A quick look at [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF] shows that for fixed number of points n, this yields an error which is proportional to the number of cubes in the partition of our domain. Since a Whitney partition is made of infinitely many such cubes (or in any case a very large number of them, see Lemma 5.1) this would lead to a very bad estimate. This is due to the fact that a function with rapid oscillations typically has large L 2 norm but small H -1 norm. To capture this, we prove in Lemma 2.1 a finer estimate for functions with "small" support. This gives an error term which, up to logarithm, does not depend on the mesh-size of our partition (see (4.10)). In order to deal with Hölder continuous densities instead of Lipschitz ones as in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF], the second idea is to replace the Knothe map used in [9, Lemma 1] by a transport map constructed via heat flow interpolation in the spirit of the Dacorogna-Moser construction (see Lemma 2.5). This might be of independent interest and may be related to similar constructions in the literature, e.g. [START_REF] Kim | A generalization of caffarelli's contraction theorem via (reverse) heat flow[END_REF][START_REF] Neeman | Lipschitz changes of variables via heat flow[END_REF] .

Further questions. We mention here some open problems related to our results:

(1) In [START_REF] Ambrosio | On the optimal map in the 2-dimensional random matching problem[END_REF] and [START_REF] Goldman | A fluctuation result for the displacement in the optimal matching problem[END_REF] quantitative rates of convergence for the optimal transport map have been obtained. It may be worth exploring whether similar quantitative rates still hold in the general setting addressed in this article. [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF] We expect that our results may be extended to the case of less regular domains and densities, with a similar asymptotic cost. However, the validity of (1.5) must require some condition on the support of ρ, such as connectedness [9, Remark 2], which is in contrast with the universal upper bounds obtained by Talagrand [START_REF] Talagrand | Matching random samples in many dimensions[END_REF] for the case p = 1. (3) The case of unbounded domains may be treated with similar techniques, possibly leading to sharper limit results for a wide class densities, e.g. Gaussian ones [START_REF] Ledoux | On optimal matching of Gaussian samples[END_REF][START_REF] Ledoux | On optimal matching of Gaussian samples II[END_REF][START_REF] Talagrand | Scaling and non-standard matching theorems[END_REF], where existence of a limit is an open question. We mention in this regard that the case d = 1 already rises yields some interesting problems, see [START_REF] Bobkov | One-dimensional empirical measures, order statistics, and Kantorovich transport distances[END_REF][START_REF] Caracciolo | Anomalous scaling of the optimal cost in the one-dimensional random assignment problem[END_REF][START_REF] Berthet | Exact rate of convergence of the expected W2 distance between the empirical and true gaussian distribution[END_REF]. (4) Finally, it may be worth exploring how large the difference m -n can be in Corollary 1.3 so that (1.3) still holds true.

Structure of the paper. In Section 2 we fix our notation and provide some general results about negative Sobolev norms, classical optimal transport and its "boundary" version and finally a construction of a transport map from a Hölder density to the uniform one via heat flow on the cube. Section 3 is devoted to the extension of the PDE approach from [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF][START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF] to the "boundary" transport cost on the square. In Section 4 and Section 5 we prove respectively Theorem 1.1 and Theorem 1.2. Section 6 describes how similar techniques allow to consider more general settings, including Riemannian manifolds. Section 7 and Section 8 contain respectively the proof of Corollary 1.3 and Corollary 1.4. Two appendices include technical results about decomposition of a domain into sets with certain properties that may be well-known in the literature, but we did not find exactly stated in a version useful for our purposes.

Notation and preliminary results

We write |A| for the Lebesgue measure of a Borel set A ⊆ R d , and A f for the integral of an integrable function f on A. If a measure is absolutely continuous with respect to the Lebesgue measure, we always identify it with its density. A cube

Q ⊆ R d of side length is a set of the form Q = d i=1 (v i , v i + ) with (v i ) d i=1 ∈ R d and > 0.
We note (Q) the sidelength of a cube Q. By a partition of Ω ⊆ R d , we always mean in fact that Ω is covered up to Lebesgue negligible sets. We denote by |v| the Euclidean norm of a vector v ∈ R d . We use the letters C, c for positive constants whose value may vary from one line to the next and ω for a generic decreasing rate function with lim t→∞ ω(t) = 0. In this paper the domain Ω is fixed and we write A B if there is C > 0 depending only on Ω (and potentially on d and p if we consider the p-Wasserstein distance in R d ) such that A ≤ CB. For a function f , we use the notation ∇f for its gradient, ∇ • f for its divergence, ∆f = ∇ • ∇f for its Laplacian. The push-forward of a measure µ with respect to a map f is denoted f µ, i.e., f µ(A) = µ(f -1 (A)).

For f : Ω ⊆ R d → R k , we write f L p (Ω) = Ω |f | p 1/p , f L ∞ (Ω) = sup x∈Ω |f (x)| and f C α (Ω) = f L ∞ (Ω) + sup x,y∈Ω x =y |f (y) -f (x)| |x -y| α ,
for α ∈ (0, 1] for the α-Hölder norm of f (when α = 1 we obtain its Lipschitz norm). We also write

Lip Ω f = sup x,y∈Ω x =y |f (y) -f (x)| |x -y| .
We omit to specify Ω when it is clear from the context and write f p , f ∞ , f C α and Lip f . We collect below some general facts useful to prove our main results. For possible future reference we state and prove them in a slightly more general form, e.g., in general dimensions or for general, non quadratic, costs.

f W -1,p (Ω) = sup |∇φ| L q (Ω) ≤1 Ω f φ,
where q is the Hölder conjugate of p (i.e. 1 = 1 p + 1 q ). For p = 2 we simply write H -1 (Ω). Notice in particular that in order to have f W -1,p (Ω) < ∞ we must have Ω f = 0. In this case we may also restrict the supremum to functions φ having also average zero. When it is clear from the context, we will drop the explicit dependence on Ω in the norms.

The heat semi-group with null Neumann boundary conditions on Ω is well-defined as the symmetric Markov semi-group of operators (P t ) t≥0 arising as the L 2 gradient flow of the Dirichlet energy ∇f 2 2 on the Sobolev space f ∈ H 1 (Ω). It is well-known [START_REF] Varopoulos | Analysis and geometry on groups[END_REF] that the validity of L 2 -Poincaré-Wirtinger inequality on Ω is equivalent to a spectral gap, i.e., for some constant c > 0, for every f ∈ L 2 (Ω) with Ω f = 0,

P t f 2 ≤ e -ct f 2 , for t ≥ 0.
(2.1)

In particular, one has the integral representation

∆ -1 f = ∞ 0 P t f dt
for the solution to the elliptic PDE, ∆u = f in Ω with null Neumann boundary conditions and also, for the negative Sobolev norm,

f 2 H -1 = Ω f ∆ -1 f = ∞ 0 Ω f P t f dt = ∞ 0 Ω (P t/2 f ) 2 dt. (2.2)
Sobolev inequality instead is equivalent to ultracontractivity: for every t ∈ [0, 1] and

f ∈ L 1 with Ω f = 0, P t f ∞ t -d/2 f 1 . (2.3)
In terms of the symmetric heat kernel p t (x, y) = p t (y, x) (defined by the equality

P t f (x) = Ω p t (x, y)f (y)dy), it reads p t ∞ t -d/2
. This inequality can be easily combined with (2.1) to obtain the stronger inequality

P t f ∞ t -d/2 e -ct f 1 .
(2.4)

The following lemma will be used to bound the negative Sobolev norm of functions supported on subsets A ⊆ Ω.

Lemma 2.1. Let d ≥ 2, Ω ⊆ R d be a bounded connected Lipschitz domain. For every f ∈ L ∞ (Ω) with Ω f = 0, f H -1      f L 1 log f 1 f ∞ + 1 if d = 2, f 1 2 + 1 d 1 f 1 2 -1 d ∞ if d > 2.
Proof. By (2.2)

f 2 H -1 = ∞ 0 Ω f P t f dt.
By (2.4) and

P t f ∞ ≤ f ∞ , for every t > 0, Ω f P t f ≤ f 1 P t f ∞ f 1 min( f ∞ , t -d 2 exp(-ct) f 1 )
.

By integration we obtain for every t 0 > 0,

f 2 H -1 f 1 t 0 f ∞ + f 1 | log t 0 | if d = 2, t 1-d/2 0 if d > 2 .
Optimizing in t 0 by setting

t 0 = f d/2 1 f -d/2
∞ concludes the proof.

2.2. Optimal transport. We introduce some notation for the Wasserstein distance and recall few simple properties that will be used throughout. Proofs can be found in any of the monographs [START_REF] Villani | Topics in optimal transportation[END_REF][START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]. Let p ≥ 1 and µ, λ be positive Borel measures with finite p-th moments and equal mass µ(R d ) = λ(R d ). The Wasserstein distance of order p between µ and λ is defined as

W p (µ, λ) = min π∈C(µ,λ) R d ×R d |x -y| p dπ(x, y) 1 p
, where C(µ, λ) is the set of couplings between µ and λ. For a Borel subset Ω ⊆ R d , we also write

W Ω,p (µ, λ) = W p (µ ¬ Ω, λ ¬ Ω), which implicitly assumes that µ(Ω) = ν(Ω).
Let us recall that W p is a distance, in particular the triangle inequality holds:

W p (µ, ν) ≤ W p (µ, λ) + W p (ν, λ).
We will also use the classical sub-additivity inequality

W p p k µ k , k λ k ≤ k W p p (µ k , λ k ), (2.5)
for a finite set of positive measures µ k , λ k . A simple combination of the two properties above and Young inequality yield a geometric subadditivity property [22, Lemma 3.1], that we report here for the reader's convenience: there exists a constant C = C(p) > 0 such that, for every Borel partition (Ω k ) k∈N of Ω, and every ε ∈ (0, 1),

W p Ω,p µ, µ(Ω) λ(Ω) λ ≤ (1+ε) k W p Ω k ,p µ, µ(Ω k ) λ(Ω k ) λ + C ε p-1 W p Ω,p k µ(Ω k ) λ(Ω k ) χ Ω k λ, µ(Ω) λ(Ω) λ .
(2.6) If we assume that Ω is sufficiently regular, then we may use the Benamou-Brenier formula in a similar fashion as in [START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF]Lemma 3.4] (see also [START_REF] Peyre | Comparison between W2 distance and Ḣ-1 norm, and Localization of Wasserstein distance[END_REF]Corollary 3]) to bound from above the Wasserstein distance by the negative Sobolev norm. Notice that we use the fact that the Euclidean distance is bounded from above by the geodesic distance in Ω (they coincide if Ω is convex). Lemma 2.2. Assume that Ω ⊆ R d is a bounded connected open set with Lipschitz boundary. If µ and λ are measures on Ω with µ(Ω) = λ(Ω), absolutely continuous with respect to the Lebesgue measure and inf Ω λ > 0, then, for every p ≥ 1,

W p Ω,p (µ, λ) 1 inf Ω λ p-1 µ -λ p W -1,p (Ω) . (2.7)
To deal with lower bounds, we rely on a "boundary" variant of the optimal transport, introduced in [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF], but also independently studied in the setting of Euclidean bipartite matching problems, see [START_REF] Barthe | Combinatorial optimization over two random point sets[END_REF][START_REF] Dereich | Constructive quantization: approximation by empirical measures[END_REF]. Here, one is allowed to transport any amount of mass from and to the boundary of an open set Ω ⊆ R d . We write, for finite positive measures µ and ν with finite p-th moment,

W b p Ω,p (µ, λ) = min π∈Cb Ω (µ,λ) Ω×Ω |x -y| p dπ(x, y) 1 p , where Cb Ω (µ, λ) is the set of positive measures π on Ω × Ω such that π 1 ¬ Ω = µ ¬ Ω and π 2 ¬ Ω = λ ¬ Ω (π 1 , π 2
denoting respectively the first and second marginals of π). The triangle inequality holds,

W b Ω,p (µ, ν) ≤ W b Ω,p (µ, λ) + W b Ω,p (ν, λ).
A geometric superadditivity property also holds: for every disjoint family

(Ω k ) k∈N of open subsets of Ω, W b p Ω,p (µ, ν) ≥ k W b p Ω k ,p (µ, ν).
The intuition behind this property is quite clear. Given any plan on Ω, one should be able to suitably restrict it on each Ω k by stopping the transport at each boundary. To prove it rigorously, it is sufficient to argue in the case of discrete measures

µ = i µ i δ x i , ν = i ν i δ y i . The general case follows by approximation. Given π ∈ Cb Ω,p (µ, ν), we define (π k ) k∈N with π k ∈ Cb Ω k ,p (µ, ν) such that Ω×Ω |x -y| p dπ(x, y) ≥ k Ω k ×Ω k |x -y| p dπ k (x, y). (2.8)
Again by approximation, we may also assume that π = π δ (x ,y ) is discrete, with (x , y ) ∈ Ω × Ω. The following algorithm can be used to define the sequence (π k ) k . Set initially π k = 0 for every k. For every , consider the pair (x , y

). If x , y ∈ Ω k , then add π δ (x ,y ) to π k . If x ∈ Ω k and y ∈ Ω j with k = j, let t -= inf {t ∈ [0, 1] : (1 -t)x + ty / ∈ Ω k } t + = inf {t ∈ [0, 1] : tx + (1 -t)y / ∈ Ω j } , and define y -= (1 -t -)x + t -y ∈ ∂Ω k , x + = t + x + (1 -t + )y ∈ ∂Ω j . Then, add respectively π δ (x ,y -) to π k and π δ (x + ,y ) to π j . If x ∈ Ω k for some k and y ∈ Ω \ j Ω j ,
set t -and y -as the previous case and add π δ

(x ,y -) to π k . Similarly, if x ∈ Ω \ j Ω j and y ∈ Ω k then add π δ (x + ,y) to π k . In all the other cases, i.e., if x , y ∈ Ω \ k Ω k , do nothing. It is not difficult to check that (π k ) k thus defined indeed satisfies π k ∈ Cb Ω k and (2.8) holds (using in particular that p ≥ 1 to argue that |x -y | p ≥ |x -y -| p +|x + -y | p ).
Arguing similarly as in the proof of (2.6) we obtain a symmetric inequality: there exists a constant C(p) > 0 such that, for every Borel partition (Ω k ) k∈N of Ω, and every ε ∈ (0, 1),

W b p Ω,p µ, µ(Ω) λ(Ω) λ ≥(1 -ε) k W b p Ω k ,p µ, µ(Ω k ) λ(Ω k ) λ - C ε p-1 W b p Ω,p k µ(Ω k ) λ(Ω k ) χ Ω k λ, µ(Ω) λ(Ω) λ .
(2.9)

To obtain lower bounds, we will also rely on the dual formulation of W b Ω,p that reads

W b p Ω,p (µ, ν) = sup Ω f dµ + Ω gdν : f (x) + g(y) ≤ |x -y| p for x, y ∈ Ω, f = g = 0 on ∂Ω .
(2.10)

In fact, only the inequality ≥ will be used, which is immediate to check.

Remark 2.3. Many of the above properties, in particular (2.6) and (2.9), hold for the Wasserstein distance defined over any length metric space, in particular on a Riemannian manifold.

2.3.

A transport map via heat flow. The following result yields a Lipschitz map on the cube (0, 1) d transporting a given Hölder probability density to the uniform one. The main point here is that the Lipschitz constant of the map is very close to 1. This construction, possibly interesting on its own, allows us to avoid the use of general boundary regularity theory for the optimal transport map with respect to the squared distance cost [START_REF] Caffarelli | Boundary regularity of maps with convex potentials-II[END_REF], as it does not seem to be applicable to the case of the cube (it requires strictly convex and C 2 boundary). Notice that for Hölder continuous densities, counterexamples can be constructed if we work with the Knothe map as in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF].

Proposition 2.4. For d ≥ 1, α ∈ (0, 1], there exists C > 0 depending on d and α only such that the following holds: for any ρ : (0, 1) d → (0, ∞) with

(0,1) d ρ = 1 and ρ -1 C α ≤ 1/2, there exists T : (0, 1) d → (0, 1) d such that T ρ = 1, with Lip T, Lip T -1 ≤ 1 + C ρ -1 C α and T (∂(0, 1) d ) = ∂(0, 1) d .
Proof. Using the Neumann heat semi-group, we define ρ t = P t ρ, so that, for t ≥ 0, the weak formulation of the heat equation reads, for every for every f ∈ H 1 ((0, 1) d ),

(0,1) d f ρ t = (0,1) d f ρ - t 0 (0,1) d
∇f ∇ρ t dt, for every t ≥ 0.

We notice first that the assumption gives inf (0,1) d ρ ≥ 1/2 so that, for every t ≥ 0, inf (0,1) d ρ t ≥ 1/2 as well, hence

ρ -1 t ∞ 1.
By standard heat kernel estimates (one has in fact an explicit representation of p t (x, y), see [2, Appendix A] for the case d = 2) we have, for every t > 0,

∇ρ t ∞ t -1 2 + α 2 e -ct ρ -1 C α and ∇ 2 ρ t ∞ t -1+ α 2 e -ct ρ -1 C α .
We thus define the time-dependent vector field (b t ) t>0 , on (0, 1)

d , b t (x) = -∇ log ρ t (x) = - ∇ρ t (x) ρ t (x) , with ∇b t (x) = - ∇ 2 ρ t (x) ρ t (x) + ∇ρ t (x) ⊗ ∇ρ t (x) ρ 2 t (x)
.

Using the previous estimates, it follows that

b t ∞ t -1 2 + α 2 e -ct ρ -1 C α and ∇b t ∞ t -1+ α 2 e -ct ρ -1 C α .
Since the right-hand sides are integrable functions with respect to t ∈ (0, ∞), standard arguments yield that the associated flow X(t, x), i.e., the solution to

X(t, x) = x + t 0 b s (X(s, x))ds for t ≥ 0, x ∈ (0, 1) d
is well-defined and Lipschitz continuous:

Lip X(t, •) ≤ exp t 0 ∇b s ∞ ds ≤ exp Ct α 2 e -ct ρ -1 C α . (2.11)
Because of the homogeneous Neumann boundary condition, b t is tangent to the boundary ∂(0, 1) d , hence X(t, x) ∈ (0, 1) d for every t ≥ 0. Moreover, lim t→∞ X(t, x) ∈ (0, 1) d exists for every x ∈ (0, 1) d , since, for s ≤ t,

|X(t, x) -X(s, x)| ≤ t s b r ∞ dr s α 2 e -cs ρ -1 C α ,
which yields that t → X(t, x) is Cauchy as t → ∞. We thus define T (x) = lim t→∞ X(t, x).

Letting t → ∞ in (2.11), we have

Lip T ≤ exp (C ρ -1 C α ) ≤ 1 + C ρ -1 C α .
To prove that T is invertible, we simply notice that for every t > 0 the inverse map of x → X(t, x) is given by y → Y t (t, y), Y t being the flow of the "backward" in time vector field

(b t s ) s∈[0,t] , b t s (y) = -b (t-s) (y).
Writing the analogue of (2.11) for y → Y t (t, y) and letting t → ∞ yields the claimed bound on Lip T -1 .

To conclude, we need to argue that T ρ = 1. This follows quite classically from the fact that ρt = X(t, •) ρ and ρ t both solve the same continuity equation for the vector field (b t ) t>0 . Since b t is locally Lipschitz continuous, with Lipschitz norm that is integrable in time, uniqueness holds for the Cauchy problem (this follows for example by a slight adaptation of [START_REF] Santambrogio | Optimal transport for applied mathematicians[END_REF]Theorem 4.4]). Therefore, ρ t = ρt and T ρ = lim t→∞ ρt = lim t→∞ ρ t = 1.

In general, if Ω ⊆ R d is open and T : Ω → Ω is Lipschitz, then for any pair of measures µ, λ with µ(Ω) = λ(Ω), W Ω,p (T (µ ¬ Ω), T (λ ¬ Ω)) ≤ (Lip T )W Ω,p (µ, λ), (2.12)
since any coupling π between µ ¬ Ω and λ ¬ Ω induces the coupling (T, T ) π. For the "boundary" optimal transport a similar inequality holds provided that T : Ω → Ω is such that T (Ω) ⊆ Ω and T (∂Ω) ⊆ ∂Ω:

W b Ω,p (T (µ ¬ Ω), T (λ ¬ Ω)) ≤ (Lip T )W b Ω,p (µ, λ). (2.13)
Combining these observations with Proposition 2.4, we obtain the following result for the bipartite matching problem, extending [9, Lemma 1] to the case of Hölder continuous densities.

Lemma 2.5. Let d ≥ 1, α ∈ (0, 1], p ≥ 1 and let ρ be a uniformly positive and α-Hölder continuous function on a domain Ω ⊆ R d . Then, there exists a constant c = C(p, ρ) > 0 such that, for every cube Q ⊆ Ω with side length r < c, the following holds:

(1) if

(X i ) n i=1 are i.i.d. on Q with common density ρ/ρ(Q) and (X r i ) n i=1 are i.i.d. uniformly distributed on [0, r] d , then E W p * 1 n n i=1 δ X i , ρ ρ(Q) -E W p * 1 n n i=1 δ X r i , 1 |Q| r α E W p * 1 n n i=1 δ X r i , 1 |Q| , (2.14)
where

W * ∈ {W Q,p , W b Q,p };
(2) with the same notation, if (X i ) n i=1 , (Y j ) m j=1 are (jointly) i.i.d. on Q with common density ρ/ρ(Q) and (X r i ) n i=1 , (Y r j ) m j=1 are (jointly) i.i.d. and uniformly distributed on [0, r] d , then

E   W p *   1 n n i=1 δ X i , 1 m m j=1 δ Y j     -E   W p *   1 n n i=1 δ X r i , 1 m m j=1 δ Y r j     r α E   W p *   1 n n i=1 δ X r i , 1 m m j=1 δ Y r j     . (2.15)
Proof. The proof is identical to the proof of [9, Lemma 1] but we include it for the reader's convenience. We only prove that for r small enough,

(1 -Cr α )E W p * 1 n n i=1 δ X r i , 1 |Q| ≤ E W p * 1 n n i=1 δ X i , ρ ρ(Q)
since the other inequalities may be proven in the same way. Without loss of generality we assume Q = (0, r) d . Let also ρ = min Ω ρ. We define ρ r (x) = ρ(rx)r d /ρ(Q) for x ∈ (0, 1) d , so that (0,1) d ρ r = 1, and for every x, y ∈ (0, 1) d ,

ρ r (x) -ρ r (y) ≤ ρ C α ρ r α |x -y| α , thus ρ r -1 C α ≤ 1/2 if r is sufficiently small. We define S : Q → Q as S(x) = rT (x/r),
where T is the map provided by Proposition 2.4. It holds Lip S = Lip T and Lip

S -1 = Lip T -1 , S ρ/ρ(Q) = 1/|Q| is uniform, hence the variables (S(X i )) n i=1 have the same law as (X r i ) n i=1 . Moreover, S 1 n n i=1 δ X i = 1 n n i=1 δ S(X i )
has the same law as 1 n n i=1 δ X r i . Since S(∂Q) = ∂Q and S(Q) = Q, it follows from (2.12) and (2.13) that

E W p * 1 n n i=1 δ X r i , 1 |Q| ≤ (Lip T ) p E W p * 1 n n i=1 δ X i , ρ ρ(Q) ,
This proves the claim provided r is small enough so that (Lip T ) p ≤ (1 -Cr α ) -1 .

Asymptotic equivalence of the usual and boundary costs for uniform points on the square

In this section we show how to modify the proofs from [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF], using also ideas from [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF], to obtain the analogue for W b of the limits

lim n→∞ n log n E W 2 2 1 n n i=1 δ X i , 1 = 1 4π , (3.1)
and

lim n→∞ n log n E W 2 2 1 n n i=1 δ X i , 1 n n i=1 δ Y i = 1 2π ,
where (X i , Y i ) ∞ i=1 are i.i.d. uniformly distributed on (0, 1) d . We deal first with the simpler case of matching to the reference measure. Proposition 3.1. Let (X i ) ∞ i=1 be i.i.d. uniformly distributed on (0, 1) 2 . Then

lim n→∞ n log n E W b 2 (0,1) 2 ,2 1 n n i=1 δ X i , 1 = 1 4π .
Proof. To simplify the notation, we write W b instead of W b (0,1) 2 ,2 and W instead of

W (0,1) 2 ,2 . Since W b ≤ W , it is sufficient to prove that lim inf n→∞ n log n E W b 2 1 n n i=1 δ X i , 1 ≥ 1 4π . (3.2)
We do it in several steps.

Step 1 (Regularization). Write µ n = 1 n n i=1 δ X i , and for t ≥ 0,

µ t n (x) = P t µ n = 1 n n i=1 p t (x, X i ),
where we recall that (P t ) t≥0 denotes the Neumann heat semi-group on (0, 1) d with kernel p t (x, y). The triangle inequality for W b and the inequality W b ≤ W give

W b 2 µ t n , 1 ≤ W b(µ n , 1) + W (µ n , µ t n ) 2 ≤ W b 2 (µ n , 1) + W (µ n , µ t n ) W (µ n , µ t n ) + 2W (µ n , 1
) . We choose t = (log n) 4 /n, so that, the refined contractivity estimate [2, Theorem 5.2] gives

E W 2 (µ n , µ t n ) log log n n .
Together with the upper bound

E W 2 (µ n , 1) log n n ,
which follows from (3.1), we obtain that for n large enough,

E W b 2 (µ n , 1) ≥ E W b 2 µ t n , 1 -C √ log n log log n n . (3.3)
Thus, it is sufficient to prove (3.2) with µ t n instead of µ n .

Step 2 (PDE and energy estimates). Let f t n be the solution to -∆f t n = µ t n -1 with null Neumann boundary conditions, i.e., since P s µ t n = µ n,t+s ,

f t n = ∞ t (µ s n -1)ds.
Recall that, as proved in [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF] and [2, Lemma 3.14],

E (0,1) 2 ∇f t n 2 = log n 4πn + O 1 n (3.4)
and

E (0,1) 2 ∇f t n log n n 2 . (3.5)
In addition to these gradient estimates, we will also need the upper bounds

E (0,1) 2 f t n 2 1 n , and
E (0,1) 2 f t n 4 1 n 2 , (3.6)
which can be proved in a similar way. Indeed, for every x ∈ (0, 1) 2 ,

µ s n (x) -1 = 1 n n i=1 (p s (x, X i ) -1) = 1 n n i=1 ξ i
is a sum of centered i.i.d. random variables (ξ i ) n i=1 . For s ∈ (0, 1), we bound from above

E ξ 2 i = (0,1) 2 p 2 s (x, y) -1 = p 2s (x, x) -1 s -1
where we used the ultracontractivity (2.3) for the heat kernel and similarly

E ξ 4 i (0,1) 2 p 4 s (x, y) + 1 s -2 (0,1) 2 p 2 s (x, y) + 1 s -3 .
This yields immediately that

(0,1) 2 E (µ s n (x) -1) 2 1 ns .
Using Rosenthal's inequality [START_REF] Rosenthal | On the subspaces of L p (p > 2) spanned by sequences of independent random variables[END_REF] (see also [START_REF] Ledoux | On optimal matching of Gaussian samples[END_REF]) we also get

(0,1) 2 E (µ s n (x) -1) 4 1 n 3 s 3 + 1 n 2 s 2 .
For s > 1, we use (2.4), and p ∈ {2, 4},

µ s n -1 p ≤ µ s n -1 ∞ e -cs/2 s -d/2 µ s/2 n -1 1 ≤ e -cs/2 s -d/2 µ s/2 n -1 p e -cs/2 µ 1/2 n -1 p .
Taking expectation and using also the upper bounds above for s = 1/2, it follows that for some constant c > 0,

(0,1) 2 E [(µ s n (x) -1) p ]
e -cs n p/2 . For p = 2, we obtain

∞ t E (µ s n (x) -1) 2 1/2 1 t 1 √ ns ds + ∞ 1 e -cs/2 √ n ds 1 √ n ,
and, for p = 4,

∞ t E (µ s n (x) -1) 4 1/4 1 t 1 n 3 s 3 + 1 n 2 s 2 1/4 ds + ∞ 1 e -cs/4 √ n ds 1 √ n .
We conclude, for p = 2 that

E (0,1) 2 |f t n | 2 = ∞ t ∞ t E (0,1) 2 (µ n,s 1 -1)(µ n,s 2 -1) ds 1 ds 2 ≤ ∞ t ∞ t E (0,1) 2 (µ n,s 1 (x) -1) 2 1/2 E (0,1) 2 (µ n,s 2 (x) -1) 2 1/2 ds 1 ds 2 , ≤   ∞ t E (0,1) 2 (µ s n (x) -1) 2 1/2 ds   2 1 n ,
and similarly for p = 4,

E (0,1) 2 |f t n | 4 ≤   ∞ t E (0,1) 2 (µ s n (x) -1) 4 1/4 ds   4 1 n 2 .
Step 3 (Dual potential). For η ∈ (0, 1) we introduce a smooth cut-off function

χ η : (0, 1) 2 → [0, 1] such that χ η (x) = 1 on [η, 1 -η] 2 , χ η (x) = 0 on (0, 1) 2 \ [η/2, 1 -η/2] 2 , with ∇χ η ∞ η -1 , ∇ 2 χ η ∞ η -2 .
The function f n = f t n χ η can be extended by periodicity to a function on the flat torus T 2 . By [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]Corollary 3.3] applied to T 2 , there exists g n ∈ C(T 2 ) such that for x, y ∈ T 2 ,

f n (x) + g n (y) ≤ d T 2 (x, y) 2 2
, and

T 2 (f n + g n ) ≥ -e ∆fn ∞ T 2 |∇f n | 2 .
We may naturally interpret g n as being defined on (0, 1) 2 , so that using d T 2 (x, y) ≤ |x-y| 2 , we have for x, y ∈ (0, 1) 2 ,

f n (x) + g n (y) ≤ |x -y| 2 2 .
To prove that g n (x) = 0 on ∂(0, 1) 2 , we recall that [4, Remark 3.4] g n is also the unique viscosity solution to the Hamilton-Jacobi equation ∂ t u + 1 2 |∇u| 2 = 0 with initial condition -f n , hence, it coincides with the Hopf-Lax solution

g n (x) = inf y∈(0,1) 2 -f n (y) + d T 2 (x, y) 2 2 .
From this representation, we see at once that, if

f n ∞ ≤ η 2 8 (3.7) then g n (x) = 0 on ∂(0, 1) 2 , because that f n (y) = f t n (y)χ η (y) = 0 if x ∈ ∂(0, 1) 2 and d T 2 (x, y) < η/2.
Let us also notice that, since ∇f n = ∇f t n χ η + f t n ∇χ η , we have by Young and Cauchy-Schwarz

(0,1) 2 ∇f n -∇f t n 2 (0,1) 2 |∇f t n | 2 (1 -χ η ) 2 + (0,1) 2 f t n ∇χ η 2 (0,1) 2 |∇f t n | 4 1/2 η 1/2 + (0,1) 2 |f t n | 2 η -2 ,
which in expectation gives by (3.5) and (3.6)

E (0,1) 2 ∇f n -∇f t n 2 log n n η 1/2 + 1 n η -2 . (3.8)
Using (3.4) and Young inequality we find for every ε > 0, lim sup

n→∞ n log n E (0,1) 2 |∇f n | 2 - 1 4π ε + C ε η 1/2 . Choosing ε = η 1/4 we conclude lim sup n→∞ n log n E (0,1) 2 |∇f n | 2 - 1 4π η 1/4 . (3.9)
Arguing similarly, we have

n log n 2 E |∇f n | 4 1 + 1 η 4 (log n) 2 .
(3.10)

Let us finally argue that

f t n ∞ + ∇f t n ∞ + ∆f t n ∞ ∇ 2 f t n ∞ . (3.11) Since the estimate ∆f t n ∞ ∇ 2 f t n ∞ is clear and f t n ∞ ∇f t n ∞ follows from (0,1) 2 f t n = 0, we are left with the proof of ∇f t n ∞ ∇ 2 f t n ∞ . Fix x 2 ∈ (0, 1
) and consider the function φ(x 1 ) = f t n (x 1 , x 2 ). We have φ (x 1 ) = ∇f t n (x 1 , x 2 ) • e 1 (here (e 1 , e 2 ) is the canonical basis of R 2 ). By the Neumann boundary conditions,

φ (0) = φ (1) = 0, so that φ ∞ ≤ φ ∞ ∇ 2 f t n ∞ and thus ∇f t n • e 1 ∞ ∇ 2 f t n ∞ .
Exchanging the roles of x 1 and x 2 concludes the proof of the claim.

Step 4 (Conclusion) In the event A n = ∇ 2 f t n ∞ < 1/ log n we use (f n , g n ), as defined in the previous step, in the duality (2.10). The event A n has large probability, i.e.

P (A c n ) n -k , (3.12)
for every k > 0 (with an implicit constant depending on k only, see [START_REF] Ambrosio | Finer estimates on the 2-dimensional matching problem[END_REF]), hence it is sufficient to bound from below the expectation of W b 2 (µ t n , 1) on A n (using the trivial bound W b 2 (µ t n , 1) ≤ 2 on A c n ). If A n occurs, we have by (3.11),

f t n ∞ + ∇f t n ∞ + ∆f t n ∞ 1 log n ,
so that (3.7) holds if n is sufficiently large (for fixed η) and

∆f n ∞ ∇ 2 (f t n χ η ) ∞ 1 η 2 log n = ω η (n),
where for fixed η, lim n→∞ ω η (n) = 0. Thus,

1 2 W b 2 (µ t n , 1) ≥ (0,1) 2 (f n + g n ) + f n (µ t n -1)
≥ -e ωη(n)

(0,1) 2 |∇f n | 2 2 - (0,1) 2 f n ∆f t n = 1 - e ωη(n) 2 (0,1) 2 |∇f n | 2 + (0,1) 2 ∇f n • ∇f t n -∇f n .
Taking expectations, we bound the second term in the right-hand side using (3.8) and (3.9), lim sup

n→∞ n log n E (0,1) 2 ∇f n • ∇f t n -∇f n η 1/4 .
For the first term, we have lim inf

n→∞ n log n E (2 -e ωη(n) )χ An (0,1) 2 |∇f n | 2 ≥ lim inf n→∞ n log n E (0,1) 2 |∇f n | 2 -lim sup n→∞ n log n E χ A c n (0,1) 2 |∇f n | 2 (3.9) = 1 4π -Cη 1/4 ,
where we used in the last line that lim sup

n→∞ n log n E χ A c n (0,1) 2 |∇f n | 2 ≤ lim sup n→∞ P(A c n ) 1/2 n log n E (0,1) 2 |∇f n | 4 1/2 (3.12)&(3.10) = 0.
Letting η → 0 we finally obtain the thesis.

Next, we deal with the analogue result for the bipartite matching. In fact, we even consider the case of different number of points. Let us point out that we only prove a lower bound since this is what we will later use in the proof of Theorem 1.2. The corresponding upper bound may be obtained as a corollary of that theorem or more elementarily by arguing as in (5.2). Proposition 3.2. Let (X i , Y i ) ∞ i=1 be i.i.d uniformly distributed in (0, 1) 2 , then there exists a rate function 2 ω such that for every n ≤ m,

E W b 2 (0,1) 2 ,2 1 n n i=1 δ X i , 1 m m i=1 δ Y i ≥ log n 4πn 1 + n m -ω(n) .
Proof. The proof is very similar to the lower bound in the previous Proposition 3.1 so we only sketch it. Define µ n = 1 n n i=1 δ X i and λ m = 1 m m i=1 δ Y i , and write W for W (0,1) 2 ,2 and similarly for W b. We denote by µ t n and λ t m the Neumann heat-kernel regularizations of µ n and λ m . We define f t n as the solution to -∆f t n = µ t n -1 and f t m the solution of -∆f t m = λ t m -1. We finally set f t n,m = f t n -f t m , and choose t = (log n) 4 /n. Notice that if t m = (log m) 4 /m then t m ≤ t.

Repeating Step 1 of the proof of Proposition 3.1 we have, by the triangle inequality and

W b ≤ W , W b 2 µ t n , λ t m -W b 2 (µ n , λ m ) ≤ W (µ t n , µ n ) + W b(µ n , λ m ) + W (λ t m , λ m ) 2 -W b 2 (µ n , λ m ) W 2 (µ t n , µ n ) + W 2 (λ t m , λ m ) + W (µ n , λ m ) W (µ t n , µ n ) + W (λ t m , λ m ) .
The refined contractivity estimate [2, Theorem 5.2] gives

E W 2 (µ t n , µ n ) log log n n ,
and moreover, since t = (log n) 4 /n = γ/m with γ = (log n) 4 (m/n) log m, applying again [2, Theorem 5.2] we have

E W 2 (λ t m , λ m ) log γ m log log n + log m n m log log n n ,
where we used that m ≥ n and thus n m log m n is bounded. Using also that

E W 2 (µ n , λ m ) E W 2 (µ n , 1) + E W 2 (λ m , 1) log n n + log m m log n n
we find the analogously to (3.3) that for some constant C > 0,

E W b 2 (µ n , λ m ) ≥ E W b 2 (µ t n , λ t m ) -C √ log n log log n n .
Turning now to Steps 2-4, we see that the proof can be repeated verbatim using f t n,m instead of f t n and the triangle inequality to obtain the various estimates from the corresponding ones on f t n and f t m . The only additional ingredient is that using integration by parts and independence, we have the following orthogonality property

E (0,1) 2 ∇f t n • ∇f t m = -E (0,1) 2 f t n ∆f t m = - (0,1) 2 E f t n E λ t m = 0.
2 recall that a rate function is a generic decreasing function ω such that limt→∞ ω(t) = 0.

Therefore, appealing once again to [2, Lemma 3.14],

E (0,1) 2 |∇f t n,m | 2 = E (0,1) 2 |∇f t n | 2 + E (0,1) 2 |∇f t m | 2 + 2E (0,1) 2 ∇f t n • ∇f t m = E (0,1) 2 |∇f t n | 2 + E (0,1) 2 |∇f t m | 2 = | log t| 4πn + | log t| 4πm + O 1 n = log n 4πn 1 + n m + ω(n) .
Remark 3.3. Arguing as in [4, Theorem 5.2], we can also provide a lower bound for the Wasserstein distance of order 1:

lim inf n→∞ 1 √ log n E W b (0,1) 2 ,1 1 n n i=1 δ X i , 1 n n i=1 δ Y i > 0.
The proof is exactly as in [START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]Theorem 5.2]. The only difference lies in the choice of the Lipschitz function in the dual description of the Wasserstein distance. The function φ from [4, Theorem 5.2] must be replaced by φχ η where χ η is defined in the proof above.

Proof of Theorem 1.1

We recall that for a given Lipschitz and connected domain Ω and a Hölder continuous and uniformly strictly positive probability density ρ on Ω, we consider (X i ) ∞ i=1 i.i.d. random variables with common distribution ρ. Letting

µ n = 1 n n i=1 δ X i ,
we want to prove that

lim n→∞ n log n E W 2 2 (µ n , ρ) = |Ω| 4π . (4.1)
Without loss of generality, we assume that |Ω| = 1. We also omit to specify p = 2 and simply write

W Q = W Q,2 and W b Q = W b Q,2 for Q ⊆ Ω.
Let us introduce some further notation: for r > 0 we consider (X r i ) ∞ i=1 , i.i.d. uniformly distributed random variables in [0, r] 2 . We then define the average cost functions

F r (n) = E W 2 [0,r] 2 1 n n i=1 δ X r i , 1 r 2 , F b r (n) = E W b 2 [0,r] 2 1 n n i=1 δ X r i , 1 r 2
By scaling, (3.1) and Proposition 3.1, we have

F r (n) = r 2 F 1 (n) ≤ r 2 log n 4πn (1 + ω(n)) , and F b r (n) = r 2 F b 1 (n) ≥ r 2 log n 4πn (1 -ω(n)) . ( 4

.2)

Step 1 (Whitney decomposition). We consider {Q k } k a Whitney decomposition of Ω, see [START_REF] Stein | Singular integrals and differentiability properties of functions (PMS-30)[END_REF]Chapter 6]. For δ = δ(n) > 0 to be chosen below we set

U δ = {Q k : diam(Q k ) ≥ δ} and V δ = {Q k : diam(Q k ) < δ}.
For r > 0 dyadic and small enough so that Lemma 2.5 applies with r diam(Ω), we split each Q k into r -2 sub-cubes. We let V r δ be the union of these cubes. Up to relabeling we have

U δ ∪ V r δ = {Ω k } k . We finally define κ k = µ n (Ω k )/ρ(Ω k ).
For ε > 0, we estimate by the subadditivity inequality (2.5)

E W 2 Ω (µ n , ρ) ≤ (1 + ε) k E W 2 Ω k (µ n , κ k ρ) + C ε E W 2 k κ k ρχ Ω k , ρ (4.3)
and the superadditivity inequality (2.9)

E W b 2 Ω (µ n , ρ) ≥ (1 -ε) k E W b 2 Ω k (µ n , κ k ρ) - C ε E W b 2 Ω k κ k ρχ Ω k , ρ . (4.4)
In both expressions, we recognize in the right-hand sides a sum of "local" contributions and an additional "global" term. We consider these contributions separately in the next two steps and show in particular that the global term does not contribute in the limit.

Step 2 (Local term). We let N k = nµ n (Ω k ) be the number of points X i which belong to the cube Ω k . This is a random variable with binomial law and parameters (n, ρ(Ω k )). We decompose the sum in the right-hand side of (4.3) as

k E W 2 Ω k (µ n , κ k ρ) = U δ E W 2 Ω k (µ n , κ k ρ) + V r δ E W 2 Ω k (µ n , κ k ρ) .
For the first term we use the naive estimate

E W 2 Ω k (µ n , κ k ρ) diam 2 (Ω k )E [µ n (Ω k )] ≤ δ 2 ρ(Ω k ),
to get

U δ E W 2 Ω k (µ n , κ k ρ) δ 2 U δ ρ(Ω k ) δ 2 |d(•, ∂Ω) δ| δ 3 . (4.5)
As for the second term, we use (2.14) from Lemma 2.5 to infer that for every

Ω k ∈ V r δ , E W 2 Ω k (µ n , κ k ρ) = E E W 2 Ω k (µ n , κ k ρ)|N k ≤ (1 + Cr α )E N k n F √ |Ω k | (N k ) (4.2) ≤ (1 + Cr α ) |Ω k | 4πn E [log(N k )(1 + ω(N k ))] .
Using the concentration properties of the binomial random variable N k and the fact that δ = δ(n) satisfies lim n→∞ nδ 2 = ∞, we see that

E [log(N k )(1 + ω(N k ))] = log(nρ(Ω k ))(1 + ω(nr 2 δ 2 )) (4.6)
and we find

E W 2 Ω k (µ n , κ k ρ) ≤ (1 + Cr α ) log(nρ(Ω k )) 4πn |Ω k |(1 + ω(nr 2 δ 2 )). (4.7)
We now claim that

k |Ω k || log |Ω k || | log r|. (4.8) Indeed, k |Ω k || log |Ω k || = U δ |Ω k || log |Ω k || + V r δ |Ω k || log |Ω k || = U δ |Q k || log |Q k || + V δ |Q k || log |Q k | + log r| k |Q k || log |Q k || + V δ |Q k || log r| 1 + | log r| | log r|,
where we used that (notice that the finiteness of the integral below is ensured, for instance, by the finiteness of the Minkowski content of Ω, in turn ensured by Lipschitz regularity)

k |Q k || log |Q k || Ω | log d(•, ∂Ω)| < ∞.
Thus, from (4.7) and (4.8),

k E W 2 Ω k (µ n , κ k ρ) ≤ (1 + Cr α ) log n n 1 4π + C | log r| n + log n n ω(nr 2 δ 2 ) + δ 3 . (4.9)
For the analogue term in (4.4), we simply discard the terms with Ω k ∈ U δ and for the remaining ones use again (2.14) from Lemma 2.5 and the function F b r instead of F r to obtain

k E W b 2 Ω k (µ n , κ k ρ) ≥ (1 -Cr α ) log n n 1 -Cδ 4π -C | log r| n + log n n ω(nr 2 δ 2 ) .
Step 3 (Global term). We now turn to the second term in the right-hand side of (4.3), for which we show that

E W 2 Ω k κ k ρχ Ω k , ρ | log r| n . (4.10)
Notice that the inequality W b ≤ W gives an inequality also for the corresponding term in (4.4). We first use Lemma 2.2 to deduce

W 2 Ω k κ k ρχ Ω k , ρ k (κ k -1)ρχ Ω k 2 H -1 = k (κ k -1) (ρχ Ω k -ρ(Ω k )χ Ω ) 2 H -1 = k (κ k -1)f k 2 H -1
, where

f k = ρχ Ω k -ρ(Ω k )χ Ω .
Taking expectation and expanding the squares we have

k (κ k -1)f k 2 H -1 = k E[(κ k -1) 2 ] f k 2 H -1 + j =k E[(κ j -1)(κ k -1)] f j , f k H -1 .
For the first sum, we use that

E[(κ k -1) 2 ] = 1 ρ 2 (Ω k ) E (µ n (Ω k ) -ρ(Ω k )) 2 ≤ 1 nρ(Ω k ) 1 n|Ω k | and, by Lemma 2.1 with f = f k , since f k ∞ 1 and f k 1 ∼ ρ(Ω k ), f k 2 H -1 ρ(Ω k ) 2 | log ρ(Ω k )| |Ω k | 2 |log |Ω k || (4.11) to get k E[(κ k -1) 2 ] f k 2 H -1 1 n k |Ω k || log |Ω k || 1 n k |Ω k || log |Ω k || (4.8) | log r| n .
For the second sum, we use that if j = k,

E[(κ j -1)(κ k -1)] = E [µ n (Ω j )µ n (Ω k )] ρ(Ω j )ρ(Ω k ) -1 = - 1 n
together with (4.11) and Cauchy-Schwarz inequality to conclude

j =k E[(κ j -1)(κ k -1)] f j , f k H -1 1 n j =k f j H -1 f k H -1 1 n j =k |Ω j || log |Ω j || 1/2 |Ω k || log |Ω k || 1/2 1 n k |Ω k || log |Ω k || 1/2 2 | log r| n ,
where we argued as above to bound

k |Ω k || log |Ω k || 1/2 | log r| 1/2
. This proves (4.10).

Step 4 (Conclusion). Putting together (4.9) and (4.10), we conclude that

E[W 2 (µ n , ρ)] ≤ (1 + ε)(1 + Cr α ) log n n 1 4π + C log n n ω(nr 2 δ 2 ) + δ 3 + 1 ε | log r| n
and similarly

E[W b 2 (µ n , ρ)] ≥ (1 -ε)(1 -Cr α ) log n n 1 -Cδ 4π - log n n ω(nr 2 δ 2 ) + 1 ε | log r| n . We now choose δ = δ(n) = n -β with β ∈ ( 1 3 , 1
2 ), so that the condition lim n→∞ nδ 2 = ∞ holds but also lim n→∞ nδ 3 / log n = 0. This yields

(1 + ε)(1 + Cr α ) 1 4π ≥ lim sup n→∞ n log n E[W 2 (µ n , ρ)] ≥ lim inf n→∞ n log n E[W b 2 (µ n , ρ)] ≥ (1 -ε)(1 -Cr α ) 1 4π .
Letting finally r → 0 and ε → 0 we conclude the proof of (4.1).

Proof of Theorem 1.2

For (X i , Y i ) ∞ i=1 i.i.d. random variables with common distribution ρ in Ω, we write

µ n = 1 n n i=1 δ X i and λ m = 1 m m j=1 δ Y j .
We now prove that

lim n,m→∞ m/n→q n log n E W 2 2 (µ n , λ m ) = |Ω| 4π 1 + 1 q .
As above, we may assume by scaling that |Ω| = 1 and we will omit to specify p = 2, simply writing

W Q = W Q,2 and W b Q = W b Q,2 for Q ⊆ Ω.
For n ≤ m, we will mostly focus on the lower bound,

E W b 2 Ω (µ n , λ m ) ≥ log n 4πn 1 + n m -ω(n) . (5.1)
Indeed, we first show how the upper bound,

E W 2 Ω (µ n , λ m ) ≤ log n 4πn 1 + n m + ω(n) , (5.2)
can be then quickly obtained as a consequence of [4, Proposition 4.9]3 . Letting T µn (respectively T λm ) be the optimal transport maps between ρ and µ n (respectively λ m ), by independence we have

E W 2 Ω (µ n , λ m ) ≤ E Ω |T µn -T λm | 2 ρ = E W 2 Ω (µ n , ρ) + E W 2 Ω (λ m , ρ) -2 Ω E [(T µn -x)] • E (T λm -x) ρ. (5.3)
We start with the case n = m. In this case µ n and λ m have the same law, hence the last term above becomes

Ω E [(T µn -x)] • E (T λm -x) ρ = (0,1) 2 |E [(T µn -x)]| 2 ρ.
By (5.1) and Theorem 1.1 (recall (4.1)) we get lim sup

n→∞ n log n Ω |E [(T µn -x)]| 2 ρ ≤ lim n→∞ n log n E W 2 Ω (µ n , ρ) - 1 2 E W 2 Ω (µ n , λ n ) = 0.
Therefore,

E W 2 Ω (µ n , ρ) + Ω |E [(T µn -x)]| 2 ρ ≤ log n 4πn (1 + ω(n)) .
We now turn to the general case n ≤ m. Using the inequality -2ab ≤ a 2 + b 2 for the last term in (5.3), we obtain (here we use that ω(m) ≤ ω(n) and log m = log m n + log n)

E W 2 Ω (µ n , λ m ) ≤ log n 4πn (1 + ω(n)) + log m 4πm (1 + ω(m)) = log n 4πn 1 + n m + log m log n -1 n m + n m log m log n ω(m) + ω(n) ≤ log n 4πn 1 + n m + 1 log n log m n m n (1 + ω(n)) + n m + 1 ω(n) ≤ log n 4πn 1 + n m + ω(n) .
This proves (5.2).

We thus focus on the proof of (5.1). For this we follow the same steps as in Theorem 1.1. However, we need to be more careful when estimating the "global" term. Indeed, we cannot apply directly Lemma 2.2 since the measure λ m is singular. To fix this issue, we first slightly modify the Whitney decomposition to avoid cubes whose measure is too small. This guarantees that with overwhelming probability, every element of the partition contains many points so that we may argue as in [START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF]Proposition 5.2]. Notice that of course we could have used Lemma 5.1 also in the proof of Theorem 1.1.

To simplify the exposition, we postpone the proof of the following geometric result to Appendix A. Lemma 5.1. Let Ω ⊂ R d be a bounded open set with Lipschitz boundary, let {Q k } k be a Whitney decomposition of Ω. For every δ > 0 small enough (depending on Ω), letting

V δ = {Q k : diam(Q k ) ≥ δ}, there exists a family U δ = {Ω k } k of disjoint open sets such that diam(Ω k ) δ, |Ω k | ∼ δ d and U δ ∪ V δ is a partition of Ω.
Step 1 (Whitney-type decomposition and reduction to a good event). For δ = δ(n) > 0 to be specified below, let U δ and V δ be given by Lemma 5.1. For r > 0 dyadic we divide each cube Q k ∈ V δ in r -2 equal sub-cubes of sidelength r (Q k ). We let V r δ be their collection and relabel

U δ ∪ V r δ = {Ω k } k . We set θ = 1 √ log n (5.4)
and for every k such that λ m (Ω k ) = 0

κ k = µ n (Ω k ) λ m (Ω k ) .
If instead λ m (Ω k ) = 0 we arbitrarily set κ k = 0. Define the event

A = ∀Ω k ∈ V r δ ∪ U δ , 1 - µ n (Ω k ) ρ(Ω k ) + 1 - λ m (Ω k ) ρ(Ω k ) ≤ θ 4 . Notice that if |1 -µn(Ω k ) ρ(Ω k ) | + |1 -λm(Ω k ) ρ(Ω k ) | ≤ θ 4 then |1 -κ k | ≤ θ 2 . (5.5)
We claim that if

δ = 1 n β , for some β ∈ (0, 1/2), r ∼ 1 log(n) , then P(A c ) exp -c n 1-2β (log n) 3 , (5.6)
so that we can restrict ourselves to the event A in the following steps. By a union bound, to prove the claim, we bound

P(A c ) ≤ k P 1 - µ n (Ω k ) ρ(Ω k ) ≥ θ + k P 1 - λ m (Ω k ) ρ(Ω k ) ≥ θ .
We focus only on the terms involving µ n (those with λ m are analogous, recalling that n ≤ m). For every k, the number of points nµ n (Ω k ) is a binomial random variable with parameters (n, ρ(Ω k )). Hence, by Chernoff bounds

P 1 - µ n (Ω k ) ρ(Ω k ) ≥ θ ≤ exp - nρ(Ω k ) 2 log n .
Summing this over k and using that by definition of a Whitney partition, for every j ∈ N such that 2 -j ≥ δ,

#{Ω k ∈ V r δ : (Ω k ) ∈ r(2 -(j+1) , 2 -j ]} r -2 2 j , we find k P 1 - µ n (Ω k ) ρ(Ω k ) ≥ θ ≤ V r δ exp - nρ(Ω k ) 2 log n + U δ exp - nρ(Ω k ) 2 log n 2 -j ≥δ r -2 2 j exp -c nr 2 2 -2j log n + δ -1 exp -c nδ 2 log n r -2 δ -1 exp -c nr 2 δ 2 log n + δ -1 exp -c nδ 2 log n r -2 δ -1 exp -c nr 2 δ 2 log n .
Up to replacing the constant c > 0 with a smaller one, to control the terms r -2 δ -1 n β (log n) 2 , we obtain (5.6). We end this step applying (2.9) and using W b ≤ W to get for every ε > 0,

W b 2 Ω (µ n , λ m ) ≥ (1 -ε) k W b 2 Ω k (µ n , κ k λ m ) - C ε W 2 Ω k κ k χ Ω k λ m , λ m ≥ (1 -ε) V r δ W b 2 Ω k (µ n , κ k λ m ) - C ε W 2 Ω k κ k χ Ω k λ m , λ m .
Step 2 (Local term). For r > 0 let (X r i , Y r i ) be i.i.d uniformly distributed random variables in [0, r] 2 and for n ≤ m let

F b r (n, m) = E W b 2 [0,r] 2 1 n n i=1 δ X r i , 1 m m i=1
δ Y r i so that by scaling and Proposition 3.2,

F b r (n, m) = r 2 F b 1 (n, m) ≥ r 2 log n 4πn 1 + n m -ω(n) . Define N k = nµ n (Ω k ) (respectively M k = mλ m (Ω k ))
to be the number of points X i (respectively Y i ) in Ω k so that in particular κ k = N k /M k . For every given k, the random variables N k and M k are independent binomial random variables with parameters respectively (n, ρ(Ω k )) and (m, ρ(Ω k )). We then define the random variables

N * k = min(N k , M k ) and M * k = max(N k , M k ). For every fixed Ω k ∈ V r
δ , recalling that Ω k is a cube of sidelength at most of order r, and using (2.15), we have

E W b 2 Ω k (µ n , κ k λ m ) ≥ (1 -Cr α )E N * k n F b √ |Ω k | (N * k , M * k ) ≥ (1 -Cr α ) 1 4πn |Ω k |E log N * k 1 + N * k M * k -ω(N * k ) ≥ (1 -Cr α ) log(nρ(Ω k )) 4πn |Ω k | 1 + n m -ω(nr 2 δ 2 ) ,
where in the last line one can argue as for (4.6). Summing over k and using (4.8) we find

E   V r δ W b 2 Ω k (µ n , κ k λ m )   ≥ (1-Cr α ) log n n 1 -Cδ 4π 1 + n m -ω(nr 2 δ 2 ) -C | log r| n . (5.7)
Step 3 (Global term). We claim that for δ

= n -β with β ∈ 1 3 , 1 2 , n log n E W 2 Ω k κ k χ Ω k λ m , λ m χ A log log n √ log n . ( 5.8)
For this we argue along the lines of [START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF]Proposition 5.2]. We define

θ k = 1 -κ k + θ
and recall that if (5.5) holds (which is the case if A occurs) then 3 2 θ ≥ θ k ≥ 1 2 θ. We now use triangle inequality to write

W 2 Ω k κ k χ Ω k λ m , λ m W 2 Ω k κ k χ Ω k λ m , k (1 -θ k )λ m + θ λ m (Ω k ) ρ(Ω k ) ρ χ Ω k + W 2 Ω k (1 -θ k )λ m + θ λ m (Ω k ) ρ(Ω k ) ρ χ Ω k , k (1 -θ k )λ m + θ k λ m (Ω k ) ρ(Ω k ) ρ χ Ω k + W 2 Ω k (1 -θ k )λ m + θ k λ m (Ω k ) ρ(Ω k ) ρ χ Ω k , λ m .
The first and last terms are estimated in a similar way so we only estimate the first one.

We use the fact that 1 -θ k = κ k -θ and subadditivity (2.5) to bound from above

W 2 Ω k κ k χ Ω k λ m , k (κ k -θ)λ m + θ λ m (Ω k ) ρ(Ω k ) ρ χ Ω k ≤ k W 2 Ω k ((κ k -θ) + θ) λ m , (κ k -θ)λ m + θ λ m (Ω k ) ρ(Ω k ) ρ (2.5) ≤ θ k W 2 Ω k λ m , λ m (Ω k ) ρ(Ω k ) ρ .
By (4.5) and (4.7), we get

E W 2 Ω k κ k χ Ω k λ m , k (κ k -θ)λ m + θ λ m (Ω k ) ρ(Ω k ) ρ χ Ω k θ log m m + δ 3 .
For the middle term, we use again subadditivity together with Lemma 2.2 and the fact that λ m (Ω k )/ρ(Ω k ) ∼ 1 in the event A, to infer

W 2 Ω k (1 -θ k )λ m + θ λ m (Ω k ) ρ(Ω k ) ρ χ Ω k , k (1 -θ k )λ m + θ k λ m (Ω k ) ρ(Ω k ) ρ χ Ω k ≤ W 2 Ω k θ λ m (Ω k ) ρ(Ω k ) χ Ω k ρ, k θ k λ m (Ω k ) ρ(Ω k ) χ Ω k ρ 1 θ k (θ -θ k ) λ m (Ω k ) ρ(Ω k ) χ Ω k ρ 2 H -1 = 1 θ k (κ k -1) λ m (Ω k ) ρ(Ω k ) χ Ω k ρ 2 H -1 1 θ k µ n (Ω k ) ρ(Ω k ) -1 χ Ω k ρ 2 H -1 + 1 θ k λ m (Ω k ) ρ(Ω k ) -1 χ Ω k ρ 2 H -1
, where in the last line we used that (κ k -1)

λm(Ω k ) ρ(Ω k ) = ( µn(Ω k ) ρ(Ω k ) -1)-( λm(Ω k ) ρ(Ω k ) -1
) and triangle inequality in H -1 . By (4.10) we obtain

E   i (1 -κ k ) λ m (Ω k ) ρ(Ω k ) χ Ω k ρ 2 H -1   | log r| n .
This proves

E W 2 Ω i κ k χ Ω k λ m , λ m χ A θ log n n + δ 3 + 1 θ | log r| n ,
which recalling the choice (5.4) of θ, δ and r concludes the proof of (5.8).

Step 4 (Conclusion). Putting together (5.7) and (5.8) and using (5.6), we find that

n log n E W b 2 Ω (µ n , λ m ) ≥ (1 -ε)(1 -Cr α ) 1 -Cδ 4π 1 + n m -ω(nr 2 δ 2 ) - C ε log log n √ log n + exp -c n 1-2β (log n) 3 ≥ 1 4π 1 + m n + ω(n) ,
where in the last line we chose ε = ε(n) → 0 as n → ∞ with ε log log n/ √ log n.

Extension to general settings

In this section we show that the techniques developed above lead to extensions of Theorem 1.1 and Theorem 1.2 to more general settings, including Riemannian manifolds. Let us give the following definition (see also [START_REF] Trillos | On the rate of convergence of empirical measures in ∞transportation distance[END_REF]Definition 3.1]). Definition 6.1. We say that a metric measure space (Ω, d, m) with m(Ω) = 1 is welldecomposable if for every ε > 0 there exist a finite family, that we call ε-decomposition,

(U k , Ω k , T k ) k , such that, for every k, (1)
U k ⊆ Ω is open and m(U k ∩ U k ) = 0 for k = k , m( k U k ) = 1, (2)
T k : U k → Ω k , is invertible with T k (U k ) = Ω k and T k (∂U k ) = ∂Ω k and Lip T k , Lip T -1 k ≤ (1 + ), (4) (T k ) (m ¬ U k (3)
) has Hölder continuous density, uniformly positive and bounded.

If (Ω, d, m) is well-decomposable we define

|Ω| = inf ε>0 k |Ω k | : (U k , Ω k , T k ) k is an ε-decomposition of Ω .
Remark 6.2. Notice that the quantity inside the brackets is a decreasing function of ε and that if Ω is a Riemannian manifold it coincides with its (Riemannian) volume.

In addition to such decomposability assumption, we need to assume the validity of the inequality:

W 2 2 (f m, m) f -1 2 L 2 (m)
, for every probability density f . (6.1)

This estimate is analog to the weaker form of the combination of (2.7) with Lemma 2.1 used in [START_REF] Benedetto | Euclidean random matching in 2d for non-constant densities[END_REF][START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF] and is intimately connected to the validity of an L

q ∈ [1, ∞], lim n,m→∞ m/n→q n log n E   W 2 2   1 n n i=1 δ X i , 1 m m j=1 δ Y j     = |Ω| 4π 1 + 1 q , (6.2)
and

lim n→∞ n log n E W 2 2 1 n n i=1 δ X i , m = |Ω| 4π . (6.3)
Examples of metric measure spaces for which the result apply include smooth domains in two-dimensional Riemannian manifolds with boundary, as shown in Appendix B. But also other, less regular cases, may be included, e.g. polyhedral surfaces.

Proof of Theorem 6.3. As usual, we simply write W = W 2 and W b = W b 2 . We write

µ n = 1 n n i=1 δ X i , λ m = 1 m m j=1 δ Y j . For ε > 0, let (U k , Ω k , T k ) k be an ε-decomposition and let κ k = µ n (U k )/λ m (U k ).
In both cases, after an application of the sub-additivity and super-additivity inequalities (2.6), (2.9) we are reduced to estimate separately a finite sum of "local" terms, one for each U k , and the global terms (using that W b ≤ W ) that are respectively

E W 2 2 k κ k χ U k λ m , λ m and E W 2 2 k µ n (U k ) m(U k ) χ U k m, m . ( 6
E W 2 (µ n ¬ U k , κ k λ m ¬ U k ) ≤ (1 + ε) 2 E   N W 2   1 N N i=1 δ X k,i , 1 M M j=1 δ Y k,j     , and
E W b 2 U i (µ n ¬ U k , κ k λ m ¬ U k ) ≥ (1 + ε) -2 E   N W b 2 Ω i   1 N N i=1 δ X k,i , 1 M M j=1 δ Y k,j     ,
where, for every k, we consider i.i.d. random variables (X k,i , Y k,i ) ∞ i=1 with common law ρ k given by (T k ) m normalized to a probability measure on Ω k and N , M are further independent random variables with binomial laws of parameters respectively (n, m(U k )), (m, m(U k )). For each k, by conditioning upon N , M , and by (5.1), (5.2) with (X k,i , Y k,i ) ∞ i=1 , we obtain that, letting

N * = min {N, M }, M * = max {N, M } E   N W 2   1 N N i=1 δ X k,i , 1 M M j=1 δ Y k,j     ≤ E N |Ω k | log N * 4πN * 1 + N * M * + ω k (N * ) ,
and similarly

E   N W b 2 Ω k   1 N N i=1 δ X k,i , 1 M M j=1 δ Y k,j     ≥ E N |Ω k | log N * 4πN * 1 + N * M * -ω k (N * ) . As n → ∞ with m/n → q ∈ [1, ∞], we have N/n → m(U k ), M/m → m(U k ) and N * /N → 1, M * /N → q,
with the usual concentration inequalities which eventually give

lim n→∞ n log n E N |Ω k | log N * 4πN * 1 + N * M * -ω k (N * ) = |Ω k | 4π 1 + 1 q .
A similar argument gives

lim n→∞ n log n E W 2 U k µ n (U k ) m(U k ) m, m = |Ω k | 4π . (6.5)
Global term. We consider the second term in (6.4), for which we apply (6.1) obtaining

E W 2 2 k µ n (U k ) m(U k ) χ U k m, m k E µ n (U k ) m(U k ) -1 2 m(U k ) C(ε) n .
After multiplying by n/ log n and letting first n → ∞ and then ε → 0, the validity of (6.2) is thus settled. To bound (6.4), we let κ = min k κ k -ε, so that, as n → ∞, κ → 1 -ε and using union and Chernoff bounds we have the crude estimate (but sufficient for our purposes)

P (∀k, |κ -κ k | > 2 ) C(ε) n . (6.6)
We bound from above, using the subadditivity inequality (2.5) and the triangle inequality

E W 2 2 k κ k χ U k λ m , λ m ≤ E W 2 2 k (κ k -κ)χ U k λ m , (1 -κ)λ m ≤ E W 2 2 k (κ k -κ)χ U k λ m , k (κ k -κ) λ m (U k ) m(U k ) χ U k m + E W 2 2 k (κ k -κ) λ m (U k ) m(U k ) χ U k m, (1 -κ)m + E W 2 2 ((1 -κ)m, (1 -κ)λ m ) . (6.7)
For the first term, we use again subadditivity (2.5),

E W 2 2 k (κ k -κ)χ U k λ m , k (κ k -κ) λ m (U k ) m(U k ) χ U k m ≤ k E W 2 U k (κ k -κ)λ m , (κ k -κ) λ m (U k ) m(U k ) m ≤ k C(ε) diam(E) 2 n 4 + εE W 2 U k λ m , λ m (U k ) m(U k ) m C(ε) diam(E) 2 n + ε k |Ω k | log m m ,
having used (6.6) and (6.5) with λ m instead of µ n . For the second term in (6.7), we use (6.1),

E W 2 2 k (κ k -κ) λ m (U k ) m(U k ) χ U k m, (1 -κ)m k E (κ k -κ) λ m (U k ) m(U k ) -(1 -κ) 2 m(U k ) C(ε) n ,
having used the variance bounds

E λ m (U k ) m(U k ) -1 2 + E (κ k -1) 2 C(ε) n .
For the third term in (6.7), we use that 1 -κ > 2ε with probability smaller than C(ε)/n and then, on the complementary event, the already settled (6.3), to obtain

E W 2 2 ((1 -κ)m, (1 -κ)λ m ) C(ε) n + ε m log m .
After collecting all these bounds, we multiply by n/ log n and let first n → ∞ and then ε → 0 to conclude.

Proof of Corollary 1.3

We recall that with our usual notation, we want to prove that for any sequence m = m(n) ≥ n with lim n→∞ (m -n)/ log n = 0, it holds

lim n→∞ 1 log n E min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 = |Ω| 2π . The inclusion S n ⊆ S n,m yields immediately min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 ≤ min σ∈Sn n i=1 |X i -Y σ(i) | 2 ,
so that it always holds, for any sequence m = m(n) ≥ n,

lim sup n→∞ 1 log n E min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 ≤ lim n→∞ 1 log n E min σ∈Sn n i=1 |X i -Y σ(i) | 2 = |Ω| 2π .
To show the converse inequality, given σ ∈ S n,m we induce a matching between (X i ) m i=1 and (Y j ) m j=1 by pairing the points (

X i ) m i=n+1 to those in (Y j ) m j=1 \ Y σ(i) n i=1 . Since Ω is bounded, we obtain the inequality min σ∈Sm m i=1 |X i -Y σ(i) | 2 ≤ min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 + diam(Ω) 2 (m -n).
Taking expectation and using the assumption (m-n)/ log n → 0, which in particular gives log n/ log m → 1 yields

|Ω| 2π = lim m→∞ 1 log m E min σ∈Sm m i=1 |X i -Y σ(i) | 2 ≤ lim inf n→∞ 1 log n E min σ∈Sn,m n i=1 |X i -Y σ(i) | 2 .
8. Proof of Corollary 1.4

Given (x i ) n i=1 , (y i ) n i=1 ⊆ R 2 , we introduce the notation To prove the converse inequality, let τ ∈ S n be a minimizer for C 2 TSP ((X i ) n i=1 ) and let σ ∈ S n be an optimal matching between (X i ) n i=1 and (Y i ) n i=1 , with respect to the quadratic cost. We let τ = σ • τ ∈ S n , so that

C 2 b TSP ((X i ) n i=1 , (Y i ) n i=1 ) ≤ n i=1 |X τ (i) -Y σ(τ (i)) | 2 + |Y σ(τ (i)) -X τ (i+1) | 2 .
For every ε > 0, using the inequality |a + b| 2 ≤ (1 + ε)|a| 2 + (1 + ε -1 )|b| 2 , we bound from above, for every i ∈ {1, . . . , n},

|Y σ(τ (i)) -X τ (i+1) | 2 = |(Y σ(τ (i)) -X τ (i) ) + (X τ (i) -X τ (i+1) )| 2 ≤ (1 + ε)|Y σ(τ (i)) -X τ (i) | 2 + (1 + ε -1 )|X τ (i) -X τ (i+1) | 2 .
Summing upon i gives The thesis follows letting ε → 0.

C 2 b TSP ((X i ) n i=1 , (Y i ) n i=1 ) ≤(2 + ε) min σ∈Sn n i=1 |X i -Y σ(i) | 2 + (1 + ε -1 )C
To prove (8.2), we rely on the known (deterministic) bound sup

(x i ) n i=1 ⊆Ω C 2 TSP ((x i ) n i=1 ) 1,
where the implicit constant depends on Ω only. We briefly recall here a proof, for the reader convenience, based on the space-filling curve heuristic [START_REF] Steele | Probability theory and combinatorial optimization[END_REF]Section 2.6]. Without loss of generality, we consider the case Ω = Q a square. Consider a Peano curve ψ : [0, 1] → Q, that is surjective and 1/2-Hölder continuous. Let t i ∈ [0, 1] be such that ψ(t i ) = x i and choose τ ∈ S n such that (t τ (i) ) n i=1 are in increasing order. Then,

C 2 TSP ((x i ) n i=1 ) ≤ n i=1 |ψ(t τ (i) ) -ψ(t τ (i+1) )| 2 ≤ ψ C 1/2 n i=1 |t τ (i) -t τ (i+1) | ≤ 2 ψ C 1/2 .
Remark 8.1. The argument above is in fact rather general and it applies on general metric spaces provided that the cost of the quadratic bipartite matching problem is asymptotically larger than that of the travelling salesperson problem. Also, already in the Euclidean setting, it seems possible to adapt it to the case p = 2, but existence of the limit for the asymptotic cost of the bipartite matching problem is not known. Finally, exactly as in [START_REF] Capelli | Exact value for the average optimal cost of the bipartite traveling salesman and two-factor problems in two dimensions[END_REF], the argument applies as well for the random bipartite quadratic 2-factor problem, i.e., for the relaxation where the single tour is replaced by a collection of disjoint ones: indeed, since the cost becomes smaller, it is sufficient to notice that the lower bound (8.1) still holds.

tubular neighbourhood of γ(S 1 ), i.e., a map Γ : S 1 × (-r 0 , r 0 ) → M , Γ(t, r) = exp γ(t) (r γ(t) ⊥ ), where γ(t) ⊥ denotes the inward unit normal to ∂Ω at γ(t). By compactness, if r 0 is sufficiently small, such a parametrization exist, is smooth (in particular Lipschitz) and Γ(S 1 × (0, r 0 )) ⊆ Ω while Γ(S 1 × (-r 0 , 0)) ⊆ Ω c .

Let us also assume that r 0 is sufficiently small so that exp p is well-defined and invertible on a ball B p (r 0 ) for every p ∈ Ω and with Lipschitz constant smaller than 1 + ε (together with its inverse). Then, normal coordinates T p = exp -1 p are well defined on the image exp p (B p (r 0 )).

Consider a finite mesh {t i } i ⊆ S 1 with d(t i , t i+1 ) < r 0 / Lip(Γ) and let p i = γ(t i ) and q i = Γ(t i , r 0 /2) ∈ Ω. Then, d(q i , q i+1 ) ≤ Lip(Γ)|t i -t i+1 | < r 0 /4. On the other side, d(q i , ∂Ω) ≥ r 0 /2, so that B(q i , r 0 /4) ∩ ∂Ω = ∅ and in particular the geodesic segment connecting q i to q i+1 does not intersect ∂Ω. We define U i to be the rectangle-like region with boundaries given by the geodesic segments connecting p i to q i , then q i to q i+1 and q i+1 and p i , together with the piece of boundary between p i and p i+1 (i.e., Γ(0, s) with s ∈ [t i , t i+1 ]. We also notice that the boundary is piecewise C 1 (with non-tangential intersections of different pieces), hence when transformed by T p i it is Lipschitz regular.

This settles the tubular neighbourhood of ∂Ω. To complete the construction we then consider a finite set of points {p j } ⊆ Ω \ Γ(S 1 , (0, r 0 /2)) such that for every q ∈ Ω there exist p j with d(p j , q) < r 0 and for every p j there are distinct points p k , p such that max {d(p j , p k ), d(p k , p ), d(p j , p )} < r 0 /2.

We connect with geodesic segments all points p j , p k with d(p j , p k ) < r 0 /2 as well as p j , q + i with d(p j , q + i ) < r 0 /2. Notice that none of these segments intersects ∂Ω. As in the proof of [24, Theorem 2.3.A.1], up to enlarging the family to include intersections of these segments, we obtain a decomposition into geodesic polygons. To conclude the construction of the ε-decomposition it is sufficient to add these polygons to the sets U i obtained above.

In the case of non-empty ∂M , the only difference is that the notion of tubular neighbourhood must take into account points p = γ(t) ∈ Ω ∩ ∂M , hence Γ(t, r) is defined only for r ∈ [0, r 0 ).

Finally, inequality (6.1) follows from the validity of an L 2 -Poincaré-Wirtinger inequality on Ω, adapting the proof of [START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF]Lemma 3.4] to the Riemannian setting. In turn, the validity of Poincaré-Wirtinger inequality seems to be folklore on smooth connected compact Riemannian manifolds and a full proof, also in presence of a smooth boundary, may be given by gluing local inequalities [23,Section 10.1].

π

  ij |x i -y j | p = W p p

2. 1 .

 1 Negative Sobolev norms. Let Ω ⊆ R d be a bounded connected open set with Lipschitz boundary. We denote the negative Sobolev norm by

  Ω k ⊆ R 2 isopen, bounded and connected, with Lipschitz boundary,

2 b

 2 i) -x τ (i+1) | 2for the costs of the Euclidean travelling salesperson problem andC 2 b TSP ((x i ) n i=1 , (y i ) n i=1 ) = min τ,τ ∈Sn n i=1 |x τ (i) -y τ (i) | 2 + |y τ (i) -x τ (i+1) | 2for the cost of its bipartite variant (writing τ (n + 1) = τ (1) for permutations τ ∈ S n ). The proof of the lower bound follows straightforwardly taking expectation in the inequalityC 2 b TSP ((X i ) n i=1 , (Y i ) n i=1 ) ≥ 2 min σ∈Sn n i=1 |X i -Y σ(i) | 2 . TSP ((X i ) n i=1 , (Y i ) n i=1 ) ≥ |Ω| π .

  2 -Poincaré-Wirtinger inequality (see e.g.[START_REF] Goldman | Convergence of asymptotic costs for random euclidean matching problems[END_REF] Lemma 3.4]).

	Theorem 6.3. Let (M, d) be a length space and Ω ⊆ M be open with diam(Ω) < ∞. Let
	m be a probability measure on M with m(Ω) = 1, such that (6.1) holds and (Ω, d, m) is
	well-decomposable and consider (X i , Y i ) ∞ i=1 be i.i.d. with common density m. Then, for
	every

  .[START_REF] Ambrosio | A PDE approach to a 2-dimensional matching problem[END_REF]

	Local terms. Using (2.12) with T -1 k	and (2.12) with T k and taking expectation give
	respectively

in fact the statement there wrongly contains an equality that should be instead an inequality ≤

A. thanks the support of the PRIN 2017 project "Gradient Flows, Optimal Transport and Metric Measure Structures". M.G. was partially supported by the project ANR-18-CE40-0013 SHAPO financed by the French Agence Nationale de la Recherche (ANR). D.T. was partially supported by Gnampa project 2020 "Problemi di ottimizzazione con vincoli via trasporto ottimo e incertezza".

Appendix A. Proof of Lemma 5.1

Let us recall that given a bounded open set Ω with Lipschitz boundary, and {Q k } k a Whitney decomposition of Ω we want to construct for every δ > 0 small enough (depending on Ω), a partition of U δ ∪ V δ of Ω such that V δ = {Q k : diam(Q k ) ≥ δ}, and

We start by constructing a partition of the set

2 . Such a family can be for instance constructed by choosing any starting point x 1 ∈ ∂Ω and then setting

i=1 {x i }) as long as ∂Ω is not covered by ∪ k-1 i=1 B δ (x i ). We then set

Notice that since {|x -

In words this means that we consider in U δ either cubes which are at distance greater than √ dδ from Ω c but which are not in V δ or we combine the sets Ω k with the cubes which intersect both A δ and its complement. The choice to which Ω k we associate Q i is arbitrary (as long as they intersect). The sets Ω k satisfy both diam(Ω k ) δ and |Ω k | ∼ δ d since for every j such that d(Q j , Ω c ) < √ dδ and Q j ∩ A c δ = ∅, we have (Q j ) ∼ δ and thus the number of such cubes which can intersect Ω k is uniformly bounded with respect to δ.

Appendix B. Decomposition of Riemannian manifolds

We show that Theorem 6.3 applies to compact connected smooth Riemannian manifolds (M, g), possibly with boundary. In fact, we argue in the case of bounded connected domains Ω ⊆ M with smooth boundary (so that M itself does not need to be compact). Notice that Ω with the restriction of the metric g is also a compact connected smooth Riemannian manifolds with boundary, but the induced distance would then be larger (and different, in general) than the restriction of the one on M , exactly as in the Euclidean case.

Lemma B.1. Let (M, g) be a two-dimensional Riemannian manifold with (possibly empty) boundary ∂M and let Ω ⊆ M be open, bounded, connected, with smooth boundary. Let also ρ be a Hölder continuous probability density on Ω (with respect to the Riemannian volume) uniformly strictly positive and bounded from above. Then, (Ω, d, ρ), where d is the restriction of the Riemannian distance on M , is well-decomposable and (6.1) holds.

Proof. To show that Ω is well-decomposable, we modify the well-known argument to obtain a decomposition of M into geodesic triangles, see e.g. [24, Theorem 2.3.A.1], to take into account also the presence of the boundary ∂Ω. To simplify the exposition, assume first that ∂M = ∅. Then, by compactness, ∂Ω is the union of a finite family of closed simple C 1 curves, and to simplify again let us assume that there is only one such curve γ : S 1 → M parametrized so that γ(t) = 0 for every t ∈ S 1 . Consider also a parametrization of the