Let $(M^n,g)$ be a closed smooth Riemannian manifold. Consider a set $\mathcal B_{\epsilon}$ which consists of a maximal number of points in $M$ with pairwise distance no smaller than $\epsilon$.
We define the measure $$ m_{\epsilon}:=\frac{\sum_{x\in \mathcal B_{\epsilon}} \delta_{x}}{|B_{\epsilon}|}. $$ Can we prove that as $\epsilon \to 0$, $m_\epsilon$ converges weakly to the $n$-dimensional Hausdorff measure up to a rescaling?
I assume that by maximal you mean with respect to inclusion. Then the answer is no. Consider the following counterexample on the real line:
Let $\mathcal{B}_\epsilon := \epsilon\mathbb{Z}$ and $\widetilde{\mathcal{B}}_\epsilon := 3/2 \epsilon \mathbb{Z}$. Both sets are maximal with respect to inclusion for any $\epsilon$, but for the corresponding measures you will get $m_\epsilon \rightharpoonup \frac{1}{2} \mathcal{H}^1$ and $\widetilde{m}_\epsilon \rightharpoonup \frac{1}{3} \mathcal{H}^1$. So if you have a third sequence $\widehat{\mathcal{B}}_\epsilon$ that alternates between the two, the corresponding measures will not converge.
What you should be able to get though by compactness is a weakly converging subsequence to something of the form $a(x)d\mathcal{H}^n$, with constants $0<c_1 \le a(x) \le c_2 <\infty$ corresponding to the optimal and "worst" sphere packing in the respective dimension.
Edit: Now that the question has been clarified to mean maximal number¹ of points, the answer is yes. We can prove this by calculating upper/lower bounds that result in the constants $c_1$, $c_2$ above coinciding.
Define the maximal packing density in the limit by $$(*)\quad c := \limsup_{R \to \infty} \frac{|B_1|}{|B_R|}\max\{\#A: A \subset B_R \wedge x,y \in A \Rightarrow |x-y| \geq 1 \}.$$
Obviously $c$ is a well defined number between $0$ and $1$. Now I am going to assume that $M^n$ is flat² and assume that $0 < \epsilon \ll \delta \ll 1$. Then I claim that for any ball $B_\delta$ of radius $\delta$, we have $$ c \left(\frac{\delta}{\delta-\epsilon}\right)^n \stackrel{(i)}\leq \frac{\#(\mathcal{B}_\epsilon \cap B_\delta)}{|B_\delta||B_\epsilon|}=\frac{m_\epsilon(B_\delta)}{|B_\delta|}\stackrel{(ii)}\leq c\left(\frac{\delta + \epsilon}{\delta}\right)^n.$$
Assume that (i) is not the case. Then we could find an $R>0$, a corresponding $A$ in $(*)$ and a Ball $B_{(\delta-\epsilon)/\epsilon} \subset B_R$ such that $\# (B_{(\delta-\epsilon)/\epsilon} \cap A) > \# (\mathcal{B}_\epsilon \cap B_\delta)$. This is a contradiction, because we could replace the points in one set with those in the other, rescaling them with a factor of $\epsilon$. (Because we pick a radius $\delta-\epsilon$, the distances to the other points are no problem.) Similarly (ii) needs to hold, otherwise we could replace points in $(*)$.
But then for any converging subsequence $m_{\epsilon_k} \rightharpoonup m$ in the sense of measures, we know that $\frac{m(B_\delta)}{|B_\delta|} =c$, which implies $m = c \mathcal{H}^n$. As any subsequence has a (locally) converging subsequence by compactness and the limit is always the same, this in turn implies that $m_\epsilon \rightharpoonup c \mathcal{H}^n$.
¹Which either means that $M^n$ is compact or could be understood as something like "If we remove N points from $\mathcal{B}_\epsilon$, we can add at most N points again." for the proof to work. With the latter I am not entirely sure that such a set needs to exist.
²If it is not, the argument still works, but one needs to increase/decrease $\epsilon$ by a $\delta$-dependent factor when comparing with the flat case, as the distances change slightly. The resulting error should vanish for $\delta \to 0$, which is the limit we need to identify $m=c \mathcal{H}^n$ anyway.
