[step:Expand the entropy to second order]
Define the scalar function $\Phi: \mathbb{R} \to \mathbb{R}$ by setting $\Phi(0) := 0$ and, for $s \neq 0$,
\begin{align*}
\Phi(s) := s^2\log(s^2).
\end{align*}
We use the following little-o notation: for a real-valued function $r$ defined near $0$, the expression $r(a)=o(a^2)$ as $a \to 0$ means $r(a)/a^2 \to 0$, and $r(a)=o(1)$ as $a \to 0$ means $r(a) \to 0$.
For $s$ near $0$, Taylor expansion at $1$ gives
\begin{align*}
\Phi(1+s) = 2s + 3s^2 + s^2\rho(s)
\end{align*}
where $\rho: (-\delta,\delta) \to \mathbb{R}$ is a function satisfying $\rho(s) \to 0$ as $s \to 0$ for some $\delta>0$. Let $M := \sup_{x\in\mathbb{R}^n}|u(x)|$. For $|\varepsilon|M<\delta$, the remainder satisfies
\begin{align*}
\sup_{x\in\mathbb{R}^n}|\rho(\varepsilon u(x))| \to 0
\end{align*}
as $\varepsilon \to 0$. Since $\nu$ is a probability measure and $u$ is bounded, this uniform Taylor remainder gives
\begin{align*}
\int_{\mathbb{R}^n} \varepsilon^2u(x)^2\rho(\varepsilon u(x)) \, d\nu(x)=o(\varepsilon^2).
\end{align*}
Therefore, integrating the expansion with $s=\varepsilon u(x)$ gives
\begin{align*}
\int_{\mathbb{R}^n} g_\varepsilon^2 \log(g_\varepsilon^2) \, d\nu(x)
= 2\varepsilon \int_{\mathbb{R}^n} u \, d\nu(x)
+ 3\varepsilon^2 \int_{\mathbb{R}^n} u^2 \, d\nu(x)
+ o(\varepsilon^2).
\end{align*}
The mean-zero condition cancels the linear term, so
\begin{align*}
\int_{\mathbb{R}^n} g_\varepsilon^2 \log(g_\varepsilon^2) \, d\nu(x)
= 3\varepsilon^2 \int_{\mathbb{R}^n} u^2 \, d\nu(x)
+ o(\varepsilon^2).
\end{align*}
Next,
\begin{align*}
\int_{\mathbb{R}^n} g_\varepsilon^2 \, d\nu(x)
= \int_{\mathbb{R}^n} (1+2\varepsilon u+\varepsilon^2u^2) \, d\nu(x)
= 1+\varepsilon^2\int_{\mathbb{R}^n} u^2 \, d\nu(x).
\end{align*}
Using $\log(1+t)=t+o(t)$ as $t \to 0$, we obtain
\begin{align*}
\left(\int_{\mathbb{R}^n} g_\varepsilon^2 \, d\nu(x)\right)
\log\left(\int_{\mathbb{R}^n} g_\varepsilon^2 \, d\nu(x)\right)
= \varepsilon^2 \int_{\mathbb{R}^n} u^2 \, d\nu(x) + o(\varepsilon^2).
\end{align*}
Subtracting the two expansions gives
\begin{align*}
\operatorname{Ent}_\nu(g_\varepsilon^2)
= 2\varepsilon^2 \int_{\mathbb{R}^n} u^2 \, d\nu(x) + o(\varepsilon^2).
\end{align*}
[/step]