$\newcommand{\ket}[1]{\left|#1\right>}$
$\newcommand{\bra}[1]{\left<#1\right|}$
(1) Is this all correct so far?
Yep.
(2) If so, is it justified to use it in numerical derivations as an
eigenfunction of the Hamiltonian, i.e. can we write
$\hat{H}\ket\Phi\overset{?}{=}E_{0}\ket\Phi$ or are we restricted to
$\bra\Phi\hat{H}\ket\Phi=E_{0}$?
First, you can write $\hat{H}\ket\Phi = E\ket\Phi$ where $\Phi$ is the HF wave function but then $\hat{H}$ has to interpreted as the so-called mean-field Hamiltonian. It is the approximate Hamiltonian in which the exact potential which represents the energy of the Coulomb repulsions between the
electrons ($\hat{V}_\mathrm{ee}$) is replaced with some approximate mean-field potential ($\hat{V}_\mathrm{MF}$) which describes the model system in which electrons do not instantaneously interact with each other, but rather each and every electron interacts with the average, or mean, electric field created by all other electrons. That is basically the definition of the HF wave function, i.e. it is an eigenfunction of the mean-field Hamiltonian:
$$
\hat{H}_\mathrm{MF} \ket\Phi = E_\mathrm{MF} \ket\Phi \, ,
\quad \text{where} \quad
\hat{H}_\mathrm{MF} = \hat{T}_\mathrm{e} + \hat{V}_\mathrm{en} + \hat{V}_\mathrm{MF} \, .
$$
Secondly, despite the fact that the HF wave function is not an eigenfunction of the exact electronic Hamiltonian
$$
\hat{H}_\mathrm{e} \ket\Phi \neq E_\mathrm{e} \ket\Phi \, ,
\quad \text{where} \quad
\hat{H}_\mathrm{e} = \hat{T}_\mathrm{e} + \hat{V}_\mathrm{en} + \hat{V}_\mathrm{ee} \, ,
$$
we could evaluate it energy $\bra\Phi\hat{H}_\mathrm{e}\ket\Phi$ using Slater rules and minimize it to find an upper bound to the ground state energy
$$
\bra\Phi\hat{H}_\mathrm{e}\ket\Phi \geq E_\mathrm{e0}
$$
But the resulting Slater determinant obtained by minimizing the energy is only an approximation to the ground state wave function.
Finally, yes, can we use $\ket\Phi$ as our trial wave function despite the fact that it is not an eigenfunction of $\hat{H}_\mathrm{e}$ since it is not required by the variation method. $\ket\Phi$ is normalizable, it (presumably) satisfies same boundary conditions as the exact wave function, and it has adjustable parameters. That is basically enough to qualify as the trial wave function. And at the end of the day, starting from a known solution for a simpler problem is a very usual practice for variation method.
