3 where r$_{∗}$ is the regular Kerr spacetime ‘tortoise’ coordinate defined by dr$_{∗}$/ dr = ( r 2 + a $^{2}$) / Δ. Our code is a 2+1D first–order code for the time evolution of modes m of φ, Π given two independent functions of two variables at t = 0 as initial data, based on [9]. The convergence order is taken to be second–order in the radial and temporal directions, and sixth order in the angular direction (although some of the higher– ℓ results were obtained with tenth–order angular operators). The standard grid resolution we use is 64 , 000 grid points in the radial direction, and 64 angular grid points (which we denote as 64 K × 64), unless stated otherwise, taking the temporal step to be the largest step consistent with the Courant condition. In practice, we choose as initial data φ ( r, θ ) | t =0 = 0 and Π( r, θ ) | t =0 = f ( r ) P$_{ℓ}$ ( cos θ ), where the radial function f ( r ) is chosen to be a gaussian f ( r ) = (1 $_{/}$√$_{2}$$_{πσ}$ $^{2}$) exp[ − ( r$_{∗}$ − r$_{∗}$ $_{0}$)$^{2}$/ (2 σ $^{2}$)] with r$_{∗}$ 0 = 25 M and σ = 6 M . The outer and inner boundaries are placed at r boundary ∗ = ± 800 M , which allows us to integrate to t = 1 , 500 M at r$_{∗}$ 0 without seeing boundary reflection effects. We choose the Kerr parameter a = 0 . 995 M . We use quadrupole or octal precision floating–point arithmetic throughout according to need. We emphasize that our choice of initial data —based on the field and its conjugate momentum— is in accord with other common choices. Specifically, when we choose a single value of ℓ , our choice of initial data is a pure multipole moment, equivalent to choices made based on the field and its time derivative. There is no ambiguity in the definition of a pure multipole [7]. Since the goal of this work is to investigate the late-time decay rate of higher multipole modes of Kerr black holes a number of improvements to our time-domain Teukolsky equation code were necessary. This is because generating accurate numerical data in the context of higher multipole modes involves a number of challenges. Firstly, these simulations need to be rather long — this is because typically the observed field exhibits an exponentially decaying oscillatory behavior (“quasi- normal ringing”) in the initial part of the evolution and only much later this transitions over to a clean power-law decay. One needs to wait for these oscillations to dissipate away. For higher multipoles this phase lasts much longer, thus requiring significantly longer evolutions. Secondly, because each multipole has a different decay rate (which increases with an increase in ℓ ) at late times one ends up with numerical data in which different multipoles have widely different amplitudes (often 30 – 40 orders of magnitude apart!). For this reason, not only does the numerical solution scheme have to be high-order (to reduce the discretization errors to the required levels) but it also requires high-precision floating-point numerical computation (due to the large range of amplitudes involved). B. Code Improvements Below we list the major improvements made since the publication of [7] to our Teukolsky equation evolution code: 1. Higher-Order Algorithm: As mentioned above, we require a higher-order numerical evolution scheme to keep the discretization errors to sufficiently low levels, otherwise these can overwhelm the numerical data associated to the weakly excited higher multipole modes of interest. It turns out that it is sufficient that only the angular differentiation (i.e. θ -derivatives of the field) be implemented using a higher-order numerical stencil. The temporal and the radial direction related operations can simply stay 2nd-order and such a mixed approach yields sufficiently good results [7]. For this reason, in this work we choose the finite-difference angular differentiation operator to be 10th-order accurate and leave the rest of the numerical scheme as a standard 2nd-order Lax- Wendroff algorithm. 2. Higher-Order Projection Integrals: We also require a matching 10th-order accuracy in mode projection integrals to prevent numerical integration errors from overwhelming the higher multipole numerical data. 3. Higher-Order Numerical Precision: As pointed out above, we also require high-numerical precision in our computations. In particular, double, quadruple and octal precision may be required depending upon the details of Kerr tails simulation being attempted. This is simply to prevent roundoff error from overwhelming the weakly excited modes of interest – especially important for “up” excited higher multiple modes! We made use of quadruple precision (in some cases, octal precision) floating-point accuracy in this work. Now, very few current computer processors support full quadruple-precision (128–bit) datatype and operations. And to the best of our knowledge, no compute hardware natively supports octal-precision (256–bit) arithmetic. Thus, we implemented such high-precision computations in software using publicly available software libraries. This made our work particularly challenging computationally, because with each such increase in numerical precision, our code’s overall performance drops by an order-of-magnitude. 4. Code Parallelism and HPC Hardware: In order to perform the numerical simulations in a reasonable time- frame, we made use of message-passing (MPI) based parallelism in our code and also used sophisticated HPC