Solution (source code)

= Solution

For a nonzero trial vector use the <Rayleigh quotient> $\mathcal E(\psi)=\langle\psi,H\psi\rangle/\langle\psi,\psi\rangle$. Its first variation in any allowed direction $\eta$ is
$$
\delta\mathcal E=\frac{2\operatorname{Re}\langle\eta,(H-\mathcal E)\psi\rangle}{\|\psi\|^2}.
$$
Vanishing for all complex directions implies $(H-\mathcal E)\psi=0$; an <eigenvector> conversely makes it vanish. The ground <energy> is the infimum of this quotient for a lower-bounded <Hamiltonian>, so minimizing over a trial family provides an upper bound to it, exact when that family contains a ground <eigenvector>.

Set $a=\langle\phi,H\phi\rangle$, $b=\langle\phi,I\phi\rangle$ and $c=\langle\phi,HI\phi\rangle$. All three are real: $H,I$ are <self-adjoint> and commute. Also $I^\dagger I=I^2=1$. Expansion of the numerator and norm gives
$$
\boxed{\mathcal E(t)=\frac{(1+t^2)a+2tc}{1+t^2+2tb}.}
$$
Differentiate the quotient, retaining only nonzero trial vectors:
$$
\mathcal E'(t)=\frac{2(c-ab)(1-t^2)}{(1+t^2+2tb)^2}.
$$
Thus the stationarity equation is $\boxed{(t^2-1)(c-ab)=0}$. If $c=ab$, the quotient is identically $a$ wherever defined. Otherwise $t=\pm1$ are the two extrema, with
$$
\boxed{\mathcal E_\pm=\frac{a\pm c}{1\pm b},\qquad \mathcal E_+-\mathcal E_-=\frac{2(c-ab)}{1-b^2}.}
$$
In this nonconstant case $|b|<1$, so both denominators are positive. If $c-ab>0$, the plus value is the maximum and the minus value the minimum; the assignments reverse when $c-ab<0$. One can also see global extremality by resolving the trial vector into the orthogonal $I=\pm1$ sectors: its <energy> is a weighted mean of their two fixed projected energies.

The null-vector exception matters. If $|b|=1$, then $I\phi=\pm\phi$, the covariance $c-ab$ vanishes and one of $\phi\pm I\phi$ is zero. That value has no <Rayleigh quotient> and is not an extra stationary physical state.

For particle exchange, change of integration variables shows that $I\phi(x_1,x_2)=\phi(x_2,x_1)$ is <unitary> and $I^2=1$, hence <self-adjoint>. Exchanging coordinates interchanges the two kinetic derivatives and leaves the symmetric potential unchanged, so $[I,H]=0$. The vectors $\phi\pm I\phi$ lie in its symmetric and antisymmetric sectors. Their normalized nonzero versions are respectively <bosonic> and <fermionic> spatial trial <wavefunctions>.

Because the <energy> of a general trial state is the weighted average of its sector energies, at least one nonzero symmetry <projection> has <energy> no greater than the original trial. Thus minimizing can be performed in a definite exchange-symmetry sector. This is a variational approximation within the trial subspace; it does not assert that projecting an arbitrary trial vector produces an exact <eigenvector> of $H$. For particles with a prescribed statistics sector one minimizes within that physical sector.