I am trying to understand a proof in "Centroidal voronoi tessellations: applications and algorithms" from Du et al. In particular, I do not understand Proposition 3.1, stating that a necessary condition for an "optimal" Voronoi tesselation its generators need be equal to the centroids of the cell they define. It begins with a first order variation of the function F, then we should divide by $\epsilon$ and construct the limes. Unfortunately, I do not understand how we get this way from the formula for the first order variation to $z_j$.
Copied from the text, the step I do not understand is going from: $F(z_j + \epsilon v) - F(z_j) = \int_{y \in Y_j} p(y) |y - z_j - \epsilon v|^2 - |y - z_j|^2 dy$ to $z_j = \frac{\int_{y \in Y_j} yp(y)dy}{\int_{y \in Y_j} p(y) dy}$.
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} } $$