Very interesting question, and it kept me up despite daylight saving time cheating me of one hour of sleep last night... A good reference is Albright, Burdett and Whangbo, Orbital Interactions in Chemistry 2nd ed. pp 282ff. which explains this in much greater detail than I can. (In general, that book is fantastic.) I will try my best to summarise what they have written - please point out any mistakes that I may have made. I think there is also much more information out there in the literature (a search for "Jahn–Teller cyclobutadiene" should give lots of good results), although it is of course not always easy to read. $\require{begingroup}\begingroup\newcommand{\ket}[1]{|#1\rangle}$
1. What is the ground state of undistorted $\ce{C4H4}$?
Based on simple Huckel theory as you have shown above, one would expect the ground state to be a triplet state with the two highest electrons in different orbitals. This is simply because of exchange energy (Hund's first rule).
However, it turns out that the singlet state with the two electrons in different orbitals is more stable. That explains the lack of signal in ESR experiments.
(Note that this orbital picture is a simplification, as I will describe later.)
In short, this is due to dynamic spin polarisation; the idea is that there are low-lying excited singlet states which mix with the ground state to stabilise it. This is the basic idea underlying so-called configuration interaction, which mixes in excited molecular states into the wavefunction in order to obtain a better estimate of the ground-state wavefunction. If you are interested in the theoretical aspect of this, I think Szabo and Ostlund's Modern Quantum Chemistry should cover it. The topic of dynamic spin polarisation is treated in depth in pp 285ff. of Albright et al.
Your question, "how to predict the singlet state a priori"? is a good one. Albright et al. write that (p 289)
In general, high-level configuration interaction calculations are needed to determine the relative stabilities of the singlet and triplet states of a diradical.
So, it does not seem to be something that you can figure out a priori, and it is most certainly not something you can figure out by just using Huckel theory (which is an incredibly simplistic treatment).
2. What is the degeneracy of the ground state?
You also wrote that "the singlet state will be doubly degenerate". However, in the case where both electrons are in different orbitals, there is only one possible spatial wavefunction, due to the quantum mechanical requirement of indistinguishability. (As we will see in the group theoretical treatment, there is actually no doubly degenerate case, not even when the electrons are paired in the same orbital.)
If I label the two SOMOs above as $\ket{a}$ and $\ket{b}$ then the combination $\ket{ab}$ is not a valid spatial wavefunction. Neither is $\ket{ba}$ a valid wavefunction.
You have to take linear combinations of these two states, such that your spatial wavefunction is symmetric with respect to interchange of the two electrons. In this case the appropriate combination is
$$2^{-1/2}(\ket{ab} + \ket{ba}) \qquad \mathrm{{}^1B_{2g}}$$
There are two more possible symmetric wavefunctions:
$$2^{-1/2}(\ket{aa} + \ket{bb}) \qquad \mathrm{{}^1A_{1g}}$$
$$2^{-1/2}(\ket{aa} - \ket{bb}) \qquad \mathrm{{}^1B_{1g}}$$
(Simply based on the requirement for indistinguishability, one might expect $\ket{aa}$ and $\ket{bb}$ to be admissible states. However, this is not true. If anybody is interested to know more, please feel free to ask a new question.) The orbital picture is actually woefully inadequate at dealing with quantum mechanical descriptions of bonding. However, Voter and Goddard tried to make it work out, by adding and subtracting electron configurations:1
Ignore the energy ordering in this diagram. It is based on Hartree–Fock calculations, which do not include the effects of excited state mixing (CI is a "post-Hartree–Fock method"), and therefore the triplet state is predicted to be the lowest in energy. Instead, just note the form of the three singlet states; they are the same as what I have written above, excluding the normalisation factor.
The point of this section is that, this singlet state is (spatially) singly degenerate. The orbital diagrams that organic chemists are used to are inaccurate; the actual electronic terms are in fact linear combinations of electronic configurations. So, using orbital diagrams to predict degeneracy can fail!
3. Why is there a Jahn–Teller (JT) distortion in a singly degenerate state?
You are right in that for a JT distortion - a first-order JT distortion, to be precise - to occur, the ground state must be degenerate. This criterion is often rephrased in terms of "asymmetric occupancy of degenerate orbitals", but this is just a simplification for students who have not studied molecular symmetry and electronic states yet. Jahn and Teller made the meaning of their theorem very clear in their original paper;2 already in the first paragraph they write that the stability of a polyatomic molecule is not possible when "its electronic state has orbital degeneracy, i.e. degeneracy not arising from the spin [...] unless the molecule is a linear one". There is no mention of "occupancy of orbitals".
So, this is not a first-order JT distortion, which requires degeneracy of the ground state. Instead, it is a second-order JT distortion. Therefore, most of the other discussion so far on the topic has unfortunately been a little bit off-track. (Not that I knew any better.)
Albright et al. describe this as a "pseudo-JT effect" (p 136). For a better understanding I would recommend reading the section on JT distortions (pp 134ff.), it is very thorough, but you would need to have some understanding of perturbation theory in quantum mechanics. Please note that there is a typo on p 136, as is pointed out in the comments to this answer. (Alternatively, there is also a discussion of second-order JT effects in a paper by Pearson.3) The idea is that, the distortion reduces the symmetry of the molecule and therefore allows ground and excited states to mix. This leads to a stabilisation of the ground state and hence an electronic driving force for the distortion.
In slightly more detail, we can use a group theoretical treatment to do this. The second-order correction to the energy is given by
$$\sum_{j\neq i}\frac{\langle i |(\partial H/\partial q)| j \rangle}{E_i^{(0)} - E_j^{(0)}}$$
where $\ket{i}$ is the ground state (which in this case is described by $2^{-1/2}(\ket{aa} - \ket{bb})$) and $\{\ket{j}\}$ are the excited states. $E_i^{(0)}$ is the unperturbed energy of the ground state $\ket{i}$ and likewise for $E_j^{(0)}$. $q$ is a vibrational coordinate of the molecule.
Note that the denominator in this term, $E_i^{(0)} - E_j^{(0)}$, is negative. So, this correction to the energy is always negative, i.e. stabilising. (There is another second-order term which isn't relevant to the discussion here.) The idea is that in order to have a significant stabilisation upon distortion, two criteria must be fulfilled:
- The denominator must be small, i.e. the excited state must be low-lying.
- The numerator must be non-zero, i.e. there must be appropriate symmetry such that $\Gamma_i \otimes \Gamma_j \otimes \Gamma_H \otimes \Gamma_q$ contains the totally symmetric irreducible representation. Since $H$ transforms as the TSIR, this is equivalent to the criterion that $\Gamma_q = \Gamma_i \otimes \Gamma_j$.
In cyclobutadiene, the two SOMOs in $D_\mathrm{4h}$ symmetry transform as $\mathrm{E_u}$. To find the irreps of the resultant states we take
$$\mathrm{E_u \otimes E_u = A_{1g} + [A_{2g}] + B_{1g} + B_{2g}}$$
(Note that there is no spatially degenerate term in this direct product, so there cannot be a first-order JT distortion!) The $\mathrm{A_{2g}}$ state in square brackets corresponds to the triplet case and we can ignore it. In our case (presented without further proof, because I'm not sure why myself) the ground state is $\mathrm{B_{1g}}$. The two low-lying singlet excited states transform as $\mathrm{A_{1g} \oplus \mathrm{B_{2g}}}$. Therefore if there exists a vibrational mode transforming as either $\mathrm{B_{1g}} \otimes \mathrm{A_{1g}} = \mathrm{B_{1g}}$ or $\mathrm{B_{1g}} \otimes \mathrm{B_{2g} = \mathrm{A_{2g}}}$, then the numerator will be nonzero and we can expect this vibrational mode to lead to a distortion with concomitant stabilisation.
In our case, it happens that there is a vibrational mode with symmetry $\mathrm{B_{1g}}$, which leads to the distortion from a square to a rectangle. This therefore allows for a mixing of the ground $\mathrm{B_{1g}}$ state with the excited $\mathrm{A_{1g}}$ state upon distortion, leading to electronic stabilisation.
Another way of looking at it is that upon lowering of the symmetry from square ($D_\mathrm{4h}$) to rectangular ($D_\mathrm{2h}$), both the $\mathrm{B_{1g}}$ ground state and $\mathrm{A_{1g}}$ excited state adopt the same symmetry $\mathrm{A_g}$. You can prove this by looking at a descent in symmetry table. Therefore, mixing between these two states will be allowed in the new geometry.
Because we are talking about the mixing of states here, the orbital picture is not quite enough to describe what is going on, sadly. The best representation we can get is probably this:
However, Nakamura et al. have produced potential energy curves for the distortion of cyclobutadiene from square to rectangular and rhomboidal geometry.4 Here is the relevant one (rectangular geometry):
Solid lines indicate singlet states, and dashed lines triplet states. Note that this diagram depicts what I have described above: the singlet ground state ($\mathrm{{}^1B_{1g}}$), as well as the first excited singlet state ($\mathrm{{}^1A_{1g}}$), both have the same irrep upon distortion. The ground state is stabilised upon the lowering of symmetry, precisely due to this mixing.
References
(1) Voter, A. F.; Goddard, W. A. The generalized resonating valence bond description of cyclobutadiene. J. Am. Chem. Soc. 1986, 108 (11), 2830–2837. DOI: 10.1021/ja00271a008.
(2) Jahn, H. A.; Teller, E. Stability of Polyatomic Molecules in Degenerate Electronic States. I. Orbital Degeneracy. Proc. R. Soc. A 1937, 161 (905), 220–235. DOI: 10.1098/rspa.1937.0142.
(3) Pearson, R. G. The second-order Jahn–Teller effect. J. Mol. Struct.: THEOCHEM 1983, 103, 25–34. DOI: 10.1016/0166-1280(83)85006-4.
(4) Nakamura, K.; Osamura, Y.; Iwata, S. Second-order Jahn–Teller effect of cyclobutadiene in low-lying states. An MCSCF study. Chem. Phys. 1989, 136 (1), 67–77. DOI: 10.1016/0301-0104(89)80129-6. $\endgroup$