# Finding eigenvalues of a differential operator

I am trying to get the eigenvalues of the following differential operator

$$L\psi(r) = -f \partial_r (f \partial_r \psi(r)) + V \psi(r)$$

which must satisfy (obviously)

$$L \psi(r) = \omega^2 \psi(r)$$

where I want to acquire both the real and imaginary part of $$\omega$$. To my problem, we have

$$f = 1 - \frac{2M}{r}$$ and $$V = f \left( \frac{l (l-1)}{r^2} + \frac{2 (1-S^2) M}{r^3} \right)$$ with $$M = 1, l = 2, S=2$$. The boundary conditions are $$\psi(2M) = 0, \psi(\inf) = 0$$. (To whomever may care, I am getting the real and imaginary oscillations of a Schwarzschild black hole. It is well studied in the literature, but I need to recover the result).

I tried three different ways to do it:

1) Using NDEigensystem:

f = 1 - 2*(M/r);
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3);
M = 1; l = 2; S = 2;
\[ScriptCapitalL] = f*D[f*D[\[Chi][r], r], r] + V*\[Chi][r];
\[ScriptCapitalB] = DirichletCo - ndition[\[Chi][r] == 0, True];
boundarydistance = 10;
{ev, ef} =
NDEigensystem[{\[ScriptCapitalL], \[ScriptCapitalB]}, \[Chi][
r], {r, 2*M, boundarydistance}, 3];
Print["The eigenvalues are = ", ev]
Print["The real part is = ", Re[Sqrt[ev]]]
Print["The imaginary part is = ", Im[Sqrt[ev]]]
Plot[Evaluate[ef], {r, 0, boundarydistance}]


THE PROBLEM:It works perfectly, for a condition which is not at infinity. When I try to set boundarydistance = :inf:, it gives me error.

2) NDSolve with "Shooting Method":

f = 1 - 2*(M/r);
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3);
M = 1; l = 2; S = 2;
\[ScriptCapitalL] = (-f)*D[f*D[\[Chi][r], r], r] + V*\[Chi][r];
NDSolve[{\[ScriptCapitalL] == \[Lambda]*\[Chi][r], \[Chi][2*M] == 0, \[Chi][Infinity] == 0}, \[Chi][r], {r, 0, 10},
Method -> {"Shooting"}]
Plot[%, {r, 0, 10}]


THE PROBLEM: It gives the following error: NDSolve::ndsv: Cannot find starting value for the variable [Chi]^[Prime]. Besides this, how to I recover the eigenvalue? It seems like, to use NDSolve, I need to put a value for it.

Needs["PacletManager"]
PacletInstall["CompoundMatrixMethod",
"Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"];
Needs["CompoundMatrixMethod"]
f = 1 - 2*(M/r);
V = f*(l*((l - 1)/r^2) + (2*(1 - S^2)*M)/r^3);
M = 1; l = 2; S = 2;
sys = ToMatrixSystem[f*D[f*D[\[Chi][r], r], r] - V*\[Chi][r] == \[Omega]*\[Chi][r],
{\[Chi][3*M] == 0, \[Chi][Infinity] == 0}, \[Chi], {r, 3, 10}, \[Omega]];
Plot[Evans[\[Omega], sys], {\[Omega], 0, 15}]


THE PROBLEM: Well, it actually runs. But I need to start from somewhere far from zero (due to the singularity) and I cannot really understand the answer which is:

{{{0, 1}, {(12 - 10 r + 2 r^2 + \[FormalLambda] r^4)/((-2 + r)^2 r^2), (4 r - 2 r^2)/((-2 + r)^2 r^2)}}, {}, {}, {r, 2.5, 10}}


and a plot, which does not resembles me anything that I would like to get.

EDIT:

Since I am not converging to an answer, it may help to add some more details. The real equation that I am trying to solve is $$\frac{d^2 \psi(r)}{dr_*^2} + (\omega^2 - V(r)) \psi(r) =0$$ where the (tortoise) coordinate is defined via $$\frac{dr}{dr_*} = f(r)$$ with the definitions of both $$V(r)$$ and $$f(r)$$ are given above. If you use these two equations together, you will get in the first equation of this page.

I thought that working with a single differential equation, although a bit more complicate, might be easier. However, I have no idea, even so, how I would work with the two differential equations separately, since the problem with this tortoise coordinate is that it cannot give me a function $$r(r_*)$$, because it is a transcendental equation.

• For my package (your method 3), the output you have quoted is the result of putting the equation into matrix form. I have a newer version in development which is much better at coping with the singularity. Do you know what the eigenvalue should be? – KraZug Jul 2 at 16:04
• A proper way to deal with this type of problems is the exterior complex scaling method. For applications in atomic and molecular physics have a look at this review, J. Phys. B 37 R137 (2004), and for applications in general relativity in this one Mathematical Research Letters 4 (1997): 103-122. Quite important, under complex scaling the resonance eigenfunction can be rendered square integrable. – yarchik Jul 2 at 16:15
• You have different signs in the equation $L\psi (r)=-f...$ and in the code \[ScriptCapitalL] = f*... – Alex Trounev Jul 2 at 16:48
• @KraZug Yes, I do know what the eigenvalues should be. However, I gave some more thinking, and this might give a different value due to the coordinate system I am using. Even so, that is not enough to explain my problemas. – Edison Cesar Jul 3 at 12:54
• @AlexTrounev It may have a different sign because I started to try several different things. However, even so, this was not what really made anything change. – Edison Cesar Jul 3 at 12:55