Integral of one error function and Gaussian
Mathematica does not evaluate the integral. Surprisingly, there is an exact result. By changing variables, $z-\mu=x$, your integral of $\operatorname{erf} \times \text{Gaussian}$ may be written as
$$
I(b)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(ax+b) \exp \left(-cx^2 \right)
$$
Since $\partial_t \operatorname{erf}(t)= \frac{2}{\sqrt{\pi}}e^{-t^2}$, differentiating with respect to $b$ produces
$$
\partial_b I(b)=\int\limits_{-\infty}^\infty dx \ \frac{2}{\sqrt{\pi}} \exp\left( -(ax+b)^2-cx^2 \right)
$$
This is just a Gaussian integral over $x$, the result is
$$
\partial_b I(b)= \frac{2}{\sqrt{a^2+c}} \exp \left( -\frac{b^2 c}{a^2 +c} \right)
$$
Now we integrate to recover $I(b)$
$$
I(b)-I(0)=\int\limits_0^b \ dt \frac{2}{\sqrt{a^2+c}} \exp \left( - t^2\frac{ c}{a^2 +c} \right)
$$
$I(0)=\int dx \ ( \text{odd function})=0$, and the integral on the right is simply the definition of $\operatorname{erf}$, so we have
$$
I(b)=\sqrt{\frac{\pi}{c}} \operatorname{erf} \left( b \sqrt{\frac{c}{a^2+c}} \right)
$$
$$
\tag*{$\blacksquare$}
$$
Integral of two error functions and Gaussian
We may use a similar method to find that
$$
\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x) \operatorname{erf}(a_2 x) e^{-cx^2} = \frac{2}{\sqrt{c\pi}} \tan^{-1}\left( \frac{a_1 a_2}{ \sqrt{c(a_1^2+a_2^2+c)}} \right)
$$
However, once we introduce $b\neq 0$, we meet the more complicated integral in an intermediate step
$$
\int\limits_0^b dt \ \operatorname{erf}(\alpha t + \beta) e^{-\gamma t^2}
$$
Update! Using the result from this post, we can make progress. Let
$$
J(\beta)=\int\limits_{-\infty}^\infty dx \ \operatorname{erfc}(\alpha x+\beta) \operatorname{erf}(a x+b) \exp \left(-cx^2 \right)
$$
Then $J(\infty)=0$. Differentiating
$$
\partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a x+b) \exp \left(-cx^2 - (\alpha x +\beta)^2 \right)
$$
Which we can evaluate using our first result
$$
\partial_{\beta} J(\beta)=-\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}} \operatorname{erf}\left( \left(b-\frac{a \alpha \beta}{c+\alpha^2} \right) \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}} \right) \exp \left( -\beta^2 \left( 1+\frac{\alpha^2}{c+\alpha^2}\right) \right)
$$
With $\Lambda=\frac{2}{\sqrt{\pi}} \sqrt{\frac{\pi}{c+\alpha^2}}$, $A=-\frac{a \alpha }{c+\alpha^2} \sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$, $B=b\sqrt{\frac{c+\alpha^2}{c+\alpha^2+a^2}}$, and $C=\left( 1-\frac{\alpha^2}{c+\alpha^2}\right)$. Now we integrate
$$
J(\infty)-J(\beta)=-\int\limits_\beta^\infty \ dt \Lambda \operatorname{erf}(At+B)e^{-Ct^2}
$$
$$
J(\beta)=\int\limits_\beta^\infty dt \ \Lambda \operatorname{erf}(At+B)e^{-Ct^2}= \frac{\Lambda}{\sqrt{2C}} \int\limits_{\sqrt{2C}\beta}^\infty du \ \operatorname{erf}(Au/\sqrt{2C}+B)e^{-u^2/2}
$$
Which may be written in terms of the 'Generalized Owen-T' function
$$
J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}}T\left( \sqrt{2C}\beta, A/\sqrt{C},\sqrt{2}B\right)
$$
Using identity ii we have
$$
J(\beta)=\frac{2 \sqrt{\pi} \Lambda}{\sqrt{C}} \left[ \frac{1}{2\pi} \left( \tan^{-1}(A/\sqrt{C}) -\tan^{-1}(A/\sqrt{C} + B/\beta\sqrt{C}) \\ -\tan^{-1}(B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} )) \right) +\frac{1}{4} \operatorname{erf}(B (1+A^2/C)^{-1/2}) \\ +T(\beta \sqrt{2C},(A\beta^{-1}+B)/\beta \sqrt{C})+T(\sqrt{2}B (1+A^2/C)^{-1/2}, B^{-1}(\beta \sqrt{C} +AB/\sqrt{C} + A^2 \beta /\sqrt{C} ) \right]
$$
Where $T$ is the Owen T function. I've checked it numerically, I hope there are no typos here.
Approximate integral of two or more error functions with Gaussian
We can approximate the integral using Laplace's method. I'll demonstrate with two error functions, but in principle you could have more. Let
$$
I(c)=\int\limits_{-\infty}^\infty dx \ \operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) e^{-cx^2}
$$
Expand the product of $\operatorname{erf}$s around the Gaussian peak: $x=0$
$$
\operatorname{erf}(a_1 x+ b_1) \operatorname{erf}(a_2 x+ b_2) = \operatorname{erf}( b_1) \operatorname{erf}( b_2) + x \left( \text{junk} \right) + x^2 \frac{2 e^{-b_1^2-b_2^2}}{\pi} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) + \mathcal{O}(x^3)
$$
Term by term integration produces the asymptotic series. The $(\text{junk})$ and all $x^{\text{odd}}$ terms integrate to zero. Noting that
$$
\int\limits_{-\infty}^\infty dx \ x^{2n} e^{-cx^2} = c^{-n-1/2} \Gamma(n+1/2)
$$
For fixed $a_i, \ b_i$ we have
$$
I(c) \sim \sqrt{\frac{\pi}{c}} \operatorname{erf}( b_1) \operatorname{erf}( b_2) \\ +\frac{1}{c^{3/2}} \frac{ e^{-b_1^2-b_2^2}}{\sqrt{\pi}} \left( 2 a_1 a_2 -a_2^2 b_2 e^{b_1^2} \sqrt{\pi} \operatorname{erf}(b_1) - a_1^2b_1 e^{b_2^2} \sqrt{\pi} \operatorname{erf}(b_2) \right) \ \ , \ \ c \to \infty
$$
Which is the two term approximation.
Best Answer
$\newcommand{\d}{\,\mathrm{d}}$Take $z$ to be in $(0,\sqrt{N})^2$ and $0<\sigma$. By the multivariate integral substitution rule we have: $$\int_{[0,\sqrt{N})^2}\exp(-\sigma^{-2}\|x-z\|^2)\d\lambda(x)\\=\sigma^2\int_{[-z_1/\sigma,(\sqrt{N}-z_1)/\sigma)\times[-z_2/\sigma,(\sqrt{N}-z_2)/\sigma)}\exp(-\|x\|^2)\d\lambda(x)$$
The difference between this and: $$\sigma^2\int_{\Bbb R^2}\exp(-\|x\|^2)\d\lambda(x)=\pi\sigma^2$$Is exactly: $$\psi(\sigma,N,z):=\sigma^2\int_{\Bbb R^2\setminus[-z_1/\sigma,(\sqrt{N}-z_1)/\sigma)\times[-z_2/\sigma,(\sqrt{N}-z_2)/\sigma)}\exp(-\|x\|^2)\d\lambda(x)$$We can bound: $$\begin{align}0&<\psi(\sigma,N,z)\\&\le\sigma^2\int_{\{|u|\ge\mu/\sigma\}\times\{|v|\ge\nu/\sigma\}}\exp(-\|(u,v)\|^2)\d\lambda(u,v)\\&\le4\sigma^2\int_{(\kappa/\sigma,\infty)^2}\exp(-\|(u,v)\|^2)\d\lambda(u,v)\end{align}$$Where $0<\mu:=\min(z_1,\sqrt{N}-z_1),\,0<\nu:=\min(z_2,\sqrt{N}-z_2)$ and $\kappa:=\min(\mu,\nu)>0$. We can continue to bound: $$\begin{align}\psi(\sigma,N,z)&\le4\sigma^2\exp(-2\kappa^2/\sigma^2)\int_{(\kappa/\sigma,\infty)^2}\exp(-(u^2-\kappa^2/\sigma^2+v^2-\kappa^2/\sigma^2))\d\lambda(u,v)\\&=4\sigma^2\exp(-2\kappa^2/\sigma^2)\left(\int_{\kappa/\sigma}^\infty\exp(-(x^2-\kappa^2/\sigma^2))\d x\right)^2\end{align}$$
If we substitute $y^2=x^2-\kappa^2/\sigma^2$ then $0\le\frac{\d x}{\d y}=\frac{y}{x}=\sqrt{1-\frac{\kappa^2}{\sigma^2x^2}}\le1$ so I can continue to bound: $$0<\psi(\sigma,N,z)\le4\sigma^2\exp(-2\kappa^2/\sigma^2)\left(\int_0^\infty\exp(-y^2)\d y\right)^2=\pi\sigma^2\cdot\exp(-2\kappa^2/\sigma^2)$$
Now - if $z$ is "pretty close" to the middle of the region, say we've fixed a constant $0<\alpha<1$ and $\alpha\sqrt{N}\le\min(z_1,z_2,\sqrt{N}-z_1,\sqrt{N}-z_2)$ and if $\sigma:=\delta\sqrt{N}$ for some defined $\delta>0$ then you've got: $$0<\left|\pi\sigma^2-\int_D\cdots\right|\le\pi\sigma^2\cdot\exp(-2\alpha^2/\delta^2)$$If $\delta\ll\alpha$ the right hand side is extremely small indeed. What's more, we're actually saying the right hand side is small even if you divide through by $\sigma^2$ (which is also very small). So the error clearly decays very fast, and the paper's heuristic is about as justified as a heuristic can be.