TeukolskyRadial Functionality for s=1/2

[maorby]

The TeukolskyRadial function currently does not support s=1/2.
Is there a way to enable this functionality?
How can I implement this myself?
It appears that the MST method is suitable for s=1/2.

[barrywardell]

I would expect most things will just work with s=1/2 since there is no inherent assumption on s in most of the equations used. It's also likely that the HeunC method will work for s=1/2.

We had a similar request for the SpinWeightedSpheroidalHarmonics package (see issue 22 for the implementation and issue 32 for some subtleties that arose). In that case, it was just a matter of removing the check that s is an integer, so I suggest starting with that and checking if everything just works.

[maorby]

I've reviewed the relevant packages and found no restrictions indicating that s must be an integer, even though it's consistently labeled as s_Integer. Also I've been searching for potential issues that could cause the code to fail without reporting errors but found none.

Is there anything else I should check?

When I run TeukolskyRadial with s=1/2, it seems like the function executes correctly, but it doesn't know how to export the results.

[barrywardell]

The _Integer part is what is restricting s to be an integer. The first step is to remove the integer restriction (leaving just s_) and check if it works.

[maorby]

I removed the _Integer part from the s_, l_ and m_ in the following files: TeukolskyRadial.m, NumericalIntegration.m, MST.m, RenormalizedAngularMomentum.m, and it works well.
How can I ensure that the results obtained by the package are correct?

[barrywardell]

Great. A good check would be to see if the solutions produced satisfy the Teukolsky equation. If they do and if they also satisfy the expected boundary conditions (see the notebook in Tests/Correctness in the repository) then I think we can be confident that they are correct.

[maorby]

I’ve completed the tests that you suggested. It seems that the solutions indeed satisfy the Teukolsky equation and the expected boundary conditions.

I also tried to do the test the solutions for the case of static modes but the test failed for some cases, and I did not find any reference article that the test file relies on, so I'm not sure how to make the necessary adjustments to the case of half-integer spin.

[barrywardell]

I'm glad to hear that the non-static case is working. We should enable that in the next release of the package.

In the static case we don't have the usual same idea of "in" and "up" solutions as we do in the non-static case. This is true also for integer spin and isn't an issue that is specific to the half-integer case. We do, however, still have a second order ODE which admits two linearly independent homogeneous solutions and we just need to make a choice of which basis of solutions to use. The convention that we decided on is to have solutions determined by asymptotic behaviour similar to the non-static case, in particular towards infinity the asymptotic behaviours are $r^{-l-s-1}$ and $r^{l-s}$ while towards the horizon they are $\Delta^{i \tau/2}$ and $\Delta^{-i \tau/2-s}$ where $\tau = - a m /\sqrt{1 - a^2}$. We define "in" and "up" solutions similar to the non-static case, by choosing the "up" solution such that the coefficients towards infinity are 0 and 1 and similarly the "in" solutions such that the solutions towards the horizon are 0 and 1. This is tested in the notebook Tests/Correctness/StaticAmplitudes.nb and there is further detail in issue #26.

[maorby]

I've completed the tests for the static mode as well. However, for some cases in the test of the IN solution near the horizon (with parameters a = 0, τ = 0), I had to multiply the IN solution by a constant factor (2, -2, -8, e, depending on the mode) in order to satisfy the boundary condition. Do you have any idea where these factors might come from?

Another problem I've encountered is when testing the s = 5/2 case: the IN solution appears to be non-continuous at infinity (specifically around r = 5000). Up to that point, the solution satisfies the boundary condition, but beyond this point, it does not. Any insights would be appreciated!

[maorby]

The problem I’m describing here does not directly relate to the main issue, but it highlights an issue with the solution to the Teukolsky equation for a neutrino field (spin s=1/2) that I obtained using the TeukolskyRadial function.

I've encountered a problem while calculating the greybody factors for a black hole in the case of a neutrino field (s = 1/2).

To compute these factors, I solve the Teukolsky equation using the Teukolsky package for each set of parameters (s, l, m, a, and ω) and then use the reflection and incidence amplitudes in the following expression:

$$\Gamma=1−\Bigg|\frac{TeukolskyRadial[s,l,m,a,\omega ][Reflection]\cdot TeukolskyRadial[−s,l,m,a,\omega ][Reflection]}{TeukolskyRadial[s,l,m,a,\omega ][Incidence]\cdot TeukolskyRadial[−s,l,m,a,\omega ][Incidence]}\Bigg|$$

For electromagnetic and gravitational fields (s = 1, 2), the calculation of these factors matches the low-frequency approximation (as derived by Don Page, 1975) and approaches 1 at high frequencies, as expected.

However, for the neutrino field (s = 1/2), the results do not align with the low-frequency approximation from the same article. Additionally, I’m getting negative values for the greybody factors, which should remain positive across all frequencies to ensure a positive emission spectrum. It's worth mentioning that at high frequencies, these factors still approach 1 as expected.

All the changes I made in the Teukolsky package amount to removing the _Integer restriction from s, l, and m in the following files:

TeukolskyRadial.m
NumericalIntegration.m
MST.m
RenormalizedAngularMomentum.m

Any insights would be great!