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.

3) Using the magical package in one of the answers in here:

   "Site" -> "http://raw.githubusercontent.com/paclets/Repository/master"]; 
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.


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.

  • 2
    $\begingroup$ 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? $\endgroup$ – KraZug Jul 2 '19 at 16:04
  • $\begingroup$ 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. $\endgroup$ – yarchik Jul 2 '19 at 16:15
  • $\begingroup$ You have different signs in the equation $L\psi (r)=-f...$ and in the code \[ScriptCapitalL] = f*... $\endgroup$ – Alex Trounev Jul 2 '19 at 16:48
  • $\begingroup$ @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. $\endgroup$ – Edison Cesar Jul 3 '19 at 12:54
  • $\begingroup$ @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. $\endgroup$ – Edison Cesar Jul 3 '19 at 12:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.