# How to plot the spectrum of infinite dimensional operators?

For a finite dimensional operator it is easy to find eigenvalues in Mathematica. However, I found this example $$H=p^2 + x^2 (ix)^\epsilon$$ in this article. The spectrum of $$H$$ is plotted in figure 1 and looks like following:

How can one obtain this figure using Mathematica?

• It is explained how to compute this spectrum in C. M. Bender and S. Boettcher, Phys. Rev. Lett. 80, 5243 (1998). Just a side remark, you cannot simply write an operator and ask for a spectrum. There are always 2 ingredients needed: operator and boundary conditions (or domain). Thus, your question, as formulated, makes no sense. Commented Oct 29, 2022 at 21:11

This restricts to a finite interval $$[-r,r]$$ and uses simple discretization with $$2n+1$$ points:

(* discretization parameters *)
r=10;
n=500;

(* matrix *)
p2=-n^2/r^2*SparseArray[{
Band[{1,1}]->-2,Band[{2,1}]->1,Band[{1,2}]->1},{2*n+1,2*n+1}]//N;
x=Range[1,n]*r/n//N;
x2[eps_]:=With[{r=x^2*(I*x)^eps},
DiagonalMatrix[Join[Conjugate[r]//Reverse,{0},r]//N]];

ev[eps_?NumericQ]:=Map[{eps,#}&,Chop[Eigenvalues[p2+x2[eps],-10]]];
evs=Table[ev[eps],{eps,-0.7,2,0.05}];

(* plot *)
ListLinePlot[Transpose[evs,{2,1}],PlotMarkers->{"\[FilledDiamond]", 10},
PlotRange->{{-0.7,2},{0,19}}]

Eigenvalues for $$\epsilon = -0.5$$ as an example:

ev[-0.5][[;;,2]]
(*
{13.5532 +4.00504 I,
11.3696 -3.01092 I,
11.3696 +3.01092 I,
9.09025 +1.99448 I,
9.09025 -1.99448 I,
6.65531 +0.951423 I,
6.65531 -0.951423 I,
4.42176,
3.19568,
1.08691}
*)

Note that there seem to be three real eigenvalues, other eigenvalues are complex. The plot above only shows the three real ones.

• why $n^2/r^2$ ? Commented Dec 4, 2022 at 17:59
• See the formula for the 2nd derivative as a limit. The step size $h$ in that formula is $r/n$ in this case, therefore $1/h^2$ is $n^2/r^2$. Commented Dec 4, 2022 at 20:36