[guided]The point of this test function is to move one discrete derivative from the test function onto the coefficient-gradient product $a_{ij}\partial_{x_j}u$. Fix a coordinate direction $e_k$ and a scale $h$ with $0<|h|<\rho$. We use
\begin{align*}
\phi_h:=-D_{k,-h}(\zeta^2D_{k,h}u).
\end{align*}
The product $\zeta^2D_{k,h}u$ is defined on $K$, where the translate $x+he_k$ lies in $U'$, and then extended by zero to $U'$. Since $\operatorname{supp}(\zeta^2D_{k,h}u)\subset K$ and, for $0<|h|<\rho$, the sets $K-he_k$ and $K+he_k$ are contained in $U'$, this zero extension belongs to $H^1_0(U')$. With the convention $D_{k,-h}\psi(x)=(\psi(x-he_k)-\psi(x))/(-h)$, the backward difference $D_{k,-h}(\zeta^2D_{k,h}u)$ can be nonzero only when $x\in K$ or $x-he_k\in K$; hence its support is contained in $K\cup(K+he_k)$, a compact subset of $U'$. Differentiating does not enlarge support, so the same translated-support condition holds for $\partial_{x_i}D_{k,-h}(\zeta^2D_{k,h}u)$.
The weak equation gives
\begin{align*}
\sum_{i,j=1}^n \int_{U'} a_{ij}(x)\partial_{x_j}u(x)\partial_{x_i}\phi_h(x)\,d\mathcal{L}^n(x)=\int_{U'} f(x)\phi_h(x)\,d\mathcal{L}^n(x).
\end{align*}
Since weak derivatives commute with difference quotients on these compactly supported translated functions, we have
\begin{align*}
\partial_{x_i}\phi_h=-D_{k,-h}\partial_{x_i}(\zeta^2D_{k,h}u).
\end{align*}
The discrete integration-by-parts identity applies because $\partial_{x_i}(\zeta^2D_{k,h}u)$ and its required translate have compact support in $U'$. It says that a backward difference can be moved from one factor to the other as a forward difference, with a minus sign. The minus sign in the definition of $\phi_h$ was chosen so that the leading term has the positive coercive sign. Thus
\begin{align*}
\sum_{i,j=1}^n \int_{U'} D_{k,h}\!\left(a_{ij}\partial_{x_j}u\right)(x)\partial_{x_i}(\zeta^2D_{k,h}u)(x)\,d\mathcal{L}^n(x)=-\int_{U'} f(x)D_{k,-h}(\zeta^2D_{k,h}u)(x)\,d\mathcal{L}^n(x).
\end{align*}
Now expand the discrete derivative of the product. With $w_h=D_{k,h}u$,
\begin{align*}
D_{k,h}(a_{ij}\partial_{x_j}u)=(\tau_{k,h}a_{ij})\partial_{x_j}w_h+(D_{k,h}a_{ij})\partial_{x_j}u.
\end{align*}
The first summand is the desired one: it contains $\partial_{x_j}w_h$, hence second-order information about $u$. Since
\begin{align*}
\partial_{x_i}(\zeta^2w_h)=\zeta^2\partial_{x_i}w_h+2\zeta(\partial_{x_i}\zeta)w_h,
\end{align*}
the highest-order part is
\begin{align*}
\sum_{i,j=1}^n \int_{U'} \zeta(x)^2(\tau_{k,h}a_{ij})(x)\partial_{x_j}w_h(x)\partial_{x_i}w_h(x)\,d\mathcal{L}^n(x).
\end{align*}
For each $x$ in the integration region, the coefficient matrix is evaluated at $x+he_k$. The ellipticity hypothesis is valid at that point, so with $\xi=\nabla w_h(x)$ it yields
\begin{align*}
\sum_{i,j=1}^n(\tau_{k,h}a_{ij})(x)\partial_{x_j}w_h(x)\partial_{x_i}w_h(x)\ge \theta|\nabla w_h(x)|^2.
\end{align*}
Multiplying by $\zeta(x)^2$ and integrating gives the coercive lower bound
\begin{align*}
\theta\int_{U'} \zeta(x)^2|\nabla w_h(x)|^2\,d\mathcal{L}^n(x).
\end{align*}[/guided]