$\renewcommand{\Re}{\operatorname{Re}}$ $\newcommand{\erf}{\operatorname{erf}}$ $\newcommand{\dag}{\dagger}$ $\newcommand{\const}{\mathrm{const}}$
Consider Example 3 of Appendix M: Schrödinger operator in 2D with a constant magnetic and no electric field \begin{equation} L=\frac{1}{2} (-i\partial_x -\frac{1}{2}B y)^2 + \frac{1}{2} (-i\partial_y +\frac{1}{2}B y)^2 \label{equ-N.1} \end{equation} with $B>0$ (or $B<0$) has a pure point spectrum. For simplicity assume that $B>0$. We apply a gauge transformation which for Schrödinger operator means multiplying it from the left and right by $e^{i\hbar^{-1} \phi (x)}$ and $e^{-i\hbar^{-1} \phi (x)}$ respectively with a real-valued $\phi$ (which is an unitary transformation) and replaces $-i\hbar\nabla$ by $-i\hbar\nabla - (\nabla\phi)$ (which is equivalent to changing vector potential $\mathbf{A}(x)$ by $\nabla \phi$ which in turn does not change $\nabla \times \mathbf{A}$. Taking $\hbar=1$ and $\phi= \frac{1}{2}B xy$ we arrive to \begin{equation*} L'=\frac{1}{2} (-i\partial_x -B y)^2 + \frac{1}{2} (-i\partial_y )^2; \end{equation*} then making Fourier transform by $x\mapsto \xi$ we get \begin{equation*} L''=\frac{1}{2} (-\xi -B y)^2 + \frac{1}{2} (-i\partial_y )^2; \end{equation*} and plugging $y=B^{-\frac{1}{2}} (y_{\textsf{new}} -B^{-1}\xi)$ we get \begin{equation*} \frac{1}{2} B (-\partial_y^2 +y^2) \end{equation*} which is a harmonic oscillator multiplied by $B$ and in virtue of Appendix K its spectrum consists of eigenvalues $E_n = |B|(n+\frac{1}{2})$, $n=0,1,2,\ldots$ which are called Landau levels.
However there is a "hidden variable" $\xi$, so eigenfunctions Hermite functions of $y$ but multiplied by arbitrary functions $C(\xi)$ rather than by constants which implies that these eigenvalues have constant multiplicities.
Consider Example 8 of Appendix M: Schrödinger operator with periodic potential in $\mathbb{R}^d$: $V(x+a)=V(x)$ for all $a\in \Gamma$ where $\Gamma$ is a lattice of periods.
Let us decompose function $u(\mathbf{x})$ into $n$-dimensional Fourier integral \begin{equation*} u(x)= \iiint _{\mathbb{R}^n} e^{i \mathbf{k}\cdot \mathbf{x}} \hat{u}(\mathbf{k})\,d^n\mathbf{k}, \end{equation*} then replace this integral by a sum of integrals over dual elementary cell $\Omega^*$ shifted by $\mathbf{n}\in \Gamma^*$ \begin{equation*} \sum_{ \mathbf{n}\in \Gamma^* }\iiint _{\Omega^* +\mathbf{n}} e^{i \mathbf{k}\cdot \mathbf{x}} \hat{u}(\mathbf{k})\,d^n\mathbf{k}, \end{equation*} then change variable $\mathbf{k}=\mathbf{k}_{\textsf{new}}+\mathbf{n}$ \begin{equation*} \iiint _{\Omega^* } e^{i (\mathbf{k}+\mathbf{n})\cdot \mathbf{x}} \Bigl( \underbrace{\sum_{ \mathbf{n}\in \Gamma^* } e^{i \mathbf{n}\cdot \mathbf{x}} \hat{u}(\mathbf{k}+\mathbf{n})}_{=U(\mathbf{k},\mathbf{x})}\Bigr)d^n\mathbf{k}, \end{equation*} we observe that $U(\mathbf{k},\mathbf{x})$ is quasiperiodic with quasimomentum $\mathbf{k}$.
In advanced Real Analysis it would be a decomposition of our Hilbert space $\mathsf{H}=L^2(\mathbb{R}^n)$ into direct integral of Hilbert spaces $\mathsf{H}(\mathbf{k})$ of such functions, and our operator is acting in each of those spaces separately, with a spectrum $\sigma (L(\mathbf{k}))=\{E_n (\mathbf{k}): n=1,2,\ldots\}$. This implies that $L$ has a band spectrum: it consists of spectral bands $\sigma_n:=[\min _{k\in \Omega^*} E_n(k) ,\max _{k\in \Omega^*} E_n(k)]$: \begin{equation*} \sigma(L) =\bigcup_{n=1}^\infty \sigma_n; \end{equation*} these spectral bands can overlap. On can prove that $E_n (\mathbf{k})$ really depend on $\mathbf{k}$ and are not taking the same value on some set of non–zero measure (another notion from Real Analysis) which implies that the spectrum $\sigma(L)$ is continuos.