This is just Bloch's theorem--- there is a parameter which is the mixing phase between adjacent wells, which is the lattice analog of momentum, and for a large translation invariant lattice of wells, this momentum is a conserved quantity, so you have an eigenstate of H for each value of this phase.
Define the translation operator T to be the operator that moves an infinite line of wells one step to the left. This is a symmetry, and it is a unitary operator, so the eigenvalue of T is a pure phase $e^{i\phi}$. As $\phi$ varies continuously, you parametrize the band. The theory is the same as the band states in a regular conductor.
Perhaps a good way to understand this is through the example of the $H_2$ molecule. Before delving into the molecule, however, let us examine the Hydrogen atom more closely. The $1s$ state of the Hydrogen atom is filled by one electron, but it can take two at most because of spin. Now, the other Hydrogen atom also has one atom in the $1s$ state and can again accommodate at most two.
If we were to bring the two atoms closer together, and then ultimately form a molecule, do you think that the number of states changes? Of course not! The molecule again must be able to accommodate a total of four electrons (two from each atom). (We are ignoring higher energy states for now).
When we think about the Hydrogen molecule, there is a very strong constraint on what the wavefunction must look like. This is because the Hydrogen molecule possesses parity or reflection symmetry. This means that when I switch the atoms, I must either get a $+$ sign or a $-$ sign on the wavefunction. (See the figure below.)
(source: tu-bs.de)
To be more explicit, let's represent the $1s$ states on the Hydrogen atoms 1 and 2 when they are far apart by $|1s_{1}\rangle$ and $|1s_{2}\rangle$. When we bring the atoms close together, the reflection symmetry constraint imposes a severe constraint. Let's see how we can combine the wavefunctions of the individual atoms to satisfy this constraint. It turns out that the following combinations work quite well:
\begin{equation}
|\psi_B\rangle = |1s_1\rangle + |1s_2\rangle
\end{equation}
\begin{equation}
|\psi_A\rangle = |1s_1\rangle - |1s_2\rangle
\end{equation}
where $B$ and $A$ denote bonding and anti-bonding respectively. Now look what happens when I switch the places of the two molecules by applying the reflection operation, $R$:
\begin{equation}
R|\psi_B\rangle = |1s_2\rangle + |1s_1\rangle = +|\psi_B\rangle
\end{equation}
\begin{equation}
R|\psi_A\rangle = |1s_2\rangle - |1s_1\rangle = -|\psi_A\rangle
\end{equation}
This set of wavefunctions therefore satisfies the constraint that we set by imposing reflection symmetry. Now it turns out that because the reflection symmetry, the wavefunctions must look something like this:
Again, each orbital can accommodate two electrons for a total of four, exactly the number of states we started out with. As a rule of thumb, the energy of the orbital with the greater curvature is usually the more energetic one, which is the anti-bonding one in this case. Therefore, the bonding orbital will be filled by the two electrons.
I think that the way you are trying to visualize the problem is flawed. You should not think of the standing wave as a sine wave extending through space infinitely. It is really cut off by the potential of the hydrogen atom and looks like the red and green curves above for the individual atoms. When you bring them closer together and they mix you get the purple and blue curves. Hope this helps!
Best Answer
Molecular orbitals are simply used to construct $N$-electron basis states (Slater determinants).
First, let us have a one-particle Hilbert space $\cal{H}$ of dimension $M$, with $M$ linearly independent vectors $\{\phi_p\}_{p=1}^M$ forming a basis on it. In non-relativistic one-electron QM, $\cal{H}=L^2(\mathbb{R}^3)\otimes\mathbb{C}^2$ (scattering states not considered), so in principle, $M\rightarrow\infty$ should be taken; we nevertheless keep $M$ finite to properly see the dimensions (also, $M$ is obviously finite in all numerical calculations). For the sake of simplicity, let us further take the basis to be orthonormal: $$\langle\phi_p|\phi_q\rangle=\delta_{pq} \ , $$ although this is not necessary. Since this is a basis, any one-particle function can be expanded in it as $$ |\psi\rangle=\sum_{p=1}^Mc_p|\phi_p\rangle \ . $$
The important thing to realize is that constructing the antisymmetrized products of $N$ such vectors in all ${M}\choose{N}$ possible ways will leave you with an orthonormal basis on the antisymmetric subspace of the $N$-particle Hilbert space. For example, for $N=2$, the antisymmetrized products $$ |\Phi_{pq}\rangle=\frac{1}{\sqrt{2}}\left[ |\phi_p\rangle\otimes|\phi_q\rangle-|\phi_q\rangle\otimes|\phi_p\rangle \right] $$ for $p,q=1,...,M$, $p<q$, form a basis in the antisymmetric subspace of the two-particle Hilbert space $\cal{H}\otimes\cal{H}$ (which is the relevant subspace for two electrons). Any two-electron function can be expanded in this ${{M}\choose{2}}=M(M-1)/2$ -dimensional basis: $$ |\Psi\rangle=\sum_{\substack{p,q=1 \\ p<q}}^Mc_{pq}|\Phi_{pq}\rangle \ , $$ orthonormality being understood as $$ \langle\Phi_{rs}|\Phi_{pq}\rangle=\delta_{pr}\delta_{qs}-\delta_{ps}\delta_{qr} \ . $$ The idea is the same for the general $N$-electron case, the antisymmetrized products (Slater determinants) $$ |\Phi_{p_1p_2...p_N}\rangle= \frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}} \left[|\phi_{p_1}\rangle\otimes|\phi_{p_2}\rangle\otimes...\otimes|\phi_{p_N}\rangle\right] $$ can be used in $$ |\Psi\rangle=\sum_{\substack{p_1,p_2,...,p_N=1 \\ p_1<p_2<...<p_N}}^Mc_{p_1p_2...p_N}|\Phi_{p_1p_2...p_N}\rangle \ . $$ Such expansions are used extensively in the many-body theory. In the most basic version of Full-CI, these ${M}\choose{N}$ Slater determinants are used to construct the matrix representation of the Hamiltonian, and the energies and coefficients $c_{p_1p_2...p_N}$ of the eigenstates are obtained from the diagonalization of this matrix.
Note that, apart from orthonormality, I assumed nothing about the one-particle basis. These orbitals are usually found from e.g. Hartree-Fock or Kohn-Sham calculations. On the level of theory, LCAO does not have much to do with any of this, it is just often (almost always) useful to expand molecular orbitals over atom-centered functions: $$ |\phi_p\rangle=\sum_{\mu=1}^Md_{\mu p}|\chi_\mu\rangle \ . $$ This is practical from a computational/interpretational point of view, as it is often easier to see e.g. bonding/antibonding orbitals in this way. But the expansion over atomic functions is by no means necessary, you can do quantum chemistry in e.g. a plane wave basis, too.
Update
I used bra-ket notations above in order to properly show the tensor product structure. The (probably) better known form of these functions is recovered by projection with the formal coordinate eigenstates $|\vec{r}\rangle$. The one-electron wave functions (orbitals) then read $$ \phi_{p}(\vec{r})= \langle\vec{r}|\phi_p\rangle= \left( \begin{matrix} \phi_{p\uparrow}(\vec{r}) \\ \phi_{p\downarrow}(\vec{r}) \end{matrix} \right) \ , $$ having two components due to their spinor character (this is the $\mathbb{C}^2$ part of ${\cal{H}}$). The components $\phi_{p\uparrow}(\vec{r})$ and $\phi_{p\downarrow}(\vec{r})$ are different in the general case; one usually wants orbitals to be eigenfunctions of $\hat{s}_z$, which forces one of the components to be zero.
The $N$-electron Slater determinants are obtained similarly after a projection with $$ |\vec{r}_1,\vec{r}_2,...,\vec{r}_N\rangle=|\vec{r}_1\rangle\otimes|\vec{r}_2\rangle\otimes...\otimes|\vec{r}_N\rangle \ , $$ leading to $$ \Phi_{p_1p_2...p_N}(\vec{r}_1,\vec{r}_2,...,\vec{r}_N)= \frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}} \left[\phi_{p_1}(\vec{r}_1)\otimes\phi_{p_2}(\vec{r}_2)\otimes...\otimes\phi_{p_N}(\vec{r}_N)\right] \ , $$ which shows that non-relativistic $N$-electron states have $2^N$ spinor components. Note that ${\cal{P}}$ permutes the orbital indices, not the coordinate indices.
When doing formal manipulations with wave functions, it is often useful to introduce a combined spatial+spin coordinate eigenstate $|x\rangle=|\vec{r},\sigma\rangle$ for $\sigma=\uparrow,\downarrow$, which implements a further projection onto the upper or lower spinor component: $$ \langle\vec{r},\uparrow\hspace{-0.1cm}|\phi_p\rangle=\phi_{p\uparrow}(\vec{r}) \ \ \ , \ \ \ \langle\vec{r},\downarrow\hspace{-0.1cm}|\phi_p\rangle=\phi_{p\downarrow}(\vec{r}) \ . $$ Projecting $|\Phi_{p_1p_2...p_N}\rangle$ with $$ |x_1,x_2,...,x_N\rangle=|x_1\rangle\otimes|x_2\rangle\otimes...\otimes|x_N\rangle $$ will similarly pick out a single one from the $2^N$ spinor components: $$ \Phi_{p_1p_2...p_N}(x_1,x_2,...,x_N)= \frac{1}{\sqrt{N!}}\sum_{{\cal{P}}\in S_N}(-1)^{\pi_{\cal{P}}}{\cal{P}} \left[\phi_{p_1}(x_1)\phi_{p_2}(x_2)...\phi_{p_N}(x_N)\right] \ . $$ In this way, it is often easier to talk about the antisymmetry of the electronic wave function under simultaneous interchanges in coordinate and spin space: $$ \Psi(x_1,x_2)=-\Psi(x_2,x_1) \ . $$