Expanding the norm squared yields $\int_{y\in V_j} \rho({\bf y})\left( \left|{\bf y} -{\bf z}_j - \epsilon {\bf v}\right|^2 - \left|{\bf y} -{\bf z}_j\right|^2 \right) d{\bf y}$
\begin{align}
& = \int_{y\in V_j} \rho({\bf y})\left( \left|{\bf y-{\bf z}_j}\right|^2 + \epsilon^2 \left|{\bf v}\right|^2 -2\epsilon \left({\bf y-{\bf z}_j}\right)\cdot{\bf v} - \left|{\bf y-{\bf z}_j}\right|^2\right) d{\bf y} \\
& = \int_{y\in V_j} \rho({\bf y}) \left( \epsilon^2 \left|{\bf v}\right|^2 -2\epsilon \left({\bf y-{\bf z}_j}\right)\cdot{\bf v}\right) d{\bf y}.
\end{align}
When dividing by $\epsilon$ and taking the limit, the $\epsilon^2$ term drops out giving,
\begin{align}
\lim_{\epsilon \rightarrow 0} \frac{{\mathcal F}({\bf z}_j + \epsilon {\bf v}) - {\mathcal F}({\bf z}_j) }{\epsilon} & = -2 \int_{y\in V_j} \rho({\bf y}) \left({\bf y-{\bf z}_j}\right)\cdot{\bf v} d{\bf y}\\
& = -2{\bf v} \cdot \int_{y\in V_j} \rho({\bf y}) \left({\bf y}-{\bf z}_j\right) d{\bf y}.
\end{align}
If ${\bf z}_j$ are minimizers of this functional, then this derivative must be zero for any possible vector ${\bf v}$. Thus
$$
{\bf 0} = -2 \int_{y\in V_j} \rho({\bf y}) \left({\bf y-{\bf z}_j}\right) d{\bf y}.
$$
Rearranging terms gives,
\begin{align}
\int_{y\in V_j} \rho({\bf y}){\bf y} d{\bf y} & = \int_{y\in V_j} \rho({\bf y}){\bf z}_j d{\bf y} \\
& = \left(\int_{y\in V_j} \rho({\bf y}) d{\bf y} \right){\bf z}_j.
\end{align}
Dividing by $\int_{y\in V_j} \rho({\bf y}) d{\bf y}$ gives the expression reported in the paper,
$$
{\bf z}_j = \frac{ \int_{y\in V_j} \rho({\bf y}){\bf y} d{\bf y} } { \int_{y\in V_j} \rho({\bf y}) d{\bf y} }
$$