Random assignment problems on ${2d}$ manifolds

We consider the assignment problem between two sets of $N$ random points on a smooth, two-dimensional manifold $\Omega$ of unit area. It is known that the average cost scales as $E_{\Omega}(N)\sim\frac{1}{2\pi}\ln N$ with a correction that is at most of order $\sqrt{\ln N\ln\ln N}$. In this paper, we show that, within the linearization approximation of the field-theoretical formulation of the problem, the first $\Omega$-dependent correction is on the constant term, and can be exactly computed from the spectrum of the Laplace--Beltrami operator on $\Omega$. We perform the explicit calculation of this constant for various families of surfaces, and compare our predictions with extensive numerics.


Introduction
The Euclidean assignment problem is a transportation problem between a set X of N "red" points and a set Y of N "blue" points. Both sets are supposed to be on a given n-dimensional Riemannian manifold Ω. A transportation map is a bijective map T : X → Y, that is, a pairing among the red and blue points. A transportation cost wxy is given for each pair (x, y) ∈ X × Y, and the cost of the map T is the expression Date: May 14, 2021.
We will also define E Ω [X, Y] as the minimum cost among the possible transportation maps. Given a probability distribution for X and a probability distribution for Y, we use E Ω (N ) as a shortcut for E[E Ω [X, Y]]. For definiteness, we will assume in this paper that the points of X and Y are uniformely distributed.
Within this geometric framework, it is natural to choose for wxy a function of the distance d(x, y) between x ∈ X and y ∈ Y, that is wxy = f (d(x, y)). We expect that, if the function f has some natural (monotonicity, smoothness,. . . ) properties, the large-N behaviour of E Ω (N ) (with the volume of Ω kept constant) is dominated by the behaviour of f near zero. In turn, this suggests to concentrate only on power-law functions, f (d) = d p for some p > 0, as any other detail of the function is either trivially rescaled, or washed out in the limit. The two cases most studied in the literature are p = 1 and p = 2, where a number of (different) useful extra features emerge. This paper makes no exception, and in fact we will only consider here the case p = 2, that is, we will set once and for all wxy = d 2 (x, y).
It is a longstanding question to understand the asymptotic behaviour, for large N , of E Ω (N ), and, when n ≥ 2, the results are very partial for any manifold Ω, including the conceptually simplest ones (like the unit hypercube, or the unit hypertorus), and any value of p, including the special cases p = 1 and p = 2. In particular, in the two-dimensional case for p = 2 (see [1,2,3,4]), it has been proved that, as long as Ω has unit volume, the leading term is Ω-independent: The main goal of the present paper is to show that, at least in the case n = 2, p = 2, the leading behaviour in N of E Ω (N ) which is Ω-dependent is a constant, which can be calculated exactly. More precisely, we are not able to establish a full perturbative expansion for E Ω (N ), up to corrections o(1), for any Ω. Nonetheless, for all pairs (Ω, Ω ), we predict that where {λi} i≥1 is the set of eigenvalues of the Laplace-Beltrami operator on Ω that are different from zero (if Ω is a manifold with a boundary ∂Ω of perimeter of order 1, it is the set of eigenvalues of the Laplace-Beltrami operator, with Neumann boundary conditions). That is, all terms in an asymptotic expansion of E Ω (N ) which do not decrease with N must be 'universal'.
One can notice that K Ω is a regularization of the trace of the Laplace-Beltrami operator. Another equivalent regularization is the so-called Robin mass R Ω , see for instance [5,6]. In particular Eq. (1.2) can be equivalently written as (1.4) E Ω (N ) − E Ω (N ) = 2(R Ω − R Ω ) + o (1) The definition of R Ω is postponed to Section 3, while its relation with K Ω is described in Section 4. Let us put this result in context, by summarising (part of) the state of the art for this problem. In contrast with the transportation problem for continuous measures, in this case the candidate optimal transportation maps are just the N ! permutations, that is, the possible bijections between two sets of cardinality N , and in particular they are a finite set. For one given choice of the N 2 weights wxy, the computational problem of finding one optimal map T , 1 and the associated cost E Ω [X, Y] is a well-studied problem, which turns out to be in the polynomial class [7,8,9]. Thus, the associated computational problem can be quickly solved.
This fact is in striking contrast with the problem in Probability Theory, of understanding the asymptotic of the average cost, on various domains Ω and statistical ensembles for the red and blue point processes. This random version of the problem has attracted much attention both in Mathematics and Physics. In the Physics community, the interest has come from the analogy with 'spin glasses' in Statistical Mechanics, and a seminal contribution was given in the eighties by Orland [10], Mézard and Parisi [ 11], that considered the problem "in infinite dimension", by introducing the so-called "random-link" approximation. This version of the problem was addressed using (non-rigorous) statistical physics techniques such as the replica theory and the cavity method [12]. Their original results were later put on rigorous ground [13,14,15]. The extension to finite-dimension of the random-link results is, however, quite challenging. A first attempt was carried on by Mézard and Parisi [16,17] that showed how, for n > 2, the random-link result can be used as a zero-order approximation for the finite-dimensional solution, adding perturbatively a series of corrections. In the same years, a remarkable result was obtained by Ajtai and coworkers [18] for n = 2: they proved that, if the problem is considered on the unit square Ω ≡ R := [0, 1] 2 , then E Ω (N ) = O(ln N ). 2 Recently, the forementioned result has been refined. In particular, by means of non-rigorous arguments, in Refs. [1,2] it was claimed that, on the unit square R, where c PP R (N ) = o(ln N ) (the factor 2 is for later convenience). This result has been later rigorously proved by Ambrosio and coworkers [3] and extended to any 2-dimensional compact closed manifold Ω [4], where it is shown that (1.6) E Ω (N ) = 1 2π ln N + 2c PP Ω (N ).
The latter paper also proves rigorous bounds for c PP Ω , namely that c PP Ω (N ) = O( √ ln N ln ln N ). It has been recently conjectured that Eq. (1.6) holds also in the case of points generated from non-uniform densities [19].
In this paper we further investigate the problem of the estimation of c Ω (N ). Extending the arguments given in [2], we argue that, on a generic two-dimensional manifold of unit area, the correction c Ω in Eq. (1.6) can be written as  Analogous claims and results hold for other variants of the problem, most notably when one set of points is still sampled with the Poisson random process, and the other set is either a deterministic regular grid (we investigate here the cases of square (S), triangular (T) or "Fibonacci" (F) [20,21] grids), or, in the variant of the problem where T is the transportation between a discrete and a continuous measure, the uniform measure (U) on Ω. In these three new cases, the factor 2 in equation (1.6) disappears, and we have the similar structure Let us stress again that the functions c * are 'universal', in the sense that they do not depend on the choice of manifold Ω (but they do depend on the choice of local randomness, e.g. among P, S, T, F, U), while the geometric correction K Ω depends on the choice of manifold, but is 'universal' in a different sense, as it is independent of the choice of local randomness (provided that the extra factor 2 in the P case is taken into account). Just as well as equation (1.2), such a decomposition is not at all granted a priori, and is somewhat surprising.
We also give numerical estimates of the associated values of c * , 3 under the hypothesis that they are indeed constant, namely The paper is organized as follows. In Section 2 we define the random matching problems we are interested in. In Section 3 we present our functional approach for the derivation of the scaling of the optimal cost, including the finite-size corrections given in Eq. (1.7). For simplicity, we concentrate only on the Poisson-Poisson case. In Section 5 we apply our theory to different domains, giving an explicit computation of K Ω for all of them. In Section 6 we compare our predictions with numerical results obtained solving a large number of instances of the problem on the domains under investigation. Finally, in Section 7 we give our conclusions.

The random assignment problem
Let us consider a connected, two-dimensional smooth Riemannian manifold Ω having finite volume and, if with a non-empty boundary, finite perimeter, with metric g. Given a system of local coordinates (x 1 , x 2 ) around a point p ∈ Ω, g = ij gij(p) dx i ⊗ dx j , and given two elements v, w in the tangent bundle in p, we will denote by v, w p := ij gij(p)viwj. For the sake of generality, we will perform our analysis in the slightly more subtle case of ∂Ω = ∅ (the arguments below can be easily adapted to the case ∂Ω = ∅). We will denote dσ the Riemannian measure for Ω, and we will assume the measure of Ω to be equal to 1. 4 Moreover, we will denote by δp(x) dσ(x) the unit measure concentrated in p ∈ Ω, that is, given a test function ϕ(x), Suppose now that two sets of points are given on Ω, namely a set of N points X := {Xi} N i=1 ⊂ Ω, that we call the red points, and a set of N points Y := {Yi} N i=1 ⊂ Ω, that we call the blue points. The assignment problem consists in assigning each red point to one blue point, in such a way that the resulting map T is a bijection, and a certain total cost function is minimized. As motivated in the introduction, the cost function is the sum of the costs for each pair (Xi, Yj) such that T (Xi) = Yj, and the cost of a pair (x, y) is the square of the Riemann distance between x and y of the selected pairs. In formulas, we have to find the optimal bijection T * : X → Y such that and d(x, y) is the Riemann distance between the points x and y, i.e., the infimum of the lengths of the curves that join the two points. Note that each feasible T corresponds to a permutation π of N elements, so that T (Xi) = Y π(i) , and searching for the optimal map is equivalent to searching for the optimal permutation. If we introduce the two atomic measures the optimal cost minT EN [T |X, Y] coincides with the 2-Wasserstein distance (squared) between the two empirical measures in Eq. (2.3), of which we shall recall here the definition (see e.g. [22]): given two probability measures ν1 and ν2 on Ω, their 2-Wasserstein distance is where the infimum is taken over the set Γ(ν1, ν2) of all the joint probability distributions J with first and second marginal given by ν1 and ν2, respectively. It is well known (see for example [3,23,24,25]) that, in our setting, the set of the optimal joint probability distributions J is a convex polytope, called Birkhoff polytope, whose extreme points are all and only the permutations π which are optimal within the probability distributions of the form J(x, y) = i δX i (x)δY π(i) (y). Accordingly, the set of optimal maps T : Ω → Ω pushing ν1 to ν2, i.e., those realising the infimum in the expression coincides with the set of maps of the form T (Xi) = Y π(i) , with π optimal in the sense above. (The situation is much simpler when ν X is absolutely continuous with respect to the Riemannian measure, as in this case the optimal transport map T would be unique. For a more complete discussion see also [24, ch. 9]). The distance in Eq. (2.5) corresponds, up to a multiplicative constant, to the cost in Eq. (2.2) when ν1 ≡ ν X and ν2 ≡ ν Y . Therefore In the following, we will consider various statistical ensembles of pairs (X, Y). At this point, many choices are possible. To be definite, we will always choose X and Y to be independent of each other. We will also choose Y to be always what we shall call, with abuse of notation, the "Poisson random process on Ω of size N ", that is, the Yj's are i.i.d., uniformly chosen on Ω (w.r.t. the measure dσ). 5

Poisson (P)
: Also X is given by a Poisson random process on Ω of size N . 5 The 'genuine' Poisson random process on Ω is defined by an intensity, not by a size. The process of intensity N dσ produces configurations Y in which the number of points is a Poissonian random variable of average N . However, in our context, the large-N convergence of the local properties of the fixed-size Poisson process to the ones of the genuine Poisson process is fast enough to justify our abuse of language.

Uniform (U):
Together with the Poisson-Poisson, the most general and interesting case is the Uniform-Poisson case, in which the cost is the distance beetwen the Poisson random process and the uniform measure dσ: . As a discrete approximation of this case, we can introduce various grid-Poisson assignment problem (GP), interesting by themselves: Square grid (S): when Ω is a flat a × b rectangular domain (possibly up to identification of the boundaries, e.g. as in a torus), and there exists a value k such that ka, kb ∈ N and k 2 ab = N , then a natural choice is to fix X to be the square grid of spacing 1/k. In the case of a torus, we can imagine identifying the horizontal sides of the fundamental rectangular region with a shift s. In this case the grid has no local defects when also ks ∈ N, and the modular parameter of the resulting surface is τ = (s + ib)/a, so that the set of points in the moduli space which can be realised by a grid with cardinality between N and N + N 1 2 + becomes dense everywhere in the limit of large N . Triangular grid (T): analogous to the square grid, in the case in which Ω is a flat hexagon (possibly up to identification of the boundaries, e.g. as in a torus), with sizes (a, b, c, a, b, c) in cyclic order. Of course, this includes as special cases the regular triangle and hexagon, and the rhombus of angle π/3. Now we require that there exists a value k such that ka, kb, kc ∈ N, and √ 3 2 k 2 (ab + bc + ca) = N , and the natural choice is to fix X to be the triangular grid of spacing 1/k. In the case of the torus, calling ω = e 2πi/3 , the associated modular parameter is τ = ωb+ω 2 c a+ω 2 c , so that also in this case the whole moduli space can be accessed by increasing N . Fibonacci grid (F): This is a subtle construction, adapted to the case of a sphere, see Fig. 1b, and based on stereographic projection from a Fibonacci spiral on the plane (from which the name), described in [20,21]. The local aspect of this grid around one given point is somewhat intermediate between the one of a square and of a triangular lattice, with variations depending on the spherical coordinates of the point, and on the precise value of N . We will not enter in the detail of this construction, and the reader is referred to the forementioned papers.
In Appendix A we give more details about the relation between grid-Poisson assignment problems and the UP problem, with an estimation of the convergence rate of the optimal cost in the former to the optimal cost in the latter. As anticipated, we are interested in the study of the asymptotic behaviour in N of the average optimal transportation cost, for which we will adopt the general notation where the average E [·] is taken over the pertinent statistical ensemble for the point processes.

Main conjecture
As we said above, the study of E Ω (N ) in any dimension n = 1 or ∞ (i.e., in the random-link model) seems rather difficult. A possible approach to the study of the asymptotic behaviour for large N is based on the fact that, in this limit, T (x) is expected to be 'close' to the identity map (more precisely, from [18] we expect that d(Xi, Yj) = O( (ln N )/N ) for pairs of points which are paired by an optimal matching), and an expansion in this small parameter, at the first non-trivial order, might still capture the relevant features of the solution of the problem. In [1,26,2] this approach has been applied to the study of the problem when Ω is the square, or the torus with modular parameter τ = i. The analysis leads to a result whose interpretation requires a regularization that takes into account the finite-N effects and avoid divergences, as we will see below.
In the approach in [1,2,26], a close analogy naturally emerges between the evaluation of the average optimal cost in the assignment problem and the evaluation of the electrostatic energy of 2N particles, N of each charge sign, pinned in random positions on Ω. This is a result a posteriori of the theory, as the obvious analogy just doesn't hold as is (in the electrostatic problem, the energy is the sum of N (2N − 1) pair contributions, not just N , which scale logarithmically with the distance of the pair, instead that quadratically). In a sense, the proposed linearization follows the opposite track of the suggestion by Born and Infeld [27] of a non-linear version of electrodynamics in order to solve the problem of divergencies. Similar ideas have been proposed recently by Brenier for fluid motion [28,29].
3.1. Linearization. Let us review the arguments of [2], in their natural generalisation to a Riemannian manifold. We start by introducing, for each map T : Ω → Ω, the cost Note that, at this stage, T is not a transportation map. First, because we haven't still imposed the fact that the push-forward of ν X is ν Y , and second, as specific to our transportation problem dealing with atomic measures, the 'true' optimal transportation map is only defined on the support of ν X , which is not the whole Ω. We start by solving the first issue. An equivalent formulation of the constraint is that, for any function φ : Ω → R, we must have Again, it would be enough to consider functions φ with the same support of ν Y , a fact of a certain relevance as it implies that, by expanding φ over the appropriate basis of functions, we have only N − 1 independent constraints, instead that infinitely many, as it would be the case if ν Y (x) were absolutely continuous with respect to the Lebesgue measure. The idea is now to write down a Lagrangian that combines the cost expression in Eq. (2.2) with the condition in Eq. (3.2) as where φ plays the role of a Lagrange multiplier. The optimal map T * satisfies the Euler-Lagrange equations obtained from the Lagrangian above (which turn out to be nonlinear).
We shall now use the idea that, for N → +∞, we expect T (x) → x for any x ∈ Ω, due to the fact that the matched pairs become infinitesimally close under the scaling in which |Ω| is kept fixed. Then, there exists a vector field µ(x) on Ω such that, at the leading order The direction of the field is the one of the geodesic curve realising the distance of T (x) from x. Pictorially We shall now introduce which is another perturbative parameter (when averaging over our statistical ensemble, monomials E[δν(x1)δν(x2) · · · δν(x k )] have a definite scaling with N , and high powers are suppressed). The Lagrangian is approximated, in this limit, by its quadratic version, We have also used the fact that, if the integrand is smooth enough, we can neglect the discrepancy between dν X and dσ (while still treating more carefully δν(x)). Extremizing the new Lagrangian, and using that if ∂Ω = ∅ the field µ is tangent to ∂Ω, we obtain the non-homogeneous linear equations which are the linearization of the original Euler-Lagrange equations in the fields µ and φ. In local coordinates: where ∂i ≡ ∂ x i , the tensor g ij is the inverse of g and |g| = det g. The two equations imply the Poisson equation to be solved with Neumann boundary conditions, since the flux of µ(x) at the boundary is zero. Here −∆ is the Laplace-Beltrami operator on Ω, i.e., in local coordinates 3.2. The divergence of the cost and the problem of regularization. The functional approach above tells us that µ = −∇φ, where −∆φ = δν. We can use the fact that 6 to write down an expression, valid for N 1, for a quantity (x) that we shall call the cost density Here we have introduced the Green function G(x, y) of −∆ on the orthogonal complement of the locally constant functions. The Green function is a symmetric function that satisfies the equations where ∂nG(x, y)| y∈∂Ω is the normal derivative in x with respect to the boundary ∂Ω of the domain. The equations above identify a unique Green function up to an additive constant: we will fix this constant adopting the convention The obtained results have, however, a fundamental problem. The quantity in Eq. (3.13) is divergent for any x ∈ Ω. The responsibility for this fact comes from several sources, one of which is having treated the field µ(x) as a continuous field, instead that a collection of N vectors, one per each point Xi ∈ X. This gives locally, in coordinates on the tangent space in Xi, a field whereμ(x) is such thatμ(Xi) is a finite quantity (here we have used the diagonal expression for the Green function G, see below Eq. (3.17)). Such an approximation could still be used through a delicate Cesàro limit, if we had to perform integrals in which µ appears linearly. However our cost density is quadratic in these fields, and locally, at a formal level, for δ 1, that is, the appropriate result, which depends on the positions of the points and is finite at finite N , is shifted by a fixed but divergent quantity. Yet again, we observe a perfect analogy with 2-dimensional electrostatics, namely with the classical problem of the field self-energy for a distribution of point charges.
Analytically, we see this feature emerging from our result by observing that for d(x, y) → 0, the Green function behaves as [5,6] with a logarithmic divergence. We perform therefore a regularization of the logarithmic divergence, along the same lines of the classical treatment of electrostatics. Let us introduce Ω δ (x) = Ω \ B δ (x), where B δ (x) = {y ∈ Ω : d(x, y) < δ} is the ball of radius 0 < δ 1 centered in x. We can introduce a regularized expression and a corresponding "regularized cost" where the second integral runs over the border ∂B δ (y) := {x ∈ Ω : d(x, y) = δ} of B δ (y), dλ(x) is the line element of ∂B δ (y) in x, and n is the outward normal to B δ . By Eq. (3.14) and Eq. (3.14c) the first integral is infinitesimally small for δ → 0. Therefore (3.20) E For 0 < δ 1, the inner integral can be estimated using the expression in Eq. (3.17), so that We finally get The integral of m(x) is sometimes called Robin mass [5,6]. Now, we suppose that the regularization by the parameter δ acts in the same way on all geometries. Under this assumption we can compare two different geometries, Ω and Ω , obtaining the following conjecture.

Conjecture.
Let Ω, Ω be two regular two-dimensional manifolds, then In other words, the differences of the average cost among different manifolds, in the large-N limit, are expected to be regularization-independent, and, in addition to this, can be expressed in terms of the Robin's masses of the Laplace-Beltrami Green function on Ω and Ω . The analytic evaluation of these differences will be the main object of our investigation, starting from Section 5. The remaining of this section is instead devoted to a further justification of the assumption at the basis of equation (3.24).
One problem at this point is that our regularization parameter δ does not have a clear relation with the perturbative parameter N −1 . In order to better understand what is the microscopic mechanism beyond the regularization, we observe that equation (3.22) can be formally written for δ → 0 as where the operator −∆ −1 is the inverse Laplace-Beltrami operator on Ω (with Neumann boundary conditions, if the boundary exists) As said above, a logarithmic divergence appears for δ → 0 and both sides of Eq. (3.25) are infinite. By the Weyl law on the asymptotics of the eigenvalue counting function N Ω (λ) for the Laplace-Beltrami operator [30] we know that, for a 2-dimensional manifold with unit volume, and under Neumann boundary conditions, the leading behaviour of N Ω (λ) for large λ is 7 Furthermore, the eigenfunction f λ associated to a given value of λ 'looks locally' like a plane wave with wavelength 1/ √ λ. This has two consequences at the level of our approximations when passing from the complete Lagrangian, equation expansion of φ(T (x)) = φ(x + µ(x)) around x, in the basis of the eigenfunctions {f λ }, is perturbative in the parameter E(µ(x)) √ λ, which we expect, from [18], to be of order λ ln N/N . Second, if our basis is orthonormal for the measure dσ, that is, (f λ , fρ) dσ = δ λρ (assuming for simplicity of notation that the spectrum is non-degenerate), under the measure dν X we get instead The result of this analysis is that, if we decompose our fields in the basis of eigenfunctions of −∆, we can neglect the corrections coming from the further terms in the Taylor expansion, and the discreteness of the measure, only for those eigenfunctions with λ N (up to possible factors (ln N ) γ in the scaling). Conversely, in the regime λ N some unknown mechanism comes into play, and we expect that its effect is to dump the sum tr ∆ −1 appearing in (3.25), possibly at a scale λ N . In [1], this unknown dumping mechanism is supposed to be encoded in a cut-off function F (λ/N ). Now, as this function is related to the local expansion of the fields µ and φ at high frequencies, and as the relation between eigenvalue λ and local wavelength is universal, the function F (λ/N ) must have one of the two flavours of universality: it shall not depend on the manifold Ω, while in general it must depend on the type of problem (among Poisson, various grids and uniform), that is, in a natural generalisation of the treatment of [1] to our setting, we should have some unknown functions F •P (λ/N ), with • being one among P, S, T, F or U, and no dependence on Ω. Going on with our analysis of the PP case (the reasoning can be repeated for all other cases similarly) and using F (λ/N ) for F PP (λ/N ) for brevity, within the assumptions of [1] we should interpret the correspondence (3.25) above as where Λ(Ω) is the set of nonzero eigenvalues of the Laplace-Beltrami operator on Ω. Following the analysis already performed in [1], this gives for any domain of unit measure, for some constant c Ω depending on the cut-off, which cannot be determined if F is not known. We recall that, as anticipated in the introduction, the leading term in Eq. (3.30) is the correct asymptotic cost, as rigorously proved in Refs. [3,4,33], the presence of a logarithm being known since the eighties [18]. Note that there is no guarantee that the cut-off function scales exactly as F (λ/N ), as the mechanism beyond the dumping of the high-wavelength contributions, and the amount of this dumping, are not under control. It may well be, for example, that the function has the form F λ N (ln N ) γ , which would give a variant of (3.30) in which instead of the constant term it will appear an universal term O(γ ln ln N ).
All these arguments lead us to reformulate our conjecture as follows.

Conjecture (Alternative formulation).
Let Ω be a regular two-dimensional manifold, then Moreover, for Ω, Ω different regular manifolds,

Remark 1.
It is important to remark that, from the simulations made in Section 4, we have evidence that the term c * (N ) can be chosen as a constant, i.e. that formula (3.30) is compatible with our numerical results. Obviously it is not possibile to deduce that c * (N ) = c * from numerical simulation only. Anyway, in case, we would get (3.31) and c Ω = R Ω + c * .

Remark 2.
We notice that starting from Eq. (3.29) and comparing two manifolds, we get the interesting fact The combination of the two integrals above can be rewritten as The universality of Weyl law implies that the factor N Ω (λ) − N Ω (λ) grows no faster than √ λ ln λ, so that, even in absence of the function F (that is, in the limit of N large), the integral is convergent at infinity (and near zero is regularised by the spectral gap). This allows us to predict (3.35) lim

Different regularization procedures
The expression (3.29) is, annoyingly, a diverging expression depending on Ω. A way of studying this expression is by introducing a regularization parameter for these contributions, and then deducing an evaluation of (3.35) from a singular expansion in around zero.
One standard way to perform this programme is the so-called zeta regularization [34]. Let us introduce the generating function  For reasons that will appear clearer below, we will call Kronecker's mass the constant K Ω . Despite the fact that there seems to be no reason a priori to believe that K Ω and R Ω are related, it has been proved by Morpurgo that R Ω − K Ω is a universal constant (that is, it does not depend on Ω), given by [5,36,37,38] (4.4) where γ E = 0.57721 . . . is the Euler-Mascheroni constant. In particular, this universality result is crucial in checking a posteriori that our two predictions (3.24) and (4.3), obtained by two different analyses, are consistent, and also implies that our Conjecture is equivalent to the statement of E. (3.35). The computation of the Kronecker's mass is often easier than the Robin's mass, as we will show below. For a few manifolds Ω, both computations, of R Ω and K Ω , can be performed with relatively small effort, and we will do this, for pedagogical reasons, in order to illustrate with an example the forementioned general result.
Another way of performing our programme is to consider the regularized sum for p > 0. This corresponds to a specific choice of function F (λ/N ), provided that is identified with γ /N p , for γ some constant. Also in this case we have, universally, and this leads to the prediction The analogue of the Morpurgo theorem reads in this case as can be evinced by comparing the two regularizations for the diverging integral 1 4π ∞ 1 dx x (which, besides the fact that it is an integral rather than a sum, it has all the appropriate asymptotics properties implied by the Weyl law). Namely, for this choice we have Another regularization in the same spirit is through the regularized sums with θ(x) is the Heaviside step function, which, yet again, corresponds to a specific choice of function F (λ/N ), provided that is identified with γ/N , for γ some constant. Also in this case we have a universal asymptotics

Examples
To verify our ansatz, we will compute the Kronecker's mass and the Robin's mass for different Ω. We will compare our analytic results with numerical simulations in Section 6. We will start considering flat manifolds having g(x) = I, and consider manifolds with uniform curvature starting from Section 5.5.
5.1. The unit rectangle. Let us start by considering the problem on the rectangle. We call R(ρ) the rectangle [0, √ ρ] × [0, 1/ √ ρ], and we consider the Laplace-Beltrami operator with Neumann boundary conditions. The eigenfunctions of −∆ on R(ρ) are given by The corresponding eigenvalues are We proceed computing the Kronecker mass using the regularized function Here we have adopted τ = iρ, in compliance with standard notation for modular forms, and have introduced the lattice zeta function ζτ (s) defined in D. This calculation is readily performed thanks to a remarkable result due to Kronecker (and known as first limit formula of Kronecker), reported in D, equation (D. 7), and that we repeat here: where η(τ ) is the Dedekind η function. Kronecker's formula allows us to immediately obtain 8 that for ρ = 1 (unit square) simplifies to (5.6) We will see in the following that the first limit formula of Kronecker will allow us to extract the Kronecker's mass for many types of flat domains: this explains our choice of 'Kronecker mass' for denoting K Ω . We can give a slightly more compact form to the function in (5.5), shifted by its minimal value above: We shall make a remark on this expression. In the limit in which the rectangle is very elongated, we get 1 8 We recall here that √ ρη(iρ) = η( i /ρ). 1 The thumb rule in performing these limits is that, for ρ → 0, +∞, that is the average cost for the Poisson-Poisson one-dimensional assignment problem on the segment of unit length [39]. This is not by accident. Indeed, in any rectangular domain we can evaluate the average energy of the permutation in which the k-th red point counting from the left is matched to the k-th blue point counting from the left. This configuration is optimal w.h.p. in the limit ρ → +∞, and would be optimal, at any ρ, if the vertical coordinates of all the points were equal. On the other side, a worst case is when all the vertical coordinates of red points are zero, and all the vertical coordinates of blue points are 1/ √ ρ, so that, calling E [0,1] (N ) the average energy for the 1-dimensional problem on the [0, 1] segment, we get which, by substituting our scaling ansatz, gives When we take a limit N → ∞, ρ → ∞ on a direction ρ √ N , we thus get which is consistent with (5.8). In the following we will encounter various other domains which allow a consistency check with a 1-dimensional limit. We will reach similar conclusions, without entering in the details of the estimates, as this is done by minor modifications of the reasonings presented here.

The flat torus.
We shall now consider the problem on the flat torus T(τ ). To describe the corresponding manifold, let us first consider the lattice of points on R 2 (5.12) Λ = ω n, n ∈ Z 2 generated by the matrix corresponding to the base vectors (5.14) ω1 := 0 , ω2 := s h .
In such lattice it is possible to define fundamental parallelograms D(ω), containing no further lattice points in its interior or boundary. A fundamental parallelogram is given for example by Fig. 2a. We will also use a shortcut adapted to rectangles, A torus T is defined as a quotient between the complex plane and a lattice Λ, T := R 2 /Λ. In other words, each point x ∈ D is identified with the set of points {x + ω n, n ∈ Z 2 }, the distance between two points in D being the minimum distance between the elements of their equivalence classes. It is well known that two matrices ω and ω identify the same lattice Λ and the same torus T (although not the same fundamental domain D) if and only if (ω) −1 ω ∈ SL(2, Z). For each ω, we introduce the half-period ratio Given a lattice Λ generated by ω, it is possible to associate to it a dual lattice Λ * generated by ω * , such that ω * ω T = I, identity matrix, i.e., Each torus T = R 2 /Λ is then naturally associated to a dual torus given by T * := R 2 /Λ * . In the following, we will restrict, without loss of generality, to the case in which the fundamental parallelograms have unit area, choosing such that ρ ∈ R + and σ ∈ R, and we will denote the corresponding torus by T(τ ), where τ := σ + iρ is the half-period ratio. The Kronecker's mass. Due to the periodicity conditions, the eigenfunctions of −∆ on T(τ ) have the form for all k * = ω * k ∈ Λ * , k = n −m ∈ Z 2 . The corresponding eigenvalue is (5.21) λ (n,m) = |2πk * | 2 = (2π) 2 |n + τ m| 2 ρ = (2π) 2 |n + τ m| 2 (τ ) .
We can compute now the Kronecker mass using the regularized function and removing the pole in s → 1, as discussed in Section 3. This calculation is readily performed, again thanks to the first limit formula of Kronecker, equation (D.7), which allows us to immediately obtain In Fig. 2b we present a contour-plot of the related expression (τ )|η(τ )| 4 in the complex plane τ , confined to the canonical fundamental region of the moduli space. In particular, this function diverges for τ → 0 and has minimum at τ = exp(iπ( 1 2 ± 1 6 )). This implies that, among all unit tori equipped with the flat metric, the "hexagonal" one, that is the one for which τ is a sixth root of unity, is the one in which the average cost of the Euclidean Random Assignment Problem is minimised. More strikingly, as deduced from results in [5] which are in turn based on the results in [35], the hexagonal torus is minimal also among unit surfaces with non-uniform metric. Example: the rectangular and rhomboidal tori. We shall call "rectangular torus" a torus in which the fundamental parallelogram is a rectangle. This case corresponds to τ = iρ, with ρ > 0 real. Our formula specialises to which is invariant under the map ρ → ρ −1 , as it should. In the region ρ ∈ (0, 1] the lowest value is achieved at ρ = 1 (see also Fig. 2b), where (5.25) We can give a slightly more compact form to this function, shifted by its minimal value: Similarly, we shall call "rhomboidal torus" a torus in which the fundamental parallelogram is a rhombus. This case corresponds to τ = e iθ , with 0 < θ ≤ π /2, and our formula specialises to that is, again shifting by the value for the standard torus, As was the case for the rectangle, the expression in equation (5.26), in the limit in which the torus is very "thin and long", becomes This happens to be the average cost for the Poisson-Poisson one-dimensional assignment problem on the circle of unit length [39], as was to be expected, by a reasoning analogous to the one presented for the case of the rectangle. The Robin mass. Let us now evaluate, for the generic flat torus T(τ ), the Robin mass R T (τ ). Calling z = z(x, y) = (x1 − y1) + i(x2 − y2), the Green's function on the torus is given in this case by [40] (5.30) where θ1(z; τ ) is an elliptic θ function. The Robin mass is obtained from It is immediately seen that, in agreement with the Morpurgo theorem, equation (4.4) is satisfied.

Other boundary conditions on the unit rectangle.
The unit rectangle and the rectangular torus are obtained starting from the fundamental domain D(ρ) in equation (5.16), and assuming respectively open (i.e., Neumann for the field φ) and periodic boundary conditions. Other choices of boundary conditions are possible, which correspond to other classical surfaces, with or without boundary. Each choice leads to a different spectrum of the Laplacian, which in turn implies a different Kronecker mass (and, according to our theory, a different finite-size correction to the optimal cost of the assignment problem).  The corresponding eigenvalues are therefore Repeating the same type of calculations performed for the rectangle (that is, expressing the regularised sum as a combination of ζτ (s) (for some τ 's) and ζ(2s)), we obtain We also remark that which is the cost density for the one-dimensional assignment problem with open boundary conditions (i.e., on the unit segment), while which is the density of cost for the one-dimensional assignment problem with periodic boundary conditions (i.e., on the unit circle), again, as was to be expected. The nontrivial solution of the equation K C (ρ) = K C is ρ = 0.625352 . . . , while the minimum value of the mass occurs for ρ = 0.793439 . . . Remark that the constants above do not appear in the study of [35], because in our context, in presence of a boundary, we should impose Neumann boundary conditions (while the authors of [35] only analyse the case of Dirichlet boundary conditions).
The parity constraint implements the twisted identification of the strip. The corresponding eigenvalues are Repeating now the usual arguments we get We also remark that which is the average cost for the problem on the segment of length 1 2 (and the fact that the length is not 1 is related to the fact that the twisted boundary conditions are effectively 'folding in two' the segment), while Proceeding as above, one can obtain so that, in particular, We also remark that both one-dimensional limits coincide with the corresponding constructions for the Möbius strip, and indeed the limits of the analytical expressions are the same, as Here K K (ρ) = K K for ρ = 1.09673 . . . , whereas the minimum is obtained at ρ = 1.04689 . . . . The calculation proceeds as in the other cases, giving which is symmetric for ρ ↔ 1 ρ , as it must be. In particular Now, both one-dimensional limits produce a domain corresponding to the segment of length 1 2 , and indeed  Figure 5. Absolute difference of average optimal costs for the assignment on the cylinder C(ρ), on the Möbius strip M(ρ) and on the Klein bottle K(ρ) with the corresponding costs for ρ = 1. The numerical results, represented by the dots, are compared with the analytical prediction obtained from Kronecker's masses.

The disc and the cone.
Up to now, we have mostly solved the problem using the zeta regularization of the Laplacian, and relying on Kronecker's first limit formula. Only for the case of the torus, we have also performed the calculation of the Robin mass, and verified the prediction of the Morpurgo theorem. In this section we will give the results for a geometry Ω in which the calculation of the Robin mass is done with relatively small effort, as the Green function can be guessed through the method of images, while the calculation of the Kronecker mass would require a sum over maxima and minima of Bessel functions. 2 Let us introduce the notation Dp(r) for the circular sector of radius r and angle 2π p , see Fig. 6a, The unit area condition implies 2πr 2 = p. We considered the case p ∈ N, and we choose periodic boundary conditions in the angular direction: this is equivalent to say that we identify the two radii of the sector, obtaining in this way a cone of height r 1 − p −2 , see Fig. 6b. This surface is interesting, as it is the first example in our list of a surface with singular curvature, the conical singularity being at the vertex of the cone. We have argued in the introduction that, because of the scaling ∼ √ N −1 ln N of the field µ, we expect that the same theory applies to the flat Euclidean space and to curved manifolds, as long as the curvature is non-singular. The case of surfaces with a finite number of conical singularities would require a different (although feasible) argument, and the verification of our theory on this family of surfaces (as well as the surfaces treated in Section 5.5) is an important validation of our predictions.
The Robin's mass for this case is obtained in Appendix B, and it is equal to where φ(z) is the digamma function. In particular, for α → 2π we recover the case of the unit disc D ≡ D1: The Kronecker's mass is readily obtained using equation (4.4). 5.5. The unit sphere and the real projective sphere. An example of transportation problem on the surface of the sphere S 2 has been considered in [42], where the problem of transporting a uniform mass distribution into a set of random points on S 2 is analyzed. As an example of applications of our approach to non-flat manifolds, here we consider the problem in our usual setting, i.e., a transportation between two atomic measures of random points. As in the previous cases, the information on the finitesize corrections is related to the spectrum of the Laplace-Beltrami operator on the manifold. It is well known that the eigenfunctions of −∆ on the surface of a sphere of radius r are the spherical harmonics Y l,m (θ, φ) with l ∈ N and m ∈ Z with −l ≤ m ≤ l. The corresponding eigenvalues are that is, the eigenvalue l(l+1) r 2 has multiplicity 2l + 1, for the range of integers −l ≤ m ≤ l. We fix the surface area of the sphere to 1 (taking r = (4π) − 1 /2 ) and we proceed using the zeta regularization, computing In this case, after some algebra, we are led to use 'just' the version of the zeta regularization for the Riemann zeta function (which is much simpler than the Kronecker formula) and we obtain The Kronecker mass for the unit sphere is therefore Alternatively, we can use one of the regularizations illustrated in Section 4. A convenient one is the evaluation of W ( 1 /2) , which gives Recalling that in our case r 2 = (4π) −1 , the final result is that, in light of (4.8), allows to rederive equation (5.61).
The spherical lune. The calculation above can be extended to the spherical lune S 2 k , a surface on a sphere of radius r, 4πr 2 = k, contained by two half great circles which meet at antipodal points with dihedral angle 2π k , see Fig. 7. In Appendix C it is shown that the Kronecker's mass corresponding to this manifold is Choosing periodic boundary conditions on the two half great circles, the Kronecker's mass takes an additional − 1 2π ln 2 contribution, and we obtain that reduces to (5.61) for k = 1, as it should. The projective sphere. A variation of the problem on the (unit) sphere is the problem on the (unit) real projective sphere PS 2 , that is, the sphere in which antipodal points are identified. The eigenfunctions of the Laplace-Beltrami operator are still the spherical harmonics Y l,m (θ, φ) with l ∈ N and m ∈ Z, −l ≤ m ≤ l, but we have to restrict ourselves to eigenfunctions that are invariant under the transformation (θ, π) → (π − θ, φ + π), i.e., to even values of l. Working on the unit-area manifold accounts to have 2πr 2 = 1. We get so that the Kronecker's mass is Using a sharp cut-off instead which, again by using (4.12) and (4.14), allows to rederive equation (5.67).

Numerical results
We have numerically investigated all the cases described above to check our predictions. To solve the assignment problem we have used an implementation of the Jonker-Volgenant algorithm [8]. For each domain Ω, we have computed the expected optimal cost averaging over at least 10 4 independent instances and different sizes N of the system, 32 ≤ N ≤ 1024. In each case, we have estimated c •P Ω (N ) assuming that they are indeed constant at the leading order, i.e., c •P Ω (N ) ≡ c •P Ω , via a least square regression. Here • = {P, S, T, F, U}. Our results are given in Table 1.
For each domain we have also computed c •P * = c •P Ω − K Ω that we expect to be domain-independent. In the PP case, our best estimation of c PP * , obtained for the PP problem on the torus (see Fig. 3a) is  Table 1. Kronecker mass and finite-size corrections c PP Ω evaluated by numerical simulations of random assignments on different domains. An estimation of c PP Ω − K Ω is also given. We also give our numerical results for c •P Ω , obtained performing random assignments on the same domains but fixing one set of points on a grid. The type of adopted grid is specified in the last column. Fig. 6c and in Fig. 7b we have also considered the differences of average optimal costs in the case of the circular sector and of the spherical lune. In all investigated cases we found a perfect agreement with our predictions.

Conclusions and perspectives
In this paper we have considered the assignment problem between two sets of random points on a generic two-dimensional smooth manifold of unit area. We have showed, by means of analytical arguments and numerical simulations, that the average optimal cost can be written as if one of the two sets is fixed on a grid (square, triangular, Fibonacci) or replaced with the uniform measure. In the equations above, K Ω is a precise quantity that can be obtained by a zeta-regularization of the trace of the inverse Laplace-Beltrami operator on Ω. The contributions c •P * are instead independent on Ω and related to the 'local details' of the problem (i.e., if the assignment is between random points, or with a grid, or with the uniform measure). We have given an exact computation of K Ω for different domains, and using different procedures.
The quantity c •P * (N ) shows a weak dependence on N (if no dependence at all): it has been proven indeed that c UP * = O( √ ln N ln ln N ) [4], a bound that, because of triangular inequality, holds for all the cases that we have considered. Assuming that c •P * (N ) are constants, we have verified, within our numerical precision, their independence on Ω in all considered transportation cases (Poisson-Poisson, grid-Poisson, uniform-Poisson).
Our results reduce the computation of the (leading) finite-size correction to the optimal cost to the computation of the Ω-independent contributions c •P * (N ). These contributions are intrinsically dependent on the local nature of the problem (and therefore change if, for example, we fix on a grid one of the two sets of points) and are inherited by the regularization of the highest part of the spectrum of −∆, as discussed in Section 3. What are the properties (and possibly the exact value, if they are constant) of c •P * (N ) remains an open question. Finally, analogous results are expected to hold in higher dimension at the leading order. In particular, the analysis in [1] suggests that, for d > 2, the local properties of the problem affect the coefficient of the leading term, with a finite-size correction depending on the spectrum of the Laplacian only. This case may also be investigable with our tools, as indeed versions of the Weyl law, which we crucially use in Section 3.2, exist in generic dimension. We leave the investigation of the higher dimensional problem for future works.
Acknowledgement. E. Caglioti and G. Sicuro would like to thank Giorgio Parisi for putting them in contact. D. Benedetto and E. Caglioti thanks Gabriele Mondello and Riccardo Salvati Manni for clarifying discussions about the case of the torus. The authors are grateful to Jürg Fröhlich for his careful reading of the manuscript. A. Sportiello is partially supported by the Agence Nationale de la Recherche, Grant Number ANR-18-CE40-0033 (ANR DIMERS).
Appendix A. Comparison of GP problem with UP problem on the flat torus As commented in the main text, the very same arguments presented for the PP case in Section 3 can be repeated for the GP case and the UP case, the only difference being that (3.12) has to be replaced by The result is an overall 1 2 factor, as shown in the final formula (1.8b). There is however no guarantee that c PP Ω (N ) = c •P Ω (N ) for • = {S, T, F, U} at fixed Ω as one might naively expect from (1.8), because these quantities depends on the regularizing function F •P (λ/N ), that, even assuming that it exists, is expected to be in general different in each case. In [3] it is proved that . As discussed in Section 3, one way to numerically estimate c UP * (N ) is to perform a transportation between two sets of points, supposing that one of them (e.g., the red ones) is fixed on a grid and not random. As intuitively expected, a grid approximation provides some information on c UP * . For example, let us consider the transportation between a set Y = {Yi}i=1,...,N of random points on the (flat) torus T and a set of hN = L 2 points X = {Xi} i=1,...,hN fixed on a square grid. This is not an assignment problem because the cardinality of the two sets is different. However, if h ∈ N the transport can be obtained as an assignment "replicating" each point in Y h times, see Fig. 1c, so that in the optimal configuration each point of the original set Y will be associated to h grid points. By classical convexity properties of the squared Wasserstein distance, it can be proved that, given hN = L 2 and considering a squared L × L grid on the unit flat 2-torus, where we denoted by a generic GP the correction in the discrete transportation problem (in particular, c GP T (N ) = c SP T (N ) for h = 1). By means of the GP problem, c UP * (N ) can be estimated for each N in the limit h → +∞.
To prove Eq. (A.4), first we notice that if a probability measure ν1 is the convex combination of measures νη, that is, if we can write where γ(η) is a probability measure, then, for any measure ν2 Let us assume now that the following empirical measure is given on the torus T, δX i concentrated on hN = L 2 points on the torus, L ∈ N, and let us assume that such points X := {Xi} i=1,...,hN lie on a square grid of step L −1 . Let us denote also with Y = {Yi}i=1,...N a set of N points sampled from the uniform distribution, and with ν Y their empirical density. Noticing that the uniform measure on the torus is a convex combination of grid measures in the torus, we get, for any N , , where we label by GP the quantities corresponding to the grid-Poisson transportation problem. Now let us estimate E UP T (N ) from below. In the optimal solution for W 2 2 (σ, ν Y ) the point x is joined, for almost all x, with one of the points in Y. Let us denote with Y J(x) ∈ Y the point in Y to which x is associated. Moreover we can associate almost every point x to the closest grid point, let us denote it by X I(x) ∈ X. Using the functions I(x) and J(x) we can build a coupling between ν Y and ν X defined by (A.9) Pi,j := dx δ i,I(x) δ j,J(x) .
In fact, it is easy to check that ∀i u Pi,u = N −1 and, ∀j, u Pu,j = (hN ) −1 . Therefore, being W 2 2 (ν X , ν Y ) the minimum on all the possible coupling between ν X and ν Y we get Now, taking the average on the locations of the Y points, and noticing that E[x − Y J(x) ] = 0, we get The points As expected, lim h→+∞ c GP T = c UP T .
Appendix B. Robin's mass for the circular sector In this appendix we compute the Robin mass of the Laplace-Beltrami Green function on the circular sector Dp of angle α = 2π p , with p ∈ N, with periodic boundary conditions in the angular direction. We will work on the sector Dp(r) defined in Eq. (5.54) and we will then restrict ourselves to the unit area case. Let us start considering the functions where R θ is the rotation matrix of an angle θ around the origin. The function A (p) (x, y) is such that A (p) (Rαx, y) = A (p) (x, y). In the circular sector Dp(r) we have n + 1 (nκ + r) s (nκ + r + 1) s If k = 2κ + 1, then 2m = 0 mod k iff m = 0 mod k: repeating the usual arguments, the same result is obtained, showing that the Kronecker's mass in the case of Neumann conditions differs from the periodic boundary conditions case by an overall 1 2π ln 2 constant.
Known For a complete discussion of the results sketched here, see for example [44].