BlackHolePerturbationToolkit / Teukolsky

A Mathematica package for computing solutions to the Teukolsky equation.
http://bhptoolkit.org/Teukolsky/
MIT License
19 stars 3 forks source link

Static modes not working #26

Closed barrywardell closed 1 year ago

barrywardell commented 1 year ago

Since the recent changes to support more general sources the solutions for the static modes are no longer working. I have traced this down to the fact that the Wronskian is now computed using the "In" incidence coefficient, but we do not currently compute asymptotic amplitudes for the static modes. We could certainly still compute the Wronskian this way if we knew the amplitudes of the solutions. Is the asymptotic behaviour of the static modes know in closed form? Should we still call them "Incidence", "Reflection" and/or "Transmission" amplitudes in this case?

znasipak commented 1 year ago

Static modes are a bit trickier. Unlike the non-static modes, there is no standard convention for normalizing the static solutions, which then leads to ambiguity in the asymptotic amplitudes. You can write everything down in closed form, since the solutions are just hypergeometric functions. I usually just evaluate the Wronskian using the Teukolsky solutions and their derivatives. Then you don't have to worry about static vs non-static modes.

barrywardell commented 1 year ago

Getting back to this as it's the main issue blocking the next release:

There are also some more fundamental questions:

barrywardell commented 1 year ago

PS. We could do as Zach says and just directly evaluate the Wronskian, but that would introduce 4 extra evaluations of the radial function so it would be better to avoid it if possible.

Chriskavn commented 1 year ago

The following works for the Wronskian as far as I can tell:

-(2^(-l - s) E^(-(1/2) I \[Kappa] \[Tau] (1 + (2 Log[\[Kappa]])/( 1 + \[Kappa]))) \[Pi] \[Kappa]^(-l - s) (s + I \[Tau]) Csc[\[Pi] (s + I \[Tau])] Gamma[2 + 2 l] If[ s >= 1 && \[Tau] == 0, 1, Gamma[1 - s - I \[Tau]]])/(Gamma[ 1 + l - s] Gamma[1 + l - I \[Tau]] Gamma[1 - s - I \[Tau]] Gamma[ 1 + s + I \[Tau]])

Chriskavn commented 1 year ago

Answering you last point Barry, we don't have a notion of ingoing and outgoing solutions, but we do have linearly independent solutions which we have chosen. At the horizon, the two possible behaviours go like R-==(r-(1+\kappa))^{-s-I \tau/2} and R+== (r-(1+\kappa))^{I \tau/2}. Our choice of Rin gives Rin=Ain*R- at the horizon, for some Ain we can derive, and Rup=Bin R- + Bout R+ at the horizon, for some Bin and Bout. (Note our choice of Rin seems to be irregular for s>0, which is strange, am I missing something?). As such the wronskian is then given by 2^(1 + s) Ain Bout \kappa^(1 + s) (s + I \tau) .

znasipak commented 1 year ago

Rin isn't always regular but \Delta^s Rin should be. I'll double check to see how the static solutions are implemented at the moment. Yeah the notion of `"in" and "up" makes less sense in this case. We are just choosing the independent solutions that are regular and smooth at the horizon and infinity. The solutions are relatively simple, so it straightforward to derive the Wronskian for these solutions. I have those written down somewhere. I'll find them and post them soon. The only hiccup is that you have to separate the cases ma == 0 and ma != 0 and we have to decide what normalization convention we want to use. I propose we normalize them so that

\Delta \bar{R}^\pm_{slm0) = R^\pm_{-slm0}
znasipak commented 1 year ago

Here is my proposal for the case where $ma = 0$, which I expect to be the most common scenario for the static solutions. We define $$R^\mathrm{in}_{slm0} = (2\kappa)^{-s} \Gamma(1-s+1) [(-x)(1-x)]^{-s/2} P_l^s(1-2x),$$ and $$R^\mathrm{up}_{slm0} = 2 (-1)^{s+1} (2\kappa)^{-s} \Gamma^{-1}(1+s+1) [(-x)(1-x)]^{-s/2} Q_l^s(1-2x),$$ where $x = (1 + \kappa - r)/(2\kappa)$. Then the rescaled Wronskian is simply $$\Delta^{s+1} W[R^\mathrm{in}_{slm0}, R^\mathrm{up}_{slm0}] = 2\kappa.$$ Additionally you get the very simple relation between the opposite spin-weight solutions, $$\Delta^{s} R^\mathrm{in/up}_{slm0} = R^\mathrm{in/up}_{-slm0}.$$

When coding this up in Mathematica, you have to be careful to pick type 3 of LegendreP and LegendreQ.

znasipak commented 1 year ago

For the $ma \neq 0$ case, we can instead define $$R^\mathrm{in}_{slm0} = (2\kappa)^{-s} \Gamma(l+1-s) (-x)^{-s-i\tau/2} (1-x)^{-i\tau/2} \mathbf{F}(l+1-i\tau, -l-i\tau; 1-s-i\tau, x),$$ and $$R^\mathrm{up}_{slm0} = -(2\kappa)^{-s} \Gamma(l+1-i\tau) (-x)^{-l-1-s+i\tau/2} (1-x)^{-i\tau/2} \mathbf{F}(1+l-i\tau, l+1+s; 2l+2, x^{-1}),$$ where $\tau = -ma/\kappa$ and $\mathbf{F}$ denotes the regularized Gauss hypergeometric function, $$\mathbf{F}(a,b;c;z) = \frac{F(a,b;c;z)}{\Gamma(c)} .$$ Then we actually retain the same rescaled Wronskian as the $ma=0$ case, $$\Delta^{s+1} W[R^\mathrm{in}_{slm0}, R^\mathrm{up}_{slm0}] = 2\kappa.$$ This formulation is also nice, because in the limit $ma \rightarrow 0$, $R^\mathrm{in/up}_{slm0}$ solutions agree with those above. So we could just use these solutions for all $ma$ cases.

barrywardell commented 1 year ago

How does this compare to what's currently implemented?

barrywardell commented 1 year ago

The current static solutions are described in issue #10.

znasipak commented 1 year ago

They currently differ in their normalization. Excluding all of the pieces that depend on $x$ in the above expressions, their amplitudes are $$A^\mathrm{in} = (2\kappa)^{-s} \Gamma(l+1-s), \qquad A^\mathrm{up} = - (2\kappa)^{-s} \Gamma(l+1-i\tau).$$ The current implementation has $$A^\mathrm{in} = (2\kappa)^{-2s} e^{-i\kappa \tau/2(1 + \log\kappa-\log(1+\kappa)}, \qquad A^\mathrm{up} = - (2\kappa)^{-s-l-1}\Gamma(2l+2),$$ except when $s < 1$ or $\tau\neq 0$, we have the slight modification $$A^\mathrm{in} = (2\kappa)^{-2s} e^{-i\kappa \tau/2(1 + \log\kappa-\log(1+\kappa)}\Gamma(1-s-i\tau).$$ I think we currently have these normalizations so that the static solutions somewhat match the non-static solutions in the limit $\omega \rightarrow 0$. There are some regions of parameter space where this kind of works. However, no matter how small you make $\omega$ there will always be some $r\gg 1$ where you are definitively in the wave-zone and the static mode is going to look fundamentally different to the non-static mode. So I'm not sure if this way of normalizing the solutions is quite the right result. I'm leaning towards dealing with static solutions as distinctly different and normalizing them in a way that makes inversions and Wronskians simple for calculations later down the road. But I understand that might not be everyone's preference.

barrywardell commented 1 year ago

There's also at least one more difference: the in solution has a different power in the (-x): s-i tau/2 instead of -s-i tau/2. Is this a possible typo in Zach's expression above?

znasipak commented 1 year ago

That's a typo in my expression above. I've edited it to the correct expression.

barrywardell commented 1 year ago

The up solutions also look quite different, at least on the face of it. Are they actually the same up to the normalisation?

znasipak commented 1 year ago

Not really. The up solutions are difficult to match because the differences between the static and non-static modes are worst at infinity, since it is a regular singular point for one equation and irregular singularity for the other. Its a bad place to compare the modes. If you want things to look nice in the limit $\omega \rightarrow 0$, then you would want to do it near the horizon, but obviously this is not straightforward for the up solutions. This is why the static up solutions have a very straightforward normalization of $$R^\mathrm{up}(r\rightarrow\infty) \sim r^{-l-s-1}$$

barrywardell commented 1 year ago

Sorry, I wasn't clear. What I meant is that the up solution above appears to differ from the one currently implemented. Plotting the two suggest it's just an overall factor, although I'm not sure if one or the other does better in corner cases?

barrywardell commented 1 year ago

In terms of normalisation, we previously agreed a general principle that the homogeneous solutions should be normalised to have a value of 1 for the transmission coefficients. The static modes are fundamentally different, but I think we could do something analogous by normalising such that the in solution is $R^{\text{in}}=A^{\text{in}} (r-(1+\kappa))^{-s-I \tau/2}$ with $A^{\text{in}} = 1$, and similarly for the up solution $R^{\text{up}}=A^{\text{up}} r^{-l-s-1}$ with $A^{\text{up}} = 1$.

znasipak commented 1 year ago

Oh I see. It is just a different permutation of the hypergeometric function based on the Kummer relations, $(-x)^{-l-1-s+i\tau/2} (1-x)^{-i\tau/2} \mathbf{F}(1+l-i\tau, l+1+s; 2l+2, x^{-1})$ $= (-x)^{-s-i\tau/2} (1-x)^{-l-1+i\tau/2} \mathbf{F}(1+l-i\tau, l+1-s; 2l+2, (1-x)^{-1}) $

barrywardell commented 1 year ago

Ok, so it seems we can keep the current implementation but should consider picking a better normalisation.

znasipak commented 1 year ago

Then we want for all $ma$ and $s$, $$R^\mathrm{up}_{slm0} = -(2\kappa)^{-s-l-1} (-x)^{-l-1-s+i\tau/2} (1-x)^{-i\tau/2} {F}(1+l-i\tau, l+1+s; 2l+2, x^{-1}).$$ For $ma \neq 0$ or $s<1$ $$R^\mathrm{in}_{slm0} = (2\kappa)^{-s-i\tau/2} \Gamma(1-s-i\tau) (-x)^{-s-i\tau/2} (1-x)^{-i\tau/2} \mathbf{F}(l+1-i\tau, -l-i\tau; 1-s-i\tau, x)$$ or $$R^\mathrm{in}_{slm0} = (2\kappa)^{-s-i\tau/2}(-x)^{-s-i\tau/2} (1-x)^{-i\tau/2} {F}(l+1-i\tau, -l-i\tau; 1-s-i\tau, x).$$ The solutions with $ma = 0$ and $s \geq 1$ have a different asymptotic behavior of $$R^\mathrm{in}_{slm0} = A^\mathrm{in} (r-(1+\kappa))^{-i\tau/2}.$$

barrywardell commented 1 year ago

This has now been fixed in pull request #37 and some additional commits immediately afterwards.