The claim follows from the following two propositions:
Proposition 1: Every $\phi$-Hermitian matrix $H$ with all eigenvalues outside the nonpositive real numbers, admits a $\phi$-Hermitian square root $\sqrt{H}$ such that $\sqrt{H} H = H \sqrt{H}$.
Proof. Let $\chi(H)$ be the real representation of $H$. Let $P^{-1}\chi(H) P = J_{m_1}(\lambda_1)\oplus\dotsb\oplus J_{m_k}(\lambda_k)$ be the Jordan form of $\chi(H)$. We know that none of the $\lambda_i$ are nonpositive real numbers. We may therefore use Hermite interpolation to construct a polynomial $p \in \mathbb C[x]$ satisfying $p(J_{m_i}(\lambda_i))^2 = J_{m_i}(\lambda_i)$ and $p(J_{m_i}(\overline{\lambda_i}))^2 = J_{m_i}(\overline{\lambda_i})$ for every $i$. Importantly, observe that all of the coefficients of $p$ are in $\mathbb R$. So we let $\sqrt{H} := p(H)$. Since the coefficients of $p$ are real, we get:
- $p(H)^\phi = p(H^\phi) = p(H)$.
- $\chi(p(H)^2) = p(\chi(H))^2 = \chi(H)$, from which we cancel $\chi$ to obtain $p(H)^2 = H$.
- $p(H) H = H p(H)$
and so we are done.
Proposition 2: Proposition 1 holds for all $H$ which are $\phi$-Hermitian and non-singular, and with no other preconditions on eigenvalues.
Proof. It only remains to consider the case when $H$ has a negative real eigenvalue. We can perturb $H$ to $H + \varepsilon$, preserving its $\phi$-Hermitianess, but making all its eigenvalues be outside $\mathbb R^-$. Let $K$ be a square root of $H + \varepsilon$. By submultiplicativity of the operator norm, we have $\lVert K \rVert \geq \sqrt{\lVert H + \varepsilon \rVert}$. Similarly, by observing $K^{-2} = H^{-1}$, we obtain $\lVert K \rVert \leq \frac 1 {\sqrt{\lVert H + \varepsilon \rVert}}$. So we obtain that the square roots of small perturbations of $H$ have matrix norms within a finite band. By Bolzano-Weierstrass, we may produce a convergent sequence converging to a square root of $H$.
