main source of error), a fully adaptive Newton-Galerkin scheme is obtained; see also the related works on monotone nonlinear problems [@ChaillouSuri:07; @El-AlaouiErnVohralik:11; @CongreveWihler:15], or the articles [@ChaillouSuri:06; @Han:94] dealing with modelling errors in linearized models. Our numerical results in [@AmreinWihler:15] reveal that sensible solutions can be found even in the singularly perturbed case, and that our scheme is reliable in the sense that the adaptive Newton-Galerkin method is able to converge to solutions which preserve the qualitative structure of (reasonably chosen) initial guesses (see also [@AmreinWihler:14; @ScWi11]).

*$hp$-FEM and $hp$-adaptivity:* The aim of the present paper is to extend our work in [@AmreinWihler:15] from a low-order to an $hp$-adaptive Newton-Galerkin approach. For simplicity of exposition, we focus on the 1d case, however, we mention that an extension to the multi-d case is possible as well. Exploiting the $hp$-framework may potentially lead to *exponential rates of convergence* in the numerical approximation. For the purpose of $hp$-adaptivity, we carry out an ${\epsilon}$-robust $hp$-version *a posteriori* residual analysis for the FEM discretizations of the Newton linearizations. We remark that using $hp$-adaptivity, in contrast to the adaptive low-order methodology developed in the previous works [@El-AlaouiErnVohralik:11; @AmreinWihler:14; @AmreinWihler:15], introduces an additional aspect of adaptivity that 