I doubt any answer will be satisfactory. My opinion is that we are still very far from a mathematical justification. If we accept the mathematical foundations of quantum mechanics, and if we make the approximation that the nucleus of the atom is just one heavy thing with $N$ positive charges, then the motion of the $N$ electrons is governed by a linear equation (Schrödinger) in ${\mathbb R}^{3N}$. The unknown is a function $\psi(r^1,\ldots,r^N,t)$ with the property (Pauli exclusion) that it has full skew-symmetry. For instance,
$$\psi(r^2,r^1,\ldots,r^N,t)=-\psi(r^1,r^2,\ldots,r^N,t).$$
In practice, we look for steady states $e^{i\omega t}\phi(r^1,r^2,\ldots,r^N)$. Then $\omega$ is the energy level.
Because of the very large space dimension, one cannot perform reliable calculations on computer, when $N$ is larger than a few units. One attempt to simplify the problem has been to postulate that $\phi$ is a Slatter determinant, which means that
$$\phi(r^1,r^2,\ldots,r^N)=\|a_i(r^j)\|_{1\le i,j\le N}.$$
The unknown is then an $N$-tuple of functions $a_i$ over ${\mathbb R}^3$. Of course, we do not expect that steady states be really Slater determinants; after all, the Schrödinger equation does not preserve the class of Slater determinants. Thus there is a price to pay, which is to replace the Schrödinger equation by an other one, obtained by an averaging process (Hartree--Fock model). The drawback is that the new equation is non-linear. Such approximate states have been studied by P.-L. Lions & I. Catto in the 90's.
Update. Suppose $N=2$ only. If we think to $\phi$ as a finite-dimensional object instead of an $L^2$-function, then it is nothing but a skew-symmetric matrix $A$. Approximation à la Slater consists in writing $A\sim XY^T-YX^T$, where $X$ and $Y$ are vectors. In other words, one approximate $A$ by a rank-two skew-symmetric matrix. The approximation must be in terms of the Hilbert-Schmidt norm (also named Frobenius, Schur): this norm is natural because of the requirement $\|\phi\|_{L^2}=N$. If $\pm a_1,\ldots,\pm a_m$ are the pairs of eigenvalues of $A$, with $0\le a_1\le\ldots\le a_m$, then the best Slater approximation $B$ satisfies $\|B\|^2=2a_m^2$, $\|A-B\|^2=2(a_1^2+\cdots+a_{m-1}^2)$. Not that good. Imagine how much worse it can be if $N$ is larger than $2$.
You are completely correct in your analysis of the structure one obtains by considering the
Schrodinger equation for Z electrons in a central potential due to a nucleus of charge Ze
when the Coulomb interaction between electrons as well as relativistic effects such as spin-orbit coupling are ignored. In such a model the periods would indeed be 2, 8, 18 etc. for
exactly the reasons you have described. In physics jargon the energy in this model depends only on the principal quantum number $n \in {\mathbb Z}, n>0$ and the allowed $\ell$ values
are $\ell\le n-1$. Orbitals with $\ell=0,1,2,3,4, \cdots$ are labelled by $s,p,d,f \cdots$ for
historical reasons. Thus at $n=1$ one can have one or two states in the $(1s)$ orbital (accounting for spin), at $n=2$ one has the $(2s)$ and $(2p)$ orbitals with $2(1+3)=8$
states. And at $n=3$ one should have $18$ states by filling the $(3s),(3p),(3d)$ orbitals.
But in the real world this is not what happens. This simple model does not give a correct description of atoms in the real world once you get past Argon. I believe the main effect leading to the breakdown is the Coulomb interaction between the electrons.
So there is no simple mathematical model based say just on the representation theory of $SU(2)$ and simple solutions to the Schrodinger equation which will account for the structure of the periodic table past Argon. However one could ask whether including Coulomb interactions between electrons does gives a model which correctly reproduces the next few rows of the periodic table past Argon. I am not an expert on this, but since I doubt there are physical chemists on MO I'll just give my rough sense of things.
To approach this problem with some level of rigor probably requires difficult numerical work
and my impression is that this is beyond the current state of the art. However there are
rough models which try to approximate what is going on by assuming that the interactions between electrons can be replaced by a spherically symmetric potential which is no longer
of the $1/r$ form. This leaves the shell structure as is, but can change the ordering of
which shells are filled first. In such a model instead of filling the $(3d)$ shell after
Argon one starts to fill the $(4s)$ and $(3d)$ shells in a somewhat complicated order. Eventually one fills the $(4s),(3d),(4p)$ shells and this leads to the line of the periodic
table starting at K and ending at Krypton.
Added note: There is one nice piece of mathematics associated with this problem that I should have mentioned, even if it doesn't by itself explain the detailed structure of the periodic table. When Coulomb interactions between electrons and relativistic effects are ignored the energy levels of the Schrodinger equation with a central $1/r$ potential depend only the quantum number $n$, but not on the quantum number $\ell$ which determines the representation $V_\ell$ of $SO(3)$ referred to above. When screening is included this is no longer the case and the energies depend on both $n$ and $\ell$. Why is this?
With a $1/r$ central potential there is an additional vector $\vec D$ which commutes with the Hamiltonian. Classically this vector is the Runge-Lenz vector and its conservation explains why the perihelion of elliptical orbits in a $1/r$ potential do not precess. Quantum mechanically the
commutation relations of the operators $\vec D$ along with the angular momentum operators $\vec L$ are those of the Lie algebra of $SO(4)$ (for bound states with negative energy). There are two Casimir invariants, one vanishes and the other is proportional to the energy.
As a result the energy spectrum depends only on $n$ and can be computed using group theory without ever solving the Schrodinger equation explicitly. Perturbations due to screening, that is from some averaged effect of the Coulomb interactions between electrons,
change the $1/r$ potential to some more general function of $r$ and break the symmetry
generated by the Runge-Lenz vector $\vec D$. As a result the energy levels depend on both $n$ and $\ell$.
Best Answer
I think you can find more in Lieb and Seiringer's book "The Stability of Matter in Quantum Mechanics", or see also Freeman Dyson http://www.webofstories.com/play/4415 and the book review http://arxiv.org/abs/1111.0170.