# NDEigensystem with FEM: parasitic solutions

I am trying to find the eigensystem of a system of 3 coupled ODEs. Analytically, the system spectrum should have a gap at [-2M,2M] (except for the degenerate state with E=0, strictly at the middle of the gap, which is of no interest). However, NDEigensystem with Finite Elements Method gives unwanted solutions inside the gap. Why is it so? How can I tune the procedure to get rid of the parasitic solutions? The MWE is below. Thank you!

Needs["NDSolveFEM"];
M = 0.01;
B = 50;
A = 5;
R = 100;

{vals, funs} = NDEigensystem[({
{-2 B u''[x] - (2 B)/x u'[x] + (2 B/x^2 + 2 M) u[x] -
2 A ( v'[x] + 2/x v[x])},
{A ( u'[x] - 1/x u[x]) - A ( z'[x] + 3/x z[x])},
{2 A ( v'[x] - 2/x v[x]) +
2 B z''[x] + (2 B)/x z'[x] + (-18 B/x^2 - 2 M) z[x]}
}), {u[x], v[x], z[x]}, {x, 0, R}, 1,
Method -> {"SpatialDiscretization" -> {"FiniteElement", \
{"MeshOptions" -> {"MaxCellMeasure" -> 0.001}}},
"Eigensystem" -> {"Arnoldi", "BasisSize" -> 100,
"MaxIterations" -> 10000, "Shift" -> -M}}];

vals
{-0.00974081}


upd: I suspect the parasitic states arise because we are trying to solve the problem within a finite region, mesh being also of finite coarseness. If I take a finer mesh, the spatial region of $$x \in (0,d)$$ where parasitic solutions do exist shrinks, i.e. $$d$$ decreases, but never disappears completely. It is vital for me however to get rid of these solution in the spectrum.

upd2. I tried to use the graded mesh,

mesh = ToGradedMesh[{Line[{{0}, {R}}], <|"Alignment" -> "Left",
"MinimalDistance" -> 0.00001, "ElementCount" -> 1000|>}];


But the problem persists.

upd3. The numerical scheme is stabilized by detuning the band dispersion slightly, in addition to the graded mesh. This gives rise to the appearance of v''[x] in the second equation. The unwanted solutions in the gap (half of the gap in this case) disappear.

• Could you give a link to the paper where this model analyzed with a statement about gap? May 8 at 3:30
• A few comments: Are you sure that {u[x], v[x], z[x]} is the right order? It matters. What are the right-hand sides of your ODEs? Do you really want to use the default NeumannValue boundary conditions, especially since v enters only as first order in derivatives? The ODEs are singular at x = 0. The eigenfunction solution for v is very noisy. May 8 at 15:24
• Which of the many sets of equations in the article you cited are you trying to solve? May 8 at 16:54
• I am glad that you have solved your problem. Please do look at the eigenfunctions, though, especially near x = 0 to see that they are well-behaved. You also should perform a series expansion of the ODEs at x = 0, which will uniquely determine the boundary conditions there. For the ODEs in your questions, the lowest order nonzero terms are u'[0] and z'''[0]. May 9 at 13:57

The symmetric model considered in the post is numerically unstable. The numerical scheme is stabilized by, physically speaking, detuning the dispersion of two electron bands: $$\pm M \pm B q^2 \to \pm M \pm B_i q^2$$, $$B_1 \ne B_2$$. Which is by the way much more relevant to experiment. This leads to the appearance of the second derivative $$(B_1-B_2)v''(x)$$ in the second equation. This stabilizes the solution, so that the unwanted energy levels disappear from the gap. In an asymmetric model the solutions do arise in half of the gap (at $$[0,2M]$$ when $$B_1>B_2$$), but at $$(-2M,0)$$ there should be a gap in the spectrum. And there is!

ClearAll["Global*"];
Needs["NDSolveFEM"];
M=0.01;
B1=120;
B2=20;
A=5;
R=100;
`