I am studying the Hartree equation for N-particles for the first time and things are not clear to me. Given the density matrix
$$\gamma_{N, t}^{(n)} (x_1,..,n_n; y_1,...,y_n) = \begin{cases} \int \gamma_{N,t}(x_1,...,x_n,x_{n+1},...,x_N; y_1,...,y_n,x_{n+1},...,x_N) dx_{n+1}...d_{x_N},\,\,\,\,\, 1 \leq n \le N\\ 0,\,\,\,\,\, o.w. \end{cases}$$
by which we define the Wigner transform
$$w_N(x;v):= \frac{1}{(3 \pi)^{3 N}} \int e^{- i v. y} \gamma_N \left( x + \frac{y}{2}, x- \frac{y}{2}\right) dy.$$
and the rescaled Wigner transform
$$w_{N,\varepsilon}(x;v)=\varepsilon^{-3N}w_N(x;v):= \frac{1}{(3 \pi)^{3 N}} \int e^{- i v. y} \gamma_N \left( x + \frac{\varepsilon y}{2}, x- \frac{\varepsilon y}{2}\right) dy.$$
The author says that using the Heisenberg equation: $$i\varepsilon \partial_t \gamma_{N,t} =[H_N, \gamma_{N,t}],\,\,\,[A,B]=AB-BA,$$
And $H_N$ is the Hamiltonian of the N particle system, we can (easily) derive the time evolution of Wigner transform, namely
$$\partial_t W_N (t; x,v) +\sum_{j=1}^N v_j . \triangle_{x_j} W_N (t;x,v)=\frac{- i \varepsilon^2}{(2 \pi)^{3N}} \int \left[ U\left(x+ \frac{\varepsilon y}{2}\right)- U\left(x- \frac{\varepsilon y}{2} \right) \right] e^{i y (u-v)} W_N(t; x,v) du dy.$$
I can not find how the author gets this result. Could you please enlighten me?