There is a problem in this estimate. In Part 2, We can't derive the second inequality from the first inequality in Part 2. The main reason is the last term on the RHS in the first inequality is in $\|\cdot\|_1$ not $\|\cdot\|_0$.
How can I fix this proof or are there any stability result controlling $\tau \sum_{n=1}^k\|\partial_\tau u_h^n\|_0^2$ where $f\in L^2(0,T; H^{-1}(\Omega))$.
For your convenience, you can focus on Part 2. Also I have shown other parts in the proof for your reference.
Here is the problem description and the proof
Let's consider a simple parabolic equation (in variational form).
Assume $f\in L^2(0,T; H^{-1}(\Omega))$. Find $u\in H^1(\Omega )$ such that $u(\cdot,0)=u_0$ and $u(x,t)=g(x,t)$ on $\partial \Omega \times (0,T)$ and \begin{equation} (u_t,v)+a(u,v)=(f,v) ~~\forall v\in H^1_0(\Omega) \end{equation} for almost all $t\in (0,T)$ where $a(\cdot,\cdot)$ and $f(\cdot)$ are given by $$ a(w,v)=\int_{\Omega} a(x) \nabla w\cdot \nabla v ~dx ,~~~f(v)=\int_{\Omega} f(x)v(x)~dx $$ (Actually, I only concern about the case that $a(x)\equiv 1)$
Then the finite element approximation can be formulated as follows:
Find $u_h^n\in V_0^h $ such that $u^0_h=\prod_h u_0$ and for $n=0,1,...,M:$ \begin{equation}\label{1} (\partial_\tau u_h^n , v_h)+a(u^n_h,v_h)=(f^n,v_h)~~\forall v_h\in V_0^h \end{equation}
Stability estimates We establish some stability estimates for the finite element solution $\{u_h^n\}$.
Part 1: Let $v_h=\tau u_h^n$ in approximation formulation, we have \begin{equation}\label{dstability2} \frac{1}{2}\|u^n_h\|_0^2 -\frac{1}{2}\|u^{n-1}_h\|_0^2+\tau\|u_h^n\|_1^2\leq \tau \|f^n\|_{-1}\|u_h^n\|_1\leq \frac{1}{2}\tau\|f^n\|_{-1}^2+\frac{1}{2}\tau\|u_h^n\|_{1}^2 \end{equation}
Summing up over $n=1,2,...,k\leq M,$ \begin{equation} \|u_h^k\|_0^2+\sum_{n=1}^k\tau\|u_h^n\|_1^2\leq \tau \sum_{n=1}^k\|f^n\|_{-1}^2+\|u_h^0\|_0^2= \tau \sum_{n=1}^k\|f^n\|_{-1}^2+\|\prod\nolimits_h u_0\|_0^2 \end{equation}
Part 2:
Now taking $v_h=\tau \partial_\tau u_h^n$ in approximation formulation, we get \begin{equation} \tau \|\partial_\tau u_h^n\|_0^2+\frac{1}{2}\|u_h^n\|_1^2-\frac{1}{2}\|u_h^{n-1}\|_1^2\leq \tau\langle f^n,\partial_\tau u_h^n\rangle \leq \frac{1}{2}\tau\|f^n\|_{-1}^2+\frac{1}{2}\tau\|\partial_\tau u_h^n\|_1^2 \end{equation}
Summing up over $n=1,2,...,k\leq M$. we derive \begin{equation} \tau \sum_{n=1}^k\|\partial_\tau u_h^n\|_0^2+\|u_h^k\|_1^2\leq \|u_h^0\|_1^2+\tau\sum_{n=1}^k\|f^n\|_{-1}^2 \end{equation}
Part 3:
Next, we will estimate the dual norm of $\partial_\tau u_h^n$. For this, we introduce a $L^2-$projection $Q_h: L^2(\Omega)\to V_0^h$ by \begin{equation} (Q_h w, v_h)=(w,v_h) ~~~\forall v_h\in V_0^h. \end{equation}
Then we can show \begin{equation} \|w-Q_h w\|_0\leq Ch |w|_1, ~~~~~\|Q_h w\|_s\leq C\|w\|_s ~~~\text{ for } s=0,1. \end{equation}
Now for any $\phi \in H^1_0(\Omega)$, we derive \begin{equation} (\partial_\tau u_h^n, \phi)=(\partial u_h^n, \phi-Q_h \phi)- (\nabla u_h^n, \nabla Q_h \phi)+ \langle f^n, Q_h \phi\rangle. \end{equation}
\begin{equation} \begin{aligned} |(\partial_\tau u_h^n, \phi)|&=\|\partial u_h^n\|_0\| \phi-Q_h \phi\|_0+\|u_h^n\|_1\|\nabla Q_h \phi\|_1+ \|f^n\|_{-1}\| Q_h \phi \|_1\\ & \leq C(h\|\partial_\tau u_h^n\|_0+\|u_h^n\|_1+\|f^n\|_{-1})\|\phi\|_1. \end{aligned} \end{equation}
This implies that \begin{equation} \tau \sum_{n=1}^M\|\partial_\tau u_h^n\|^2_{H^{-1}(\Omega)} \leq C \tau \sum_{n=1}^M \{h^2\|\partial_\tau u_h^n\|_0^2+\|u_h^n\|_1^2+\|f^n\|_{-1}^2\}. \end{equation}