Short maths version at end
Consider the self-adjoint operator (the occupation number operator that forms the Hamiltonian of the quantum harmonic oscillator) $$N = a^\dagger a$$ Where $a$ and $a^\dagger$ are the "creation" and" annihilation" operators satisfying $[a, a^\dagger ]=1$. It can be shown that the eigenvalues of this operator are $n \in \mathbb{N} $ and the associated eigenstates are denoted $|n\rangle$.
Now in physics we appeal to the spectral theorem, so that the eigenstates of a self-adjoint bounded operator are complete in the Hilbert space so that, for example, the identity admits a decomposition $$I = \sum_{n=0}^\infty |n\rangle \langle n|, $$ or an arbitrary state can be written as a linear combination of the eigenstates $$|a\rangle = \sum c_n |n\rangle $$ for coefficients $c_n$.
However as far as I can tell, the number operator is not bounded - there is no real number $K$ such that $$| \, N|n\rangle \, | < K|\, |n\rangle\, | ~~ \forall n$$ where the vertical lines indicate the norm on the Hilbert space. This is because the eigenvalues of the operator have no upper bound and the states are normalised to unity.
So why does the spectral theorem apply? Is there a generalisation to unbounded operators with nice (which?) properties? Or am I mistaken and the operator is bounded after all?
Mathematicians' version:
Suppose there is a linear self-adjoint operator $N: H \rightarrow H$ mapping a vector space over $\mathbb C$ to itself with eigenvectors $v_n$ satisfying $N v_n = n v_n$ for n a natural number. It seems N is unbounded yet a version of the spectral theorem holds so that the eigenvectors are complete.
The number operator is indeed unbounded; but there is a spectral theorem for general unbounded self-adjoint operators from the 1930's going back to Stone and von Neumann. The core statement reads as follows (compare for example Thm.10.4 in the book "Quantum Theory for Mathematicians" (2013) by B. Hall):
Edit: The only difference between the spectral theorem for bounded and the one for unbounded self-adjoint operators is the sense in which the above integral has to be understood. While (1) for bounded $A$ holds in the usual norm sense, for unbounded $A$ one has to consider a sequence of bounded Borel functions $h_n(x)$ with $h_n(x)\to x$ as $n\to\infty$ for each $x$ and $|h_n(x)|\leq|x|$ for all $x$ and $n$. Then for all $\psi$ in the domain of $A$, eq.(1) changes into $$ \lim_{n\to\infty}\int_{\sigma(A)} h_n(\lambda)\,d\mu^A(\lambda)\psi=A\psi\,.\tag{2} $$ For more details on this topic, refer to "Methods of Modern Mathematical Physics. I: Functional Analysis" (1980) by Reed & Simon (Theorem VII.2 & VIII.5)
Now the number operator is a special case because it is a discrete operator, i.e. the spectrum of $N$ consists only of eigenvalues which do not accumulate anywhere. Thus its spectral measure---loosely speaking---corresponds to an orthonormal basis of $\mathcal H$, and not "just" a resolution of the identity in terms of an integral $I=\int 1\,d\mu^A$, cf. also here.