Is there a way to prove that a bound state wavefunction can always be chosen real for an arbitrary potential in Quantum Mechanics?

As we can prove many things that always (at least in introductory quantum mechanical problems) apply using an arbitrary potential (like that $E>V_{\rm min}$ or else the solutions are non-normalizable and superpositions of them can't produce normalizable wavefunctions), is there a way to generally prove for an arbitrary potential that bound states always correspond to real functions?

Answer

I) Yes, the time independent Schrödinger equation (TISE) of the form $$\left(-\frac{\hbar^2}{2m} {\bf \nabla}^2 +V({\bf r}) -E \right) \psi({\bf r}) ~=~0 $$ is $\mathbb{C}$-linear and invariant$^1$ under complex conjugation. So if the wave function $\psi$ is a solution with finite square-norm, then so will $\psi^{\ast}$, $$\frac{\psi+\psi^{\ast}}{2}\quad\text{and}\quad \frac{\psi-\psi^{\ast}}{2i}$$

be. The two latter are real solutions, and at least one of them has non-zero square-norm if $\psi$ has non-zero square-norm. Hence we can always choose a normalizable solution to be real. See also Problem 2.1b in Griffiths, Intro to QM, and this related Phys.SE post.

II) Note that the same argument does not apply to scattering states in the continuous spectrum, since the boundary conditions at infinity of $\psi^{\ast}$ could be off.

Also it does not apply to the time dependent Schrödinger equation (TDSE).

--

$^1$ Here we have implicitly used that the eigenvalue $E\in\mathbb{R}$ is real because the Hamiltonian $H$ is self-adjoint. Note that self-adjointness alone is not enough. E.g. $$H~=~\frac{1}{2m}\left(\frac{\hbar}{i} {\bf \nabla}-q {\bf A}({\bf r}) \right)^2 +V({\bf r})$$ is self-adjoint, but the corresponding TISE is not invariant under complex conjugation.