Can Hartree-Fock mix more orbitals than there are electrons? I'm trying to understand Hartree-Fock well enough to write my own implementation. A point that I seem to be coming to that I'm not sure is correct: the method only produces one-electron orbitals that are linear combinations of the starting orbitals?

That is, if I have $n$ electrons, I start them off in some $n$ orbitals (say, atomic orbitals). Is it correct that the final result will only mix those $n$?

It seems very natural that I might want to include other orbital shapes -- maybe I want the first $m \gg n$ atomic orbitals in my molecule, and let all of those contribute to building lower-energy one-electron orbitals. As best s I understand it, standard HF does not allow for this, because at each step of diagonalizing the Fock matrix it's just using combinations of the existing orbitals -- with no room for the excited states.

Is my understanding correct? Are there other extensions of HF which allow broader mixings?

One possibility that I could imagine would be starting off with $m$ electrons, so that I have $m$ orbitals available. Then I run HF, and at the end I only use the $n$ lowest orbitals. This runs into the problem that I will still be computing Coloumb terms from the $(m-n)$ extra "fake" electrons, which will disturb the shape of my $n$ lowest orbitals. So perhaps the correct approach would be to use $m$ electrons mixing, but when I compute the Fock matrix $F$, I only include $J$ and $K$ terms from the $n$ lowest modes (but letting these affect the energies of all $n$ modes). But that seems messy, because now $F$ isn't symmetric... in the end this seems like a risky path to go down.

Help appreciated, thank you.
First some general considerations

Using more orbitals than electrons is not just possible, but almost essential for Hartree-Fock. Although in practice one can do calculations in a so called minimal basis set, this approach is lacking a lot of accuracy. Usually one adds additional basis functions, which leads to basis sets classified as double zeta (triple, quadruple and so on are also possible). 

For open shell systems even a minimal basis set includes more atomic orbitals than electrons are available. For example if you want to calculate a single carbon atom, you would have one basis functions for each of the $1s$, $2s$ and all three $2p$ orbitals. This gives a total of 5 atomic spatial orbitals (=10 spin orbitals), but only 6 electrons. (This is not really a good example, as occupying just one of the three $2p$ orbitals breaks the spherical symmetry of the atom, but I am leaving the scope of your question ...)

Virtual orbitals

Since all molecular orbitals (MOs) are expressed in this basis of atomic orbitals (AOs), it means we will get as many MOs as we put AOs in. If ordered by their energy, we first get a set of occupied MOs until all electrons are distributed. After that comes a set of virtual (or unoccupied) MOs. Those are still soluations of the Roothaan-Hall equations, i.e. they are eigenfunctions of the Fock-Operator, but they do not contribute Hartree-Fock energy.

However, virtual orbitals are important in correlation methods going beyond Hartree-Fock, e.g. Coupled Cluster or Pertubation Theory. Here the virtual orbitals become partially occupied and thus influence the ground state energy.

Equations

The ground state energy in case of closed-shell HF reads

\begin{equation}
  E_0 = 2\sum_a h_{aa} + \sum_a\sum_b 2J_{ab} - K_{ab}
\end{equation}

where the sum only goes over the occupied spatial orbitals, which equals half the number of electrons.
And for the Fock operator we have
\begin{equation}
  f(1) = h(1) + \sum_a 2J_a(1)-K_a(1)
\end{equation}
with the same restrictions to the summation. Since the Fock operator is applied to all MOs, this means the effect (coulomb repulsion) of the occupied orbitals is considered in each virtual orbital, but the virtual orbitals do not affect the occupied ones.

(Equations are taken from the book by Szabo and Ostlund, p. 134)