Beyond the educational value, Bent's rule has had significant impact on how we understand wave functions of molecules. It's simplicity is almost unique in the world of quantum chemistry; almost everybody can immediately understand it and apply it. Together with the VSEPR model it can lead to quite accurate predictions on the back of an envelope.
Molecular Orbital Theory vs. Valence Bond Theory (MO vs. VB)
In their respective infinite treatment these theories are congruent. For more on this topic please read What is actually the difference between valence bond theory and molecular orbital theory?
Important for this post is that both theories are approximations for the wave function of a molecule, hence they do not describe bonding beyond the stationary point they are applied to, and they operate within the clamped nuclei (Born-Oppenheimer) approximation.
In MO theory it is common to express the molecular orbitals as linear combinations of atomic orbitals. Those MO are usually chosen to be orthonormal. The MO themselves are delocalised, but there are algorithms to transform them into more intuitive VB type orbitals. This will provide us with a Lewis-like bonding picture. Famously used in this regard is Natural Bond Orbital (NBO) theory.
Bent's rule
Unfortunately I cannot find my copy of Bent's paper, so this answer will be based purely on the definition of Bent's rule.
The IUPAC gold book defines Bent's rule (DOI: 10.1351/goldbook.BT07000):
In a molecule, smaller bond angles are formed between electronegative ligands since the central atom, to which the ligands are attached, tends to direct bonding hybrid orbitals of greater p character towards its more electronegative substituents.
Water is a good example to demonstrate Bent's rule and its interplay with the VSEPR model. In $\ce{H2O}$ you expect oxygen to follow the octet rule, hence there will be four electron pairs around it, and it predicts a generally tetrahedral structure of the molecule. The angle for these shapes is about $\newcommand{\degree}{^\circ}109\degree$. Since hydrogen is more electronegative than no ligand, Bent's rule predicts a smaller angle than that.
Probably due to its simplicity, Bent's rule can be applied within the MO and VB bonding pictures.
Coulson's Theorem
Coulson's theorem is the formal theory which applies to Bent's rule (see Wikipedia).
We choose a set of orthonormal AO, as would be consistent with the LCAO approximation.
$$
\langle\chi_i|\chi_j\rangle = S_{ij} = \delta_{ij}
\begin{cases}
0, & i \neq j\\
1, & i = j
\end{cases}\tag{1}\label{orthonormalAO}
$$
We can transform these orbitals into sets of orthonormal hybrid orbitals.
\begin{align}
\varphi_a &= \sum_i \lambda_i \chi_i &
\langle\varphi_a|\varphi_b\rangle = S'_{ab} &= \delta_{ab}
\begin{cases}
0, & a \neq b\\
1, & b = b
\end{cases}\tag{2}\label{orthonormalHybridAO}
\end{align}
Let's look at a special case, where we choose atomic orbitals to be $\chi_i$ with $i \in \mathrm{s}, \mathrm{p}_x, \mathrm{p}_y, \mathrm{p}_z$. We can construct following hybrid orbitals (also compare Mathematical form of four hybrid orbitals):
\begin{align}\tag{3}\label{orthonormalSP3}
\varphi_{xyz}
&= \lambda_\mathrm{s}\chi_\mathrm{s}
+ \lambda_{\mathrm{p}_x} \chi_{\mathrm{p}_x}
+ \lambda_{\mathrm{p}_y} \chi_{\mathrm{p}_y}
+ \lambda_{\mathrm{p}_z} \chi_{\mathrm{p}_z}\\
\varphi_x
&= \lambda_\mathrm{s}\chi_\mathrm{s}
+ \lambda_{\mathrm{p}_x} \chi_{\mathrm{p}_x}
- \lambda_{\mathrm{p}_y} \chi_{\mathrm{p}_y}
- \lambda_{\mathrm{p}_z} \chi_{\mathrm{p}_z}\\
\varphi_y
&= \lambda_\mathrm{s}\chi_\mathrm{s}
- \lambda_{\mathrm{p}_x} \chi_{\mathrm{p}_x}
+ \lambda_{\mathrm{p}_y} \chi_{\mathrm{p}_y}
- \lambda_{\mathrm{p}_z} \chi_{\mathrm{p}_z}\\
\varphi_z
&= \lambda_\mathrm{s}\chi_\mathrm{s}
- \lambda_{\mathrm{p}_x} \chi_{\mathrm{p}_x}
- \lambda_{\mathrm{p}_y} \chi_{\mathrm{p}_y}
+ \lambda_{\mathrm{p}_z} \chi_{\mathrm{p}_z}
\end{align}
If all $\lambda_i$ are chosen to be $\frac{1}{2}$, then we would come to the famous $\mathrm{sp}^3$ hybrid orbitals.
Let's generalise this a bit and drop the normalisation:
\begin{align}\tag{4}\label{hybrid}
\varphi_a
&= N\cdot(\chi_\mathrm{s} + \lambda_{\mathrm{p}_k} \chi_{\mathrm{p}_k}) &
\leadsto
\varphi_a
&= \chi_\mathrm{s} + \lambda_{\mathrm{p}_k} \chi'_{\mathrm{p}_k}
\end{align}
From \eqref{orthonormalHybridAO} and \eqref{hybrid}, and ignoring the trivial case:
\begin{align}\tag5
\delta_{ab}
&= \langle\varphi_a|\varphi_b\rangle\\
0
&= \langle
\chi_\mathrm{s} + \lambda_{\mathrm{p}_k} \chi'_{\mathrm{p}_k} |
\chi_\mathrm{s} + \lambda_{\mathrm{p}_l} \chi'_{\mathrm{p}_l}
\rangle\\
&= \langle\chi_\mathrm{s}|\chi_\mathrm{s}\rangle
+ \lambda_{\mathrm{p}_k} \langle\chi_\mathrm{s}|\chi'_{\mathrm{p}_k}\rangle
+ \lambda_{\mathrm{p}_l} \langle\chi_\mathrm{s}|\chi'_{\mathrm{p}_l}\rangle
+ \lambda_{\mathrm{p}_k}\lambda_{\mathrm{p}_l}
\langle\chi'_{\mathrm{p}_k}|\chi'_{\mathrm{p}_l}\rangle
\end{align}
Given \eqref{orthonormalAO} we can simplify this. Note that $\chi'$ is actually a linear combination of $\mathrm{p}$ orbitals, and therefore they don't necessarily have to be orthogonal. The angle between two such orbitals is given by the inner product, i.e. $\langle\chi'_{\mathrm{p}_k}|\chi'_{\mathrm{p}_l}\rangle = \cos \theta_{kl}$.
\begin{align}\tag6
0 &= 1 + 0 + 0 + \lambda_{\mathrm{p}_k}\lambda_{\mathrm{p}_l} \cos \theta_{kl}
\end{align}
Therefore Coulson's Theorem is
\begin{align}\tag{7}\label{Coulson}
\cos \theta_{kl} &= \frac{-1}{\lambda_{\mathrm{p}_k}\lambda_{\mathrm{p}_l}}
\end{align}
If you call $\lambda$ the hybridisation index, then you see that the angle between two hybrid orbitals is dependent on the this index. In other words, the angle between two ligands around a central atom determines the $\mathrm{p}$ character of the hybrid orbitals that ideally describe this bond.
What you can also derive from this formula is, given a constant angle $\theta_{kl}$, if you use a hybrid-orbital with increased $\mathrm{p}$ character, the other hybrid-orbital must have less $\mathrm{p}$ character, i.e. $\lambda_{\mathrm{p}_k} < \lambda_{\mathrm{p}_l}$.
This is in principle all that Bent's rule requires in the framework of LCAO-MO theory to be consistent; the argument about the electronegativity of a ligand is based on observation of molecular structure.
I hope the above answers the titular question sufficiently. In the second part I'll try to show where the original arguments fail.
First of all, you are basing your argumentation on even more handwavy mathematics then I did above. You are also not considering the actual molecular structure and the implications of that in terms of the LCAO-MO approximation of the wave function.
The orbitals contributing to the O-H bonds are the lower two in energy.
This is not quite correct. In LCAO-MO (or MO in general) theory all orbitals contribute to bonding.
Ignoring orbital mixing for now, these are formed from linear combinations of the hydrogen s orbitals with the oxygen 2s (giving 2a1) and 2py orbitals (giving 1b2). The lone pairs are a pure px orbital (1b1) and a combination of pz with a small contribution from the hydrogen s orbitals (3a1).
The shown MO scheme does not ignore orbital mixing, and neither do you in your argumentation. If the oxygen's $\mathrm{s}$ and $\mathrm{p}$ orbitals were far enough apart to not mix, then only one sort of these orbitals could interact with the hydrogen orbitals to lead to bonding. This would either lead to a situation where there is no oxygen $\mathrm{s}$-character in the bonds, or all of it.
This basically follows the argument that for the heavier homologues of water the oxygen valence $\mathrm{s}$-orbital becomes the lone pair.
Based on the contributing oxygen atomic orbitals, that gives 100% p lone pairs and 50% p/50% s for the bonds, the opposite of the result predicted by Bent's rule, which is that the lone pairs should have more s character.
I don't have an optimal LCAO-MO solution at hand, but I have the next best thing: a HF/STO-3G wave function with $C_\mathrm{2v}$ symmetry, $\angle(\ce{HOH}) = 100.0\degree$, and $d(\ce{OH}) = \pu{98.9 pm}$. That is off by plenty, but it's Hartree-Fock, so that was expected. The general tenor wouldn't change, if you were to run it on the experimental structure. This is the output of the wave function:
Atomic contributions to Alpha molecular orbitals:
Alpha occ 1 OE=-20.252 is O1-s=1.0006
Alpha occ 2 OE=-1.258 is O1-s=0.7899 H2-s=0.0928 H3-s=0.0928 O1-p=0.0245
Alpha occ 3 OE=-0.594 is O1-p=0.5365 H3-s=0.2317 H2-s=0.2317
Alpha occ 4 OE=-0.460 is O1-p=0.6808 O1-s=0.1328 H2-s=0.0932 H3-s=0.0932
Alpha occ 5 OE=-0.393 is O1-p=1.0000
Alpha vir 6 OE=0.582 is H3-s=0.3144 H2-s=0.3144 O1-p=0.2947 O1-s=0.0766
Alpha vir 7 OE=0.693 is O1-p=0.4635 H2-s=0.2683 H3-s=0.2683
Unfortunately this doesn't tell us much, apart from that there is obviously one $\mathrm{p}$ lone pair, which was expected due to symmetry. (In MO theory though, this would still count as a bonding π orbital).
From the above table, one would derive that the majority of bonding comes from MO 3, which is a nearly pure oxygen $\mathrm{p}$ orbital.
I've run an NBO6 analysis on this calculation, which is basically a unitary transformation of the MO, and arrived at the following description:
(Occupancy) Bond orbital / Coefficients / Hybrids
------------------ Lewis ------------------------------------------------------
1. (2.00000) CR ( 1) O 1 s(100.00%)
1.0000 0.0000 0.0000 0.0000 0.0000
2. (2.00000) LP ( 1) O 1 s( 71.29%)p 0.40( 28.71%)
0.0000 0.8443 0.0000 0.0000 0.5358
3. (2.00000) LP ( 2) O 1 s( 0.00%)p 1.00(100.00%)
0.0000 0.0000 1.0000 0.0000 0.0000
4. (1.99926) BD ( 1) O 1- H 2
( 59.24%) 0.7697* O 1 s( 14.36%)p 5.97( 85.64%)
0.0000 0.3789 0.0000 -0.7071 -0.5970
( 40.76%) 0.6384* H 2 s(100.00%)
1.0000
5. (1.99926) BD ( 1) O 1- H 3
( 59.24%) 0.7697* O 1 s( 14.36%)p 5.97( 85.64%)
0.0000 0.3789 0.0000 0.7071 -0.5970
( 40.76%) 0.6384* H 3 s(100.00%)
1.0000
---------------- non-Lewis ----------------------------------------------------
6. (0.00074) BD*( 1) O 1- H 2
( 40.76%) 0.6384* O 1 s( 14.36%)p 5.97( 85.64%)
0.0000 -0.3789 0.0000 0.7071 0.5970
( 59.24%) -0.7697* H 2 s(100.00%)
-1.0000
7. (0.00074) BD*( 1) O 1- H 3
( 40.76%) 0.6384* O 1 s( 14.36%)p 5.97( 85.64%)
0.0000 -0.3789 0.0000 -0.7071 0.5970
( 59.24%) -0.7697* H 3 s(100.00%)
-1.0000
This results in two $\mathrm{sp}^6$ bonding orbitals, one pure $\mathrm{p}$ lone pair and a $\mathrm{sp}^{0.4}$ lone pair. From this we would conclude that more $\mathrm{p}$ character is directed towards the hydrogen ligands.
If we properly account for orbital mixing, there is no effect on the 1b1 lone pair orbital, [...]
That is correct, this is due to symmetry.
[...] while the 3a1 lone pair does in fact get more s character, [...]
In your argumentation 3a1 would have been a pure $\mathrm{p}$ lone pair, with mixing it does gain $\mathrm{s}$ character, which is obviously more than none.
[...] but this increased s character is not expected to exceed 50%, since that would just mean a switch of label between the 2a1 and 3a1 orbitals. 50% s in the 3a1 represents the maximum mixed case.
The part with the switch is only true if the energy of these orbitals also switches and the contribution of the hydrogen were equal in both orbitals.
Does this argument also imply that an $\mathrm{sp}$ orbital was the maximum mixed orbital?
I'm very sorry, but I cannot follow the logic as I do not understand what is meant with maximum mixed case. In LCAO you are able to continuously mix $\mathrm{s}$ and $\mathrm{p}$ orbitals, if that actually produces a sensible approximation to the wave function is a completely different topic.
That gives a final result of a maximum of 25% s in the lone pairs (50% in the 3a1 and 0% in 2b1) and 25% in the O-H bonds (2a1 and 1b2), still contradicting Bent's Rule, which predicts that the lone pairs should have more s than the bonds, not the same amount.
This would be exactly the result which you were to expect from tetrahedral coordination: two $\mathrm{sp}^3$ bonds, one $\mathrm{sp}$ and one $\mathrm{p}$ lone pair.
[One a side note, we know that H2O does not represent the maximum mixing case, as SH2 has a smaller bond angle consistent with even greater orbital mixing, but that's beside the point here.]
I am really interested in that side note. I still want to understand what is meant by the maximum mixing case, and especially how a smaller bond angle is consistent with more mixing.