## 2 Fractional-in-Space Parabolic Equation
### Adaptive Forward-Backward Euler Discretization in Time  
The local refinement discussed above implies a more specific regularity of the solution in time that is a priori known. We do not need any such assumptions in this paper. The developed adaptive time stepping algorithm estimates the current error and automatically adjusts the next step size based on it. The procedure is as follows. The approximate solution is calculated using the stable backward Euler scheme (3), while a forward Euler scheme is used to estimate the size of the next time step. This type of methods are widely used in standard (local) diffusion problems. For example, adaptive time stepping for coupled multiscale problems is studied in [23; 24], where the solution dynamics arise from nonlinear transitions of material properties. A parabolic equation with a jumping right-hand side is considered in [25], which is closest to the test problem we used in our numerical experiments. The time step adaptation algorithm explored in this paper is largely in the spirit of [25]. The great challenge and novelty of the results presented here are due to the non-locality of fractional diffusion. Let us consider the semi-discrete Cauchy problem  
\[\frac{\partial\mathbf{U}}{\partial t}+K\mathbf{U}=\mathbf{F} \tag{5}\]  
corresponding to the FEM discretization in space of (1), where the vectors \(\mathbf{U}\) and \(\mathbf{F}\) represent the nodal functions of the unknown solution \(u(x,t)\) and the right-hand side \(f(x,t)\) on \(\mathcal{T}\). The error \(\mathbf{z}^{i+1}\) of the backward Euler scheme for (5) can be estimated by  
\[||\mathbf{z}^{i+1}||\leq\sum_{i=0}^{j}\tau_{i+1}||\Psi_{i+1}||, \tag{6}\]  
where the truncation error \(\Psi_{i+1}\) has the form  
\[\Psi_{i+1}=\mathbf{F}^{i+1}-\frac{\mathbf{U}^{i+1}-\mathbf{U}^{i}}{\tau_{i+1} }-K\mathbf{U}^{i+1}.\]  
Following the notations introduced in the fully discrete case, the superscript \(j\) means that the values are taken in the node \(t_{j}\in\omega\). Now, let \(\delta>0\) be the parameter defining the desired accuracy. We are looking for such a time step \(\tau_{i+1}\) that the inequality  
\[\Psi_{i+1}\leq\tau_{i+1}\delta \tag{7}\]holds. If (7) is satisfied, from (6), we conclude that the error accumulates linearly and is majorized by  
\[||\mathbf{z}^{i+1}||\leq\delta t_{i+1}.\]  
Unfortunately, it is not possible to apply (6) and (7) directly. This is because, to calculate the truncation error, we need the exact solution \(\mathbf{U}=\mathbf{U}(t)\), which is of course unknown. In the algorithm for the time step selection presented below, we use an approximate truncation error. For this purpose, an auxiliary forward Euler step is performed. Thus, the numerical solution of the parabolic problem is calculated by applying as the main discretization tool the unconditionally stable forward Euler scheme. Importantly, the additionally included inverse-free forward scheme is significantly cheaper from a computational point of view. It should be noted that this scheme is unstable, which is not a problem, however, because we only apply it locally to estimate \(\tau_{i+1}\). Algorithm for Time Step Selection  
The adaptive size of \(\tau_{i+1}\) is computed in the following three steps:  
1. Prognostic solution. Calculation of \(\tilde{n}^{i+1}\) in \(\tilde{I}_{i+1}=t_{i}+\tilde{\tau}_{i+1}\) using the forward Euler scheme (8). 2. Prognostic error. The prognostic solution \(\tilde{n}^{i+1}\) is used to compute the prognostic truncation error \(\tilde{\Psi}_{i+1}\) given by (9). 3. Step size selection. The time step \(\tau_{i+1}\) is selected based on the prognostic error and the parameter \(\delta\) that defines the desired accuracy.