physical applications and draw some conclusions.

Construction of the $\mathcal{C}$ Operator {#s2}
==========================================

We begin by summarizing the mathematical properties of the solution to the Sturm-Liouville differential equation eigenvalue problem $$-\phi_n''(x)+x^2(ix)^\epsilon\phi_n(x)=E_n\phi_n(x)
\label{e3}$$ associated with the Hamiltonian $H$ in (\[e1\]). This differential equation must be imposed on an infinite contour in the complex-$x$ plane. For large $|x|$ this contour lies in wedges placed symmetrically with respect to the imaginary-$x$ axis [@rr3]. The boundary conditions on the eigenfunctions are that $\phi(x)\to0$ exponentially rapidly as $|x|\to\infty$ along the contour. For $0\leq\epsilon<2$, the contour may lie on the real axis.

When $\epsilon\geq0$, the Hamiltonian has an unbroken $\mathcal{PT}$ symmetry. Thus, we know that the eigenfunctions $\phi_n(x)$ are simultaneously eigenstates of the $\mathcal{PT}$ operator: $$\mathcal{PT}\phi_n(x)=\lambda_n\phi_n(x).
\label{e4}$$ As we argued above, $\lambda_n$ is a pure phase and, without loss of generality, for each $n$ this phase can be absorbed into $\phi_n(x)$ by a multiplicative rescaling so that the new eigenvalue under $\mathcal{PT}$ is unity: $$\mathcal{PT}\phi_n(x)=\phi_n^*(-x)=\phi_n(x).
\label{e5}$$

It is not known rigorously yet, but there is strong evidence that when properly normalized the eigenfunctions $\phi_n(x)$ are complete. The coordinate-space statement of completeness reads $$\sum_n(-1)^n\phi_n(x)\phi_n(y)=\delta(x-y)\qquad(x,y~{\rm real}).
\label{e6}$$ This nontrivial result has been verified numerically to extremely high accuracy 