How to derive the Mott gap mathematically? From the one-band Hubbard model, $$H=-t\sum\limits_{<ij>, \sigma}c_{i\sigma}^{\dagger}c_{j\sigma}+U\sum\limits_{i}n_{i\uparrow}n_{i\downarrow},$$ we know if $U\gg t$, the energy cost of two electrons staying on the same site is too high, so there the band will be split into two bands with gap U. How to prove this mathematically from the Hubbard Hamiltonian?
 A: I think the easiest way to look at this problem would be using perturbation theory. Your hamiltonian has two pieces: the first piece is a hopping term, and the second piece is an on-site potential. You are considering the case where the on-site potential is much more important than the hopping.
Therefore, to zeroth order in perturbation theory, we ignore the hopping term, and the eigenstates of the hamiltonian are just the particle number eigenstates. The eigenvalues are just multiples of $U$, where the multiplier on $U$ is the number of doubly occupied sites.
Now we turn on hopping. In each energy level the perturbation is degenerate, and so each energy level will broaden into a band. This perturbation even mixes states in different energy levels. However, since the energy levels are so far separated, they do not mix much and the contribution to the energy is only of order $t^2/U \ll U$. Thus the energy is not perturbed enough to cause the bands to overlap and we have distinct bands. Since the broadening of the bands is of order $t$, the band spacing is still nearly $U$. 
The explanation of perturbation theory should be in a reference, like wikipedia. 
One problem with my answer is that I predicted the number of bands is equal to the number of sites instead of two. Perhaps I misinterpreted your question (or I might just be wrong).
A: I think the most insightful view of the band split occurring in the Hubbard model comes from the Hubbard III paper (you can find it here - https://royalsocietypublishing.org/doi/abs/10.1098/rspa.1964.0190). In it, he shows that the DOS to go from overlapping to split the increase of U using the equations of motion.
You might also find it useful to look into CP approximation of the Hubbard Model which views the model as a binary alloy for the opposite spin being either present or absent at a site. There are certain limitations to the model but it produces very good results. 
