Are orbitals observable physical quantities in a many-electron setting? Orbitals, both in their atomic and molecular incarnations, are immensely useful tools for analysing and understanding the electronic structure of atoms and molecules, and they provide the basis for a large part of chemistry and in particular of chemical bonds.
Every so often, however, one hears about a controversy here or there about whether they are actually physical or not, or about which type of orbital should be used, or about whether claimed measurements of orbitals are true or not. For some examples, see this, this or this page. In particular, there are technical arguments that in a many-body setting the individual orbitals become inaccessible to experiment, but these arguments are not always specified in full, and many atomic physics and quantum chemistry textbooks make only casual mention of that fact.
Is there some specific reason to distrust orbitals as 'real' physical quantities in a many-electron setting? If so, what specific arguments apply, and what do and don't they say about the observability of orbitals?
 A: Generally speaking, atomic and molecular orbitals are not physical quantities, and generally they cannot be connected directly to any physical observable. (Indirect connections, however, do exist, and they do permit a window that helps validate much of the geometry we use.)
There are several reasons for this. Some of them are relatively fuzzy: they present strong impediments to experimental observation of the orbitals, but there are some ways around them. For example, in general it is only the square of the wavefunction, $|\psi|^2$, that is directly accessible to experiments (but one can think of electron interference experiments that are sensitive to the phase difference of $\psi$ between different locations). Another example is the fact that in many-electron atoms the total wavefunction tends to be a strongly correlated object that's a superposition of many different configurations (but there do exist atoms whose ground state can be modelled pretty well by a single configuration).
The strongest reason, however, is that even within a single configuration $-$ that is to say, an electronic configuration that's described by a single Slater determinant, the simplest possible many-electron wavefunction that's compatible with electron indistinguishability $-$ the orbitals are not recoverable from the many-body wavefunction, and there are many different sets of orbitals that lead to the same many-body wavefunction. This means that the orbitals, while remaining crucial tools for our understanding of electronic structure, are generally on the side of mathematical tools and not on the side of physical objects.

OK, so let's turn away from fuzzy handwaving and into the hard math that's the actual precise statement that matters. Suppose that I'm given $n$ single-electron orbitals $\psi_j(\mathbf r)$, and their corresponding $n$-electron wavefunction built via a Slater determinant,
\begin{align}
\Psi(\mathbf r_1,\ldots,\mathbf r_n)
& = 
\det
\begin{pmatrix}
\psi_1(\mathbf r_1) & \ldots & \psi_1(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\psi_n(\mathbf r_1) & \ldots & \psi_n(\mathbf r_n)
\end{pmatrix}.
\end{align}

Claim
If I change the $\psi_j$ for linear combinations of them,
  $$\psi_i'(\mathbf r)=\sum_{j=1}^{n} a_{ij}\psi_j(\mathbf r),$$
  then the $n$-electron Slater determinant
  $$
\Psi'(\mathbf r_1,\ldots,\mathbf r_n)
=
\det
\begin{pmatrix}
\psi_1'(\mathbf r_1) & \ldots & \psi_1'(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\psi_n'(\mathbf r_1) & \ldots & \psi_n'(\mathbf r_n)
\end{pmatrix},
$$
  is proportional to the initial determinant,
  $$\Psi'(\mathbf r_1,\ldots,\mathbf r_n)=\det(a)\Psi(\mathbf r_1,\ldots,\mathbf r_n).$$
  This implies that both many-body wavefunctions are equal under the (very lax!) requirement that $\det(a)=1$.

The proof of this claim is a straightforward calculation. Putting in the rotated orbitals yields
\begin{align}
\Psi'(\mathbf r_1,\ldots,\mathbf r_n)
&=
\det
\begin{pmatrix}
\psi_1'(\mathbf r_1) & \cdots & \psi_1'(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\psi_n'(\mathbf r_1) & \cdots & \psi_n'(\mathbf r_n)
\end{pmatrix}
\\&=
\det
\begin{pmatrix}
\sum_{i}a_{1i}\psi_{i}(\mathbf r_1) & \cdots & \sum_{i}a_{1i}\psi_{i}(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\sum_{i}a_{ni}\psi_{i}(\mathbf r_1) & \cdots & \sum_{i}a_{ni}\psi_{i}(\mathbf r_n)
\end{pmatrix},
\end{align}
which can be recognized as the following matrix product:
\begin{align}
\Psi'(\mathbf r_1,\ldots,\mathbf r_n)
&=
\det\left(
\begin{pmatrix}
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn} \\
\end{pmatrix}
\begin{pmatrix}
\psi_1(\mathbf r_1) & \cdots & \psi_1(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\psi_n(\mathbf r_1) & \cdots & \psi_n(\mathbf r_n)
\end{pmatrix}
\right).
\end{align}
The determinant then factorizes as usual, giving
\begin{align}
\Psi'(\mathbf r_1,\ldots,\mathbf r_n)
&=
\det
\begin{pmatrix}
a_{11} & \cdots & a_{1n} \\
\vdots & \ddots & \vdots \\
a_{n1} & \cdots & a_{nn} \\
\end{pmatrix}
\det
\begin{pmatrix}
\psi_1(\mathbf r_1) & \cdots & \psi_1(\mathbf r_n)\\
\vdots & \ddots & \vdots \\
\psi_n(\mathbf r_1) & \cdots & \psi_n(\mathbf r_n)
\end{pmatrix}
\\
\\&=\det(a)\Psi(\mathbf r_1,\ldots,\mathbf r_n),
\end{align}
thereby proving the claim.

Disclaimers
The calculation above makes a very precise point about the measurability of orbitals in a multi-electron context. Specifically, saying things like

the lithium atom has two electrons in $\psi_{1s}$ orbitals and one electron in a $\psi_{2s}$ orbital

is exactly as meaningful as saying

the lithium atom has one electron in a $\psi_{1s}$ orbital, one in the $\psi_{1s}+\psi_{2s}$ orbital, and one in the $\psi_{1s}-\psi_{2s}$ orbital,

since both will produce the same global many-electron wavefunction. This does not detract in any way from the usefulness of the usual $\psi_{n\ell}$ orbitals as a way of understanding the electronic structure of atoms, and they are indeed the best tools for the job, but it does mean that they are at heart tools and that there are always alternatives which are equally valid from an ontology and measurability standpoint.
However, there are indeed situations where quantities that are very close to orbitals become accessible to experiments and indeed get measured and reported, so it's worth going over some of those to see what they mean.
The most obvious is the work of Stodolna et al. [Phys. Rev. Lett. 110, 213001 (2013)], which measures the nodal structure of hydrogenic orbitals (good APS Physics summary here; discussed previously in this question and this one). These are measurements in hydrogen, which has a single electron, so the multi-electron effect discussed here does not apply. These experiments show that, once you have a valid, accessible one-electron wavefunction in your system, it is indeed susceptible to measurement.
Somewhat more surprisingly, recent work has claimed to measure molecular orbitals in a many-electron setting, such as Nature 432, 867 (2004) or Nature Phys. 7, 822 (2011). These experiments are surprising at first glance, but if you look carefully it turns out that they measure the Dyson orbitals of the relevant molecules: this is essentially the overlap
$$
\psi^\mathrm{D}=⟨\Phi^{(n-1)}|\Psi^{(n)}⟩
$$
between the $n$-electron ground state $\Psi^{(n)}$ of the neutral molecule and the relevant $(n-1)$-electron eigenstate $\Phi^{(n-1)}$ of the cation that gets populated. (For more details see J. Chem. Phys. 126, 114306 (2007) or Phys. Rev. Lett. 97, 123003 (2006).) This is a legitimate, experimentally accessible one-electron wavefunction, and it is perfectly measurable.
A: The general answer is that when there is electron correlation the picture of each electron occupying an orbital is no longer adequate. In this case a single Slater determinant is no longer sufficient.
The Hartree-Fock or selfconsistent field approach to atomic and molecular problems approximates the many electron wave function by a single Slater determinant. Although a Slater determinant is invariant under an orthonormal transformation of its orbitals, as remarked above, the eigenfunctions and eigenvalues of the Hartree-Fock operator have special meaning. These can be used to estimate energies and other properties of excited states through the Hellmann-Feynman theorem.
https://en.wikipedia.org/wiki/Hartree–Fock_method
https://en.wikipedia.org/wiki/Hellmann–Feynman_theorem
