When solving the Hydrogen Atom using the Schrödinger's equation, we find the wave function represented in the position bases, and it turns out that the wave function has parameters $n$, $\ell$, $m$, and $s$.
What I don't understand is: Why do we know that the energy eigenbasis is $\lvert n\ell ms \rangle$? I mean, is it because the wave function in the position representation has parameters $n$, $\ell$, $m$, and $s$?
The energy actually does not depend on $\ell,m_\ell$ and in fact on any of the quantum numbers except for $n$.
Thus, a linear combination of states with the same $n$ but different $\ell,m_\ell$ compatible with that $n$: \begin{align} \vert\psi_n\rangle = \sum_{\ell,m_\ell} c_{\ell,m_\ell} \vert{n \ell m_\ell}\rangle \tag{1} \end{align} would also be a state with energy $E_n$. The state (1) is not necessarily a state of "good" $\ell$ or $m_\ell$. However, it is known that commuting operators have common eigenvectors, so the $\vert n\ell m_\ell\rangle$ kets are in fact common eigenvectors of $\hat H$, $\hat L^2$ and $\hat L_z$. The common eigenfunctions are \begin{align} \langle r\theta\varphi\vert n\ell m_\ell\rangle = R_{n\ell}(r) Y_{\ell}^m(\theta,\varphi) \end{align} and are obtained from the usual separation of variables in the Schrödinger equation.
If we find the quantum numbers $n\ell m_\ell$ are not sufficient to completely label states, we have to find one or more additional operators that will commute with the $3$ above. $\hat S^2$ and $\hat S_z$ are such operators and, as far as we know, the set $n,\ell,m_\ell,S,m_s$ is enough to completely label the states. The extension of (1) to include spin is would be \begin{align} \vert\psi_n\rangle = \sum_{\ell,m_\ell,m_s} c_{\ell,m_\ell,m_s} \vert{n \ell m_\ell;S=\textstyle\frac{1}{2},m_s}\rangle\, , \end{align} although the resulting state $\vert \psi_n\rangle$ would not necessarily be an eigenstate of $\hat L^2$, $\hat L_z$, or $\hat S_z$.
Thus far, the set of $n\ell m_\ell S m_s$ has proven enough to uniquely label all the possible states so there's no need to drum up additional labels.
When one solves the Schrodinger equation for the hydrogen atom, we solve
$$\hat{H}\left|\psi\right>=E\left|\psi\right>$$
The above equation defined $\left|\psi\right>$ to be an eigenstate. We can separate the equation into $\left|\psi\right>=R(r)Y(\phi,\theta)$ and then solve the radial part $R(r)$ and angular part $Y(\phi,\theta)$ seperately.
When we solve the angular part of this equation, we find that it reduces to a Legendre equation that has the variables $m$ and $\ell$ in it. This equation is can be solved and normalized if and only if $m$ and $\ell$ are defined as: $0\leq\ell\leq n-1$, and $-\ell\leq m \leq \ell$.
The value $n$ comes from solving the radial part of this equation (and it shares the same $\ell$ as the angular part, which gives us its relation).
Unfortunately, my class did not solve this with spin states, but presumably those again are conditions of solving the above equation.
If you've ever solved the particle in a box problem, you'll also see that the energy there is quantized and parametrized by a single integer value. In this case the problem is more complciated and as it turns out a set of $4$ integers is enough to describe any state.
In short, we know $\left|n\ell m s\right>$ is the energy eigenbasis because we solved the eigenvalue equation.
