Elliptic regularization of the heat equation This is from PDE Evans, 2nd edition: Chapter 8, Exercise 3:

The elliptic regularization of the heat equation is the PDE $$ u_t - \Delta u -\epsilon u_{tt}=0 \quad \text{in }U_T, \tag{$*$}$$ where $\epsilon > 0$ and $U_T = U \times (0,t]$. Show that $(*)$ is the Euler-Lagrange equation corresponding to an energy functional $I_\epsilon[w] := \iint_{U_T} L_\epsilon(Dw,w_t,w,x,t) \, dx \, dt$.
(Hint: Look for a Lagrangian with an exponential term involving $t$.)

In my previous question with a different EL equation, I had to find a Lagrangian in terms of only three variables (that is, find $L=L(Du,u,x)$).
Now for this problem, I have to find a Lagrangian that is based with five variables and epsilon (that is, find $L=L_\epsilon(Du,u_t,u,x,t)$).
When the hint says "look for a Lagrangian with an exponential term involving $t$, do they mean something like $$e^{-t}(|Du|^2 +\epsilon u_t)$$ or similar?
Now, I got confused again because, unlike my previous question, for this one I am asked to work with the time variable $t$ sin addition to the spatial variable $x$. If the exponential term involving $t$ is indeed $e^{-t}$ times something, then I'm assuming upon differentiation with respect to $t$, product rule will be used to find our way back to the given EL equation.
If any answerers have any ideas, a good hint works best for me.
 A: Edit2:Let's start with the following energy functional
$$ I(u) =\int_{U \times (0,T]} \frac{e^{- t/\epsilon}}{2} ( | D u|^2 + \epsilon | u_t | ^2 ) dx dt $$
Recall that the Euler-Lagrange Equation is just an critical solution to the first variation $\delta I$. The first variation is defined as
$$ \delta I = \lim_{h \to 0}\frac{ I( u + h \xi ) - I (u) }{h} $$
and a critical solution is $u_*$ s.t. $\delta I(u_*) =0$ for all test functions $\xi$.This amounts to finding the first order term in $h$. The expansion computation goes as follows: 
$$ |D ( u + h \xi ) |^2 = | Du + h D \xi |^2 = (Du + h D \xi) \cdot (D u + h  D \xi ) =|Du|^2 + 2h Du \cdot D \xi + h^2 |D\xi|^2 $$
Then similarly
$$ |\frac{d}{dt} ( u + h \xi ) |^2 = |u_t + h \xi _t |^2 = |u_t|^2 + 2h u_t \cdot \xi_t + h^2 |D \xi|^2 $$
Thus plugging these two expressions into the variation formula we see
$$ \delta I = \int_{U \times (0,T]} e^{- t /\epsilon} (  Du \cdot D \xi + \epsilon \cdot u_t \cdot \xi_t) dx dt $$
Recall the integration by parts formula for space:
$$ \int_U D u \cdot  D \xi dx = \int_{ \partial U} \underbrace{\xi Du \cdot n dS}_{\xi \equiv 0} - \int_U \xi \Delta u d x $$ 
where $n$ is the normal on the boundary. Now due to the product structure of the measure and Fubini's Theorem we have
$$ \int_{U\times (0,T] }e^{- t /\epsilon} (  Du \cdot D \xi) dx dt = \int_{(0,T]} e^{-t / \epsilon } \left ( \int_U Du \cdot D \xi dx \right ) dt =-\int_{(0,T]} e^{-t / \epsilon } \left ( \int_U \xi  \Delta u dx \right ) dt  = -\int_{U\times (0,T] }e^{-t / \epsilon }\xi  \Delta u dx  dt$$
Now recall integration by parts in one variable 
$$ \int _a^b f g'dt = f g \Big |^b_a - \int_a^b f' g dt $$
Here we have $f = \epsilon e^{- t / \epsilon } u_t $ and $g = \xi$. Thus we see a similar thing happens to the 2nd term.
$$ \int_{(0,T]} \epsilon e^{- t / \epsilon} u_t \xi _t dt = -\int_{(0,T]} \epsilon \xi \frac{d}{dt} (e^{- t / \epsilon} u_t) dt  = \int_{(0,T]} e^{-t / \epsilon} \xi (u_t - \epsilon u_{tt}) dt $$
Thus we see
$$ \delta I = \int_{U \times (0,T] } \xi e^{- t / \epsilon } ( - \Delta u +u_t - \epsilon u_{tt} ) dx dt $$
Now since we need $\delta I =0$ for all $\xi$, this implies that 
$$- \Delta u +u_t - \epsilon u_{tt} = 0$$
since the exponential is never zero.
Original Post: Well, my go to approach is to take the inner product against the solution. i.e. if $\mathcal{L}(u) =u_t - \Delta u - \epsilon u_{tt}$, then we look at
$$(\mathcal{L}(u),u) = \int_{U_t} \mathcal{L}(u) u dx dt = \int_{U_t} (u_t u - u\Delta u -  \epsilon u_{tt} u)dx dt $$
(space and time!).
If we integrate by parts in space on the middle guy and time on the last guy ( noticing the derivative of $t$ on the magnitude on the first guy) we see it should look something like:
$$ (u_t u - u\Delta u -  \epsilon u_{tt} u) \approx \frac{1}{2}\frac{d}{dt} |u|^2 + \frac{1}{2} |Du|^2 + \frac{1}{2} \epsilon | u_t|^2 $$
See what happens if you try 
$$ L(u) =  e^{-c(\epsilon)t}\left (\frac{1}{2}|Du|^2 + \frac{1}{2} \epsilon | u_t|^2 \right ) $$
The $e^{-c(\epsilon)t}$ will give us a product rule situation without affecting the derivative term.
Hint: When we perform a variation, we're solving
$$ \delta I = \lim_{h \to 0} \frac{ I( u + h \xi ) - I ( u) }{h}$$
for any suitable test function $\xi$(notably dies on boundary). This amounts to finding first order terms in $h$ from the perturbation. i.e.
$$L( u + h \xi ) = e^{- c( \epsilon) t } \left ( \frac{1}{2} | D ( u + h \xi )|^2 + \frac{1}{2} \epsilon | u_t + h \xi_t |^2 \right) = L(u) + h [e^{ -c(\epsilon)t} \left( D u \cdot D \xi + \epsilon u_t \cdot \xi_t\right) ] + \mathcal{O}(h^2)$$
Thus we see
$$ \delta I = \int_{U_t} e^{ -c(\epsilon)t} \left( D u \cdot D \xi + \epsilon u_t \cdot \xi_t \right) dx dt $$
Now perform the integration by parts I mentioned in the motivation and match $c( \epsilon)$ to the equation you want. You should end up with
$$\delta I = \int_{U_t} e^{-c(\epsilon)t} (- \Delta u + \epsilon c(\epsilon) u_t - \epsilon u_{tt} )\cdot \xi dx dt $$
Thus we see $u$ is critical if 
$$- \Delta u + \epsilon c(\epsilon) u_t - \epsilon u_{tt} =0$$
by the fundamental lemma of variational calculus.
A: I can't comment, so I am posting  a solution.
This is based in 5 variables, but the variable $t$ is, like more one variable in $x$, this is, we can think $x \in \mathbb{R}^{n+1}$.
So the Euler-Lagrange equation has the same form that when we work with only 3 variables. 
If you consider $L(p,q,z,x,t) = ( \frac{1}{2} |p|^2 + \frac{1}{2} \epsilon q^2) e^{-\frac{1}{\epsilon} t}$, I think you'll have the answer.
P.S: If I am wrong, say me please. I'm studying this chapter in this moment too.
