The point probability of a random variable I'm trying to understand the following integral:
$\begin{align*}
\mathcal P(draw \mid \boldsymbol \mu, \boldsymbol \Sigma) 
&= 
\lim_{{ \epsilon}\rightarrow 0}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z
\\&= 
\mathcal N \left(
0;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) 
\end{align*}$
Which appears in J. Moser, “The Math Behind TrueSkill”.
To clarify, this is the integration of a multivariate Gaussian distribution. The vector $\mathbf z$ between the bounds $-\epsilon\mathbf 1$ and $-\epsilon\mathbf 1$, where $\mathbf 1$ is just a vector of same dimension as $\mathbf z$ filled with $1$ meaning the elements of $\mathbf z$ individually all are bounded from $-\epsilon$ to $\epsilon$.
Clearly the normal distribution is being written as $\mathcal N(variable; mean, variance)$ which is more commonly in my experience seen as $\mathcal N(variable \mid mean, variance)$.
This integral of course is of a Normal PDF, and thus is the CDF, or calculating the probability that all the elements of $\mathbf z$ fall between $-\epsilon$ and $\epsilon$. I'm good with that. We could of course also be tempted to write this integral as:
$\Phi\left(
\mathbf \epsilon\mathbf 1;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) -
\Phi\left(
\mathbf -\epsilon\mathbf 1;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)$
If we lean on fairly ordinary CDF notation, but would be mistaken as is clearly spelled out for the 2D case here:
http://orfe.princeton.edu/~rvdb/245/Fall14/lectures/chap3.pdf
The details of the mean and variance aren't immediately relevant to my question (I understand what they all are, whence they stem and what they mean). 
The problem I am grappling with is how to understand the limit. As $\epsilon$ tends to $0$ I expect the integral too to diminish to zero, as the $\mathcal P(-\epsilon)$ and $\mathcal P(\epsilon)$ converge toward the same value of $\mathcal P(0)$
But the stated result is not $0$, it is $\mathcal P(0)$.
I'm confused on two fronts I think. Firstly how it is that this limit does not yield $0$, but secondly how the integral (a CDF) becomes a point probability (PDF).
And lastly I guess, central to my confusion is how $\mathcal P(x)$ is not equal to $0$ for any given x for a continuous random variable. It translates to the probability of getting exactly one specific value in the real number space and one might intuit that to be effectively $0$.
To wit, I conclude that there is something I am not understanding about the interpretation of a PDF and the probability of a given point in the real number space.
 A: The first quoted equation is certainly wrong -- that limit is zero, as you expected. However, apparently the source article actually uses the following "normalized match quality":
$$\begin{align} 
q_{\text draw}(\mathbf{ \mu,\Sigma,\beta,A}):&=\lim_{{ \epsilon}\rightarrow 0}{
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z
\over
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf 0,
\mathbf  \beta^2 \mathbf A^\intercal\mathbf A
\right) d\mathbf z
}\tag{1a}\\[2ex]
&=\lim_{{ \epsilon}\rightarrow 0}{{1\over (2\epsilon)^d}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z
\over {1\over (2\epsilon)^d}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf 0,
\mathbf  \beta^2 \mathbf A^\intercal\mathbf A
\right) d\mathbf z
}\tag{1b}\\[2ex]
&={\lim_{{ \epsilon}\rightarrow 0}{1\over (2\epsilon)^d}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z
\over \lim_{{ \epsilon}\rightarrow 0}{1\over (2\epsilon)^d}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf 0,
\mathbf  \beta^2 \mathbf A^\intercal\mathbf A
\right) d\mathbf z
}\tag{1c}\\[2ex]
&={
\mathcal N \left(
\mathbf 0;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) 
\over
\mathcal N \left(
\mathbf 0;
\mathbf 0,
\mathbf  \beta^2 \mathbf A^\intercal\mathbf A
\right)\tag{1d}
}\end{align}
$$
where $\mathcal N$ is the PDF of a $d$-dimensional multivariate Normal distribution.
Note that (1c) is because the limit of a ratio is the ratio of limits (the denominator being nonzero), and (1d) is a consequence of the Mean Value Theorem for Multiple Integrals, one version of which is as follows (see $^*$ below for links to proofs):

If $f:S\to\mathbb{R}$ is a continuous Riemann integrable function on a connected Jordan measurable set $S\subset\mathbb{R}^d$, then there exists $\mathbf z_0\in S$ such that $$\int_{S}f\mathbf{\,dz}=f(\mathbf{z_0})\,\int_{S}\mathbf{\,dz} $$

This applies to the numerators of (1), for example, with $f(\mathbf{z})=\mathcal N \left(
\mathbf{z}\ ;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right),$ and with $S=S_\epsilon$ (now depending on $\epsilon$) a $d$-dimensional cube of edge-length $2\epsilon$ centered on the origin $\mathbf{0}$. For all $\epsilon>0$ as $\epsilon\to 0,$ the theorem ensures that there exists a $\mathbf z_\epsilon\in S_\epsilon$ such that
$$\begin{align}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z
&= 
\mathcal N \left(
\mathbf z_\epsilon;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)\,\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1} d\mathbf z
\\[2ex] &= 
\mathcal N \left(
\mathbf z_\epsilon;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)\,(2\epsilon)^d .\end{align}$$
Then, since $S_\epsilon\to\{\mathbf{0}\}$ (and hence $\mathbf z_\epsilon\to\mathbf 0$) as $\epsilon\to 0,$
$$\lim_{{ \epsilon}\rightarrow 0}{1\over (2\epsilon)^d}
\int_{-\epsilon \mathbf 1}^{\epsilon \mathbf 1}
\mathcal N \left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) d\mathbf z\\ =\lim_{{ \epsilon}\rightarrow 0}{1\over (2\epsilon)^d}\ \mathcal N \left(
\mathbf z_\epsilon;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)\,(2\epsilon)^d\\
=\lim_{{ \epsilon}\rightarrow 0} \mathcal N \left(
\mathbf z_\epsilon;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)\\
=\mathcal N \left(
\mathbf 0;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) $$
and similarly for the denominators in (1). 
Here's a Wolfram Alpha link for the one-dimensional version. 

$^*$ Here are some links to proofs of various versions of the theorem (I used Apostol's): 
Mathematical Analysis (Apostol, Theorem 14.16, p.401)
Is there a mean-value theorem for volume integrals?
The Mean Value Theorem for Multiple Integrals

NB: The following is too long for a comment: You wrote ...

We could of course also write this integral as:$$\Phi\left(
\mathbf \epsilon\mathbf 1;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right) -
\Phi\left(
\mathbf -\epsilon\mathbf 1;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)$$ If we lean on fairly ordinary CDF notation. 

To see that's not correct, see here for a nice graphical demonstration in the two-dimensional case. E.g., if we have $\mathbf{z}=(z_1,z_2)$, and define
$$F(z_1,z_2):=\Phi\left(
\mathbf z;
\mathbf A^\intercal\boldsymbol \mu,
\mathbf A^\intercal (\beta^2\mathbf I+\boldsymbol \Sigma)\mathbf A
\right)$$
 then the probability of the square $(-\epsilon,\epsilon)\times(-\epsilon,\epsilon)$ is not generally $F(\epsilon,\epsilon)-F(-\epsilon,-\epsilon)$, but rather $$F(\epsilon,\epsilon)+F(-\epsilon,-\epsilon)-F(-\epsilon,\epsilon)-F(\epsilon,-\epsilon).$$ 
