\documentclass[11pt,letterpaper]{amsart} \usepackage[margin=1in]{geometry} \usepackage[T1]{fontenc} \usepackage{lmodern,microtype,amsmath,amssymb,amsthm,mathtools} \usepackage{enumitem,booktabs} \usepackage[hidelinks]{hyperref} \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{lemma}[theorem]{Lemma} \newtheorem{proposition}[theorem]{Proposition} \newtheorem{corollary}[theorem]{Corollary} \theoremstyle{definition} \newtheorem{definition}[theorem]{Definition} \theoremstyle{remark} \newtheorem{remark}[theorem]{Remark} \newcommand{\N}{\mathbb N} \newcommand{\Z}{\mathbb Z} \newcommand{\R}{\mathbb R} \newcommand{\eps}{\varepsilon} \newcommand{\E}{\mathbb E} \newcommand{\Prob}{\mathbb P} \newcommand{\ind}{\mathbf 1} \DeclareMathOperator{\supp}{supp} \DeclareMathOperator{\dist}{dist} \DeclareMathOperator{\diver}{div} \DeclareMathOperator{\diam}{diam} \DeclareMathOperator{\Lip}{Lip} \title[Quartic exponential matching]{Quartic exponential matching in the unit square} \author{Will Howe} \email{howew@purdue.edu} \author{Grant Huffman} \email{huffma97@purdue.edu} \author{Jensen Kohlmeyer} \email{jenwin@purdue.edu} \author{Liam Kruer} \email{kruerl@purdue.edu} \thanks{The mathematical draft and exposition were developed and written with OpenAI Codex under the direction of the named contributors. Substantial portions of the argument and text are AI-generated.} \date{November 9, 2026} \begin{document} \begin{abstract} Take two independent samples of $n$ uniform points in the unit square. We prove that, with probability at least one half, the points can be paired so that the average of the exponential of the fourth power of the normalized pair distance is at most two. The distance normalization is $L\sqrt{\log n/n}$, with $L$ an absolute constant. One pairing therefore satisfies both coordinate inequalities in the symmetric $(4,4)$ case of Talagrand's matching problem. The argument first controls a vector field transporting smoothed point densities, and then converts that control to a pairing. The main estimate comes from a Gaussian construction whose successive partitions depend only on noise already revealed; a correction using the classical bottleneck matching preserves the required density balance. \end{abstract} \maketitle \section{Statement and notation} Imagine $n$ red points and $n$ blue points, sampled independently and uniformly in the unit square. A matching pairs every red point with one blue point, using each blue point exactly once. We seek a matching in which long pairs are rare, with a bound that controls the horizontal and vertical errors at the same time. Put $Q=[0,1]^2$, and write $|\cdot|$ for the Euclidean norm. All logarithms are natural. A permutation $\pi$ specifies the pairing of the red point $X_i$ with the blue point $Y_{\pi(i)}$. Talagrand's symmetric $(4,4)$ question \cite{TalagrandMatching} asks for one permutation satisfying the two inequalities \eqref{eq:coordinate}, one for each coordinate. We prove the Euclidean-distance bound \eqref{eq:main}; each coordinate difference is at most that distance, so it gives the requested simultaneous control. \begin{theorem}\label{thm:main} There is an absolute constant $L>0$ such that the following holds for every integer $n\ge2$. Let $X_1,\ldots,X_n,Y_1,\ldots,Y_n$ be independent, uniformly distributed points of $Q$. With probability at least $1/2$ there is a permutation $\pi$ of $\{1,\ldots,n\}$ such that \begin{equation}\label{eq:main} \frac1n\sum_{i=1}^n \exp\left(\frac{|X_i-Y_{\pi(i)}|^4}{L^4 a_n^4}\right)\le2, \qquad a_n=\sqrt{\frac{\log n}{n}}. \end{equation} Consequently the same permutation satisfies, for $j=1,2$, \begin{equation}\label{eq:coordinate} \sum_{i=1}^n \exp\left(\left(\sqrt{\frac n{\log n}} \frac{|X_i^j-Y_{\pi(i)}^j|}{L}\right)^4\right)\le2n. \end{equation} \end{theorem} \subsection*{What the estimate says} Write $d_i=|X_i-Y_{\pi(i)}|$ for the lengths of the matched pairs. Inequality \eqref{eq:main} immediately gives, for every $t\ge0$, \[ \frac1n\#\{i:d_i>La_nt\}\le 2e^{-t^4}. \] It also gives both \[ \frac1n\sum_i d_i\le La_n(\log 2)^{1/4},\qquad \max_i d_i\le La_n(\log(2n))^{1/4}. \] The first follows from Jensen's inequality, and the second by bounding each summand in \eqref{eq:main} by $2n$. Thus the same matching has small average distance, a controlled largest distance, and an explicit bound on the fraction of pairs at every intervening scale. The Leighton--Shor theorem supplies the largest-distance bound at this order, while the AKT theorem supplies the average-distance bound; see \cite{TalagrandMatching} and \cite[Theorem 4.7.1]{TalagrandBook}. The additional task here is to obtain the fourth-power exponential estimate for one matching. We treat only the symmetric $(4,4)$ case. \subsection*{How the proof is organized} The argument has four stages. The purpose of the continuous and Gaussian constructions is to build a pairing of the original finite point sets. \begin{enumerate}[leftmargin=*,itemsep=4pt] \item \emph{Replace the points by densities and reduce matching to a flow estimate.} Spread each point over a neighborhood of radius $h\asymp a_n$. With high probability the two resulting densities $\rho_{X,h}$ and $\rho_{Y,h}$ are bounded above and below. A vector field $J$ with $\operatorname{div}J=\rho_{X,h}-\rho_{Y,h}$ and zero normal component at the boundary describes how to move the excess of the first density to fill the deficit of the second. Lemma~\ref{lem:transport} shows that a fourth-power exponential bound on $J/a_n$ gives the desired bound on the pair lengths. Its proof follows the flow, obtains a fractional matching, and chooses one permutation from a convex decomposition of that matching. \item \emph{Control a Gaussian model of the density imbalance.} The main construction reveals noise one scale at a time. Before the next independent noise is revealed, it chooses the shapes of the next cells using the field already built. This is the meaning of ``predictable'' in Section~4. The area-preserving matrix in \eqref{eq:Az} contracts in the direction of the existing vector and expands in the perpendicular direction. Its identity $|A(z)z|\le1$ controls fluctuations of the squared length even when the vector is large. Section~5 turns this control into an $e^{-ct^4}$ tail up to a cutoff $R=(\log(2/h))^{1/4}$. The exact Gaussian law of the resulting source is checked in Lemma~\ref{lem:white}. \item \emph{Pass from Gaussian noise to the actual samples.} The coupling in Section~7 compares the smoothed empirical imbalance with the Gaussian source. Their difference is the divergence of a field whose supremum is $O(a_n)$, which can be absorbed at the desired scale. The result of the Gaussian construction is therefore useful for the sampled points themselves. \item \emph{Repair large field values without changing the imbalance.} The stopped estimate does not yet give a full exponential moment at all field values. Section~8 uses the classical bottleneck matching to produce a second field with the same divergence and a uniform bound. The difference of the two fields is divergence-free, so it is represented by a stream function. Truncating that function repairs the exceptional large values while retaining the required divergence. Lemma~\ref{lem:splice} gives the full exponential bound, and Lemma~\ref{lem:transport} then returns a matching of the original points. \end{enumerate} The last section assembles these statements and accounts for their exceptional events. We prove the probability threshold $1/2$ in the question; we do not use a stronger probability assertion as an input. Constants denoted by $c,C$ are positive and absolute, and may change from line to line. When a fixed failure probability is prescribed, they may also depend on that probability. An isonormal process on a real Hilbert space $H$ is a centered Gaussian linear map $W:H\to L^2(\Omega)$ with $\E W(f)W(g)=\langle f,g\rangle$. White noise on $Q$ means the isonormal process on $L^2(Q)$. Its mean-zero part is $f\mapsto W(f-\int_Q f)$. Boundary flux identities are interpreted distributionally after extension by zero; thus no boundary measure is present. \section{Compact smoothing and transport} Let $T=(\R/2\Z)^2$ with its ordinary area measure. A scalar field on $Q$ is extended evenly in both coordinates to $T$; a vector field is extended so that its normal component is odd and its tangential component even at each side. Choose a nonnegative smooth function $\eta$ of integral one, even in each coordinate, supported in $[-2,2]^2$, and bounded below by a positive constant on $[-1,1]^2$. For $00$ implies $|x-z|\le2\sqrt2h$. Indeed the circle distance between $x_j$ and $\pm z_j$ is at least $|x_j-z_j|$ for $x_j,z_j\in[0,1]$. For a sample $U=(U_i)_{i\le n}$ put $\rho_{U,h}(x)=n^{-1}\sum_i k_h(x,U_i)$. Choose a sufficiently large fixed $A$, and set \begin{equation}\label{eq:mesh} m=\lfloor (A a_n)^{-1}\rfloor,\qquad h=m^{-1}. \end{equation} We first take $n$ large enough that $m\ge4$. Then $Aa_n\le h\le2Aa_n$. The finitely many other values of $n$ will be handled at the end. \begin{lemma}\label{lem:density} With probability $1-o(1)$ both $\rho_{X,h}$ and $\rho_{Y,h}$ lie between fixed positive constants $c_0,C_0$ on $Q$. On the same event, smoothing either sample at any radius $s\ge h$, with $s<1/4$, gives a density at most an absolute constant. \end{lemma} \begin{proof} Partition $Q$ into the $m^2$ squares of side $h$. Each count has binomial law with mean $nh^2$. The exponential Markov inequality for a binomial variable gives \[ \Prob\{N_C\notin[nh^2/2,2nh^2]\}\le2e^{-cnh^2}. \] A union bound over both samples is at most $4m^2e^{-cnh^2}=o(1)$ when $A$ is sufficiently large. A point's own cell gives the lower density bound using the lower bound on $\eta$, and a bounded number of neighboring cells give the upper bound. At radius $s$, the support meets $O((s/h)^2)$ cells and the kernel is $O(s^{-2})$, proving the last assertion, including reflected cells. \end{proof} \begin{lemma}[Flux-to-matching reduction]\label{lem:transport} Suppose $c_0\le\rho_{X,h},\rho_{Y,h}\le C_0$ and $J$ is a smooth field on $Q$ with zero normal component and $\diver J=\rho_{X,h}-\rho_{Y,h}$. If \begin{equation}\label{eq:quarticflux} \int_Q \exp\bigl((|J|/b)^4\bigr)\le D, \end{equation} then a permutation satisfies \eqref{eq:main} with $L a_n$ replaced by $C\max(h,b)$, where $C$ depends only on $c_0,C_0,D$. \end{lemma} \begin{proof} Set $\rho_t=(1-t)\rho_{X,h}+t\rho_{Y,h}$ and $v_t=J/\rho_t$. The smooth vector field $v_t$ is tangent to the boundary and its flow $T_t$ stays in $Q$. The identity $\partial_t\rho_t+\diver(\rho_t v_t)=0$, together with the Jacobian formula for the flow, gives $(T_t)_\#(\rho_{X,h}\,dx)=\rho_t\,dx$. For any convex nondecreasing $\Psi$ with $\Psi(0)=0$, Jensen's inequality along trajectories gives \begin{align*} \int_Q\Psi\bigl(c_0|T_1(z)-z|/b\bigr)\rho_{X,h}(z)\,dz &\le\int_0^1\int_Q\Psi\bigl(c_0|v_t(x)|/b\bigr)\rho_t(x)\,dx\,dt\\ &\le C_0\int_Q\Psi(|J(x)|/b)\,dx. \end{align*} Choose a uniform label $i$, move $X_i$ to $z$ with density $k_h(z,X_i)$, then to $w=T_1(z)$, and finally choose label $k$ with probability $k_h(w,Y_k)/(n\rho_{Y,h}(w))$. The last label is uniform because the law of $w$ is $\rho_{Y,h}\,dw$. Both smoothing moves have length at most $2\sqrt2h$. Applying convexity to these three displacements with $\Psi(t)=e^{t^4}-1$ bounds the mean label-pair cost at scale $C\max(h,b)$ by a constant depending only on the displayed parameters. The matrix of label-pair probabilities, multiplied by $n$, is doubly stochastic. For completeness, it is a convex combination of permutation matrices: the support of a doubly stochastic matrix satisfies Hall's condition, since rows in a set of size $r$ have total mass $r$ and each neighboring column has mass at most one. Hall's matching criterion follows by the usual augmenting-path argument: if a maximum matching leaves a row unmatched, the rows reachable from it by alternating paths have fewer neighbors than rows. Thus a support permutation exists. Subtract the smallest entry on this permutation, normalize the remaining matrix, and repeat; at each nonterminal step an additional positive entry vanishes. This proves the decomposition in finitely many steps. Some component has no greater total cost than the fractional plan. Finally, if $n^{-1}\sum_i e^{t_i}\le D_1$, then $n^{-1}\sum_i e^{t_i/q}\le D_1^{1/q}$ by concavity, for $q\ge1$. Taking $q$ large enough makes the bound two and completes the proof. \end{proof} \section{A local Gaussian flux on affine Whitney cells} Let $Q_0=(-1/2,1/2)^2$, and put $d(u)=\dist(u,\partial Q_0)$. Fix a sufficiently small absolute $\eps>0$. If $D$ is a real matrix with determinant one and $q=\|D^{-1}\|$, take the maximal dyadic squares $S$ in $DQ_0$ satisfying, on their closures, \begin{equation}\label{eq:whitney} \ell(S)\le\eps q^{-1}\inf_{w\in S}d(D^{-1}w). \end{equation} Here $q\ge1$. Half-open representatives partition $DQ_0$ up to a null set, since every interior point off dyadic boundaries has arbitrarily small eligible squares. The function $d(D^{-1}w)$, extended by zero, is $q$-Lipschitz. If a square's parent is ineligible or meets the complement, this Lipschitz bound and maximality imply \begin{equation}\label{eq:whitneysize} c_\eps d(u)/q\le\ell(S)\le\eps d(u)/q, \qquad u\in C=D^{-1}S. \end{equation} The children $C$ are parallelograms. One of $D^{-1}e_1,D^{-1}e_2$ has length at least $q/\sqrt2$; choose its unit direction $e$. All children have an edge parallel to $e$, of length $l_C\asymp_\eps d(u)$ on $C$. Their diameters, and the norms of their affine linear parts as images of $Q_0$, are also comparable to $d(u)$. By reducing $\eps$, all these linear parts have norm at most a fixed $\rho<1$. No lower bound on the short edge is used. Let $P$ be orthogonal projection in $L^2(Q_0)$ onto child-constant functions, and $P_0$ projection onto constants. Choose a nonnegative smooth bump $b$ of integral one supported in a fixed central ball of $Q_0$. The Bogovski\u\i\ kernel, in dimension two, is \begin{equation}\label{eq:bogovskii} K(u,v)=(u-v)\int_1^\infty b(v+t(u-v))t\,dt. \end{equation} This is the classical integral construction; see \cite{CostabelMcIntosh}. The estimates needed here are proved below. \begin{lemma}\label{lem:local} There is an absolute constant $C$, independent of $D$, such that \begin{align} \|(P-P_0)K(u,\cdot)\|_2&\le C,\label{eq:rowbound}\\ \|(P-P_0)(K(u,\cdot)-K(u',\cdot))\|_2 &\le C\min\{1,(|u-u'|/d_C)^{1/3}\}\label{eq:increment} \end{align} for $u,u'$ in the closure of the same child $C$, where $d_C=d(c_C)$ and $c_C$ is its center. Projections of the integrable kernel mean its cell averages. For isonormal $W$ define the two components of $H(u)$ by applying $W$ to the two rows of $(P-P_0)K(u,\cdot)$. Then \[ \diver H=W\bigl((P-P_0)\,\cdot\bigr) \] with zero normal boundary condition. On each child a continuous modification satisfies, with $\gamma=1/6$, \begin{equation}\label{eq:holdertail} \Prob\left\{sup_C|H|+d_C^\gamma[H]_{C^\gamma(C)}>t\right\} \le C e^{-ct^2},\qquad t\ge1. \end{equation} All statements remain valid conditionally if $D$ is determined before the independent process $W$ is sampled. \end{lemma} \begin{proof} Directly from \eqref{eq:bogovskii}, \begin{equation}\label{eq:kernelbounds} |K(u,v)|\le C|u-v|^{-1},\quad |\partial_uK(u,v)|\le C|u-v|^{-2},\quad K(u,v)=0\quad\hbox{if }|u-v|>C d(u). \end{equation} For the first two bounds the integration is restricted to $t\le C/|u-v|$; differentiating $b$ gives one additional factor $t$. For the third, if $z=v+t(u-v)$ is in the central ball, then $u=(1-t^{-1})v+t^{-1}z$. Signed distances to the four sides are affine on this segment, so $d(u)\ge c/t$ while $|u-v|=|z-v|/t\le C/t$. We claim, for $0r$. The integral of the square of this logarithm is $O(r)$. Sum over the children to prove \eqref{eq:near}. This also controls the signed vector components. Outside the ball of radius $c d(u)$, \eqref{eq:kernelbounds} gives a bounded $L^2$ norm: the possible radii range between two constant multiples of $d(u)$. Projection contracts $L^2$. The parent average is the constant projection of the child average. This proves \eqref{eq:rowbound}, and the same argument bounds the projection of the absolute kernel. For \eqref{eq:increment} let $\delta=|u-u'|\ll d_C$ and take $r=d_C(\delta/d_C)^{2/3}$. Near the two points, \eqref{eq:near} gives an $L^2$ bound $C(r/d_C)^{1/2}$. Outside these balls the derivative bound in \eqref{eq:kernelbounds}, integrated in $L^2$ along the segment, gives $C\delta/r$. These bounds agree at $C(\delta/d_C)^{1/3}$. The remaining values of $\delta$ are covered by \eqref{eq:rowbound}. If $g\in L^2(Q_0)$ has mean zero, $Bg(u)=\int K(u,v)g(v)\,dv$ obeys $\diver Bg=g$, including after extension by zero. Indeed pairing with $\nabla\varphi$ and using $z=v+t(u-v)$ gives \[ \int Bg\cdot\nabla\varphi =\int g(v)\int b(z)\int_0^1 (z-v)\cdot\nabla\varphi(v+s(z-v))\,ds\,dz\,dv =-\int g\varphi. \] The term $\varphi(z)$ vanishes by the mean-zero condition. Contributing segments stay in $Q_0$, and the size bound makes the pairing absolutely integrable. Apply this identity to $(P-P_0)g$. The row bounds and the corresponding absolute-kernel bound justify interchanging $W$ and the $u$-integral, or equivalently first using finite sums of cell averages and then passing in $L^2$. This proves the asserted divergence identity. Finally cover a child by nets of mesh $2^{-j}d_C$ with $O(2^{2j})$ points, join points to a nearest point in the preceding net, and include nearby same-level pairs. By \eqref{eq:increment}, the probability that any level-$j$ link exceeds $t2^{-j\gamma}$ is bounded by \[ C2^{4j}\exp\{-ct^2 2^{2j(1/3-\gamma)}\}. \] The sum over $j$ is $Ce^{-ct^2}$ for $t\ge1$. Add a base value using \eqref{eq:rowbound}; summing links gives the modification and the H\"older seminorm bound. Nets can be chosen from a fixed countable dense set in square coordinates, so they are measurable in $D$. Countably many children can be treated simultaneously. Values from different children need not be identified on their boundaries. \end{proof} \section{Predictable partitions and exact white noise} Fix $0C K_*\right\}\le Ce^{-cK_*}. \end{equation} The exponent in this estimate can be made any prescribed fixed multiple of $K_*$ by enlarging $C$. Stop the pointwise calculation at the first step with $|Z_i|>R$, retaining that step. Subsequent terms in the calculation are zero. At an active step let $\widehat z_i$ be the capped state used for its old metric. For $i\ge2$ it is the cap of $Z_{i-2}$ at the center of the preceding parent; for $i=1$ it is zero. Before the crossing put $\delta_i=|Z_{i-1}-\widehat z_i|$, and put it zero at inactive or subsequently stopped steps. Both vectors in this difference lie in the radius-$R$ ball, so $\delta_i\le2R$. \begin{lemma}\label{lem:mismatch} On the good event $\sum_{i=1}^N G_i^2\le C K_*$ from \eqref{eq:Gsum}, \begin{equation}\label{eq:deltasum} \sum_i\delta_i^2\le C R^2. \end{equation} \end{lemma} \begin{proof} Projection onto the cap ball cannot increase the distance from $Z_{i-1}\in\overline B(0,R)$. The last increment at $x$ has size at most $C R G_{i-1}/\sqrt{K_*}$. For an earlier increment $j\le i-2$, both $x$ and the center used in step $i-1$ belong to the $j$th child. In the $j$th parent's coordinates their distance is at most \[ C\|B_{P_{j-1}}^{-1}B_{P_j}\|\rho^{i-2-j}. \] The first factor is comparable to the $d_C$ of that child by \eqref{eq:whitneysize}. Apply \eqref{eq:holdertail} to the oscillation, and use $\|A_{\rm old}\|\le C R$. It follows that \[ \delta_i\le\frac{C R}{\sqrt{K_*}} \left(G_{i-1}+\sum_{j\le i-2}\rho^{\gamma(i-2-j)}G_j\right). \] The first two steps obey the same formula with empty sums. The geometric sequence has bounded $\ell^1$ norm. Its convolution therefore has $\ell^2$ norm at most a constant times that of $(G_j)$, which, together with \eqref{eq:Gsum}, proves \eqref{eq:deltasum}. \end{proof} \begin{proposition}\label{prop:tail} For $1\le r\le R$ the stopped process satisfies \begin{equation}\label{eq:quartictail} \Prob\{\max_{i\le\text{crossing or }N}|Z_i|>r\} \le Ce^{-cr^4}+Ce^{-cK_*}. \end{equation} There is an absolute $\theta>0$ such that, if \begin{equation}\label{eq:Phi} \Phi_R(t)= \begin{cases} e^{\theta t^4},&0\le t\le R,\\ e^{\theta R^4}\bigl(1+4\theta R^3(t-R)\bigr),&t>R, \end{cases} \end{equation} then the continued, unstopped field satisfies \begin{equation}\label{eq:moment} \E\int_Q\Phi_R(|Z_N(x)|)\,dx\le C. \end{equation} \end{proposition} \begin{proof} At an active step the increment has the form \[ \Delta_i=K_*^{-1/2}A(\widehat z_i)U_i, \] where $U_i$ is conditionally centered Gaussian with covariance at most $CI$. At inactive steps it is zero. All terms below are zero after the crossing. Let $S_i=1+|Z_i|^2$, frozen after crossing, and put $M_i=2Z_{i-1}\cdot\Delta_i$. From \eqref{eq:Aproperties}, \[ \E(M_i^2\mid\mathcal F_{i-1}) \le\frac C{K_*}(1+R^2\delta_i^2),\qquad |\Delta_i|^2\le\frac C{K_*}(S_{i-1}+\delta_i^2)|U_i|^2. \] By Lemma~\ref{lem:mismatch}, the sum of the conditional variance bounds is at most $C(1+R^4/K_*)=C$, except on the event in \eqref{eq:Gsum}. Conditional Gaussian exponential martingales imply \begin{equation}\label{eq:marttail} \Prob\left\{\max_\ell\left|\sum_{i\le\ell}M_i\right|>u\right\} \le Ce^{-cu^2}+Ce^{-cK_*}. \end{equation} Here the exponential martingale is stopped when its predictable variance exceeds the fixed bound. Thus this argument restricts an event; it does not condition on the future event \eqref{eq:Gsum}. The conditional Gaussian-square formula also gives \[ \Prob\{\sum_i|U_i|^2>C K_*\}\le Ce^{-cK_*}. \] Furthermore, for $u\ge1$, \begin{equation}\label{eq:weighted} \Prob\left\{\sum_i\delta_i^2|U_i|^2/K_*>u, \ \sum_i\delta_i^2\le C R^2\right\} \le C\exp(-cK_*u/R^2). \end{equation} To verify this, $\delta_i$ is predictable and at most $2R$. For a small fixed $c$, the conditional logarithm of the moment generating function of $(c/R^2)\delta_i^2|U_i|^2$ is at most $C\delta_i^2/R^2$. Iteration gives a compensated exponential bound; restrict its compensator to $\sum_i\delta_i^2\le C R^2$ and apply Markov's inequality. This proves \eqref{eq:weighted} without conditioning on that future restriction. The exact identity $S_i-S_{i-1}=M_i+|\Delta_i|^2$ and discrete Gronwall give, on $\sum_i|U_i|^2\le CK_*$, \[ \max_i S_i\le C\left(1+max_\ell\left|\sum_{i\le\ell}M_i\right| +\sum_i\delta_i^2|U_i|^2/K_*\right). \] Indeed the remaining nonnegative coefficient of $S_{i-1}$ sums to a constant and the resulting Gronwall factor is its exponential. Use \eqref{eq:marttail} at $u\asymp r^2$ and \eqref{eq:weighted} at the same threshold. Since $K_*/R^2=R^2\ge r^2$, this proves \eqref{eq:quartictail}, with enlarged constants for bounded $r$. In particular the crossing probability is at most $Ce^{-cK_*}$. For the continued field, conditional orthogonality of increments and \eqref{eq:rowbound} give $\E|Z_N|^2\le C R^2$. The function \eqref{eq:Phi} is positive, nondecreasing and convex. On no crossing, integrate \eqref{eq:quartictail} up to $R$, with $\theta$ smaller than the tail constant, to bound its expectation. On crossing use the global bound \[ \Phi_R(t)\le e^{\theta R^4}(1+4\theta R^3t). \] Cauchy--Schwarz bounds its contribution by $Ce^{(\theta-c/2)K_*}(1+R^4)$. Reduce $\theta$ again. This is bounded, and averaging over the independent uniform $x$ proves \eqref{eq:moment}. \end{proof} \section{Smoothing and completion of the Gaussian source} Let $G$ be the mean-zero Green function of $\Delta$ on $T$, so that $\Delta G=\delta_0-1/4$, and put $L_h=\nabla G*\eta_h$. For $x,z\in Q$ set \begin{equation}\label{eq:F} F_h(x,z)=\sum_{g\in\{\pm1\}^2}L_h(x-gz). \end{equation} It has zero mean in $z$, zero $x$-normal component, and $\diver_xF_h(x,z)=k_h(x,z)-1$. \begin{lemma}\label{lem:green} For all $x\in Q$, \begin{align} \|F_h(x,\cdot)\|_\infty&\le C/h,& \|F_h(x,\cdot)\|_2^2&\le C\log(2/h),\label{eq:greenbasic}\\ \|\nabla_zF_h(x,\cdot)\|_2^2&\le C/h^2,& \|D_x^jF_h\|_\infty&\le C_jh^{-1-j}\quad(0\le j\le3). \label{eq:greenderiv} \end{align} If $P$ is projection onto constants of any convex partition with cell diameters at most $h$, then \begin{equation}\label{eq:resvar} \|(I-P)F_h(x,\cdot)\|_2^2\le C. \end{equation} \end{lemma} \begin{proof} Cut off $(2\pi)^{-1}\log|w|$ near zero in $T$. Its Laplacian is $\delta_0$ plus a smooth function; subtract the required mean, and solve the smooth remainder by Fourier series. This constructs $G$ and gives $|\nabla G(w)|\le C(1+\dist_T(w,0)^{-1})$ and, away from zero, $|D^j\nabla G(w)|\le C_j(1+\dist_T(w,0)^{-1-j})$. For the smoothed kernel split the convolution into distances $O(h)$ and its complement. On the former move derivatives to the bump; integration of $|w|^{-1}$ on a radius-$Ch$ disk is $O(h)$. On the latter use the unsmoothed derivatives. Thus \[ |L_h(w)|\le\frac C{h+\dist_T(w,0)},\qquad |\nabla L_h(w)|\le\frac C{(h+\dist_T(w,0))^2}, \] and the indicated derivative supremum bounds follow in the same way. Integration in dimension two proves \eqref{eq:greenbasic}--\eqref{eq:greenderiv}, with the fixed reflection sum included. For $z,z'$ in one cell, $|z-z'|\le h$. Along their segment the denominator $h+\dist_T(x-gz,0)$ changes by at most a constant factor. Consequently \[ |L_h(x-gz)-L_h(x-gz')| \le\frac{Ch}{(h+\dist_T(x-gz,0))^2}. \] Average over $z'$ in that cell, square, and integrate in $z$. Summing over cells and reflections gives a bound by a constant times $\int_T h^2(h+\dist_T(w,0))^{-4}\,dw$, which is bounded. This is \eqref{eq:resvar}. \end{proof} \begin{proposition}\label{prop:gaussian} Given a fixed $\eps_0>0$, for all sufficiently small $h$ there is a joint law of mean-zero white noise $\xi$ and a smooth zero-normal field $V_h$ such that \begin{equation}\label{eq:gaussianflux} \diver V_h=\xi(k_h(x,\cdot)-1),\qquad \int_Q\Phi_R\left(\frac{|V_h|}{C\sqrt{K_*}}\right)\le C \end{equation} with probability at least $1-\eps_0$. Here $K_*,R$ are as in \eqref{eq:KR}, and $C$ may depend on $\eps_0$. \end{proposition} \begin{proof} Reflect $Z_N$ as a vector and convolve by $\eta_h$. The distributional zero-normal identity makes reflection and convolution commute with divergence. Jensen and \eqref{eq:moment} retain the bounded expected $\Phi_R$ integral. Its source, before normalization by $\sqrt{K_*}$, is the smoothing of $\xi-\xi_{\rm res}$. Add the residual field \[ E_h(x)=W_{N+1}((I-P_N)F_h(x,\cdot)). \] All final cells have diameter at most $h$, so \eqref{eq:resvar} bounds its conditional pointwise variance by $C$. By the derivative estimates, its conditional expected squared $H^3(Q)$ norm is $O(h^{-8})$, uniformly in the random partition. The Sobolev embedding $H^3(Q)\hookrightarrow C^1(Q)$ and Markov's inequality imply, for example, a $C^1$ norm at most $h^{-6}$ except with probability $O(h^4)$. On an $h^8$-mesh, the Gaussian union bound then gives \[ \|E_h\|_\infty\le C\sqrt{K_*} \] with probability $1-o(1)$; interpolation contributes at most $h^2$. Higher Sobolev orders give smooth versions. The divergence is exactly the missing residual source and its normal component is zero. Adding a bounded normalized error preserves the estimate after a fixed rescaling: if $\|E\|_\infty\le B$, convexity gives \[ \Phi_R(|J+E|/2)\le\tfrac12\Phi_R(|J|)+\tfrac12\Phi_R(B). \] For fixed $B$, the last term is uniformly bounded. Apply this observation and Markov's inequality to the expected integral, with failure budget $\eps_0$, to prove the proposition. \end{proof} \section{The empirical--Gaussian coupling} We record all details of the transfer because the source must remain exact, and the two empirical samples must retain their independent product law. The only non-elementary probabilistic input is the following local coupling: for every $\lambda\ge1$ there is a coupling of $N_\lambda\sim\operatorname{Pois}(\lambda)$ and $Z\sim N(0,1)$ such that \begin{equation}\label{eq:poissoncoupling} \Prob\{|N_\lambda-\lambda-\sqrt\lambda Z|>t\}\le Ce^{-ct} \qquad(t\ge1). \end{equation} This is \cite[Lemma 4.5]{LachiezeMuirhead}. The error has mean zero, since both marginals have their stated means. \begin{lemma}\label{lem:coupling} For any fixed $\eps_0>0$, $h$ chosen by \eqref{eq:mesh}, and sufficiently large $n$, one can jointly realize independent uniform samples $X,Y$, mean-zero white noise $\xi$, and a smooth zero-normal error field $E_n$ so that \begin{equation}\label{eq:empiricalcoupling} \rho_{X,h}-\rho_{Y,h} =\sqrt{2/n}\,\xi(k_h(x,\cdot)-1)+\diver E_n, \qquad \|E_n\|_\infty\le C a_n \end{equation} with probability at least $1-\eps_0$. \end{lemma} \begin{proof} First treat one uniform sample $U$. Let $P_h$ be projection onto the constants of the regular $h$-grid, and define $b_C(x)=h^{-2}\int_C F_h(x,z)\,dz$. Lemma~\ref{lem:green} gives \begin{gather*} \sum_C b_C(x)=0,\quad \max_C|b_C(x)|\le C/h,\quad \sum_C|b_C(x)|^2\le Ch^{-2}\log(2/h),\\ \Lip_x b_C\le C/h^2,\qquad \|(I-P_h)F_h(x,\cdot)\|_2^2\le C. \end{gather*} We compare successive fluxes; applying $\diver$ to each always gives the corresponding exact difference of smoothed mean-zero sources. Replace $F_h(x,U_i)$ by $b_{C(U_i)}(x)$. The differences are independent, centered, bounded by $C/h$, and have variance at most $C$. Scalar Bernstein inequalities in the two coordinates give at threshold $Da_n$ for the normalized sum the exponent \[ \frac{cD^2\log n}{1+D/A}. \] A union bound over a mesh of size $n^{-8}$ succeeds for sufficiently large fixed $D$. The normalized fields have deterministic Lipschitz bound $C/h^2$, so interpolation is negligible. Now take $P\sim\operatorname{Pois}(n)$ independent of an infinite iid sequence extending $U$, and replace the first $n$ points by the first $P$ points, still dividing by $n$. On $|P-n|\le B\sqrt n$ the difference is a sum of at most $B\sqrt n$ independent centered $b_C$ variables, conditionally on $P$. The probability of the complementary size event is at most $B^{-2}$. Bernstein's exponent at $Da_n$ is at least \[ c\min\{D^2\sqrt n/B,\ Dn a_n h\} \ge c\min\{D^2\sqrt n/B,\ DA\log n\}. \] Choose $B$ for the desired failure budget and then $D$ for the same mesh argument. Subtracting the constant mode for every point means that a change of sample size creates no unaccounted constant source. The Poisson cell counts $N_C$ are independent with common mean $\lambda=nh^2$. Couple them independently to standard Gaussians $Z_C$ using \eqref{eq:poissoncoupling}, conditionally on the already sampled counts. The errors $E_C=N_C-\lambda-\sqrt\lambda Z_C$ are independent centered variables with uniform exponential tails. The error flux is $n^{-1}\sum_C E_Cb_C$. Weighted Bernstein gives exponent at least \[ c\min\left\{\frac{D^2 n^2a_n^2h^2}{\log(2/h)},\ Dn a_nh\right\} \ge c'\min\{D^2A^2,DA\}\log n. \] Also $\max_C|E_C|\le C\log n$ with probability $1-o(1)$. On this event the flux Lipschitz bound is $O(h^{-2})$, because $h^{-2}\log n/n=O(1)$. A mesh argument again bounds this error uniformly by $Ca_n$. Realize $Z_C$ as $W_Q(\ind_C/h)$, and fill the orthogonal complement by independent white noise. Then \[ \frac{\sqrt\lambda}{n}\sum_C Z_C b_C(x) =\frac1{\sqrt n}W_Q(P_hF_h(x,\cdot)). \] The omitted Gaussian field $n^{-1/2}W_Q((I-P_h)F_h(x,\cdot))$ has variance at most $C/n$. Its expected squared $H^3$ norm is $O(n^{-1}h^{-8})$. The same Sobolev, net and interpolation argument as above bounds its supremum by $Ca_n$ with probability $1-o(1)$. For example its Lipschitz norm is at most $n^5$ with probability $1-o(1)$, and the $n^{-8}$ mesh suffices. Combining these four differences proves the one-sample assertion. Use independent copies of the entire construction for $X$ and $Y$. They retain their required product law, and $\xi=(W_X-W_Y)/\sqrt2$ is standard white noise; constants are irrelevant when paired with $k_h(x,\cdot)-1$. Summing the smooth zero-normal errors proves \eqref{eq:empiricalcoupling}. \end{proof} \section{A bounded comparison field and the quartic splice} The classical Leighton--Shor theorem supplies, with probability $1-o(1)$, a matching from either sample to the same evenly spread $n$-point set whose largest edge has length at most $C(\log n)^{3/4}/\sqrt n$. Compose the two matchings. Thus, writing $K_*=\log(2/h)$ and $R=K_*^{1/4}$ as before, one obtains a matching of $X$ to $Y$ of maximum length at most $C a_nR$. The stronger classical failure estimate is $C\exp(-c(\log n)^{3/2})$ \cite[Theorem 4.7.1]{TalagrandBook}. \begin{lemma}\label{lem:bounded} On the count event of Lemma~\ref{lem:density} and this bottleneck event, there is a smooth vector-reflection-invariant field $B$ such that \begin{equation}\label{eq:B} \diver B=\rho_{X,h}-\rho_{Y,h},\qquad \|B\|_\infty\le C a_nR. \end{equation} \end{lemma} \begin{proof} Let $s_*=C a_nR$ bound the lengths in the classical matching. First convolve its straight paths at radius $s_*$ on the reflected torus. The field is \[ B_{s_*}(x)=\frac1n\sum_i\sum_g\int_0^1 g(Y_{\pi(i)}-X_i) \eta_{s_*}\bigl(x-g((1-t)X_i+tY_{\pi(i)})\bigr)\,dt. \] Differentiation in $t$ gives divergence $\rho_{X,s_*}-\rho_{Y,s_*}$. A contributing path starts within $Cs_*$ of $x$. The count bound implies that at most $Cns_*^2$ starting points contribute. Each velocity is at most $s_*$ and each kernel at most $Cs_*^{-2}$, so $\|B_{s_*}\|_\infty\le Cs_*$. For each marginal repair the smoothing radius from $s_*$ to $h$ by the dilation paths $gX_i+s(t)w$, where $w$ has density $\eta$ and $s(t)$ runs linearly between the two radii. Their time densities are the corresponding reflected smoothings, uniformly bounded by Lemma~\ref{lem:density}, and their speeds are $O(s_*)$. Their integrated fluxes are bounded by $Cs_*$ and have exactly the divergence of the endpoint density difference. Choose the orientations and subtract the $Y$ repair to obtain \eqref{eq:B}. All fields have the required vector reflection symmetry and hence zero normal component on $Q$. \end{proof} \begin{lemma}[Divergence-preserving splice]\label{lem:splice} Let $R\ge1$. Suppose smooth vector-reflection-invariant fields $J,B$ on $T$ have the same divergence and \[ \int_Q\Phi_R(|J|)\le D,\qquad \|B\|_\infty\le C_B R. \] Then a smooth vector-reflection-invariant field $J_*$ with that divergence satisfies \begin{equation}\label{eq:spliced} \int_Q e^{c_*|J_*|^4}\le C_*, \end{equation} where $c_*>0,C_*<\infty$ depend only on $D,C_B,\theta$. \end{lemma} \begin{proof} Let $M$ be the centered Hardy--Littlewood maximal operator on $T$, and set $E=\{M|J|\le\lambda\}$, where $\lambda=C_2R\ge R$ and $C_2$ is fixed. Since $J-B$ is divergence-free and has zero mean by vector reflection, its nonzero Fourier coefficients give a periodic smooth stream function $u$ with $J-B=\nabla^\perp u$ and mean zero. We recall the pointwise Sobolev estimate \[ |u(x)-u(y)|\le C\dist_T(x,y) \bigl(M|\nabla u|(x)+M|\nabla u|(y)\bigr). \] To prove it, telescope the averages over balls of radii $2^{-j}\dist_T(x,y)$ about each point, using the $L^1$ Poincar\'e inequality, and compare the two largest balls in a ball of three times the radius. Poincar\'e on a ball follows by integration on segments, and the same argument in finitely many charts treats large torus distances. On $E$, $M|\nabla u|\le\lambda+C_BR$, so $u|_E$ is $CR$-Lipschitz. The McShane formula $v(x)=\inf_{y\in E}(u(y)+CR\dist_T(x,y))$ extends it with the same Lipschitz constant. If $E$ is empty take $v=0$. The mean-zero stream function has the determinant reflection symmetry $u(gx)=\det(g)u(x)$. Averaging $v$ with these signs retains its Lipschitz constant and agreement with $u$ on the invariant set $E$. Consequently $K=B+\nabla^\perp v$ has the desired divergence, $\|K\|_\infty\le CR$, and $K=J$ almost everywhere on $E$. For the last statement use locality of weak derivatives on the set where two Lipschitz functions agree; it follows, for instance, by applying the one-dimensional absolute-continuity property on almost every parallel line and then to both coordinate derivatives. Convexity implies $\Phi_R(M|J|)\le M(\Phi_R(|J|))$. The weak $(1,1)$ maximal inequality gives \begin{equation}\label{eq:badset} |E^c|\le C D/\Phi_R(\lambda)\le C D e^{-\theta R^4}. \end{equation} Here weak type follows from the disjoint-ball covering argument: choose balls whose averages exceed the threshold, select disjoint ones by decreasing radius, and cover the original set by their fixed dilates. For $x\in E$ with $|J|\le R$, a sufficiently small $c_*$ bounds $e^{c_*|J|^4}$ by $\Phi_R(|J|)$. The part with $R<|J|\le\lambda$ has measure at most $D e^{-\theta R^4}$, and its exponential is at most $e^{c_*C_2^4R^4}$. On $E^c$ use \eqref{eq:badset} and $|K|\le CR$. Choosing $c_*$ below both relevant constant multiples of $\theta$ bounds $\int_Qe^{c_*|K|^4}$. Finally mollify $v$ with an arbitrarily narrow even torus kernel. The resulting field equals $K*\eta_\delta+B-B*\eta_\delta$. Its divergence is still $\diver B$, and the last term can have supremum at most one. Jensen and a further fixed reduction of $c_*$ preserve the bound. This proves \eqref{eq:spliced} with a smooth field. \end{proof} \section{Proof of the theorem and the indexing convention} \begin{proof}[Proof of Theorem~\ref{thm:main}] Use the mesh \eqref{eq:mesh}. Choose a small fixed failure budget $\eps_0$, for example $1/20$. The white-noise marginal in Proposition~\ref{prop:gaussian} agrees exactly with the white-noise marginal in Lemma~\ref{lem:coupling}. Join their two joint laws by conditioning on this common marginal. This is legitimate on countable isonormal coordinates; the fields are also determined by countably many coordinates. Equivalently, regular conditional probabilities on these standard Borel spaces give the usual gluing of probability laws. Since $\sqrt{K_*/n}\asymp a_n$, the field $J=\sqrt{2/n}\,V_h+E_n$ has exact divergence $\rho_{X,h}-\rho_{Y,h}$ and obeys \[ \int_Q\Phi_R\bigl(|J|/(C a_n)\bigr)\le C \] except on an event of probability at most $2\eps_0+o(1)$. The bounded normalized error is absorbed by convexity as in the proof of Proposition~\ref{prop:gaussian}. The count and classical bottleneck events have probability $1-o(1)$. On their intersection, Lemma~\ref{lem:bounded} supplies a same-divergence comparison field. Apply Lemma~\ref{lem:splice} in units of $C a_n$ to obtain a smooth zero-normal field with \[ \int_Q\exp\left((|J_*|/(C'a_n))^4\right)\le C''. \] The density bounds and Lemma~\ref{lem:transport} now give a single permutation satisfying \eqref{eq:main}. The intersection probability is larger than $1/2$ for all sufficiently large $n$. For the finitely many remaining $n\ge2$, increase $L$ so that the deterministic bound $|X_i-Y_k|\le\sqrt2$ makes every summand at most two. Finally $|X_i^j-Y_{\pi(i)}^j|\le|X_i-Y_{\pi(i)}|$ proves \eqref{eq:coordinate} for both coordinates of the same permutation. \end{proof} \begin{corollary}\label{cor:indexing} There is an absolute $L'>0$ such that, for each $N\ge2$ and independent uniform samples $(X_i)_{i=0}^N,(Y_i)_{i=0}^N$, with probability at least $1/2$ a permutation of $\{0,\ldots,N\}$ satisfies \[ \sum_{i=0}^N\exp\left(\left(\sqrt{\frac N{\log N}} \frac{|X_i^j-Y_{\pi(i)}^j|}{L'}\right)^4\right)\le2N, \qquad j=1,2. \] \end{corollary} \begin{proof} Apply Theorem~\ref{thm:main} with $n=N+1$. Multiply the displacement scale by $4^{1/4}$. Concavity of $t\mapsto t^{1/4}$ bounds the new exponential average by $2^{1/4}<4/3$. Its sum is therefore at most $(4/3)(N+1)\le2N$. The scales $\sqrt{\log(N+1)/(N+1)}$ and $\sqrt{\log N/N}$ are comparable by absolute constants for $N\ge2$. Increasing $L'$ absorbs this comparison. This also covers the literal indexing of \cite{TalagrandMatching}. \end{proof} \begin{thebibliography}{9} \bibitem{CostabelMcIntosh} M. Costabel and A. McIntosh, \emph{On Bogovski\u\i\ and regularized Poincar\'e integral operators for de Rham complexes on Lipschitz domains}, Math. Z. \textbf{265} (2010), no.~2, 297--320. \href{https://doi.org/10.1007/s00209-009-0517-8}{doi:10.1007/s00209-009-0517-8}. \bibitem{LachiezeMuirhead} R. Lachi\`eze-Rey and S. Muirhead, \emph{Asymptotics for the critical level and a strong invariance principle for high intensity shot noise fields}, arXiv:2111.09180v2 (2022), Lemma~4.5. \url{https://arxiv.org/abs/2111.09180}. \bibitem{TalagrandBook} M. Talagrand, \emph{Upper and Lower Bounds for Stochastic Processes: Decomposition Theorems}, second ed., Ergebnisse der Mathematik und ihrer Grenzgebiete, vol.~60, Springer, 2021. \href{https://doi.org/10.1007/978-3-030-82595-9}{doi:10.1007/978-3-030-82595-9}. \bibitem{TalagrandMatching} M. Talagrand, \emph{Matchings}, problem note, \href{https://michel.talagrand.net/prizes/matchings.pdf}{available online}. \end{thebibliography} \end{document}