on perturbative methods for calculating $\mathcal{C}$.

Let us consider the $\mathcal{PT}$-symmetric Hamiltonian for a harmonic oscillator perturbed by an imaginary cubic potential: $$H={\mbox{$\textstyle\frac{1}{2}$}}p^2+{\mbox{$\textstyle\frac{1}{2}$}}x^2+i\epsilon x^3 \quad(\epsilon~{\rm real}).
\label{e27a}$$

Following the procedure in Ref. [@rr11], we note that the energy eigenstates are solutions of the Schrödinger equation $$H\phi_n(x)=\left(-{\mbox{$\textstyle\frac{1}{2}$}}\frac{d^2}{dx^2}+{\mbox{$\textstyle\frac{1}{2}$}}x^2+i\epsilon x^3\right)\phi_n(x)
=E_n\phi_n(x).
\label{e27b}$$ The eigenstates and the corresponding eigenvalues may be expressed as series in powers of $\epsilon$ by perturbing around the known energy eigenstates and eigenvalues of the harmonic oscillator. To second order in perturbation theory the eigenstates have the form $$\label{e27c}
\phi_n(x)=\frac{i^n a_n}{\sqrt{2^nn!\sqrt{\pi}}}e^{-\frac{1}{2}x^2}
\Big[H_n(x)+\epsilon P_n(x)+\epsilon^2 Q_n(x)\Big]$$ with energies given by $$E_n=n+{\mbox{$\textstyle\frac{1}{2}$}}+\epsilon A_n+\epsilon^2 B_n.
\label{e27d}$$ Here, $P_n(x)$ and $Q_n(x)$ are polynomials in $x$ of degree $n+3$ and $n+6$, respectively; $a_n$ is a normalization constant to be determined. We include a factor of $ i^n$ because the unperturbed wavefunctions have the form $$\phi_n^{(0)}(x)=\frac{i^n}{\sqrt{2^nn!\sqrt{\pi}}}e^{-\frac{1}{2}x^2}H_n(x),
\nonumber$$ where $H_n(x)$ are Hermite polynomials. This ensures that the unperturbed wavefunctions are eigenstates of the $\mathcal{PT}$ operator with unit eigenvalue: $$\mathcal{PT}\phi_n^{(0)}(x)=\phi_n^{(0)}(x).$$ The wave functions are $\mathcal{PT}$-normalized according to $$\int_{-\infty}^{\infty}dx\left[\phi_n^{(0)}(x)\right]^2=(-1)^n.$$

First-Order Calculation of the Energy Eigenstates and Eigenvalues
-----------------------------------------------------------------

Note that to order $\epsilon^0$, (\[e27b\]) becomes $$\label{e27e}
\left[-{\mbox{$\textstyle\frac{1}{2}$}}\frac{d^2}{dx^2}+x\frac{d}{dx}-n \right]H_n(x)= 0$$ and to order $\epsilon^1$, we have $$\left[-{\mbox{$\textstyle\frac{1}{2}$}}\frac{d^2}{dx^2}+{\mbox{$\textstyle\frac{1}{2}$}}x^2-n-{\mbox{$\textstyle\frac{1}{2}$}}\right]\left(e^{-\frac{1}{2}x^2}
P_n(x)\right)=(A_n-ix^3)e^{-{\mbox{$\textstyle\frac{1}{2}$}}x^2}H_n(x),$$ or $$\left[-{\mbox{$\textstyle\frac{1}{2}$}}\frac{d^2}{dx^2}+x\frac{d}{dx}-n\right]P_n(x)=(A_n-ix^3)H_n(x).$$

Since 