In the LCAO method for $H_2^+$, why do we only add or subtract orbitals but never alter one of the phases of the orbitals?

Question

TLDR: Atkins' physical chemistry contains the quote "In general, from N atomic orbitals we can build N molecular orbitals". In the case of $H^+_2$, we use a combination of N=2 orbitals to create the bonding and anti-bonding orbitals as shown in eq 1. But if we allow one of the orbitals to differ by a phase $C_0$ so that $\Psi_{0}(r)=N\left[\psi_{H1s_A}(r)+C_o\psi_{H1s_B}(r)\right]$, then from two atomic orbitals we should be able to create an infinite amount of molecular orbitals. Atkins is correct. I am incorrect. But why?

TLDR Over

Using the LCAO method for the $H_2^+$ ion, we use as our trial functions $$\Psi_{\pm}(r)=N\left[\psi_{H1s_A}(r)\pm\psi_{H1s_B}(r)\right] \tag{1}$$ where $N$ is a normalization factor, and $\psi_{H1s_i}$ is the 1s orbital centered around the $i$'th nucleus. The wave function $\Psi_{+}$ results in a significant reduction in the systems energy whilst the $\Psi_{-}$ wavefunction significantly increases the systems energy. For this reason, $\Psi_{+}$ is called the bonding orbital whilst $\Psi_{-}$ is called the anti-bonding orbital.

This phenomenon is often described as playing a role in the origin of the band structure of solids. When N independent potential wells are brought near each other, the initially N-fold degenerate ground state "splits" into N energy eigenstates. This is simply an N-potential generalization to the $H^+_2$ case described above. In that case, we simply have $N=2$ and the $N=2$ energy eigenstates that arise as the nuclei approach each other are simply the bonding and anti-bonding states.

My issue is that I can't seem to think of any reason as to why we should restrict ourselves to only the two states $\Psi_{+}(r)$ and $\Psi_-(r)$ when performing the LCAO method. Surely the state $\Psi_{0}(r)=N\left[\psi_{H1s_A}(r)+C_o\psi_{H1s_B}(r)\right]$ (where $C_o$ is a complex number) should be just as valid as either of the $\Psi_{\pm}$ states? (Note that anti-symmetrization is not a requirement here since we are calculating 1 electron wave functions).

I also feel like depending on the value of $C_o$, the energy of this new state should vary continuously between the energy of the bonding state and the energy of the anti-bonding state? If this is the case, then the coupling that occurs when $N=2$ hydrogen nuclei are brought to within bonding distance of each other should result in an infinite continuum of energy eigenstates. Not merely two eigenstates. Similarly, when $N\approx 10^{23}$ atoms are brought within bonding distance of each other, we shouldn't get $N$ energy levels but rather an infinite amount of energy levels. This would then ruin the whole notion that energy bands in a solid of N atoms should each contain $N$ energy levels.

Answer 1

Your problem seems to stem from confusing the number of linear independent basis functions(which is the same number as the size of the atomic orbital basis with which we started) and the number of all possible functions that can be built using this basis. Just as in a two dimensional vector space, where you have a maximum of two independent basis functions, but an infinite number of possible vectors by adding those two basis vectors by varying the coefficients.

A way to obtain eigenfunctions is to express the Hamilton operator in a basis, for example, the atomic orbital basis. If the basis is finite, we obtain a finite matrix representation for the Hamilton operator which can be diagonalized to obtain its eigenvectors. These eigenvectors are the molecular orbitals expressed with respect to the atomic orbital basis, which means that we get for each a molecular orbital a set of coefficients that tells us how to linearly combine the original atomic orbitals to form each molecular orbital.

In the finite case, we obtain as many eigenvectors/molecular orbitals, as basis functions that we originally started with. So using a linearly independent atomic orbital basis of size N, will yield a NxN Hamiltonian matrix, which will upon diagonalization yield N eigenvectors=molecular orbitals. These molecular orbitals are not all functions that can be formed. But they are a specifically useful basis set of functions that span the N dimensional space of functions. You can add and subtract these basis functions and generate different functions that are no longer eigenfunctions of the Hamilton operator but it is trivially to calculate the energy as long as you work in this basis since you know the action of the Hamilton operator on each molecular orbital.

The complex number field is simply the field for our vector space of functions but it doesn't change any of the linear algebra or how many eigenvectors/molecular orbitals we obtain.

Now, as mentioned in the comments already, you can build an infinite number of functions/vectors with your basis simply by varying the coefficients. But these functions are not eigenfunctions/molecular orbitals. And you are correct that the energy of these functions varies.

Answer to question in comments:

You can indeed use more basis functions. Most quantum chemistry codes do that actually. The basis where we have one function per electron is the minimal basis set. You may get lower energies for your bonding orbitals with a larger basis. However, you will not neccessarily get more bonding orbitals. Most of the newly obtained molecular orbitals will simply be anti-bonding and energetically way above the "bonding" band. The newly added basis function is, loosely speaking, used to optimize the bonding orbitals but most of its contribution will end up as part of an antibonding molecular orbital as long as your initial smaller basis was not totally inadequate to begin with. They sort of get stacked upon the orbitals that you had already at higher energies. There should only be a little mixing of the new function unless your initial smaller basis was chosen very poorly. Figuratively it should look like this, where the left side uses 2 basis functions and the right side three basis functions,

The lowest molecular orbital should be almost the same in both cases and if we have a system of only two electrons, we would only occupy this one. Likewise with bands in solids. You would fill your band up to the number of available electrons. The new basis functions would mostly lead to an extension of the band at higher energies that won't be occupied. And in this sense the number of bonding orbitals is constant even if you increase the number of basis functions.