Why do bonding and anti-bonding states split in energy (mathematically)?
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_-$.