The goal here is to show that Michael's "usually Poisson" reasoning can be made rigorous. For $d>0$, let $c_d$ denote the number of other lattice points in $\mathbb{Z}^2$ within distance $d$ of the point $(n,n)$. My claim will be that for fixed $d$ the number of matching pairs of distance at most $d$ is asymptotically Poisson with mean $c_d/2$, and in particular that the probability there is at least one pair within distance $d$ converges to $1-e^{-c_d/2}$.
As pointed out by Michael in his comment, this doesn't directly extend to a guarantee of bounded expectation without some further argument showing there isn't a part of the tail that has probability tending to $0$ but expectation tending to infinity.
Call an unordered pair $\{ (x_1, y_1), (x_2, y_2) \}$ of points in our grid close if the points are within distance $d$ from each other. The number of close pairs is equal to $2n^2 c_d + O(n)$, where here (and throughout) the constant in the $O$ notation may depend on $d$ and $m$. (Equivalently, there's $4n^2 c_d+O(n)$ ordered pairs of close points. The idea is that each of the $4n^2-O(n)$ points $(x_1, y_1)$ not within distance $d$ of the boundary of the grid is the first point in $c_d$ pairs).
For each close pair $p$ of points, let $X_p$ be equal to $1$ if the pair matches, and $0$ otherwise. Then we have $$E(X_p) = \frac{1}{4n^2-1}$$ and the number of matching pairs is $$X= \sum_{p} X_p,$$ Let $\{Z_p\}$ be a family of independent random variables each equal to $1$ with probability $\frac{1}{4n^2-1}$, and let $Z=\sum Z_p$. Then $Z$ converges to the Poisson distribution as $n$ goes to infinity. To show that $X$ also converges to Poisson, it is enough to show that the moments of $X$ are asymptotically equal to the moments of $Z$, i.e. that $E(X^m-Z^m) \rightarrow 0$ for each $m$.
We have \begin{eqnarray*} E(X^m) &=& E\left((\sum_p X_p)^m\right) \\ &=& E \left( \sum_{(p_1, \dots, p_m)} X_{p_1} X_{p_2} \dots X_{p_m}\right) \\ &=& \sum_{(p_1, \dots, p_m)} E(X_{p_1} \dots X_{p_m}) \end{eqnarray*} and similarly for $Z$.
We now divide up our $m-$tuples into two classes:
Class 1: For each $i \neq j$ we either have $p_i=p_j$ or $p_i \cap p_j = \emptyset$ (i.e. no two pairs overlap in exactly one point).
Class 2: There is at least one pair $(i,j)$ for which $p_i$ and $p_j$ intersect exactly in one point.
Notice that in class $2$ it's impossible for both $p_i$ and $p_j$ to be matching pairs (If $(x_1, y_1)$ and $(x_2, y_2)$ match, it's impossible for $(x_1, y_1)$ and $(x_3, y_3)$ to match). So the contribution of Class $2$ to the expectation of $X$ is $0$. The bound we want on the expectation would follow from the following two claims:
Claim 1: For any tuple in class $1$ we have $E(X_{p_1} X_{p_2} \dots X_{p_n}) = E(Z_{p_1} Z_{p_2} \dots Z_{p_n}) (1+o(1))$. (So the total contribution of Class $1$ to $E(X)$ and $E(Z)$ is roughly the same).
Claim 2: The total contribution of class $2$ to $E(Z)$ is $o(1)$.
The intuition here is that because $m$ is fixed and the size of the grid tends to infinity, most collections of points won't overlap at all, so class $2$ isn't that significant.
For claim $1$: If $r$ is the number of distinct pairs in $(p_1, p_2, \dots, p_m)$, then we have \begin{eqnarray*} E(Z_{p_1} Z_{p_2} \dots Z_{p_m}) &=& \left(\frac{1}{4n^2-1}\right)^r \\ E(X_{p_1} X_{p_2} \dots X_{p_m}) &=& \left(\frac{1}{4n^2-1}\right)\left(\frac{1}{4n^2-3}\right) \dots \left(\frac{1}{4n^2-(2r-1)}\right)\end{eqnarray*} (Knowing certain pairs match makes it slightly more likely that other pairs match too).
In particular, we have $$E(Z_{p_1} Z_{p_2} \dots Z_{p_m}) \leq E(X_{p_1} X_{p_2} \dots X_{p_m}) \leq E(Z_{p_1} Z_{p_2} \dots Z_{p_m})\left(\frac{4n^2-1}{4n^2-2m}\right)^m$$
A sketchier argument for Claim $2$: Let $r$ be as in claim $1$. If we consider all $m$-tuples containing $r$ distinct pairs, then only an $O(\frac{1}{2n^2c_d})$ of them have any overlap between the $r$ pairs. (the denominator here coming from how there's $2n^2c_d$ close pairs to choose from. Again the constant in the $O$ notation can depend on $m$ or $d$). All tuples with $r$ pairs have the same contribution to $E(Z)$, so the contribution from pairs with overlap is only $O(\frac{1}{n^2})$ the contribution from all pairs. The claim follows since $E(Z)=O(1)$.
