I am reading the following proof:
I do not understand where the equality
$$\int_U f(x)\frac{\phi(x+he_i)-\phi(x)}{h}dx=-\int_U g_i^h(x)\phi(x+he_i)dx$$
is coming from. If we were integrating over all of $\mathbb{R}^n$, then it would follow by Fubini + an integration by parts + a change of variables, but the change of variables shifts the domain $U$, so I don't think that works in this case.
Also, in the last step, when he takes the limit, why is the convergence justified? We knwo that $\int_U g^h_i \psi dx \to \int_U g_i \psi dx$ for every $\psi\in C^1_c(V)$, but $\phi(x+he_i)$ also depends on $h$.
EDIT: I think the convergence follows by adding and substracting $\phi(x)$ and using the uniform continuity of $\phi$ along with the boundedness of $g_i^h$, but my first question remains.
The idea is to split the integral and use change of variables in one of them: \begin{align} \int_U f(x)\frac{\phi(x+he_i)-\phi(x)}{h}dx &= \frac{1}{h}\int_U f(x)\phi(x+he_i)dx - \frac{1}{h} \int_U f(x)\phi(x)dx \\ & = \frac{1}{h}\int_{U} f(x)\phi(x+he_i)dx - \frac{1}{h} \int_{U-he_i} f(x+he_i)\phi(x+he_i)dx \\ & = \frac{1}{h}\int_{U} f(x)\phi(x+he_i)dx - \frac{1}{h} \int_{U} f(x+he_i)\phi(x+he_i)dx \\ &= -\int_U g_i^h(x)\phi(x+he_i)dx \end{align}
In the second equality we used the change of variables $x \mapsto x+he_i$ and in the third, we used that if $h$ is small enough $supp \phi \subset U\cap(U-he_i)$.
For the last step, it follows as you mentioned.