Why do bonding and anti-bonding states split in energy (mathematically)? I am a physicist, and when I hear about the splitting of energy levels for MO antibonding and bonding states, I have a rough intuition that tells me this is reasonable to expect, but I can't think of (nor find) a good formal explanation.

My best shot at understanding this is borrowing from Avoided Crossing and saying that if I have a 2 atom system with a total Hamiltonian $ H = H_1 + H_2 $ that is the sum of each individual atom Hamiltonians, then I can take the Eigenstates for each individual Hamiltonian as a basis.

$$\begin{align}
H_1 | \psi_1 \rangle &amp;= E_1 | \psi_1 \rangle \\
H_2 | \psi_2 \rangle &amp;= E_2 | \psi_2 \rangle
\end{align}$$

Then the total hamiltonian in this basis becomes,
$$\begin{align}
H &amp;= \left(\begin{matrix}
E_1 &amp; \langle \psi_2 | H_1 | \psi_2 \rangle \\
\langle \psi_1 | H_2 | \psi_1 \rangle &amp; E_2
\end{matrix}\right)
\end{align}$$

In the limit of the atoms being really far apart, I would assume the off diagonal terms would vanish, and then the eigenstates are just the atomic orbitals. But if the atoms come closer together, then we get off diagonal terms which causes the splitting as described in avoided crossing.

Is this correct?
You are (essentially) correct and the off-diagonal terms are related to so-called "overlap" between atomic orbitals (AOs). It is common to talk of the AOs "overlapping" to form molecular orbitals (MOs), and as you rightly said, the off-diagonal terms (or "overlap") goes to zero in the limit of infinite bond length. The term "avoided crossing" is not commonly seen in introductory MO theory, but I think it is fair to describe it as such.

Of course there are some slight subtleties. Let's stick to a one-electron system such as $\ce{H2+}$ for now. $\psi_i$ refers to the 1s orbital on hydrogen atom $i$ and $E$ is the energy of a 1s orbital on hydrogen. The minor issue is that $\hat{H} \neq \hat{H}_1 + \hat{H}_2$ (following your notation where $\hat{H}$ is the molecular Hamiltonian and $\hat{H}_i$ is the Hamiltonian for hydrogen atom $i$), because:

$$\begin{align}
\hat{H} &amp;= -\frac{\nabla^2}{2} - \frac{1}{r_1} - \frac{1}{r_2} + \frac{1}{R} \\
\hat{H}_1 &amp;= -\frac{\nabla^2}{2} - \frac{1}{r_1} \\
\hat{H}_2 &amp;= -\frac{\nabla^2}{2} - \frac{1}{r_2} \\
\end{align}$$

(I used atomic units. $r_i$ is the distance between the electron and nucleus $i$, and $R$ is the distance between the two nuclei). So actually, it is more like

$$\hat{H} = \hat{H}_1 - \frac{1}{r_2} + \frac{1}{R}$$

and therefore the first term in the matrix

$$\begin{align}
\mathbf{H}_{11} &amp;= \left&lt;\psi_1\middle|\hat{H}_1 - \frac{1}{r_2} + \frac{1}{R}\middle|\psi_1\right&gt; \\
&amp;= E + \left&lt;\psi_1\middle|-\frac{1}{r_2}\middle|\psi_1\right&gt; + \frac{1}{R} \neq E
\end{align}$$

(since $\langle\psi_1|\hat{H}_1|\psi_1\rangle = E$, and also $\langle\psi_1|(1/R)|\psi_1\rangle = (1/R)\langle\psi_1|\psi_1\rangle = 1/R$).

For more information I refer you to Atkins' Molecular Quantum Mechanics 5th ed., pp 262–266. The terminology can be a bit confusing and also Atkins does not use atomic units (in his text he uses $j_0 = e^2/4\pi\varepsilon_0$ which is simply equal to $1$ in atomic units). Adjusting for this, the final quoted formulae for the matrix elements are:

$$\alpha = \mathbf{H}_{11} = \mathbf{H}_{22} = E - j' + \frac{1}{R}; \qquad j' = \left&lt;\psi_1\middle|\frac{1}{r_2}\middle|\psi_1\right&gt;$$

(as found earlier), and

$$\beta = \mathbf{H}_{12} = \mathbf{H}_{12} = \left(E + \frac{1}{R}\right)S - k'; \qquad S = \langle\psi_1|\psi_2\rangle;\quad k' = \left&lt;\psi_1\middle|\frac{1}{r_2}\middle|\psi_2\right&gt;$$

If you now solve the generalised eigenvalue equation $\mathbf{Hc} = E\mathbf{Sc}$ (the simple $\mathbf{Hc} = E\mathbf{c}$ only works when your basis set is orthonormal, but in this case it is not) you obtain the two eigenvalues

$$E_\pm = E + \frac{1}{R} - \frac{j' \pm k'}{1 \pm S}$$

Note that when the atoms are infinitely far apart both $S$ and $k'$ vanish, so (1) the off-diagonal element $\beta$ vanishes and (2) the "bonding" and "antibonding" orbitals are of the same energy, i.e. $E_+ = E_-$.