The strategy is identical to the [elliptic interior regularity](/theorems/95) argument — test with difference quotients — but applied at each time slice. The proof plan is: use spatial difference quotients $D_k^h u$ as building blocks for a [test function](/page/Test%20Function), derive uniform $L^2$ bounds on $D_k^h \nabla u$ independent of $h$, and invoke the difference quotient characterisation to conclude $u \in H^2$ locally in space. The time [derivative](/page/Derivative) is then recovered from the equation itself. The parameter $\delta > 0$ absorbs any lack of regularity at $t = 0$. **Step 1: Localisation.** Fix $V \subset\subset W \subset\subset U$ and $0 < \delta/2 < \delta < T$. Let $\zeta \in C_c^\infty(U)$ be a spatial cutoff with $\zeta \equiv 1$ on $V$, $\operatorname{supp}(\zeta) \subset W$, and $|\nabla \zeta| \le C/\operatorname{dist}(V, \partial W)$. Let $\eta \in C^\infty([0, T])$ be a temporal cutoff with $\eta \equiv 0$ on $[0, \delta/2]$, $\eta \equiv 1$ on $[\delta, T]$, and $0 \le \eta \le 1$. **Step 2: Difference quotient test function.** For $|h| < \operatorname{dist}(W, \partial U)$ and direction $e_k$, define the test function \begin{align*} v := -D_k^{-h}\bigl(\eta^2 \zeta^2 D_k^h u\bigr). \end{align*} This belongs to $L^2(0, T; H^1_0(U))$ because: $D_k^h u \in L^2(0,T; H^1_0)$ (translation preserves Sobolev regularity), $\eta^2 \zeta^2 D_k^h u$ has compact spatial support in $W$, and $D_k^{-h}$ translates the support by at most $|h|$ while remaining inside $U$. **Step 3: Energy identity.** Insert $v$ into the weak formulation and use the discrete [integration by parts](/theorems/210) identity $\int D_k^h f \cdot g = -\int f \cdot D_k^{-h} g$ to move $D_k^{-h}$ onto the other factors. After expanding and using the product rule for difference quotients, the principal term is \begin{align*} \int_0^T \eta^2 \int_U \zeta^2 \sum_{i,j} a_{ij} \partial_{x_j}(D_k^h u) \partial_{x_i}(D_k^h u)\, d\mathcal{L}^n\, dt, \end{align*} which, by uniform ellipticity, is bounded below by $\theta \int_0^T \eta^2 \|\zeta \nabla(D_k^h u)\|_{L^2}^2\, dt$. **Step 4: Estimation.** The remaining terms — commutator terms from $D_k^h(a_{ij})$, cross terms involving $\nabla \zeta$, the time derivative contribution, and the right-hand side — are bounded using Cauchy-Schwarz, the $L^\infty$ bounds on the coefficients and their difference quotients (using smoothness of $a_{ij}$), and Young's inequality. After absorbing the $\|\zeta \nabla(D_k^h u)\|_{L^2}^2$ term from the right side into the left, one obtains \begin{align*} \int_\delta^T \|D_k^h \nabla u\|_{L^2(V)}^2\, dt \le C\bigl(\|f\|_{L^2(U_T)}^2 + \|u\|_{L^2(U_T)}^2\bigr), \end{align*} with $C$ independent of $h$. **Step 5: Passage to $H^2$.** The uniform bound on $\|D_k^h \nabla u\|_{L^2(V)}$ (for all small $h$ and all directions $k = 1, \ldots, n$) implies, by the reverse direction of the difference quotient characterisation (for $p = 2 > 1$), that $\nabla u \in L^2(\delta, T; H^1(V))$, i.e., $u \in L^2(\delta, T; H^2(V))$. **Step 6: Recovery of $u_t$.** Since $u_t = f - Lu$ and $u \in L^2(\delta, T; H^2(V))$ with $f \in L^2(U_T)$, we get $u_t \in L^2(\delta, T; L^2(V))$. The quantitative bound follows by combining the difference quotient estimate with the equation.