set $\gamma_0=\gamma_N=0$.

We emphasize that the constants $\alpha_j$ and ${\epsilon}^2\alpha_j\gamma_j$ appearing in the error indicators $\eta_{K_j}$ from  remain bounded as $h_j,{\epsilon}\to 0$ (and $p_j\to\infty$). In addition, we note that $\nicefrac{1}{2} \min\{\beta_j,\beta_{j+1}\} \leq \gamma_j \leq \min\{\beta_j,\beta_{j+1}\}$, and that $2\sqrt{{\epsilon}^{-1} \alpha_j } \leq \beta_j \leq 3 \sqrt{{\epsilon}^{-1} \alpha_j}$.

------------------------------------------------------------------------

For *linear problems*, the (adaptive) Newton linearization is redundant. In particular, upon setting $\Delta t=1$, we infer that $\delta_{n,\Omega}=0$ in , and the residual is bounded by the local error indicators $\eta_{n,j}$ only. These quantities, in turn, can be exploited for the purpose of designing an $hp$-adaptive mesh refinement strategy, as outlined, for example, in Section \[sc:hpadapt\] below. We emphasize that this approach has been studied earlier in the report [@MelenkWihler:15], where a number of numerical tests in the context of $hp$-adaptive FEM for singularly perturbed linear problems have been performed; in particular, these experiments have underlined the robustness of the bound  as ${\epsilon}\to 0$ in the linear case. In Section \[sc:adaptivity\], we will extend this idea by proposing a suitable interplay between $hp$-mesh refinements and the (adaptive) Newton linearization from Algorithm \[al:zs\].

------------------------------------------------------------------------

Given $u_{n+1}^{(\Delta t,{{\rm hp}})}$ and $f^{\Delta t}({u^{\rm hp}}_{n+1})$ from  and , respectively, and any $v\in H^1_0(\Omega)$, there holds the 