I am using a space-time discretization to discretize an initial value problem, using P1-FEM for space discretization and BDF-2 for time discretization.
To fix the notations: $\Omega\subsetneq\mathbb{R}^3$ is a domain, $\{\mathcal{T_h}\}_{h\ge0}$ a family of regular tetrahedral triangulations of $\Omega$, $\mathcal{S}^1(\mathcal{T}_h)^3$ the space of globally continuous piecewise polynomials of degree $1$. Following the standard notations, $\{t_j\}_{j=0}^N$ is a uniform partition of the time interval $[0,T]$ with time step $k$ and, given a sequence of functions $\{u^j\}_{j=0}^N$, the backward difference operator is clearly defined as $d_tu^j:=\frac{1}{k}(u^j-u^{j-1})$.
Let us call $m$ is the solution to my problem, $m^0$ the initial condition and let $v:=\partial_t m$. It holds that $|m(t,x)|=1\,\, \forall (t,x)\in[0,T]\times\mathbb{R}^3$ and therefore $m\cdot\partial_t m=0$.
The algorithm used for time discretization is a 2-steps BDF method, therefore it produces the sequence of discrete functions approximating the solution $\{m_h^{j}\}_{j=1}^{N}\subset H^1(\Omega)$ and $\{v_h^j\}_{j=0}^{N-1}\subset H^1(\Omega)$ approximating the derivative, where
$ \bullet\quad m_h^1 $ is obtained with Backward Euler method, i.e., $m_h^1=m_h^0+kd_tm_h^1$;
$ \bullet\quad v_h^0:= d_tm_h^1$ and $v_h^j=\dot{m}_h^{j+1}\,\,\forall \,1\le j\le N-1$ is the discrete time derivative provided by BDF 2, i.e. $$\dot{m}_h^j=\frac{1}{k}\left(\frac{3}{2}m_h^{j}-2m_h^{j-1}+\frac{1}{2}m_h^{j-2}\right);$$
$ \bullet\quad m_h^j,\,\forall\, 2\le j\le N $ is obtained with BDF2, i.e., $$m_h^j=\frac{4}{3}m_h^{j-1}-\frac{1}{3}m_h^{j-2}+\frac{2}{3}k\dot{m}_h^j.$$
I am interested in the weak convergence as $h,k\to0$ of the piecewise linear and piecewise constant approximations defined as follows: $\forall\, 0\le j \le N-1$ and $t\in[t_j,t_{j+1})$ $$m_{hk}(t)=\frac{t-t_j}{k}m_h^{j+1}+\frac{t_{j+1}-t}{k}m_h^j,\\v_{hk}^-(t)=v_h^j.$$
I managed to bound $||m_{hk}||_{H^1(\Omega_T)}$ and $||v_{hk}^-||_{L^2(\Omega_T)}$ which allow to extract weakly convergent subsequences of $\{m_{hk}\}$ and $\{v_{hk}^-\}$.
I have problems in showing the convergence of derivatives, i.e., I would like to show that $$||\partial_tm_{hk}-v_{hk}^-||_{L^1(\Omega_T)}\to0\text{ when }h,k\to0.$$
Here is my try: $$||\partial_tm_{hk}-v_{hk}^-||_{L^1(\Omega_T)}=\int_{0}^T||\partial_tm_{hk}(t)-v_{hk}^-(t)||_{L^1(\Omega)}dt=\sum_{j=0}^{N-1}\int_{t_j}^{t_{j+1}}||\partial_tm_{hk}(t)-v_{hk}^-(t)||_{L^1(\Omega)}dt=\\ =\sum_{j=1}^{N-1}\int_{t_j}^{t_{j+1}}||\frac{m_h^{j+1}-m_h^j}{k}-v_h^{j}||_{L^1(\Omega)}dt=\sum_{j=1}^{N-1}\int_{t_j}^{t_{j+1}}||\frac{m_h^{j+1}-m_h^j}{k}-\dot{m}_h^{j+1}||_{L^1(\Omega)}dt=\\=k\sum_{j=1}^{N-1}||\frac{m_h^{j+1}-m_h^j}{k}-\dot{m}_h^{j+1}||_{L^1(\Omega)}=\sum_{j=1}^{N-1}||m_h^{j+1}-m_h^j-k\dot{m}_h^{j+1}||_{L^1(\Omega)},$$ where the sum goes from 1 because $v_h^0=d_tm_h^1$ and therefore it would give a null contribution and clearly $\partial_tm_{hk}(t)=d_tm_h^j$ for all $0\le j\le N-1$.
Now, I don't know how to proceed because I struggle in trying to bound $||m_h^{j+1}-m_h^j-k\dot{m}_h^{j+1}||_{L^1(\Omega)}$ in a way that goes to $0$ with $h$ and $k$. I have the feeling that I should bound it with a squared norm, otherwise I lose one order of convergence.
Is there any property of standard proof of BDF methods which guarantees that the 2-steps derivative $\dot{m}_h^j$ converges to the backward difference operator $d_tm_h^{j}$?
Hints that may help:
$\bullet$ I have a nodal orthogonality property, i.e., $\dot{m}_h^j(z)\cdot\hat{m}_h^j(z)$ for all $j\ge 2$ and for every node $z$ of the triangulation, where $\hat{m}_h^j:=2m_h^{j-1}-m_h^{j-2}.$
$\bullet$ All the results from https://arxiv.org/abs/2310.00381 hold also in my case and I have the feeling that can be used to bound properly the quantity I am interested in.
Thanks!