Define using $H(t) = \int (u_t)^2 + |\nabla u|^2 ~dx $ the standard energy. Taking the time derivative you find $$ \frac{d}{dt}H(t) = 2 \int u_t( g + \frac{2}{t} u_t)$$ Writing $\|\cdot \|$ for the $L^2$ integral, you have then $$ \frac{d}{dt} H(t) \geq - 2 \|u_t\|\cdot \|g\| + \frac{4}{t} \|u_t\|^2 \tag{A}$$ by Cauchy-Schwarz, and then by AM-GM you get $$ \frac{d}{dt} H(t) \geq - \frac{4}{t} \|u_t\|^2 - \frac{t}{4} \|g\|^2 + \frac{4}{t} \|u_t\|^2 = - \frac{t}{4} \|g\|^2 $$ So integrating you find $$ H(1) + \int_\epsilon^1 \frac{t}{4} \|g\|^2 ~dt \geq H(\epsilon) $$

So if you assume that $$\tag{*} \int_0^1 \int_{\mathbb{R}^3} |g|^2 t ~ dx~dt $$is finite, the above shows that $H(\epsilon)$ is uniformly bounded as $\epsilon \searrow 0$ and hence what you defined as the weighted energy $E(\epsilon)$ converges to zero at rate $\epsilon^2$.

If you only assume (*), then it is impossible to prove that the energy decays with any polynomial rate. This is because of the ODE examples

$$ y(t) = t^{1+\alpha}$$

with $\alpha\in (0,1)$. You find that

$$ \ddot{y} - \frac{2}{t} \dot{y} = (1+\alpha)(\alpha - 2) t^{\alpha - 1} =: g(t)$$

and this $g(t)$ satisfies $\int_0^1 t |g(t)|^2 ~dt < \infty$. So the energy ($(\dot{y})^2$) cannot be proven to go to zero as $t^\beta$ for any $\beta > 0$. By spatial truncation this ODE example can be upgraded to a bona fide solution of the wave equation you wrote down.

On the other hand, if instead of (*) you know something stronger about $g$, then some decay may be concluded.

For example: suppose you know that $\int_0^1 \|g\|^2 ~dt$ is bounded. Then starting from (A) you can instead find $$ \frac{d}{dt} H(t) \geq \left( \frac4t-1\right) \|u_t\|^2 - \|g\|^2 $$ This would give you $$ H(1) + \int_0^1 \|g\|^2 ~dt \geq H(\epsilon) + \int_{\epsilon}^1 \frac{3}t\|u_t\|^2 ~dt $$ This gives you an integrated decay of time component of the standard energy.

But this would depend on what you want to know and what you are willing to provide for (in terms of properties of $g$).

Define using $H(t) = \int (u_t)^2 + |\nabla u|^2 ~dx $ the standard energy. Taking the time derivative you find $$ \frac{d}{dt}H(t) = 2 \int u_t( g + \frac{2}{t} u_t)$$ Writing $\|\cdot \|$ for the $L^2$ integral, you have then $$ \frac{d}{dt} H(t) \geq - 2 \|u_t\|\cdot \|g\| + \frac{4}{t} \|u_t\|^2 \tag{A}$$ by Cauchy-Schwarz, and then by AM-GM you get $$ \frac{d}{dt} H(t) \geq - \frac{4}{t} \|u_t\|^2 - \frac{t}{4} \|g\|^2 + \frac{4}{t} \|u_t\|^2 = - \frac{t}{4} \|g\|^2 $$ So integrating you find $$ H(1) + \int_\epsilon^1 \frac{t}{4} \|g\|^2 ~dt \geq H(\epsilon) $$

So if you assume that $$\tag{*} \int_0^1 \int_{\mathbb{R}^3} |g|^2 t ~ dx~dt $$is finite, the above shows that $H(\epsilon)$ is uniformly bounded as $\epsilon \searrow 0$ and hence what you defined as the weighted energy $E(\epsilon)$ converges to zero at rate $\epsilon^2$.

If you only assume (*), then it is impossible to prove that the standard energy decays with any rate like $\epsilon^\beta$. This is because of the ODE examples

$$ y(t) = t^{1+\alpha}$$

with $\alpha\in (0,1)$. You find that

$$ \ddot{y} - \frac{2}{t} \dot{y} = (1+\alpha)(\alpha - 2) t^{\alpha - 1} =: g(t)$$

and this $g(t)$ satisfies $\int_0^1 t |g(t)|^2 ~dt < \infty$. So the energy ($(\dot{y})^2$) cannot be proven to go to zero as $t^\beta$ for any $\beta > 0$. By spatial truncation this ODE example can be upgraded to a bona fide solution of the wave equation you wrote down.

On the other hand, if instead of (*) you know something stronger about $g$, then some decay may be concluded.

For example: suppose you know that $\int_0^1 \|g\|^2 ~dt$ is bounded. Then starting from (A) you can instead find $$ \frac{d}{dt} H(t) \geq \left( \frac4t-1\right) \|u_t\|^2 - \|g\|^2 $$ This would give you $$ H(1) + \int_0^1 \|g\|^2 ~dt \geq H(\epsilon) + \int_{\epsilon}^1 \frac{3}t\|u_t\|^2 ~dt $$ This gives you an integrated decay of time component of the standard energy.

But this would depend on what you want to know and what you are willing to provide for (in terms of properties of $g$).