Numerical investigation of the dynamics of linear spin $s$ fields on Kerr background I. Late time tails of spin $s = \pm 1, \pm 2$ fields

TL;DR

This paper numerically studies the late-time decay of linear spin $s=\pm 1,\, \pm 2$ fields on Kerr black holes, analyzing their behavior at different characteristic locations using conformal compactification and advanced numerical methods.

Contribution

It introduces a robust numerical framework combining conformal compactification and spectral methods to analyze spin fields on Kerr backgrounds, including new decay rate measurements at multiple locations.

Findings

Determined decay rates of spin fields at horizon, infinity, and in the domain of outer communication.

Validated the numerical approach using energy and angular momentum balance relations.

Confirmed the robustness and accuracy of the method for late-time tail analysis.

Abstract

The time evolution of linear fields of spin $s = \pm 1$ and $s = \pm 2$ on Kerr black hole spacetimes are investigated by solving the homogeneous Teukolsky equation numerically. The applied numerical setup is based on a combination of conformal compactification and the hyperbolic initial value problem. The evolved basic variables are expanded in terms of spin-weighted spherical harmonics which allows us to evaluate all the angular derivatives analytically, whereas the evolution of the expansion coefficients, in the time-radial section, is determined by applying the method of lines implemented in a fourth order accurate finite differencing stencil. Concerning the initialization, in all of our investigations single mode excitations---either static or purely dynamical type initial data---are applied. Within this setup the late time tail behavior is investigated. Due to the applied conformal compactification the asymptotic decay rates are determined at three characteristic locations---in the domain of outer communication, at the event horizon and at future null infinity---simultaneously. Recently introduced new type of `energy' and `angular momentum' balance relations are also applied in order to demonstrate the feasibility and robustness of the developed numerical schema, and also to verify the proper implementation of the underlying mathematical model.

Tables

Table 1. (a) s = − 1 𝑠 1 s=-1 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|8
2	5|4	6|5	7|6	8|(7)	9|?
3	6|5	7|5	8|(5)	9|?	10|?
4	7|6	8|6	9|6	10|?	11|?
5	8|7	9|7	10|7	?|?	?|?

Table 2. (a) s = − 1 𝑠 1 s=-1 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|8
2	5|4	6|5	7|6	8|(7)	9|?
3	6|5	7|5	8|(5)	9|?	10|?
4	7|6	8|6	9|6	10|?	11|?
5	8|7	9|7	10|7	?|?	?|?

Table 3. (b) s = 1 𝑠 1 s=1 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	6|5|2	7|6|3	8| 7|4	9| 8|5	10| 9 | 6
2	6|5|2	7|6|3	8| 7|4	9| 8|5	10| 9 | 6
3	7|6|3	8|7|3	9| 8|4	10| 9|5	11| 10 |(7)
4	8|7|4	9|8|4	10| 9|4	11|10|?	12| 11 |(7)
5	9|8|5	10|9|5	11|10|5	12|11|5	13|(12)|(7)

Table 4. (c) s = − 1 𝑠 1 s=-1 , dynamic initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7| 6	8 | 7	9|(9)
2	6|5	7|5	8| 6	9 | 7	10|(8)
3	5|4	6|5	7| 6	8 | 7	9| ?
4	6|5	7|5	8|(5)	9 |(7)	10| ?
5	7|6	8|6	9| 6	(10)| ?	?| ?

Table 5. (d) s = 1 𝑠 1 s=1 , dynamic initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	6|5|2	7|6|3	8|7|4	9| 8| 5	10| 9 |6
2	7|6|3	8|7|3	9|8|4	10| 9| 5	11| 10 |6
3	6|5|2	7|6|3	8|7|4	9| 8| 5	10| (8)|6
4	7|6|3	8|7|3	9|8|4	10| 9| 5	11| 10 |6
5	8|7|4	9|8|4	10|9|4	11|10|(4)	12|(10)|6

Table 6. (a) s = − 2 𝑠 2 s=-2 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|?	9|?	10|?
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 7. (a) s = − 2 𝑠 2 s=-2 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|?	9|?	10|?
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 8. (b) s = 2 𝑠 2 s=2 , static initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	8|7|2	9| 8|3	10| 9|4	11|10|5
3	X	8|7|2	9| 8|3	10| 9|4	11|10|5
4	X	9|8|3	10| 9|3	11|10|4	12|11|5
5	X	10|9|4	11|10|4	12|11|4	13|12|5

Table 9. (c) s = − 2 𝑠 2 s=-2 , dynamic initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	8|7	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 10. (d) s = 2 𝑠 2 s=2 , dynamic initial data.

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	8|7|2	9|8|3	10| 9|4	11|10|5
3	X	9|8|3	10|9|3	11|10|4	12|11|5
4	X	8|7|2	9|8|?	10| 9|4	11|10|5
5	X	9|8|3	10|9|3	11|10|4	12|11|5

Table 11. (a) s = − 1 𝑠 1 s=-1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|?
2	5|4	6|5	7|6	8|(7)	9|?
3	6|5	7|5	8|6	9|(6)	10|?
4	7|6	8|6	9|6	10|?	11|?
5	8|7	9|7	10|7	?|?	?|?

Table 12. (a) s = − 1 𝑠 1 s=-1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|?
2	5|4	6|5	7|6	8|(7)	9|?
3	6|5	7|5	8|6	9|(6)	10|?
4	7|6	8|6	9|6	10|?	11|?
5	8|7	9|7	10|7	?|?	?|?

Table 13. (b) s = 1 𝑠 1 s=1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|5|2	6|6|3	7| 7|4	8| 8|5	9| 9 | 6
2	5|5|2	6|6|3	7| 7|4	8| 8|5	9| 9 | 6
3	6|6|3	7|7|3	8| 8|4	9| 9|5	10| 10 |6
4	7|7|4	8|8|4	9| 9|4	10|10|5	11| 11 |(6)
5	8|8|5	9|9|5	10|10|5	11|11|5	12|(12)|(6)

Table 14. (c) s = − 1 𝑠 1 s=-1 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|?
2	5|4	6|5	7|6	8|?	9|?
3	6|5	7|5	8|6	9|?	10|?
4	7|6	8|6	9|6	10|?	11|?
5	8|7	9|7	10|7	?|?	?|?

Table 15. (d) s = 1 𝑠 1 s=1 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	5|5|2	6|6|3	7|7|4	8|8|5	9|9|6
2	5|5|2	6|6|3	7|7|4	8|8|5	9|9|6
3	6|6|3	7|7|3	8|8|4	9|9|5	10|10|6
4	7|7|4	8|8|4	9|9|4	10|10|5	11|11|(6)
5	8|8|5	9|9|5	10|10|5	11|11|5	12|(12)|(6)

Table 16. (a) s = − 1 𝑠 1 s=-1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	8|6	9|(7)	10|?
4	X	8|6	9|6	10|?	10|?
5	X	9|7	10|7	?|?	?|?

Table 17. (a) s = − 1 𝑠 1 s=-1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	8|6	9|(7)	10|?
4	X	8|6	9|6	10|?	10|?
5	X	9|7	10|7	?|?	?|?

Table 18. (b) s = 1 𝑠 1 s=1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|3	8| 8|4	9| 9|5	10| 10 | 6
3	X	7|7|3	8| 8|4	9| 9|5	10| 10 | 6
4	X	8|8|4	9| 9|4	10|10|5	11| 11 |6
5	X	9|9|5	10|10|5	11|11|5	12|(12)|6

Table 19. (c) s = − 1 𝑠 1 s=-1 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	8|6	9|(7)	10|?
4	X	8|6	9|6	10|?	11|?
5	X	9|7	10|7	?|?	?|?

Table 20. (d) s = 1 𝑠 1 s=1 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|3	8|8|4	9|9|5	10|10|6
3	X	7|7|3	8|8|4	9|9|5	10|10|6
4	X	8|8|4	9|9|4	10|10|5	11|11|6
5	X	9|9|5	10|10|5	11|11|5	12|(12)|6

Table 21. (a) s = − 2 𝑠 2 s=-2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|?
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 22. (a) s = − 2 𝑠 2 s=-2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|?
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 23. (b) s = 2 𝑠 2 s=2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8| 8|3	9| 9|4	10|10|5
3	X	7|7|2	8| 8|3	9| 9|4	10|10|5
4	X	8|8|3	9| 9|3	10|10|4	11|11|5
5	X	9|9|4	10|10|4	11|11|4	12|12|5

Table 24. (c) s = − 2 𝑠 2 s=-2 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|?
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 25. (d) s = 2 𝑠 2 s=2 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8| 8|3	9| 9|4	10|10|5
3	X	7|7|2	8| 8|3	9| 9|4	10|10|5
4	X	8|8|3	9| 9|3	10|10|4	11|11|5
5	X	9|9|4	10|10|4	11|11|4	12|12|5

Table 26. (a) s = − 2 𝑠 2 s=-2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 27. (a) s = − 2 𝑠 2 s=-2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 28. (b) s = 2 𝑠 2 s=2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8|8|3	9|9|4	10|10|5
3	X	7|7|2	8|8|3	9|9|4	10|10|5
4	X	8|8|3	9|9|3	10|10|4	11|11|5
5	X	9|9|4	10|10|4	11|11|4	12|12|5

Table 29. (c) s = − 2 𝑠 2 s=-2 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	8|?	9|?	10|?
4	X	8|7	?|?	?|?	?|?
5	X	9|8	?|?	?|?	?|?

Table 30. (d) s = 2 𝑠 2 s=2 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8|8|3	9|9|4	10|10|5
3	X	7|7|2	8|8|3	9|9|4	10|10|5
4	X	8|8|3	9|9|3	10|10|4	11|11|5
5	X	9|9|4	10|10|4	11|11|4	12|12|5

Table 31. (a) s = − 1 𝑠 1 s=-1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|(8)
2	5|4	6|5	7|6	8|7	?|(8)
3	5|4	6|5	7|6	8|(7)	9|?
4	6|5	7|5	8|6	9|(7)	10|?
5	7|6	8|6	9|6	(9)|?	?|?

Table 32. (a) s = − 1 𝑠 1 s=-1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|(8)
2	5|4	6|5	7|6	8|7	?|(8)
3	5|4	6|5	7|6	8|(7)	9|?
4	6|5	7|5	8|6	9|(7)	10|?
5	7|6	8|6	9|6	(9)|?	?|?

Table 33. (b) s = 1 𝑠 1 s=1 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	5|5|2	6|6|3	7| 7|4	8| 8|5	9| 9 | 6
2	5|5|2	6|7|3	7|8|4	8|9|5	9|10| 6
3	5|5|2	6|6|3	7|(7)|4	8| ?|5	9| ? |6
4	6|6|3	7|7|3	8| 8|4	9|9|5	10| 10 |6
5	7|7|4	8|8|4	9| 9|4	10|10|(5)	11| 11 |(6)

Table 34. (c) s = − 1 𝑠 1 s=-1 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	5|4	6|5	7|6	8|7	9|(8)
2	5|4	6|5	7|6	8|7	9|(8)
3	5|4	6|5	7|6	8|(7)	9|?
4	6|5	7|5	8|6	9|?	10|?
5	7|6	8|6	9|6	10|?	?|?

Table 35. (d) s = 1 𝑠 1 s=1 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	5|5|2	6|6|3	7|7|4	8|8|5	9|9|6
2	5|5|2	6|7|3	7|8|4	8|9|5	9|10|6
3	5|5|2	6|6|3	7|?|4	8|?|5	9|?|6
4	6|6|3	7|7|3	8|8|4	9|9|5	10|10|6
5	7|7|4	8|8|4	9|9|4	10|10|5	11|11|6

Table 36. (a) s = − 1 𝑠 1 s=-1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	9|6	10|7	11|?
4	X	7|5	8|6	9|?	10|?
5	X	8|6	9|6	10|?	11|?

Table 37. (a) s = − 1 𝑠 1 s=-1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	9|6	10|7	11|?
4	X	7|5	8|6	9|?	10|?
5	X	8|6	9|6	10|?	11|?

Table 38. (b) s = 1 𝑠 1 s=1 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|3	8| 8|4	9| 9|5	10| 10 | 6
3	X	7|7|3	8| 9|4	9|10|5	10|11 | 6
4	X	7|7|3	8| 8|4	9| ?|5	10| ? |6
5	X	8|8|4	9| 9|4	10|10|5	11|(11)|6

Table 39. (c) s = − 1 𝑠 1 s=-1 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|5	8|6	9|7	10|?
3	X	7|5	?|6	?|7	?|?
4	X	7|5	8|6	9|?	10|?
5	X	8|6	9|6	10|?	11|?

Table 40. (d) s = 1 𝑠 1 s=1 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|3	8|8|4	9|9|5	10|10|6
3	X	7|7|3	8|9|4	9|10|5	10|11|6
4	X	7|7|3	8|8|4	9| ?|5	10| ?|6
5	X	8|8|4	9|9|4	10|10|5	11|?|6

Table 41. (a) s = − 2 𝑠 2 s=-2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	?|?	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 42. (a) s = − 2 𝑠 2 s=-2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	?|?	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 43. (b) s = 2 𝑠 2 s=2 , m = 1 𝑚 1 m=1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8| 8|3	9| 9|4	10|10|5
3	X	7|7|2	8| 9|3	9|10|4	10|11|5
4	X	7|7|2	8| 8|3	9| 9|4	10|10|5
5	X	8|8|3	9| 9|3	10|10|4	11|11|5

Table 44. (c) s = − 2 𝑠 2 s=-2 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	?|?	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 45. (d) s = 2 𝑠 2 s=2 , m = − 1 𝑚 1 m=-1

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8| 8|3	9| 9|4	10|10|5
3	X	7|7|2	8| 9|3	9|10|4	10|11|5
4	X	7|7|2	8| 8|3	9| 9|4	10|10|5
5	X	8|8|3	9| 9|3	10|10|4	11|11|5

Table 46. (a) s = − 2 𝑠 2 s=-2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 47. (a) s = − 2 𝑠 2 s=-2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 48. (b) s = 2 𝑠 2 s=2 , m = 2 𝑚 2 m=2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8|8|3	9|9|4	10|10|5
3	X	7|7|2	8|9|3	9|10|4	10|11|5
4	X	7|7|2	8|8|3	9|9|4	10|10|5
5	X	8|8|3	9|9|3	10|10|4	11|11|5

Table 49. (c) s = − 2 𝑠 2 s=-2 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|6	8|7	9|8	10|9
3	X	7|6	?|?	?|?	?|?
4	X	7|6	8|?	9|?	10|?
5	X	8|7	?|?	?|?	?|?

Table 50. (d) s = 2 𝑠 2 s=2 , m = − 2 𝑚 2 m=-2

$\ell$’	$\ell$=1	2	3	4	5
1	X	X	X	X	X
2	X	7|7|2	8|8|3	9|9|4	10|10|5
3	X	7|7|2	8|9|3	9|10|4	10|11|5
4	X	7|7|2	8|8|3	9|9|4	10|10|5
5	X	8|8|3	9|9|3	10|10|4	11|11|5

Equations

$${\varphi }_0=$$

$${\varphi }_1=$$

$${\varphi }_2=$$

$${\psi }^{(+1)}=$$

$${\psi }^{(-1)}=$$

$${\Psi }_2=-\frac{M}{(r-\mathrm{i}a\cos \vartheta {)}^3}.$$

$${\Psi }_0=-C_{abcd}{l}^a{m}^b{l}^c{m}^d\text{and}{\Psi }_4=-C_{abcd}{n}^a\overline{m}{}^b{n}^c\overline{m}{}^d$$

$${\psi }^{(+2)}=$$

$${\psi }^{(-2)}=$$

$$\begin{gathered} (\mathrm{d}s{)}^2=\left(1-\frac{2Mr}{\Sigma }\right)(\mathrm{d}t{)}^2+\frac{4arM}{\Sigma }{\sin }^2\vartheta \mathrm{d}t\mathrm{d}\varphi \\ -\frac{\Sigma }{\Delta }(\mathrm{d}r{)}^2-\Sigma (\mathrm{d}\vartheta {)}^2-\frac{(r^2+a^2{)}^2-a^2\Delta {\sin }^2\vartheta }{\Sigma }{\sin }^2\vartheta (\mathrm{d}\varphi {)}^2, \end{gathered}$$

$$\begin{gathered} \left[\frac{(r^2+a^2{)}^2}{\Delta }-a^2{\sin }^2\vartheta \right]\frac{{\partial }^2{\psi }^{(s)}}{\partial t^2}+\frac{4Mar}{\Delta }\frac{{\partial }^2{\psi }^{(s)}}{\partial t\partial \varphi }+\left[\frac{a^2}{\Delta }-\frac{1}{{\sin }^2\vartheta }\right]\frac{{\partial }^2{\psi }^{(s)}}{\partial {\varphi }^2}- \\ -{\Delta }^{-s}\frac{\partial }{\partial r}\left({\Delta }^{s+1}\frac{\partial {\psi }^{(s)}}{\partial r}\right)-\frac{1}{\sin \vartheta }\frac{\partial }{\partial \vartheta }\left(\sin \vartheta \frac{\partial {\psi }^{(s)}}{\partial \vartheta }\right)-2s\left[\frac{a(r-M)}{\Delta }+\frac{\mathrm{i}\cos \vartheta }{{\sin }^2\vartheta }\right]\frac{\partial {\psi }^{(s)}}{\partial \varphi }- \\ -2s\left[\frac{M(r^2-a^2)}{\Delta }-r-\mathrm{i}a\cos \vartheta \right]\frac{\partial {\psi }^{(s)}}{\partial t}+(s^2{\cot }^2\vartheta -s){\psi }^{(s)}=0. \end{gathered}$$

$$\left[\left({\nabla }^a+s{\Gamma }^a\right)\left({\nabla }_a+s{\Gamma }_a\right)-4s^2{\Psi }_2\right]{\psi }^{(s)}=0,$$

$${\Gamma }^t$$

$${\Gamma }^r$$

$${\Gamma }^{\vartheta }$$

$${\Gamma }^{\varphi }$$

$$\tau$$

$$\phi$$

$$\tau$$

$$r$$

$$\Phi^{(s)}(T,R,\vartheta,\varphi)=\big{[}\,r(R)\cdot\Delta^{s}(R)\,\big{]}\cdot\psi^{(s)}(T,R,\vartheta,\varphi)\,.$$

$$\displaystyle\partial_{TT}\Phi^{(s)}=\frac{1}{\mathscr{A}+\mathscr{B}\cdot Y_{2}^{0}}\,\Big{[}c_{RR}\cdot\partial_{RR}\Phi^{(s)}+c_{TR}\cdot\partial_{TR}\Phi^{(s)}+c_{T\varphi}\cdot\partial_{T\varphi}\Phi^{(s)}+c_{R\varphi}\cdot\partial_{R\varphi}\Phi^{(s)}$$

$$\displaystyle+c_{\vartheta\vartheta}\cdot\overline{\eth}\eth\,\Phi^{(s)}+c_{T}\cdot\partial_{T}\Phi^{(s)}+\mathrm{i}\,c_{Ty}\,Y_{1}^{0}\cdot\partial_{T}\Phi^{(s)}+c_{R}\cdot\partial_{R}\Phi^{(s)}+c_{\varphi}\cdot\partial_{\varphi}\Phi^{(s)}+c_{0}\cdot\Phi^{(s)}\Big{]}\,,$$

$$ðf$$

$$\overline{ð}f$$

$$\overline{ð}ðf={\partial }_{\vartheta \vartheta }f+\cot \vartheta {\partial }_{\vartheta }f+\frac{1}{{\sin }^2\vartheta }{\partial }_{\phi \phi }f+2\mathrm{i}s\frac{\cot \vartheta }{\sin \vartheta }{\partial }_{\phi }f+s(1-s{\cot }^2\vartheta )f.$$

$$\mathcal{L}=-({\nabla }^a-s{\Gamma }^a){\psi }^{(-s)}({\nabla }_a+s{\Gamma }_a){\psi }^{(s)}-4s^2{\Psi }_2{\psi }^{(s)}{\psi }^{(-s)}.$$

$$E^a$$

$$J^a$$

$${\int }_{\Omega }{\nabla }_a{E}^a={\int }_{\partial \Omega }{n}_a{E}^a$$

$${\int }_{\Omega }{\nabla }_a{J}^a={\int }_{\partial \Omega }{n}_a{J}^a$$

$$\begin{gathered} 0={\int }_{T=T_i}{n}_a^{(T)}{E}^a\sqrt{|h_T|}\mathrm{d}R\mathrm{d}\vartheta \mathrm{d}\phi +{\int }_{R=R_{in}}{n}_a^{(R)}{E}^a\sqrt{|h_R|}\mathrm{d}T\mathrm{d}\vartheta \mathrm{d}\phi \\ -{\int }_{T=T_f}{n}_a^{(T)}{E}^a\sqrt{|h_T|}\mathrm{d}R\mathrm{d}\vartheta \mathrm{d}\phi -{\int }_{R=R_{out}}{n}_a^{(R)}{E}^a\sqrt{|h_R|}\mathrm{d}T\mathrm{d}\vartheta \mathrm{d}\phi . \end{gathered}$$

$${\Phi }^{(s)}(T,R,\vartheta ,\phi )=\sum _{\ell =|s|}^{{\ell }_{max}}\sum _{m=-\ell }^{\ell }{\varphi }_{\ell }{}^m(T,R)\cdot {}_s{Y}_{\ell }{}^m(\vartheta ,\phi ).$$

$$(ST):\{\begin{array}{l}{\varphi }^{(s)}(R,\vartheta ,\phi )={\varphi }_{\ell }{}^m(T_0,R)\cdot {}_s{Y}_{\ell }{}^m(\vartheta ,\phi )\\ {\varphi }_T^{(s)}(R,\vartheta ,\phi )\equiv 0,\end{array}$$

$$(PD):\{\begin{array}{l}{\varphi }^{(s)}(R,\vartheta ,\phi )\equiv 0,\\ {\varphi }_T^{(s)}(R,\vartheta ,\phi )=({\varphi }_T{)}_{\ell }{}^m(T_0,R)\cdot {}_s{Y}_{\ell }{}^m(\vartheta ,\phi )\end{array}$$

$$\mathcal{B}(R)=\{\begin{array}{ll}\frac{2R}{1-R^2}\exp \left(-\frac{1}{R-c+w/2}+\frac{1}{R-c-w/2}+\frac{4}{w}\right),& \text{mbox}ifc-w/2<R<c+w/2,\\ \text{hfill}0,& \text{mbox}otherwise.\end{array}$$

$${\varphi }_{\ell }{}^m(T,R_0)\sim T^{-n},$$

Full text

1 Introduction

The study of the long-time evolution of various linear fields on given black hole backgrounds has served for decades as a preparation for the more involved study of linear and possibly nonlinear stability of the background black hole solutions themselves. In the last couple of years considerable progress has been made concerning the linear stability of the Kerr solution (see, e.g. [9, 1] and references therein). Notably, even these analytic results rest on various technical assumptions which apparently all boil down to the long-time behavior of linear spin $s$ fields satisfying the Teukolsky master equation [9, 1]. This provides immediate motivation for a thorough investigation of the asymptotic in time behavior of solutions to the Teukolsky equation on a Kerr black hole background.

This paper reports about our findings concerning numerical investigations of the dynamics of electromagnetic and gravitational perturbations on a Kerr background. Our aim was to carry out comprehensive investigations of the time evolution of linear spin $s$ fields, in particular, to study their tail behavior. This is done by solving the homogeneous Teukolsky master equation for generic linear spin $s=\pm 1$ and $s=\pm 2$ fields numerically. Since previous numerical studies focused mostly on axisymmetric configurations, we aimed to provide a more detailed study of nonaxisymmetric configurations. In doing so, the analytic framework was chosen such that it incorporates both techniques of conformal compactification and the hyperboloidal initial value problem. The time slices in the latter were chosen to be horizon penetrating, which allows us to determine the decay rates at the three characteristic locations simultaneously: at the black hole event horizon, in the domain of outer communication, and at future null infinity. In addition, the applied mathematical setup also makes it possible to use a spherical spectral representation of all the basic variables. This is so as they are expanded in terms of spin-weighted spherical harmonics based on the foliation of the Kerr background by topological two-spheres, which in practice are the Boyer-Lindquist $t=const$ and $r=const$ ellipsoids. As all angular derivatives can be given then either in terms of the “eth” and “ethbar”, $\eth$ and $\overline{\eth}$, operators or by Lie derivatives with respect to the axial symmetry of the Kerr background, all angular derivatives are evaluated analytically. This, in turn, guarantees that the pertinent multipole expansion coefficients—as elements of a large set of coupled scalar fields—get to be subject of an involved but otherwise only $(1+1)$-dimensional time-evolution problem. Notably, even the case of the most generic linear spin $s$ fields can be tested within the very same mathematical and numerical setup which, in turn, allows us to treat all the nonaxisymmetric configurations as well. Though in practice, our attention was restricted to the range $m=0,\pm 1,\pm 2$ of the azimuthal parameter, in principle, one could investigate the full spectra of nonaxisymmetric configurations within the chosen framework.

In checking and interpreting our numerical findings it was very informative and supportive to compare them to both analytic and numerical results prior to our investigations. On a Kerr background, the first systematic analytic estimates on the decay exponents were derived by Barack and Ori [2, 4]. They developed a method capable of studying the evolution of linear spin $s=0$, $s=\pm 1$, and $s=\pm 2$ in the time domain. By applying this they discovered an interesting phenomenon. Namely, for axially symmetric configurations, the decay rates at the event horizon are larger by $1$ for spin $s>0$ with respect to those relevant for $s<0$ [3]. By applying the same method, they could also study the behavior of linear perturbations while approaching the Cauchy or inner horizon of the Kerr background [17, 18]. The alternative analytic investigations by Hod were carried out in the frequency domain [12, 13, 14]. In [6, 7] Casals et al. developed a technique to further improve on the results of Hod. They presented a method to include higher order terms in the low-frequency expansion, which is necessary for computing self-force. Some complementary analytic studies can also be found in the Appendix of [11]. Note that in Ref. [11], Harms et al. also report impressive and comprehensive numerical investigations concerning the long-time evolution of axially symmetric linear spin $s$ fields on a Kerr background. It is important to emphasize that the initial data we apply are always of pure mode excitation, which allows us to study the long-time behavior of various higher mode excitations. It is of crucial importance then that the analytic and numerical investigations carried out in [14, 11] are also based on the use of single-mode-excitation-type initial data. In this respect, from among all the aforementioned excellent investigations, the ones reported in [14, 11] will be at the center of our interest. Most of the predictions made in these studies are confirmed by our investigations. Nevertheless, as there were some slight disagreements even between the predictions of these two sets of investigations, in our studies particular attention was given to clear up of the corresponding cases. Notably, our numerical results led us to the conclusion that, though only in some very special subcases, neither of the former predictions was satisfactory: i.e., they did not give the correct value for the decay exponents.

In coming up with a firm statement of the above type, it is extremely important to guarantee the self-consistency of the underlying numerical results. It is indeed critical to verify that our findings are not simply numerical artifacts of the applied method, but they are rooted in the true nature of the investigated fields. For this reason, we implemented—in addition to the conventional convergence rate checks—the recently introduced [22] new type of “energy” and “angular momentum” balance relations to verify both the proper implementation of the underlying mathematical model and the feasibility and robustness of the developed numerical schema.

Before proceeding in presenting our main results, it appears to be rewarding to recall that linear spin $s=\pm 1$ or $s=\pm 2$ fields satisfying the homogeneous Teukolsky master equations arise in a straightforward way in studying evolution of source-free electromagnetic fields in Maxwell theory or in that of linear metric perturbations on a Kerr background [20, 21]. To see this, e.g. in the electromagnetic case, recall first that in applying the Newman-Penrose formalism [16], one starts by fixing a complex null tetrad $\{l^{a},n^{a},m^{a},\overline{m}{}^{a}\}$ comprising at each point two real $l^{a}$ and $n^{a}$ and two complex $m^{a}$ and $\overline{m}{}^{a}$ null vectors such that their only nonzero inner products are $l^{a}n_{a}=-m^{a}\overline{m}_{a}=1$. We fix the remaining guage freedom by using the Kinnersley tetrad. The algebraically independent components of the Faraday tensor $F_{ab}$ can always be represented by the three complex Maxwell scalar fields

[TABLE]

Two of these, $\phi_{0}$ and $\phi_{2}$, represent the outgoing and ingoing radiations, respectively, whereas $\phi_{1}$ stores the Coulombic part of the electromagnetic field.

If $F_{ab}$ satisfies the source-free Maxwell equations on a Kerr black hole background, then some of the one-order-higher wave equations—which can be deduced for the individual Maxwell scalar components $\phi_{0}$ and $\phi_{2}$, respectively—decouple. In particular, the fields $\psi^{(+1)}$ and $\psi^{(-1)}$ defined as

[TABLE]

satisfy the homogeneous Teukolsky master equation (11) with spin $s=+1$ and $s=-1$, respectively [20]. Here, $\Psi_{2}$ is the only nonvanishing and gauge invariant Weyl-scalar component on a Kerr background which has the following form in Boyer-Lindquist coordinates using the Kinnersley tetrad:

[TABLE]

By a completely analogous argument, it can also be verified [20, 21] that in the case of linear metric perturbations, the specific contractions

[TABLE]

of the Weyl tensor $C_{abcd}$—which is also equal to the curvature tensor in the case of vacuum spacetimes—and elements of the complex null tetrad $\{l^{a},n^{a},m^{a},\overline{m}{}^{a}\}$ are distinguished like $\phi_{0}$ and $\phi_{2}$ were in the electromagnetic case. Indeed, the linear wave equations—these can be deduced from the Bianchi identity, and they are relevant for the linearly perturbed configurations—really do decouple, and the fields

[TABLE]

satisfy the homogeneous Teukolsky master equation (11) relevant for these spin $s=+2$ and $s=-2$ fields, respectively [20].

This paper is organized as follows. In Sec. 2 first the analytic framework is introduced. This part (see Sec. 2.2) is to set up the form of the wave equation to be solved once the conformal compactification and suitable regularizations of the basic variables have been carried out. Some of the details of the new type of energy and angular momentum balance relations are also discussed in Sec. 2.3. More details on the analytic and numerical setup are given in Sec. 3 containing the discussion of multipole expansions, the choice made for the initialization of time evolution, the determination of the decay rates and a summary of the results prior to ours. Sec. 4 presents all of our numerical findings. We start by discussing axially symmetric configurations in Sec. 4.2 which is followed by a detailed investigation of nonaxisymmetric configurations in Sec. 4.3. A thorough discussion of the use of the energy and angular momentum balance relations, as well as their use in verifying the convergence properties of the applied numerical implementation is given in Sec. 4.4. The discussions are completed by our final remarks in Sec. 5. The paper is closed by several appendixes providing useful details on the applied analytic and numerical settings.

2 Linear fields of spin $s$ on a Kerr background

This section provides all the details related to the applied analytic setup. In addition to the various forms of the Teukolsky master equation, the energy and angular momentum balance relations will also be discussed briefly. This section also lays down the mathematical groundwork that will be applied in later papers of this series.

2.1 Teukolsky equation

The metric of the Kerr background in Boyer-Lindquist coordinates $(t,r,\vartheta,\phi)$ can be given by the line element

[TABLE]

where $\Sigma=r^{2}+a^{2}\cos^{2}\vartheta$ and $\Delta=r^{2}-2Mr+a^{2}$, whereas $M$ and $a$ are the mass and the angular momentum per unit mass parameters of the Kerr black hole.

A linear sourceless field of spin $s$—where $s$ is integer or half-integer—on a Kerr background using the Kinnersley tetrad is known to be subject to the homogeneous Teukolsky master equation [20]

[TABLE]

A very simple covariant form of the Teukolsky master equation was introduced by Bini et al. in [5] which allows us to write (11) in the compact form

[TABLE]

where the components of the “connection vector” $\Gamma^{a}$ are [5]

[TABLE]

2.2 Regularization of the basic variables

The solutions to the homogeneous Teukolsky master equation are known to get singular either in the $\Delta\rightarrow 0$ or $r\rightarrow\infty$ limit. In order to regularize by suitable rescaling of the field variables and also to get the desired conformal compactification of the Kerr background, new coordinates are introduced. This is done by following the proposal in [19] in two succeeding steps. First, the Boyer-Lindquist coordinates $(t,r,\vartheta,\phi)$ are replaced by the ingoing Kerr coordinates $(\tau,r,\vartheta,\varphi)$, where the new time and azimuthal coordinates $\tau$ and $\varphi$ are defined via the relations

[TABLE]

Note that the $\tau=const$ hypersurfaces, as indicated in Fig.1. in [19], are horizon penetrating, whereas they all tend to spacelike infinity in the $r\rightarrow\infty$ limit.

In the second step, new coordinates $(T,R)$ replacing $(\tau,r)$ are introduced by the implicit relations

[TABLE]

The most important advantage of the application of these new coordinates is that they allow the use of a conformal compactification of the Kerr spacetime such that future null infinity ($\mathscr{I}^{+}$) gets to be represented by the $R=1$ hypersurface and also that all the $T=const$ hypersurfaces are such that they are both horizon penetrating and intersecting $\mathscr{I}^{+}$ in regular spherical cuts at $R=1$.

In order to get the desired suitably regularized basic field variables on the Kerr background $\psi^{(s)}$ is replaced by $\Phi^{(s)}$ defined as

[TABLE]

Once all the foregoing steps have been performed, the homogeneous Teukolsky master equation (12) can be seen to take the form

[TABLE]

where $Y_{1}^{0}$ and $Y_{2}^{0}$ stand for the zero spin-weight spherical harmonics with $\ell=1,2$ and $m=0$, whereas the explicit form of the involved $R$-dependent coefficients is given in Appendix A. Note also that by making use of the $\eth$ and $\overline{\eth}$ operators acting on a spin $s$ field $f$ as

[TABLE]

all the $\vartheta$ derivatives present in the Laplace-Beltrami operator can be incorporated into a single operator $\overline{\eth}\eth$ via the relation

[TABLE]

This provides us considerable simplification and enhances the accuracy of our numerical integrator significantly, as our approach is based on the use of multipole expansions of the basic variables in terms of spin-weighted spherical harmonics, and all the angular derivatives in (2.2) can be evaluated analytically. (More details on the use of spin-weighted spherical harmonics and their relations to the $\eth$ and $\overline{\eth}$ operators can be found in Appendix B.)

2.3 Conserved currents

It is known that there is no way to construct a Lagrangian out of a single spin $s$ variable $\psi^{(s)}$ and its first order derivatives such that one could get (12) as the corresponding Euler-Lagrange equation for $\psi^{(s)}$. Nevertheless, as pointed out in [22] recently, it is possible to associate a meaningful Lagrangian to a pair of spin $s$ and $-s$ fields via the relation

[TABLE]

It was also shown in [22] that by making use of some suitable infinitesimal transformations of the form $\psi^{(\pm s)}\to\psi^{(\pm s)}-\varsigma h^{a}\partial_{a}\psi^{(\pm s)}$ of the Lagrangian in (26), canonical conserved Noether currents can also be associated with a pair of spin $s$ and $-s$ fields. Note that the involved fields may be completely independent. The only requirement for the conservation of the currents is that they both satisfy their respective Teukolsky master equations. In particular, as the Lagrangian in (26) is invariant with respect to the one-parameter group of diffeomorphisms induced by the Killing vector fields $h^{a}=T^{a}=(\partial_{T})^{a}$ and $h^{a}=\varphi^{a}=(\partial_{\varphi})^{a}$ (which are also coordinate basis fields), the corresponding infinitesimal transformations endow us with the canonical energy- and angular-momentum-type currents defined as

[TABLE]

Notice that for $s=0$, (27) and (28) produce the well-known energy and angular momentum of the massless complex scalar field if \psi^{(-s)}\big{|}_{s=0}=\overline{\psi^{(s)}}\big{|}_{s=0}.

As the covariant divergence of these currents vanish, the balance relations

[TABLE]

hold. The spacetime domain of integration $\Omega$ in all of our applications is chosen to be the rectangular coordinate domain in $(T,R)$ such that it is bounded by some initial and final time slices $T=T_{i}$ and $T=T_{f}$ and by some inner and outer timelike or null cylinders given by the $R=R_{in}$ and $R=R_{out}$ hypersurfaces, respectively. In particular, denoting by $n_{a}^{(T)}=(\mathrm{d}T)_{a}/\sqrt{g^{TT}}$ and $n_{a}^{(R)}=(\mathrm{d}R)_{a}/\sqrt{-g^{RR}}$ the respective normals of the $T=const$ and $R=const$ hypersurfaces and by $h_{T}$ and $h_{R}$ the determinant of the restriction of the metric to these hypersurfaces, the energy balance relation can be given as

[TABLE]

A completely analogous balance relation can be derived for the angular momentum simply by replacing $E^{a}$ in (31) with $J^{a}$. The explicit form of the integrands $n_{a}^{(T)}E^{a}\sqrt{|h_{T}|}$, $n_{a}^{(R)}E^{a}\sqrt{|h_{R}|}$, $n_{a}^{(T)}J^{a}\sqrt{|h_{T}|}$, and $n_{a}^{(R)}J^{a}\sqrt{|h_{R}|}$ involved in the energy and angular momentum balance relations can be found in Appendix C.

3 More on the analytic and numerical setup

The discussions in this section will remain valid for any linear field of spin $s$ satisfying the homogeneous Teukolsky master equation on a Kerr black hole background.

3.1 Multipole expansions

In solving (2.2), our basic variables $\Phi^{(s)}$ are expanded in terms of spin-weight $s$ spherical harmonics ${}_{s}{Y_{\ell}}{}^{m}$ as

[TABLE]

In this way, (2.2) becomes a set of coupled $(1+1)$-dimensional linear wave equations for the expansion coefficients $\phi_{\ell}{}^{m}(T,R)$, whereas all the angular derivatives are evaluated analytically by making use of the $\eth$ and $\overline{\eth}$ operators (more details on spin-weighted spherical harmonics and the $\eth$, $\overline{\eth}$ operators can be found in Appendix B). The summation goes from $\ell=|s|$ to some $\ell=\ell_{max}$ value, which is chosen to be suitably large in order to keep the truncation error tolerably small (which, in practice, corresponds to numerical precision).

Note that in the frequency domain analysis, the eigenfunctions of the angular part of Teukolsky master equation (TME) are spin-weighted spheroidal harmonics instead of spin-weighted spherical harmonics. In the work of Casals et al. [7], the authors demonstrated that expanding in terms of spherical harmonics instead of spheroidal ones results in branch cuts in the complex Green’s function. However, in the integral, the contributions of these extra cuts cancel out, so in the end the Green’s function will be the same, and it is safe to use spin-weighted spherical harmonics as a basis. Since the spheroidal harmonics themselves can be expanded in terms of spherical harmonics, at least the decay rate of the slowest decaying modes is the same regardless of the basis of expansion.

By applying standard order reduction techniques—by introducing $(\phi_{T}){}_{\ell}{}^{m}=\partial_{T}\phi_{\ell}{}^{m}$ and $(\phi_{R}){}_{\ell}{}^{m}=\partial_{R}\phi_{\ell}{}^{m}$ as additional dependent variables—a first order strongly hyperbolic system is introduced for a vector variable that is composed of the multipole expansion coefficients. These equations were evolved in our numerical code by applying the method of lines in a fourth order Runge-Kutta integrator and also using a sixth order dissipation term for suppressing high-frequency spurious solutions [10].

3.2 The applied initial data

In solving (2.2) in any time-evolution scheme it is necessary to initialize $\Phi^{(s)}$ by specifying, on some $T=T_{0}\,(\in\mathbb{R})$ initial data surface, a pair of functions $(\phi^{(s)},\phi_{T}^{(s)})$ such that $\phi^{(s)}=\Phi^{(s)}|_{T=T_{0}}$ and $\phi_{T}^{(s)}=(\partial_{T}\Phi^{(s)})|_{T=T_{0}}$ hold. In order to be able to uncover the coupling between various modes in (32) characterized by their multipole indices $s,\ell,m$, the applied initial data will always be a single mode excitation. In addition, they will be either “static” or “purely dynamical”.

Accordingly, for some fixed values of the $s,\ell,m$ indices, a single mode excitation type of initial data $(\phi^{(s)},\phi_{T}^{(s)})$ is called static (ST) if

[TABLE]

while it will be called purely dynamical (PD) if

[TABLE]

hold on the $T=T_{0}$ initial data surface, and no summation is meant in either (33) or in (34). Because of the linearity of the TME (2.2) and the use of single mode excitations, with multipole indices $s,\ell,m$, all the excited modes in (32) will share the $s$ and $m$ values. In what follows, the $\ell$ parameter of the exciting mode will be distinguished by priming it, i.e., denoting it by $\ell^{\prime}$, while the $\ell$ parameter of the excited modes will be referred to without priming.

The $R$ dependence of the nonvanishing part of the initial data—$\phi_{\ell^{\prime}}{}^{m}(T_{0},R)$ in the static case or $(\partial_{T}\phi)_{\ell^{\prime}}{}^{m}(T_{0},R)$ in the purely dynamical case—is restricted by choosing it to be the “bump” with center $c$ and width $w$,

[TABLE]

Note that $\mathcal{B}(R)$ is a smooth function of compact support, the parameters of which were fixed in our numerical simulations as $c=0.7$ and $w=0.1$.

3.3 The late-time behavior

Likewise, in the case of spin zero fields—when monitored at certain fixed spatial locations—after an initial dynamical phase, each of the excited modes $\phi_{\ell}{}^{m}(T,R)$ go through a lasting quasinormal ringing period which is supplemented by a late-time tail behavior. In particular, this means that for sufficiently large values of $T$, the individual multipole expansion coefficients at any fixed $R=R_{0}$ spatial location decrease as

[TABLE]

where $n$ is a positive integer, the value of which may depend on all the involved parameters; i.e., in general, it has the functional form $n=n(s,\ell,m,\ell^{\prime})$.

3.3.1 Earlier analytic and numerical results on the decay rates

In studying the tail behavior the functional form of the decay rate $n=n(s,\ell,m,\ell^{\prime})$ is the center of interest. In interpreting our numerical findings, we may use as our reference frames two independent investigations prior to ours. First, we may refer to the detailed analytic studies carried out by Hod [14]. Second, the accurate numerical investigations carried out in [11] provide us important clues concerning the functional dependence of the decay rates.

Both of these investigations have some limitations in their scopes. For instance, the decay rates determined in [14] by using some Green’s-function-based analytic arguments apply only to purely dynamical initial data; i.e., no decay rates are derived therein for static initial data. It is worth also mentioning that—based on analogous analytic investigations carried out in [11] (see, in particular, Appendix B therein)—in certain (though very limited) subcases, these estimates of Hod were claimed to be imprecise. The numerical studies in [11] were also somewhat restricted, as apart from a few special cases, they were limited to axisymmetric configurations.

Despite the slight limitations in their scopes—as the aforementioned investigations are also complementary to each other—they provide us important guidance in carrying out and interpreting our numerical results. Below, we summarize the most important findings reported in [14, 11]. In particular, for purely dynamical initial data the pertinent results obtained by Green’s-function-based analytic arguments and by numerical simulations can be summarized as follows:

1. (1)

At the horizon $R=R_{+}$,

[TABLE]

where $\ell_{0}=\max\{|m|,|s|\}$ is the lowest allowed value of $\ell$ for given $s$ and $m$, and $\alpha=0$ in all cases except if $s>0$ and $m=0$ when $\alpha=1$. We find $\delta=1$ in the analytic investigations of Hod [14] but numerical results in [11] suggest that at least in the case $m=0$, the correct value is $\delta=0$. 2. (2)

At finite intermediate spatial locations with $R_{+}<R<1$, the values of $n$ can be deduced from (37) by the substitution of $\alpha=0$, i.e.,

[TABLE]

Similarly, as in (37) $\delta=1$ in [14]; however, in this case, in Appendix B of [11] Harms et al. pointed out a missed case in the argument of Hod, and as a result, $\delta=0$ when $m=0$. This correction is also consistent with the numerical results of [11]. 3. (3)

At $R=1$, representing future null infinity $\mathscr{I}^{+}$,

[TABLE]

where $\gamma=0$ in all cases except if $m=0$, $\ell^{\prime}=\ell_{0}+1$, and $\ell=\ell_{0}$ when $\gamma=1$. With vanishing $\gamma$, (39) reproduces the result of Hod [14]; nevertheless, in Appendix B of [11] the authors pointed out a missed case in the corresponding argument of [14], necessitating the inclusion of $\gamma$ in (39). This correction was also verified numerically in [11].

Although the numerical results presented in [11] are mostly for axisymmetric configurations, there are two sets of decay exponents gained from nonaxisymmetric purely dynamical initial data with $s=0$, $m=1$ and $s=-2$, $m=2$, respectively. The exponents for $s=0$, $m=1$ and $s=-2$, $m=2$ are in agreement with the predictions of [13, 14].

In the case of static initial data, we may only refer to the empirical results in [11] where, as mentioned above, the functional forms are limited to axially symmetric configurations with $m=0$. The pertinent findings in [11] can be summarized as follows:

1. (1’)

At the horizon $R=R_{+}$,

[TABLE]

where now $\alpha=0$ for $s=0$ and $\alpha=1$ otherwise. 2. (2’)

At any finite intermediate spatial location with $R_{+}<R<1$, the values of $n$ can be given as in (40), with the distinction that $\alpha=0$ in all possible cases, i.e.,

[TABLE] 3. (3’)

At $R=1$,

[TABLE]

As both the analytic and numerical setups outlined in the previous sections allow us to study nonaxisymmetric configurations, one of our main motivations in this paper is, besides verifying (40)–(42), if possible, to deduce the corresponding generalizations of these relations.

3.3.2 The local power index

As mentioned above, the excited modes following a long-lasting quasinormal ringing phase all end up in a power-law decay $\phi_{\ell}{}^{m}\sim T^{-n}$ [11, 14, 19] as depicted in Fig. 1.

The question is then how to determine the specific value of $n$. In practice, the value of the decay rate $n$—for each specific mode—is approximated by the local power index (LPI) $\mu$ determined as

[TABLE]

where, for simplicity, the multipole indices $\ell$ and $m$ of $\phi_{\ell}{}^{m}$ are suppressed. Note that, just as the decay rate $n$, the local power index $\mu$ also depends on the spatial location as well as on the multipole indices $s$, $\ell,m$, and $\ell^{\prime}$ of the involved modes [11, 19].

4 Numerical results

The discussions up to this point are valid for any linear field of spin $s$ satisfying the homogeneous Teukolsky master equation on a Kerr black hole background. In this section, considerations will be restricted to the numerical investigations of linear fields of spin $s=\pm 1,\pm 2$, and our numerical findings will be reviewed.

4.1 More specific settings

4.1.1 The input parameters

In most of our numerical simulations the parameters of the Kerr background were fixed as $M=1$ and $a=0.5$. A cutoff in the multipole expansion (32) at $\ell_{max}=8$ was found to be completely satisfactory for achieving the necessary precision. Some experimentation with doubling $\ell_{max}$ shows that the accuracy of our results is not limited by this value. This suggests that the error introduced by radial discretization and numerical arithmetic is still bigger than the error introduced by this cutoff. In many of the simulations, the use of $1024$ spatial grid points on the $[R_{+},1]$ interval was satisfactory (see Fig. 2), though, in most cases we used $2048$, and there were also cases where the use $4098$ spatial grid points was more rewarding.

As for the applied numerical precision, note that in most of our simulations the use of long double arithmetic was sufficient. Nevertheless, as for negative values of the spin parameter $s$, the LPI values were always harder to be determined, and we used the computationally more expensive quadruple precision in those cases. Unfortunately, for negative $s$ values even the use of quadruple precision could not always guarantee the required level of precision at late time which, on the other hand, is essential to draw a conclusive result on the tail behavior. In many cases, the pertinent LPIs could not even be determined despite that for positive values of the spin parameter the LPI values are already sufficiently accurate with the use of long double precision and, in general, with the use of lower spatial resolution.

4.1.2 Assembling the determined LPIs

In the succeeding subsections, our numerical findings will be reported. First, the LPIs determined are collected in tables, and then the implications of the observed phenomena are described briefly.

The tables of LPIs are given for various values of the relevant parameters. In particular, $s=\pm 1,\pm 2$ with exciting modes $\ell^{\prime}$ and excited modes $\ell$ both taking (the allowed) positive integer values from the interval $1\leq\ell^{\prime},\ell\leq 5$. The values of $m$ will also be restricted to $0,\pm 1,\pm 2$.

In advance of turning to the contents of these tables, note first that concerning the dependences of the LPIs on the spin parameter $s$, the following two cases have significantly different characters. Whenever $s<0$, it suffices to determine the values of $\mu$ at the outer horizon $R=R_{+}$ and at the future null infinity $R=1$, as the LPIs relevant for the intermediate values $R_{+}<R<1$ are exactly the same as those at the outer horizon $R=R_{+}$. As opposed to this, whenever $s>0$, the LPIs have a higher variety; i.e., the value of $\mu$ at the intermediate location $R_{+}<R<1$ differs from that at the outer horizon $R=R_{+}$ though asymptotically the LPI takes the same value for any intermediate finite locations $R_{+}<R<1$. As it is clearly indicated by the graphs in Fig. 3, for $s>0$, the values of $\mu$ have to be monitored not only at the outer horizon ($R=R_{+}$) and at future null infinity ($R=1$) but for several intermediate $R=const$ locations as well.

Remarkably, the closer the intermediate $R$ value is to $R=R_{+}$ or to $R=1$, the longer it takes to settle down to the LPI value that is relevant for pertinent intermediate locations.

In accordance with the above outlined observations, the tables will be structured as follows: As for negative $s$ (with $s=-1$ or $s=-2$), the LPI values at the horizon and intermediate finite locations are the same, and for each slot labeled by $\ell^{\prime}$ and $\ell$, only the values “$\mu_{R_{+}}$” at the horizon and “$\mu_{{\mathscr{I}}^{+}}$” at future null infinity will be indicated by separating them with a vertical line in writing “$\mu_{R_{+}}|\,\mu_{{\mathscr{I}}^{+}}$.” As for the positive $s$ (with $s=1$ or $s=2$), the LPI values at the intermediate locations, with $R_{+}<R<1$ differ from that of $\mu_{R_{+}}$ and $\mu_{{\mathscr{I}}^{+}}$, and we arrange these three values by separating them with vertical lines in writing “$\mu_{R_{+}}\,|\,\mu_{R}\,|\,\mu_{{\mathscr{I}}^{+}}$.” In practice, the LPI value indicated in the middle slot will always correspond to the one measured along the $R=0.88$ timeline.

There is various additional information indicated in the succeeding tables. Note first that in most cases the precise value of the LPIs could be determined in a clean way by making use of (43) (see also Figs. 1–3). Nevertheless, there were some cases where the limited accuracy at the very late time did not allow us to draw a clear conclusion this way. In such cases, instead of using (43), a line fitting on the log-log plots was applied in order to determine the approximate value of the LPIs relevant for the pertinent modes. Note that, as the output of this method depends, for instance, on the choice of the subintervals where the fitting is carried out, the yielded LPI values have to be taken with some caveats. To warn the reader, the corresponding LPI values are always indicated by round brackets around the pertinent numbers, i.e., by writing $(\mu)$ instead of $\mu$. We also have to admit that there were cases—especially for negative values of $s$—where neither of the above-discussed methods could be sufficiently conclusive. All these cases will be indicated by filling up the pertinent slot labeled by the $\ell^{\prime}$ and $\ell$ values with a question mark “?.”

4.2 Axially symmetric configurations

Note first that in [11] detailed numerical investigations of the axisymmetric case were carried out. This gives us the opportunity to compare the observed LPI values and, thereby, to check the performance and the reliability of our code.

The LPIs obtained for the axially symmetric cases with $|s|=\pm 1,\pm 2$ and with $m=0$ are assembled in Table 1 for $|s|=\pm 1$ and in Table 2 for $|s|=\pm 2$.

The most characteristic features which deserve notice are as follows:

The values of the LPIs at $R_{+}$ and at intermediate locations $R_{+}<R<1$ are always increased by $1$ if the value of the excited mode $\ell$ is increased by $1$. 2. 2.

The values of the LPIs at ${\mathscr{I}}^{+}$ do not necessarily increase in an analogous, strictly monotonous way, though they never decrease either while the value of $\ell$ is increased. 3. 3.

The above two observations imply that the mode with the lowest possible $\ell$ value always decays at the slowest rate, though at ${\mathscr{I}}^{+}$ there may be some other $\ell$ modes which are also decaying at the same rate. 4. 4.

For $s>0$, the value $\mu_{R_{+}}$ is always larger by $1$ than the LPI at intermediate locations $R_{+}<R<1$. 5. 5.

All observations made in the previous point are in full agreement with the analytic results and numerical findings in [14] and [11], respectively. 6. 6.

For $|s|=2$, the most significant distinction is that the first row and the first column of the tables signified by $\ell^{\prime}$, $\ell$ slots are missing, as the corresponding modes are excluded by the fact that all the spin-weighted spherical harmonics with $\ell<\max\{|s|,|m|\}$ are annihilated by the $\eth$ operator.

Note, finally, that in Table 1 there are seven numerically determined LPI values—they are indicated by boldface (in color by orange)—that differ from those which can be deduced from the rules laid down in [14, 11]. Notably, the difference is always $1$. Nevertheless, these seven numbers are all in round brackets indicating that the corresponding estimates have to be taken with some caveats.

4.3 Nonaxially symmetric configurations

Note that in the nonaxisymmetric scenarios, due to the fact that the spin-weighted spherical harmonics with $\ell<\max\{|s|,|m|\}$ are annihilated by the $\eth$ operator the first row and the first column of the tables signified by the $\ell^{\prime}$, $\ell$ slots are always missing if either $|m|=2$ or $|s|=2$. Nevertheless, for the sake of easier comparisons between the axisymmetric and nonaxisymmetric cases, we assemble the LPI values using the interval $1\leq\ell^{\prime},\ell\leq 5$, as done previously.

4.3.1 $|m|=1,2$ with static initial data

In this subsection, excitations generated by the $|m|=1,2$ modes with static initial data are considered. The relevant LPIs are collected in Tables 3–6, and the main observations concerning them are summarized as follows:

Now, even for positive values of $s$, i.e. for $s=1$ and $s=2$, the pertinent LPI values at the horizon and the LPI values at intermediate locations coincide. 2. 2.

For any fixed value of $s$, the LPIs are independent of the sign of $m$. 3. 3.

The LPI values $\mu_{R_{+}}$ at the horizon are also independent of the sign of $s$. As opposed to this, the LPI values $\mu_{{\mathscr{I}}^{+}}$ at future null infinity depend on the sign of $s$, such that the values of $\mu_{{\mathscr{I}}^{+}}$ relevant for negative $s$ are always larger by $2\,|s|$ than the pertinent LPI values $\mu_{{\mathscr{I}}^{+}}$ for positive $s$. 4. 4.

The LPI values $\mu_{R_{+}}$ and $\mu_{{\mathscr{I}}^{+}}$ relevant for the considered nonaxisymmetric case are pairwise equal to the LPI values at intermediate locations $R_{+}<R<1$ and $\mu_{{\mathscr{I}}^{+}}$ relevant for the corresponding axially symmetric configurations (discussed in the previous section).

4.3.2 $|m|=1,2$ with purely dynamical initial data

The excitations generated by purely dynamical initial data are visibly more interesting. This gets to be transparent, for instance, when $\ell^{\prime}=\ell_{0}+1$. In this case the Green’s-function-based argument of Hod [14] was corrected by results in [11], though only for $m=0$. Nevertheless, our findings indicate that corrections are needed also for $m\neq 0$. The relevant LPI values are highlighted by boldface (in color: green) characters in the Tables 7–10.

As in other cases, the LPI values at the horizon $R_{+}$ are always increased at least by $1$ while the value of $\ell$ is increased by $1$. 2. 2.

The LPI values at future null infinity ${\mathscr{I}}^{+}$ do not necessarily increase, though they never decrease while the value of $\ell$ is increased. 3. 3.

For any fixed value of $s$, the LPIs are independent of the sign of $m$. 4. 4.

The LPI values $\mu_{R_{+}}$ at the horizon are independent of the sign of the spin parameter $s$, with some exceptions when $\ell^{\prime}=\ell_{0}+1$, $\ell\geq\ell^{\prime}$, $|s|=1$, $m=2$. Unfortunately, we do not have the data to check if a similar effect is present when $|s|=2$. Nevertheless, using also Table 5 of [11], it can be seen that the effect is not present for $|s|=2$, $m=2$. 5. 5.

As above, the LPI values $\mu_{{\mathscr{I}}^{+}}$ do depend on the sign of $s$ such that the value of $\mu_{{\mathscr{I}}^{+}}$ relevant for negative $s$ is always larger by $2\,|s|$ than the LPI value $\mu_{{\mathscr{I}}^{+}}$ pertinent to the corresponding positive $s$.

Concerning the two Green’s-function-based arguments, we draw the following conclusions from our numerical findings:

The LPI values at the horizon $\mu_{R_{+}}$ are reassuring that the modifications in [11] are valid only in the case $m=0$. A notable exception is the row $\ell^{\prime}=3$ in Table 8-a. There, the data suggest that for $\ell>2$ a similar modification is needed as in the $m=0$ case. 2. 7.

For $s>0$, apart from $\ell^{\prime}=\ell_{0}+1$, the values $\mu_{R_{+}}$ and the LPI values at intermediate locations $R_{+}<R<1$ are always equal to each other. Whenever $\ell^{\prime}=\ell_{0}+1$, our LPI values $\mu_{R}$ at intermediate locations, with the exception of $\ell=\ell_{0}$, are systematically larger by $1$ than predicted by (38). Again, this suggests that a correction similar to the one presented in [11] is needed whenever $s>0$, $\ell^{\prime}=\ell_{0}+1$, and $\ell>\ell_{0}$. 3. 8.

At future null infinity our numerically determined values for $\mu_{{\mathscr{I}}^{+}}$ are consistent with the results of [14, 11].

All in all, it appears that in the very limited subcase with $\ell^{\prime}=\ell_{0}+1$, there are certain special circumstances where neither the predictions made in [14] nor those made in [11] are supported by our numerical findings. The corresponding subcase would definitely desire more thoughtful analytic inspection.

4.4 Energy and angular momentum conservation

As it was shown in [22], to any pair of spin $s$ and $-s$ solutions of the Teukolsky master equations, there always exist some conserved canonical energy- and angular-momentum-type currents. Since the corresponding vector fields $E^{a}$ and $J^{a}$ are divergence free, i.e., $\nabla_{a}E^{a}=0$ and $\nabla_{a}J^{a}=0$, and by construction are complex, they provide us two complex balance relations as formulated by (31) along with its correspondent yielded by the replacement of $E^{a}$ by $J^{a}$. All in all, the real and imaginary parts of these two complex balance relations give us four real ones, and these four together provide us—in addition to the conventional checks such as the numerical convergence rate—a very important verification of the correctness and robustness of our numerical results.

In this subsection, first these energy and angular momentum balance relations will be applied. In particular, in the panels of Fig. 4, the time dependence of the real and imaginary parts of the numerical errors

[TABLE]

are plotted against time. The applied rectangular domain of integration $\Omega$ is bounded by an initial time slice $T=T_{i}$ and by a running $T$ slice, and also by the null cylinders at the horizon $R=R_{+}$ and at future null infinity at $R=1$, respectively. Note that on analytic grounds, the integrals $\delta E$ and $\delta J$ in the balance relations in (44) should vanish identically. Accordingly, they have to be (and they are expected to stay) small if the numerical implementation is correct.

It is also important to emphasize that in each panel of Fig. 4, the aforementioned time dependences are plotted for three different resolutions, with $1024,2048,4096$ spatial grid points. In this respect, the panels in Fig. 4 are not only to demonstrate that the real and imaginary parts of the balance relations (29) and (30) hold at a satisfactorily accurate level but also to make it transparent that our numerical implementation has a fourth order convergence rate—even in the considered highly complicated expressions deduced from the numerically determined basic variables—as it should happen as a fourth order accurate Runge-Kutta time integrator is applied in the $T-R$ section. To make this transparent, the numerical values of $\delta E$ and $\delta J$ relevant for the resolutions $2048$ and $4096$ are multiplied by $16$ and $256$, respectively.

If the numerical values of $\delta E$ and $\delta J$ were exactly proportional to the fourth power of the grid spacing $\delta R$, then the lines corresponding to the three different resolutions shown in the panels of Fig. 4 would exactly coincide. The apparent deviation from this prediction, which is very small in Fig. 4(b) but larger in the other three panels, can be explained by higher order corrections to the leading $\propto(\delta R)^{4}$ power law. In particular, the assumption that the dependence of the numerical values of $\delta E$ and $\delta J$ on $\delta R$ has the functional form $c_{4}(\delta R)^{4}+c_{5}(\delta R)^{5}+\dots$, where $c_{4}$ and $c_{5}$ are suitable constants, implies that the distance between the lines for the $2048$ and $4096$ grid points should be approximately half the distance between the lines for the $1024$ and $2048$ grid points, which can indeed be seen in Fig. 4.

The conserved currents $E^{a}$ and $J^{a}$ can also be used to analyze transport in the energy- and angular-momentum-type quantities, and we shall do so in the rest of this subsection. Nevertheless, in advance of doing so it is important to emphasize that these currents, apart from the spin-weight $s=0$ case, are not defined for individual fields but only for a pair of spin-weight $s$ and $-s$ solutions to the Teukolsky master equation. This restriction applies even though the spin-weight $s$ and $-s$ solutions can be generic; i.e., they could completely be independent of each other.

In the panels in Fig. 5,

the time dependence of the real and imaginary parts of the total energy and angular momentum,

[TABLE]

where $\Sigma_{T}$ signifies the $T=const$ slices, as well as the time dependence of the integrals of the corresponding flows

[TABLE]

where the latter integrals are supposed to be evaluated at the cylinders $R_{+}\times[T_{i},T]$ at the black hole event horizon and $R_{{\mathscr{I}}^{+}}\times[T_{i},T]$ at future null infinity, are plotted, respectively. These plots are showing only the initial parts of the time evolution; nevertheless, as $99\%$ of the transport processes happens in the corresponding initial intervals, these plots are informative about the characteristic features.

Even though the spin-weight $s$ and $-s$ solutions to the Teukolsky master equation may be completely independent in our simulations, they were yielded by applying exactly the same $R$-dependent bumpy profile (35), with exactly the same parameters as described in Sec. 3.2. Probably for this reason, there is a striking similarity between the real parts of the energy- and angular-momentum-type transport processes, as depicted in Figs. 5(a) and 5(c) and between the corresponding figures (5) and (7) in [19] depicting the true energy and angular momentum transport processes relevant for a single spin-weight zero scalar field investigated in detail in [19].

In closing this subsection, note that it is obviously interesting to know what happens with the conserved currents and with their integrals under the flip $a\mapsto-a$ of the sign of the rotation parameter of the background Kerr spacetime. Remarkably, only the following two flippings are induced by this transformation. Notably, a flipping happens in the signs of the imaginary part $\mathfrak{Im}(\mathscr{E}_{T})$ of the total energy $\mathscr{E}_{T}$ and with the imaginary part $\mathfrak{Im}(\mathscr{F}_{E})$ of the integrated energy flow $\mathscr{F}_{E}$ evaluated at $R_{+}$ and at ${\mathscr{I}}^{+}$, whereas a completely analogous flipping happens in the real part $\mathfrak{Re}(\mathscr{J}_{T})$ of the total angular momentum $\mathscr{J}_{T}$ and with the real part $\mathfrak{Re}(\mathscr{F}_{J})$ of the integrated energy flow $\mathscr{F}_{J}$ evaluated at $R_{+}$ and at ${\mathscr{I}}^{+}$. The sign flips in the real part of the angular momentum integrals seem to be in accordance with intuition as under the flip $a\mapsto-a$ of the rotation parameter the background Kerr black hole gets to be counterrotating with respect to the spin $s$ and $-s$ solutions of the Teukolsky master equation.

5 Final remarks

The time evolution and the late-time behavior of spin $s$ fields with $s=\pm 1,\pm 2$ on a fixed Kerr background was examined numerically. The applied mathematical setup incorporated the techniques of conformal compactification and the hyperbolic initial value formulation of these spin $s$ fields. A new code was introduced that was developed based on the two-parameter foliation of the Kerr background by topological two-spheres as determined by the Boyer-Lindquist $t=const$ and $r=const$ level surfaces. The angular dependences along the foliating two-surfaces were treated by applying multipole expansions of the basic variables in terms of spin-weighted spherical harmonics, by which the applied method gets to be fully spectral in these angular directions. In the complementary time-radial section, the method of lines in a fourth order accurate Runge-Kutta time integration scheme was applied.

The time evolution of purely static and purely dynamical initial data were investigated. The asymptotic decay rates were evaluated at the black hole event horizon, in the domain of outer communication and at future null infinity. This was done systematically by scanning through a significantly wide range of the input and output parameters $\ell^{\prime},m,s,\ell$. The deduced decay rates were compared to those which were deduced prior to our investigation. In particular, our most important findings can be summarized as follows by referring to the results covered by [11, 12, 13, 14].

As noted earlier, in the case of static initial data, we compared our numerical results with those deduced also numerically in [11]. The only limitation in doing so was that (apart from a few exceptions) the decay exponents were determined in [11] only for axially symmetric configurations with $m=0$. Remarkably, for purely static initial data, if $|s|$ is replaced by $\ell_{0}=\max\{|m|,|s|\}$ everywhere in (40) and (41)—these were deduced in [11] for the $m=0$ case—the yielded relations [see (47) and (48) below] are automatically valid to the fully general (not necessarily axially symmetric) configurations. Accordingly, for static configurations, our findings—these are collected in Tables 1-a, 1-b, 2-a, 2-b, for axisymmetric configurations with $m=0$ and in Tables 3–6 for nonaxisymmetric configurations with $m\not=0$—can be summarized as follows:

1. (ST1)

At the horizon $R=R_{+}$

[TABLE]

where—now as it was in (37)—$\ell_{0}=\max\{|m|,|s|\}$ and $\alpha=0$ in all cases except if $s>0$ and $m=0$, when $\alpha=1$. 2. (ST2)

At any finite intermediate spatial location with $R_{+}<R<1$,

[TABLE] 3. (ST3)

Finally, at future null infinity signified by $R=1$,

[TABLE]

As for the solutions of the Teukolsky master equation starting with purely dynamical initial data, as indicated earlier, the situation gets more involved at least in certain subcases. The relevant data are collected in Tables 1-c, 1-d, 2-c, 2-d for axisymmetric and in Tables 7-10 for nonaxisymmetric initial data.

1. (PD1)

For instance, at the horizon with $R=R_{+}$, for $m\neq 0$ and $\ell^{\prime}=\ell_{0}+1$ the decay rates based on our numerical findings do not fully agree with the former predictions. Our results can be described by the following modification of (37):

[TABLE]

where $\ell_{0}=\max\{|m|,|s|\}$, $\alpha=0$ in all cases except if $s>0$ and $m=0$, when $\alpha=1$, and also $\tilde{\delta}=0$ if $m=0$ or $m=2$ with $s=-1$ and $\ell>\ell_{0}$, $\tilde{\delta}=1$ in all other cases. 2. (PD2)

At finite intermediate spatial locations with $R_{+}<R<1$, we also found that there are differences between (38) and our results if $\ell^{\prime}=\ell_{0}+1$ and $m\neq 0$. Again, this difference can be taken into account by a minor modification of $\delta$:

[TABLE]

where $\ell_{0}$ and $\alpha$ are as above, and also $\widehat{\delta}=0$ if $m=0$ or $m\neq 0$ but $\ell>\ell_{0}$ and $\widehat{\delta}=1$ otherwise. 3. (PD3)

At $R=1$ representing future null infinity $\mathscr{I}^{+}$, all of our pertinent numerical findings support the predictions of (39); i.e., the following rules apply

[TABLE]

where $\gamma=0$ in all cases except if $m=0$, $\ell^{\prime}=\ell_{0}+1$, and $\ell=\ell_{0}$, when $\gamma=1$.

The apparent dependence of the decay rate on the value of the azimuthal parameter $m$ at the horizon $R=R_{+}$ and intermediate locations $R_{+}<R<1$ indicate that there may be interesting features which will definitely desire further investigation. To clear this up, along with possibly some other interesting issues is, however, out of the scope of the present paper and they are left for future studies.

Acknowledgments

K. C. and G. Z. T., were supported in part by the NKFIH Grants No. K-115434 and No. K-116505, respectively. I. R. was supported by the POLONEZ program of the National Science Centre of Poland (under Project No. 2016/23/P/ST1/04195) which has received funding from the European Union’s Horizon 2020 research and innovation program under the Marie Skłodowska-Curie Grant Agreement No. 665778.

Appendixes

The succeeding appendixes are to provide the explicit form of the expressions mentioned and applied in the main part of the discussions of the foregoing sections. Also, we give a brief summary of the use of spin-weighted spherical harmonics.

Appendix A The homogeneous TME

As mentioned in Sec. 2.2, after regularization, the homogeneous Teukolsky master equation takes the form

[TABLE]

All the involved $R$-dependent coefficients are listed below. In particular, $\mathscr{A}$ and $\mathscr{B}$ are

[TABLE]

[TABLE]

The rest of the involved coefficients are

[TABLE]

[TABLE]

[TABLE]

Appendix B Spin-weighted spherical harmonics

One of the advantages of applying the multipole expansion

[TABLE]

in terms of spin-weight $s$ spherical harmonics ${}_{s}{Y_{\ell}}{}^{m}$ is that all the angular derivatives in (2.2) can be evaluated analytically. This is done by making use of the $\eth$ and $\overline{\eth}$ operators which act on a function $f$ of spin-weight $s$ as

[TABLE]

It is straightforward to check that for their commutator, the relation

[TABLE]

holds. It is also important to mention that all the $\vartheta$ derivatives present in the Laplace-Beltrami operator relevant for spin $s$ fields are contained by the operator $\overline{\eth}\eth$ as the relation

[TABLE]

can be seen to hold.

In addition, the spin-weighted spherical harmonics ${}_{s}{Y_{\ell}}{}^{m}$ are eigenfunctions of the operators $\overline{\eth}\eth$ and $\partial_{\varphi}$ and their spin weight is shifted by $\eth$ and $\overline{\eth}$, as they satisfy

[TABLE]

Note that as $\eth\,{}_{\ell}{Y_{\ell}}^{m}=0$ for $s=\ell$ and $\overline{\eth}\,{}_{-\ell}{Y_{\ell}}^{m}=0$ for $s=-\ell$, all the spin-weighted spherical harmonics ${}_{s}{Y_{\ell}}^{m}$ with $\ell<|s|$ vanish.

Using these relations, it is also straightforward to check that the spin-weighted spherical harmonics can be generated by the spin-raising $\eth$ and the spin-lowering $\overline{\eth}$ operators starting with the conventional (zero spin-weight) spherical harmonics $Y_{\ell}{}^{m}$ as

[TABLE]

if $0\leq s\leq\ell$ and

[TABLE]

if $0>s\geq-\ell$. Note finally that the complex conjugation acts as

[TABLE]

In virtue of Eq. (2.2), one of the coefficients involves a multiplication by $Y_{1}{}^{0}$, whereas a division by the expression $\mathscr{A}+\mathscr{B}\cdot Y_{2}^{0}$ must also be performed. Concerning multiplications, as it was explained in great detail in the appendixes of [8], once we have the expansion of scalar variables in terms of spherical harmonics, the products of these variables can easily be evaluated by making use of the Gaunt coefficients $G_{\ell_{1}}^{m_{1}}{}_{\ell_{2}}^{m_{2}}{}_{\ell_{3}}^{m_{3}}=\int Y_{\ell_{1}}^{m_{1}}Y_{\ell_{2}}^{m_{2}}Y_{\ell_{3}}^{m_{3}}\sin\vartheta\,\mathrm{d}\vartheta\mathrm{d}\varphi$. Notably, completely analogous arguments apply when our spin $s$ variables are expanded in terms of spin-weighted spherical harmonics. The corresponding Gaunt coefficients—which now acquire three additional spin indices—can be given as

[TABLE]

which can also be expressed via the Wigner-3j symbols as

[TABLE]

Then, in particular, the product $Y_{\ell_{1}}{}^{0}\cdot{}_{s}Y_{\ell_{2}}{}^{m}$ can be evaluated as

[TABLE]

As for the division by the term $\mathscr{A}+\mathscr{B}\cdot Y_{2}^{0}$, note that by following the ideas introduced in [8] (see also [15, 19]) the division by this term can also be traced back to multiplications. The essential observation made in [8] was that the Neumann series expansion

[TABLE]

can be applied to do so. In particular, by replacing $A$ with $1+x$, where $x=\mathscr{B}/\mathscr{A}\cdot Y_{2}^{0}$, and by choosing the value $k_{max}$ sufficiently large—as the term $\mathscr{B}/\mathscr{A}\cdot Y_{2}^{0}$ for any value of $\vartheta$ is much smaller than $1$—the approximate relation

[TABLE]

holds. Note that in practice numerical precision with an error tolerance of $10^{-20}$ does not require the use of more than $k_{max}=12$ terms in this series.

Appendix C The conserved currents

In order to evaluate the energy and angular momentum balance relations we need the explicit form of $\sqrt{|h_{T}|}\,n^{(T)}_{a}E^{a}$, $\sqrt{|h_{R}|}\,n^{(R)}_{a}E^{a}$, $\sqrt{|h_{T}|}\,n^{(T)}_{a}J^{a}$ and $\sqrt{|h_{R}|}\,n^{(R)}_{a}J^{a}$. The following subsections list the explicit form of the implemented expressions.

C.1 Energy density

[TABLE]

where

[TABLE]

C.2 Energy current

[TABLE]

where

[TABLE]

C.3 Angular momentum density

[TABLE]

where

[TABLE]

[TABLE]

C.4 Angular momentum current

[TABLE]

where

[TABLE]