Do solutions of the Schrödinger equation for multiple particles automatically obey spin-statistics? Consider the Hamilitonian for a general two-electron system subject to an external potential $V_\mathrm{ext}$ and an interaction potential $V_\mathrm{ee}$. In this case
$$H\psi(x, y) = -\frac{1}{2} \Delta_x \psi(x, y) - \frac{1}{2} \Delta_y \psi(x, y) + V_\mathrm{ext}(x) \psi(x, y) + V_\mathrm{ext}(y) \psi(x, y) + V_\mathrm{ee}(|x - y|)\psi(x, y)$$
where $x$ is the position of the first electron and $y$ is the position of the second electron.
It is well known that the (time-independent) Schrödinger equation $H\psi = E\psi$ does not have an analytic solution in the majority of cases but there is a solution to the equation.
Do we expect this solution to obey the proper spin-statistics for the problem? More specifically, does the differential equation $H\psi = E\psi$ have a unique minimum eigenvalue and eigenfunction pair and if so will it be symmetric or anti-symmetric?
Here are my thoughts so far, the eigenfunction with minimum eigenvalue cannot be unique as the Hamiltonian commutes with an exchange of particles and therefore if $\psi(x, y)$ is an eigenfunction then both $$\psi_\mathcal{S}(x, y) \equiv \frac{1}{2} \left[\psi(x, y) + \psi(y, x)\right]$$ and $$\psi_\mathcal{A}(x, y) \equiv \frac{1}{2} \left[ \psi(x, y) - \psi(y, x) \right]$$ will also be eigenfunctions so long as they are not identically zero.
If the ground state (eigenfunction with minimum eigenvalue) is in fact unique, it must be the case that either $\psi(x, y) = \psi_\mathcal{S}(x, y)$ or $\psi(x, y) = \psi_\mathcal{A}(x, y)$. From intuition, we would expect the latter case will never occur as heuristically bosons always have an equal or lower ground state energy compared to fermions (as both particles can occupy the lowest energy state) and the former case to occur in a "typical" case. Is there a mathematical justification for this?
In the case where $\psi = \psi_\mathcal{S}$ the Schrödinger equation has seemingly failed us as the solution does not obey the correct spin-statistics for electrons. Recalling that the Schrödinger equation is equivalent to the minimization problem
$$\text{minimize } \langle \psi, H \psi \rangle \text{ subject to } \langle \psi, \psi \rangle = 1$$
we might try to salvage things by considering the minimization problem
$$\text{minimize } \langle \psi, H \psi \rangle \text{ subject to } \langle \psi, \psi \rangle = 1 \;\text{ and }\; \psi(x, y) = -\psi(y, x)$$
Is this problem wellposed and if so can we recover a "fermionic" Schrödinger equation by looking at the Euler-Lagrange equations for this problem? If not, is the only option to move on to an equation with the correct spin-statistics built in, i.e., the Dirac equation.
 A: The spin-statistics theorem is a relativistic effect, and arises naturally in relativistic quantum mechanics. The Schrödinger equation is non-relativistic, and spin-statistics correlations must be added as an additional constraint. It’s totally possible to construct a wavefunction which is a solution to the Schrödinger  equation but has the wrong symmetry; that wavefunction just doesn’t describe the system you want.
As a commenter points out, it’s perfectly acceptable for a two-fermion wavefunction to have a symmetric spatial part, so long as the spin part is antisymmetric.  My favorite example is molecular hydrogen, where antisymmetrizing the proton wavefunction means that all of the even-$L$ rotational states are spin singlets, while all of the odd-$L$ rotational states are spin triplets.  Because nuclear spin transfer by collision is inefficient, the even-$L$ “parahydrogen” and the odd-$L$ “orthohydrogen” act like distinct materials at low temperature. They have different thermal properties, and liquifying hydrogen that was recently at room temperature releases about twice as much heat as liquifying hydrogen which has been cold for a while.
This behavior of liquid hydrogen, when discovered in the 1930s, was the original evidence that the two-proton wavefunction is antisymmetric under exchange. Some other diatomic molecules, such as nitrogen and oxygen, have nuclei which obey Bose-Einstein statistics and have different low-temperature behavior. (Note that the ortho-vs-para effect in hydrogen is magnified because of hydrogen’s small mass and low boiling temperature, so the comparison with heavier cryoliquids gets into the weeds experimentally.) Under the Schrödinger equation, both symmetries are allowed; the spin dependence must be tacked on manually to match the experimental behavior.
For an additional example of an antisymmetry being accomplished by "tacking on" an experimentally-observed degree of freedom, consider the deuteron, which is the only bound state in the two-nucleon system.  The spatial part of the wavefunction is symmetric, mostly $s$-wave with a little admixture of $d$-wave.  The spin part of the wavefunction is also symmetric, with the proton and neutron spins aligned.  The antisymmetry requirement appears if you consider the strong-interaction perspective that the proton and neutron are different "isospin projections" of the same particle, the "nucleon."  The deuteron wavefunction is antisymmetric in isospin space, which makes the total wavefunction antisymmetric.
A: The wave function of a two-electron system (i.e. a system of two identical spin $1/2$ fermions) must be antisymmetric under a simultaneous interchange of both the spatial $and$ the spin variables,
$$ \psi(\vec{x}_1,s_1; \vec{x}_2,s_2)=-\psi(\vec{x}_2,s_2;\vec{x}_1,s_1) \qquad (1).$$
This condition is required by the spin-statistics theorem proven by Wolfgang Pauli in 1940 as a consequence of relativisic quantum field theory. It says that particles with half-integer spin ($1/2, 3/2, \ldots$) are fermions and particles with spin $0,1, 2,\ldots$ are bosons. In nonrelativistic quantum mechanics, this condition has to be imposed "by hand" and the Hilbert space of $n$ identical fermions (half-integer spin) has to be chosen as  the totally antisymmetric tensor product of the corresponing one-particle Hilbert space, whereas in the case of $n$ identical bosons, the totally symmetric tensor product of the corresponding one-particle Hilbert space constitutes the suitable state space.
Coming back to to the problem of the two-electron system, it is perfectly consistent with the requirement of the antisymmetry condition $(1)$  that the wave function factorizes into a $symmetric$ spatial wavefunction $\phi(\vec{x}_1,\vec{x}_2)= \phi(\vec{x}_2,\vec{x}_1)$ and an $antisymmetric$ spin part $\chi_{s_1,s_2}=-\chi_{s_2,s_1}$, the full wave function given by  $\psi(\vec{x}_1,s_1;\vec{x}_2,s_2)= \phi(\vec{x}_1,\vec{x}_2) \chi_{s_1,s_2}$. Alternatively, the wave function may also factorize into an $antisymmetric$ spatial factor $\phi(\vec{x}_1,\vec{x}_2)=-\phi(\vec{x}_2,\vec{x}_1)$ and a $symmetric$ spin factor $\chi_{s_1,s_2}=\chi_{s_2,s_1}$.
A well known example is the two-electron system of the helium atom. Parahelium (with total spin $S=0$) corresponds to the situation where the spin part of the wave function is antisymmetric ($\sim \, |\uparrow \downarrow\rangle - |\downarrow \uparrow \rangle$) with a symmetric spatial wave function, whereas orthohelium (with total spin $S=1$) has a symmetric spin triplet $\sim \, |\uparrow \uparrow\rangle, \,  |\uparrow \downarrow \rangle+ |\downarrow \uparrow \rangle, \, |\downarrow \downarrow \rangle$ and an antisymmetric spatial wave function.
