On the quadratic random matching problem in two-dimensional domains

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{\"o}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 [29], probability and statistics [43] and even theoretical physics [32,31]. 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 [30]. The assignment problem and its common linear programming relaxation, the optimal transport problem, provide also useful tools in machine learning and data science [5,34], 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 [29] 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 [40,49]. 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 [1] 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, Talagrand [44] subsequently investigated the case of general laws, possibly not even absolutely continuous, providing a universal upper bound for p = 1, with the same rate.
On the quadratic random matching problem in two-dimensional domains 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 [7]. The classical literature focused mostly on the case p = 1 -but also on the case p = ∞, see [28,46].
A recent breakthrough in the quadratic case p = 2 was obtained by the statistical physics community, starting from the seminal work [15] and further developed in [16,39]. 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 (1.4) A rigorous mathematical proof of (1.4) was later obtained in [3]. Further simplifications and improvements of this method have been done in [2], leading to quantitative bounds for optimizers [4]. See also [20] for a justification of the linearisation ansatz of [15] down to mesoscopic scales based on a large-scale regularity theory for the Monge-Ampère equation [21]. 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 [ 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 [2,3] 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 [ δ Yi ≤ |Ω| 2π , (1.5) has been obtained in [9,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 [10], 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 EJP 27 (2022), paper 54. 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.
In the latter case, we also allow for possibly a different number of points n, m among Building on these results for Euclidean planar domains, it is possible using similar sub-additivity 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 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 [8,36,51,50]. Its random bipartite variant requires that the cycle must alternatively connect points from two given sets of points. By extending the argument from [14] 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.
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 [3,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 [7,18] and later improved in [9,22]. 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 [9] 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 [7,49]. 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 [2,3]. 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 [6] 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 [9,22] (see also [35,3,21,26]) we first rely on the estimate to reduce ourselves to an estimate in H −1 . In [9,22] (and also [21]) the right-hand side is then estimated thanks to Poincaré inequality by an L 2 norm. A quick look at [9] 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 [9], 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. [25,33].

Further questions
We mention here some open problems related to our results: 1. In [4] and [20] 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.
2. 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 EJP 27 (2022), paper 54.
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 [42] 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 [26,27,45], 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 [12,17,11].
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 [2,3] 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 ∈ 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 α ∈ (0, 1] for the α-Hölder norm of f (when α = 1 we obtain its Lipschitz norm). We also write Lip Ω f = sup 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.

Negative Sobolev norms
Let Ω ⊆ R d be a bounded connected open set with Lipschitz boundary. We denote the negative Sobolev norm by 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 [47] 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, In particular, one has the integral representation for the solution to the elliptic PDE, ∆u = f in Ω with null Neumann boundary conditions and also, for the negative Sobolev norm, Sobolev inequality instead is equivalent to ultracontractivity: for every t ∈ [0, 1] and f ∈ L 1 with Ω f = 0, In terms of the symmetric heat kernel p t (x, y) = p t (y, x) (defined by the equality P t f (x) = The following lemma will be used to bound the negative Sobolev norm of functions supported on subsets A ⊆ Ω.
By integration we obtain for every t 0 > 0, Optimizing in t 0 by setting concludes the proof.

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 [48,38].
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 is the set of couplings between µ and λ. For a Borel subset Ω ⊆ R d , we also write which implicitly assumes that µ(Ω) = ν(Ω). Let us recall that W p is a distance, in particular the triangle inequality holds: We will also use the classical sub-additivity inequality 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), If we assume that Ω is sufficiently regular, then we may use the Benamou-Brenier formula in a similar fashion as in [22,Lemma 3.4] (see also [35,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, To deal with lower bounds, we rely on a "boundary" variant of the optimal transport, introduced in [19], but also independently studied in the setting of Euclidean bipartite matching problems, see [7,18]. 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, ¬ Ω = λ ¬ Ω (π 1 , π 2 denoting respectively the first and second marginals of π).
The triangle inequality holds, A geometric superadditivity property also holds: for every disjoint family (Ω k ) k∈N of open subsets of Ω, 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 δ xi , ν = i ν i δ yi . The general case follows by approximation. Given π ∈ Cb Ω,p (µ, ν), we define (π k ) k∈N with π k ∈ Cb Ω k ,p (µ, ν) such that 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 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 EJP 27 (2022), paper 54.
ε ∈ (0, 1), To obtain lower bounds, we will also rely on the dual formulation of W b Ω,p that reads (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.

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 [13], 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 [9].
, there exists C > 0 depending on d and α only such that the following holds: for any ρ : 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 ), We notice first that the assumption gives inf (0,1) d ρ ≥ 1/2 so that, for every t ≥ 0, 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, We thus define the time-dependent vector field . EJP 27 (2022), paper 54.
Using the previous estimates, it follows that 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 is well-defined and Lipschitz continuous: Because of the homogeneous Neumann boundary condition, b t is tangent to the boundary 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 . 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 [38,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 µ(Ω) = λ(Ω), 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 (∂Ω) ⊆ ∂Ω: 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, 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 , 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 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 [3], using also ideas from [2], to obtain the analogue for W b of the limits We deal first with the simpler case of matching to the reference measure.
Proof. To simplify the notation, we write W b instead of W b (0,1) 2 ,2 and W instead of We do it in several steps.
Step 1 (Regularization). Write µ n = 1 n n i=1 δ Xi , and for t ≥ 0, 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 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, log n log log n 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 , Recall that, as proved in [3] and [2, Lemma 3.14], In addition to these gradient estimates, we will also need the upper bounds which can be proved in a similar way. Indeed, for every x ∈ (0, 1) 2 , where we used the ultracontractivity (2.3) for the heat kernel and similarly This yields immediately that ns .
Using Rosenthal's inequality [37] (see also [26]) we also get For s > 1, we use (2.4), and p ∈ {2, 4}, Taking expectation and using also the upper bounds above for s = 1/2, it follows that for some constant c > 0, On the quadratic random matching problem in two-dimensional domains and, for p = 4, We conclude, for p = 2 that and similarly for p = 4, Step 3 (Dual potential). For η ∈ (0, 1) we introduce a smooth cut-off function χ η : The function f n = f t n χ η can be extended by periodicity to a function on the flat torus T 2 . By [3, Corollary 3.3] applied to T 2 , there exists g n ∈ C(T 2 ) such that for x, y ∈ T 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 , To prove that g n (x) = 0 on ∂(0, 1) 2 , we recall that [3, 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 From this representation, we see at once that, if 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- (3.8) Using (3.4) and Young inequality we find for every ε > 0,   ∇ 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.
for every k > 0 (with an implicit constant depending on k only, see [2]), 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 where for fixed η, lim n→∞ ω η (n) = 0. Thus, Taking expectations, we bound the second term in the right-hand side using (3.8) and 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).
i.d uniformly distributed in (0, 1) 2 , then there exists a rate function 2 ω such that for every n ≤ m,

Repeating
Step 1 of the proof of Proposition 3.1 we have, by the triangle inequality 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, 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 Therefore, appealing once again to [2, Lemma 3.14],

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 ( random variables with common distribution ρ. Letting we want to prove that lim n→∞ n log n E W 2 2 (µ n , ρ) = |Ω| 4π . 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 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)) .
For ε > 0, we estimate by the subadditivity inequality (2.5) and the superadditivity inequality (2.9) 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 For the first term we use the naive estimate to get As for the second term, we use (2.14) from Lemma 2.5 to infer that for every Ω k ∈ V r δ , Using the concentration properties of the binomial random variable N k and the fact that δ = δ(n) satisfies lim n→∞ nδ 2 = ∞, we see that and we find We now claim that k |Ω k || log |Ω k || | 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 Thus, from (4.7) and (4.8), 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 Step 3 (Global term). We now turn to the second term in the right-hand side of (4.3), for which we show that 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 Taking expectation and expanding the squares we have For the first sum, we use that For the second sum, we use that if j = k, together with (4.11) and Cauchy-Schwarz inequality to conclude where we argued as above to bound k |Ω k || log |Ω k || 1/2 | log r| 1/2 . This proves (4.10).

Proof of Theorem 1.2
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, can be then quickly obtained as a consequence of [3, Proposition 4.9] 3 . Letting T µn (respectively T λm ) be the optimal transport maps between ρ and µ n (respectively λ m ), by independence we have We start with the case n = m. In this case µ n and λ m have the same law, hence the last term above becomes On the quadratic random matching problem in two-dimensional domains By (5.1) and Theorem 1.1 (recall (4.1)) we get lim sup Therefore, (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) 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 [22,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. 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 If instead λ m (Ω k ) = 0 we arbitrarily set κ k = 0. Define the event δ = 1 n β , for some β ∈ (0, 1/2), r ∼ 1 log(n) , so that we can restrict ourselves to the event A in the following steps. By a union bound, to prove the claim, we bound 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 Summing this over k and using that by definition of a Whitney partition, for every j ∈ N such that 2 −j ≥ δ, 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, 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  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 where in the last line one can argue as for (4.6). Summing over k and using (4.8) we find Step 3 (Global term). We claim that for δ = n −β with β ∈ 1 3 , 1 2 , The first and last terms are estimated in a similar way so we only estimate the first one.
On the quadratic random matching problem in two-dimensional domains By (4.5) and (4.7), we get 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 This proves which recalling the choice (5.4) of θ, δ and r concludes the proof of (5.8).

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 [46,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, has Hölder continuous density, uniformly positive and bounded.
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: for every probability density f . 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 δ Xi , λ m = 1 m m j=1 δ Yj . 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 EJP 27 (2022), paper 54.
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), and similarly with the usual concentration inequalities which eventually give Global term. We consider the second term in (6.4), for which we apply (6.1) obtaining 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) We bound from above, using the subadditivity inequality (2.5) and the triangle inequality For the first term, we use again subadditivity (2.5), having used (6.6) and (6.5) with λ m instead of µ n . For the second term in (6.7), we use 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 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 The inclusion S n ⊆ S n,m yields immediately so that it always holds, for any sequence m = m(n) ≥ n, EJP 27 (2022), paper 54.
To show the converse inequality, given σ ∈ S n,m we induce a matching between Since Ω is bounded, we obtain the inequality Taking expectation and using the assumption (m − n)/ log n → 0, which in particular gives log n/ log m → 1 yields 8 Proof of Corollary 1.4 ⊆ R 2 , we introduce the notation for the costs of the Euclidean travelling salesperson problem and for 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 inequality Using Theorem 1.2, we have then lim inf 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 ) , with respect to the quadratic cost. We let τ = σ • τ ∈ S n , so that 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}, Summing upon i gives The thesis follows letting ε → 0.
To prove (8.2), we rely on the known (deterministic) bound 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 [14], 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.

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 U δ = {Ω k } k with Ω k open, diam(Ω k ) δ and |Ω k | ∼ δ d .
We start by constructing a partition of the set A δ = {d(·, Ω c ) < √ dδ} by open sets Ω k satisfying diam( Ω k ) δ, | Ω k | ∼ δ d . For this we first consider a δ−net {x k } k of ∂Ω, that is a family satisfying ∂Ω ⊂ ∪ i B δ (x i ) and min i =j |x i − x j | ≥ δ 2 . Such a family can be for instance constructed by choosing any starting point x 1 ∈ ∂Ω and then setting x k ∈ argmax ∂Ω d(x, ∪ k−1 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 − x k | = |x − x i |} is a hyperplane and thus of Lebesgue measure zero, { Ω k } k is indeed a partition of A δ in disjoint open sets. Notice also that by definition of a Whitney partition, if Q k ∈ V δ then Q k ∩ A δ = ∅. Now on the one hand, by triangle inequality it is immediate that diam( Ω k ) δ and thus also | Ω k | δ d . On the other hand, still by triangle inequality, B δ 4 (x k ) ∩ Ω ⊂ Ω k so that by Lipschitz regularity of Ω, we also obtain | Ω k | δ d .
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 δ.

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 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 pi 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 [22,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].