I understand why all non-degenerate energy states can be chosen to be real up to an overall phase (as is highlighted here Is a non-degenerate wavefunction real or complex?).
However I've been told the argument holds in the opposite direction -- (that all energy eigenstates that can be represented as purely real are non-degenerate).
Is this true? I'm having trouble finding a proof of this.
Since the time-independent Schrödinger equation is real, is it such a surprise that the solutions should be real?
The reverse of your claim is definitely not true (at least not in general): it is perfectly possible to have real combinations of hydrogen solutions that are still degenerate. This is because energies of hydrogenoid atoms are of the form $-13.6$eV/$n^2$ and do not depend on the angular momentum $\ell$. Thus, you can take real (or complex) linear combinations such of two solutions with the same energy (i.e. same $n$) but different $\ell$ and the result will still be real (or complex). You can even arrange the angular momentum part to be real by choosing $Y_{\ell}^m(\theta,\phi)+(Y_{\ell}^m(\theta,\phi))^*$, which will again not affect the energies as they do not depend on $\ell$ or $m$. (Indeed many chemists work with "real" spherical harmonics.)
What is true is that, in a given effective potential (which would include the $\ell(\ell+1)$ term for hydrogen for instance) for fixed $\ell$, the solutions $R_{n \ell}$ and $R_{n' \ell}$ with $n'\ne n$ have different energies.
It is also true for $1d$ motion, in which case we're back to the original argument that you are solving a real differential equation so the solutions can be chosen as real. There is indeed a "non-degeneracy" theory for 1d motion in an effective potential. The proof can be found in many textbooks, such as Messiah (which I believe has the clearest proof).
Sketch of proof (if this is not plagued by typos):
Suppose $\psi\left( \xi \right)$ is solution to the Schrodinger equation describing the $n$'th \emph{bound state }of a system, having energy $E_{n},$ i.e. $\psi \left( \xi \right)$ is solution to the Schr"odinger equation
\begin{equation} \frac{d^{2}}{d\xi ^{2}}\psi \left( \xi \right) +\frac{2m}{\hbar ^{2}}\left[ E_{n}-V(\xi )\right] \psi \left( \xi \right) =0. \end{equation} Then, there is no other function, other than multiple of $\psi \left( \xi \right) ,$ that has this energy.
This is a proof by contradiction. Indeed assume $\psi \left( \xi \right) $ and $\phi \left( \xi \right) $ are both solutions, with $\psi \left( \xi \right) $ $\neq c\phi \left( \xi \right) $ for an arbitrary complex constant $c.$ Then: \begin{equation} \frac{1}{\psi \left( \xi \right) }\frac{d^{2}}{d\xi ^{2}}\psi \left( \xi \right) =-\frac{2m}{\hbar ^{2}}\left[ E_{n}-V(\xi )\right] =\frac{1}{\phi \left( \xi \right) }\frac{d^{2}}{d\xi ^{2}}\phi \left( \xi \right) , \end{equation} which can be rearranged into \begin{eqnarray} 0 &=&\phi \left( \xi \right) \frac{d^{2}}{d\xi ^{2}}\psi \left( \xi \right) -\psi \left( \xi \right) \frac{d^{2}}{d\xi ^{2}}\phi \left( \xi \right) , \nonumber \\ &=&\frac{d}{d\xi }\left( \phi \left( \xi \right) \left( \frac{d}{d\xi }\psi \left( \xi \right) \right) -\psi \left( \xi \right) \left( \frac{d}{d\xi }% \phi \left( \xi \right) \right) \right) , \end{eqnarray} so that \begin{equation} \phi \left( \xi \right) \left( \frac{d}{d\xi }\psi \left( \xi \right) \right) -\psi \left( \xi \right) \left( \frac{d}{d\xi }\phi \left( \xi \right) \right) =\gamma , \tag{1} \end{equation} where $\gamma $ is a constant.
However, for the wave functions to be normalizable, we need \begin{equation} \lim_{\xi \rightarrow \infty }\phi \left( \xi \right) =\lim_{\xi \rightarrow \infty }\psi \left( \xi \right) =0. \end{equation} Therefore, evaluating Eqn.(1) at infinity gives $\gamma =0$ and thus \begin{eqnarray} \phi \left( \xi \right) \left( \frac{d}{d\xi }\psi \left( \xi \right) \right) &=&\psi \left( \xi \right) \left( \frac{d}{d\xi }\phi \left( \xi \right) \right) , \\ \frac{1}{\psi \left( \xi \right) }\frac{d}{d\xi }\psi \left( \xi \right) &=& \frac{1}{\phi \left( \xi \right) }\frac{d}{d\xi }\phi \left( \xi \right) . \end{eqnarray} Integrating gives $\psi \left( \xi \right) $ $=c\phi \left( \xi \right) ,$ with $c$ a constant, in contraction with the assumption that $\psi \left( \xi \right) $ and $\phi \left( \xi \right) $ have the same energy but $\psi \left( \xi \right) $ $\neq c\phi \left( \xi \right) .$
This shows the solution is unique up to normalization and phase. You can always choose the phase so the function is real since the function is a solution to a real differential equation.
Even if the basis of the diagonalization of the Hamilton operator, that is a real symmetric operator, can be constructed as an othonormal system of real valued functions, its unproductive: All such states have spatial current density zero of the velocity form $$0=j = \frac{1}{2 i m} \left( \psi^* (-i\nabla \psi -\frac{e}{c} A)\ \psi -\text{cc}\right) = \Im \psi^* v \ \psi $$
For scalar theories a real basis represents a "bad" quantum basis for calulating transport and magnetic properties.
For free particles in a relativistic theories they even are a mix of classically retarded and advanced solutions, eg outgoing and incoming spherical waves and are useless in quantum field theory.
In Schrödinger theory, the degenerated subspaces of fixed angular momentum squared $$L^2= -\Delta_{\theta,\phi} = -\sin \theta^{ -1} \partial_ \theta \ \sin \theta - \sin \theta^{-2} \partial_{\phi,\phi }, \qquad (S_2, \ d\theta\wedge \sin \theta d\phi)) $$ are the essential paradigma of all quantum theories. The eigenspace is diagonalized by waves with the real eigenvalues of $L_z$ representing circular currents by a factor $e^{i m \phi }$ around the axis, while the $\theta$ waves are real standing waves between the poles with current density zero.The group theoretic analysis is displaying the deep connection of quantization of charge and mass, magnetic moments and geometric Bohr-de Broglie quantization on a compact manifold wihtout boundary.
