Numerical computation of electromagnetically sourced nonlinear tails

TL;DR

This paper numerically investigates nonlinear electromagnetic effects on black hole ringdown, revealing second-order gravitational tails that decay slower than linear ones, potentially impacting multi-messenger astrophysical observations.

Contribution

It introduces a numerical method to compute second-order gravitational tails induced by electromagnetic sources in black hole environments, highlighting their decay rates and significance.

Findings

Second-order tails decay as t^{-2l-2} at fixed position.

Tails decay as u^{-l-3} at null infinity.

Nonlinear tails decay slower than linear counterparts.

Abstract

Amazingly, recent studies indicate that nonlinear effects are of great significance for modelling black hole ringdown. Transient electromagnetic events in the astrophysical environment are typically high energetic, potentially responsible for some nonlinearities in ringdown. Motivated by the desire to understand these nonlinearities, we solve the inhomogeneous Bardeen-Press-Teukolsky equation numerically, and find second-order gravitational tails induced by an electromagnetic source. Our results suggest that the second-order tails of curvature perturbations with multipole numbers $l\geq4$ decay as $t^{-2l-2}$ at fixed spatial position and $u^{-l-3}$ in retarded-time $u$ at null infinity, slower than their linear counterparts, which can play a role in multi-messenger observations.

Tables

Table 1. Table 1: Relation between the rescaled scalars ψ ^ s , 𝒮 ^ s {}_{s}\hat{\psi},{}_{s}\hat{\mathcal{S}} and the NP scalars.

$s$	${}_s{}^{}\hat{\psi }$	${}_s{}^{}\widehat{𝒮}$
-2	${\Psi }_4/R$	$\left(2L^4/R^3\right){}_{-2}{}^{}𝒮$
-1	${\varphi }_2/R$	$\left(2L^4/R^3\right){}_{-1}{}^{}𝒮$
+1	${\varphi }_0{\left({\Psi }_2/M\right)}^{-2/3}/R$	$\left(2L^4/R^3\right){\left({\Psi }_2/M\right)}^{-2/3}{}_{+1}{}^{}𝒮$
+2	${\Psi }_0{\left({\Psi }_2/M\right)}^{-4/3}/R$	$\left(2L^4/R^3\right){\left({\Psi }_2/M\right)}^{-4/3}{}_{+2}{}^{}𝒮$

Table 2. Table 2: The position of source 𝒮 ^ [ 55 ] − 2 {}_{-2}\hat{\mathcal{S}}^{[55]} ’s peak at different time slices corresponding to Fig. 9 . The coordinate speed of the peak is much less than the outgoing characteristic c + ​ ( R = 0 ) = − 1 / 8 c_{+}(R=0)=-1/8 .

$T/M$	700	800	900	1000	1100	1200
$r/{10}^{3}M$	$1.4964$	$1.7798$	$1.9528$	$2.1523$	$2.3841$	$2.6555$
$R/{10}^{-3}\frac{L^2}{M}$	$0.6683$	$0.5619$	$0.5121$	$0.4646$	$0.4195$	$0.3766$

Equations

$${}_{-1}\mathcal{T}{\varphi }_2^{(1)}={}_{-1}\mathcal{S}=0,$$

$${}_{-2}\mathcal{T}{\Psi }_4^{(2)}={}_{-2}\mathcal{S}({\varphi }_2^{(1)};{\varphi }_2^{(1)}),$$

$${}_{-2}\mathcal{S}=d_4^{(0)}{\mathcal{R}}_d^{(2)}+d_3^{(0)}{\mathcal{R}}_h^{(2)},$$

$${\mathcal{R}}_d^{(2)}=$$

$$-(\Delta +\bar{\mu }+2\gamma -2\bar{\gamma }{)}^{(0)}{\Phi }_{20}^{(2)},$$

$${\mathcal{R}}_h^{(2)}=$$

$$-(\bar{\delta }+2\alpha +2\bar{\beta }-\bar{\tau }{)}^{(0)}{\Phi }_{22}^{(2)},$$

$$d_3$$

$$d_4$$

$$T=t-r+2M\ln \frac{r-2M}{2M}-4M\ln \frac{r}{2M},$$

$$R=\frac{L^2}{r},$$

$$\left(C_{TT}{\partial }_T^2+C_{TR}{\partial }_T{\partial }_R+C_{RR}{\partial }_R^2+{}_s{C}_T{\partial }_T+{}_s{C}_R{\partial }_R+{}_{s}C-{}_s/\Delta \right){}_s\hat{\psi }={}_s\hat{\mathcal{S}},$$

$$C_{TT}$$

$$C_{TR}$$

$$C_{RR}$$

$${}_s{C}_T$$

$${}_s{C}_R$$

$${}_{s}C$$

$${}_s\hat{\psi }(T,R,\theta ,\phi )=\sum _{l,m}{}_s{\hat{\psi }}^{[lm]}(T,R){}_s{Y}_{lm}(\theta ,\phi ).$$

$$c_{\pm }$$

$$=\frac{8M^2{R}^2-L^4\mp L^4}{16M^2\left(L^2+2MR\right)}.$$

$$\bar{\delta }{\varphi }_1^{(1)}=\frac{R}{\sqrt{2}{L}^2}{\text{dh}}^{\prime }{\varphi }_1^{(1)}=(D-\rho +2ϵ{)}^{(0)}{\varphi }_2^{(1)}.$$

$${\varphi }_1^{(1)}=\sum _{l,m}{\varphi }_1^{[lm]}(T,R)Y_{lm},$$

$${\varphi }_1^{[lm]}$$

$$=-\sqrt{\frac{2}{l(l+1)}}\frac{R^2}{L^2}\left[\left({\hat{l}}^T{\partial }_T+{\hat{l}}^R{\partial }_R\right)+M\right]{}_{-1}{\hat{\psi }}^{[lm]},$$

$$(\Delta +2\mu {)}^{(0)}{\varphi }_1^{(1)}=(\delta +2\beta {)}^{(0)}{\varphi }_2^{(1)}=\frac{R}{\sqrt{2}{L}^2}\text{dh}{\varphi }_2^{(1)},$$

$$\Delta {\varphi }_1^{[lm]}$$

$$=\frac{R^2}{L^2}\left[2R\left(\frac{{\varphi }_1^{[lm]}}{R^2}\right)+\sqrt{\frac{l(l+1)}{2}}{}_{-1}{\hat{\psi }}^{[lm]}\right],$$

$$(\mathcal{D}-2\rho {)}^{(0)}{\varphi }_1^{(1)}=(\bar{\delta }-2\alpha {)}^{(0)}{\varphi }_0^{(1)}=\frac{R}{\sqrt{2}{L}^2}{\text{dh}}^{\prime }{\varphi }_0^{(1)}.$$

$${\varphi }_0^{(1)}=\sum _{l,m}{\varphi }_0^{[lm]}(T,R){}_1{Y}_{lm}(\theta ,\phi ),$$

$${\varphi }_0^{[lm]}$$

$$=-\sqrt{\frac{2}{l(l+1)}}\frac{R^3}{L^2}\left({\hat{l}}^T{\partial }_T+{\hat{l}}^R{\partial }_R\right)\left(\frac{{\varphi }_1^{[lm]}}{R^2}\right).$$

$$\delta {\varphi }_1^{(1)}=(\Delta +\mu {)}^{(0)}{\varphi }_0^{(1)},$$

$$\Delta {\varphi }_0^{[lm]}$$

$$=\frac{R^3}{L^2}\left[R\left(\frac{{\varphi }_0^{[lm]}}{R^3}\right)+\sqrt{\frac{l(l+1)}{2}}\left(\frac{{\varphi }_1^{[lm]}}{R^2}\right)\right].$$

$${\Delta }^2{\varphi }_0^{[lm]}=\frac{R^3}{L^2}\left[2R\frac{\Delta {\varphi }_0^{[lm]}}{R^3}+\sqrt{\frac{l(l+1)}{2}}\frac{\Delta {\varphi }_1^{[lm]}}{R^2}\right].$$

$${}_{-2}\hat{\mathcal{S}}=\sum _{l_3,m_3}{}_{-2}{\hat{\mathcal{S}}}^{[l_3{m}_3]}{}_{-2}{Y}_{l_3{m}_3}.$$

Full text

I Introduction

The birth of gravitational waves (GWs) astronomy heralds a new era of gravitational physics, deepening our understanding of fundamental interactions in the strong-field regime [1, 2, 3, 4]. GWs from merging binary black holes (BHs), observed by the ground-based LIGO-Virgo-KAGRA detectors (also the forthcoming space-based missions such as LISA [5], Taiji [6], and TianQin [7]), allow us to test the validity of general relativity [8, 9, 10, 11, 12, 13]. The BH spectroscopy [14, 15, 16] is such a program concerned with the ringdown stage, where the remnant relaxes to a stationary BH.

In the BH ringdown, perturbative radiation fields relax in a set of characteristic quasi-normal modes (QNMs) and then are dominated by inverse power-law tails, which have been extensively explored in the linear BH perturbation theory (see e.g., [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35]). Recently, however, increasing studies reveal rich phenomena arising from nonlinearities in BH ringdown, such as quadratic QNMs [36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61] and nonlinear power-law tails [62, 63, 64, 65, 66, 67, 68, 69]. Surprisingly, not only can the quadratic QNM’s amplitude be even larger than that of their linear counterpart [36, 37], but the nonlinear tails can decay more slowly than the linear Price’s law [62, 63, 64, 65, 66, 67]. These results indicate that the linear BH perturbation theory alone is insufficient to model BH ringdown accurately, necessitating the nonlinear corrections.

Transient electromagnetic (EM) events in the astrophysical environment are typically high energetic, e.g., short gamma-ray bursts (GRBs) with energy released ranging from $10^{50}$ to $10^{52}$ [70, 71], supernovae events, or long GRBs releasing energy of around $10^{54}$ ergs [72, 73]. Moreover, if a neutron star has magnetic energy up to $10^{49}$ ergs and a mass of about one solar mass $M_{\odot}$, then its EM effect could surpass its second-order gravitational effect in extreme mass ratio inspiral systems with a supermassive BH $M\gtrsim 10^{6}M_{\odot}$ [58]. Such high energy renders EM sources potentially responsible for some nonlinearities in ringdown [58, 59]. Hence, multi-messenger astronomy [74, 75], which integrates GW detection and EM observation, could help the interpretation of GW signals. For example, by analyzing the QNMs of GWs generated from transient high-energetic EM waves, one can detect solitary black holes, whose number and distribution in the Milky Way encode essential information about BH formation and the existence of primordial BHs [76].

Then a natural question arises: “Can we further the understanding or observation of the electromagnetically sourced nonlinearities in BH ringdown?” Inspiringly, the analytic calculation for an ideal dipole radially free falling into a Schwarzschild BH in [58] suggests a sign of polynomial tails in gravitational perturbations sourced by electromagnetic perturbations. However, conclusive evidence still needs to be found.

To this end, we investigate second-order gravitational perturbations sourced by first-order electromagnetic perturbations under the Schwarzschild spacetime. We find that electromagnetically sourced nonlinear tails do exist by numerically solving the Bardeen-Press-Teukolsky (BPT) equations in the Newman-Penrose (NP) formalism which can be easily extended to Kerr spacetime [77, 78].

The paper is organized as follows. In Section II, we describe our methodology to study the electromagnetically sourced nonlinear tails, including hyperboloidal foliations, reconstruction of the Maxwell scalars, and calculations of quadratic source terms. Our numerical scheme, presented in Section III, is comprised of spatial discretization, analytical mesh refinement, a first-order time reduction, and time-symmetric integration. Detailed numerical results are shown in Section IV, including discussions of initial data and mode coupling, the power law of second-order tails and its astrophysical implications, and numerical checks. Finally, in Section V, we sum up our concluding remarks.

In this paper, we use the geometric units $c=8\pi G=1$ and the conventions of Chandrasekhar [79]. For example, the metric signature is $(+,-,-,-)$ and the complex conjugate of a quantity $f$ is denoted as $\bar{f}$. However, we use Greek letters as spacetime indices instead. Furthermore, the original NP spin coefficient $\pi$ is denoted by a variant pi “$\varpi$” to avoid confusion, and a variant epsilon “$\varepsilon$” denotes a perturbative parameter distinguished from the original NP spin coefficient “$\epsilon$.” Finally, an $n$th order quantity in a perturbative expansion is denoted with a trailing superscript (n), e.g., $\phi=\phi^{(0)}+\varepsilon\phi^{(1)}+\varepsilon^{2}\phi^{(2)}+\cdots$.

All codes used to create data corresponding to the findings in this manuscript are openly available [80].

II Methodology

In this Section, we present the inhomogeneous BPT equation in horizon-penetrating hyperboloidally compactified coordinates, along with the reconstruction of Maxwell scalars via the Maxwell equations, which are required for the calculation of the quadratic source term.

II.1 Newman-Penrose formalism and hyperboloidal foliations

We consider a Schwarzschild BH with mass $M$ that is perturbed by a sourceless electromagnetic field in the form of the Maxwell scalar in the NP formalism $\phi_{2}=\varepsilon\phi^{(1)}_{2}+\mathcal{O}(\varepsilon^{2})$. Because the energy-momentum tensor of the Maxwell field $T_{\mu\nu}^{\text{EM}}$ is a quadratic form of the Maxwell scalars $\phi$, the gravitational perturbation is nonlinearly sourced at the leading order $\Psi_{4}=\varepsilon^{2}\Psi_{4}^{(2)}+\mathcal{O}(\varepsilon^{3})$. Thus, the system of equations under consideration is

[TABLE]

[TABLE]

where ${}_{s}\mathcal{T}$ is the Teukolsky operator for the NP scalar with a spin-weight $s$. The source term ${}_{-2}\mathcal{S}$ reads [81]

[TABLE]

where the Ricci terms $\mathcal{R}_{d}^{(2)}$ and $\mathcal{R}_{h}^{(2)}$ are

[TABLE]

[TABLE]

and two derivatives $d_{3}$ and $d_{4}$ are defined as

[TABLE]

As $T_{\mu\nu}^{\text{EM}}$ is traceless, one can obtain the Ricci scalars $\Phi_{mn}^{(2)}=-\phi_{m}^{(1)}\bar{\phi}_{n}^{(1)}$ from the Einstein equation $R_{\mu\nu}=-T_{\mu\nu}$.

To avoid complicated boundary condition problems [82, 83] and to extract astrophysically relevant results at future null infinity, we use horizon-penetrating hyperboloidally compactified coordinates $\{T,R,\theta,\varphi\}$ in the minimal gauge [84, 85] and a rotated Kinnersley tetrad [see eqs.˜71, 72 and 73], which is regular on the horizon [86]. Specifically, the transformation from the Schwarzschild coordinates $\{t,r,\theta,\varphi\}$ to the hyperboloidal coordinates $\{T,R,\theta,\varphi\}$ is

[TABLE]

with $L$ a constant parameter which, as well as $M$, is set to 1 in this work. Future null infinity $\mathcal{I}^{+}$ is located at $R=0$, and the future event horizon $\mathcal{H}^{+}$ is located at $R_{\text{H}}=L^{2}/2M$. Note that the hyperboloidal time coordinate $T$ approaches a retarded time $u$ when $R\to 0$ and approaches an advanced time $v$ when $R\to R_{\text{H}}$.

In our coordinates and tetrad, the BPT equation reads

[TABLE]

where the equation coefficients read

[TABLE]

[TABLE]

and ${}_{s}\not{\Delta}$ is the spin-weight $s$ Laplace-Beltrami operator on the unit two-sphere. The rescaled master function ${}_{s}\hat{\psi}$ and rescaled source term ${}_{s}\hat{\mathcal{S}}$ are related to the NP scalars (see Table 1).

Due to spherical symmetry of the Schwarzschild background, we expand ${}_{s}\hat{\psi}$, as well as ${}_{s}\hat{\mathcal{S}}$, in terms of the spin-weight spherical harmonics ${}_{s}Y_{lm}(\theta,\varphi)$, the eigenfunction with an eigenvalue $-(l-s)(l+s+1)$ of ${}_{s}\not{\Delta}$,

[TABLE]

Substituting this expansion into the master equation (9) gives a series of decoupled (1+1)-dimensional partial differential equations (PDEs) for each mode ${}_{s}\hat{\psi}^{[lm]}$, where the characteristic speeds read

[TABLE]

Note that the ingoing one $c_{-}$ vanishes at the future null infinity $(R=0)$, while the outgoing one $c_{+}$ vanishes at the future event horizon $(R=L^{2}/2M)$, for which the physical quasi-normal boundary conditions are inherently satisfied.

II.2 Reconstruction of the Maxwell scalars

Before calculating the source term ${}_{s}\hat{\mathcal{S}}$, we need to solve, via the Maxwell equations, other Maxwell scalars $\phi_{1}^{(1)}$ and $\phi_{0}^{(1)}$ (also their derivatives) consistent with the solution $\phi_{2}^{(1)}$ to eq. (1). In this and next subsection, we will drop the spin coefficients that vanish in our tetrad (see A).

The Maxwell scalar $\phi^{(1)}_{1}$ is reconstructed from the following Maxwell equation

[TABLE]

Note that the $\textnormal{\dh }^{\prime}$ operator lowers the spin weight by 1. 111For example, $\textnormal{\dh }^{\prime}{}_{s}Y_{lm}=-\sqrt{(l+s)(l-s+1)}{}_{s-1}Y_{lm}$. See [87] for detailed discussion on the ð operator. That is, the $(l,m)$ mode of $\phi_{1}^{(1)}$ can only be excited by the same $(l,m)$ mode of $\phi_{2}^{(1)}$. Thus, if we expand $\phi_{1}^{(1)}$ in terms of $Y_{lm}$

[TABLE]

and substitute this into eq. (18), then we obtain

[TABLE]

where $\hat{l}^{T}=4M^{2}$ and $\hat{l}^{R}=-(L^{2}-2MR)/2$. Note that the fall-off behavior $\phi_{1}\sim R^{2}$ is important for computing the rescaled source $\hat{\mathcal{S}}$ numerically.

The directional derivative $n^{\mu}\nabla_{\mu}\phi_{1}^{(1)}=\Delta\phi_{1}^{(1)}$ is reconstructed from

[TABLE]

i.e.,

[TABLE]

where we use the ð operator’s property of raising the spin weight by 1, i.e., $\textnormal{\dh }{}_{s}Y_{lm}=\sqrt{(l-s)(l+s+1)}\,{}_{s+1}Y_{lm}.$ The fall-off behavior of $\Delta\phi_{1}$ is also $\sim R^{2}$.

The reconstruction equation of $\phi_{0}$ reads

[TABLE]

Similarly, if we expand $\phi_{0}^{(1)}$ in terms of ${}_{1}Y_{lm}$

[TABLE]

then eq. (23) reads

[TABLE]

Finally, the reconstruction equation of $\Delta\phi_{0}$ reads

[TABLE]

i.e.,

[TABLE]

Note that the second-order derivative $\Delta^{2}\phi_{0}$ is also required to compute ${}_{-2}\mathcal{S}$, which can be obtained from eq. (27)

[TABLE]

The falloff behavior of $\phi_{0}$, $\Delta\phi_{0}$ and $\Delta^{2}\phi_{0}$ are all $\sim R^{3}$.

II.3 Quadratic source

As mentioned above, the rescaled source ${}_{-2}\hat{\mathcal{S}}$ in (9) is expanded in terms of ${}_{-2}Y_{lm}$

[TABLE]

For convenience, we divide it into three parts

[TABLE]

where222A property ${}_{s}\bar{Y}_{lm}=(-1)^{m+s}{}_{-s}{Y}_{l,-m}$ is utilized in the following derivation.

[TABLE]

[TABLE]

and

[TABLE]

Here, the mode coupling is determined by the corresponding Gaunt coefficients

[TABLE]

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

[TABLE]

This completes our derivation of the inhomogeneous BPT equation with an electromagnetic source.

III Numerical scheme

In this Section, we describe our numerical scheme to solve the (1+1)-dimensional PDEs, i.e., eq. (9) with ${}_{s}\not{\Delta}$ replaced by the eigenvalue $-(l-s)(l+s+1)$.

High-precision floating-point numbers are required to obtain accurate tail decay rates. In our numerical simulations, we employ DoubleFloat64 floating-point numbers provided by the Julia package DoubleFloats.jl [88].

III.1 Spatial discretization and analytical mesh refinement

For spatial discretization, the Chebyshev pseudo-spectral method is used in the radial $R$ direction. However, growing gradients close to the future null infinity $\mathcal{I}^{+}$ occur in the late-time profile of the master function ${}_{s}\hat{\psi}^{[lm]}$ due to the differences between the decay rates at $\mathcal{I}^{+}$ and finite radii. We use analytical mesh refinement (AnMR) to solve this problem.

Specifically, $N+1$ scaled and shifted Chebyshev-Gauss-Lobatto collocation points $\{R_{i}^{\text{Cheb}}\}$ are used within $R\in[0,R_{\text{H}}]$ during the quasi-normal ringing phase,

[TABLE]

While ${}_{s}\hat{\psi}^{[lm]}(T,R)$ enters the polynomial tail decay phase, its profile behaves like a function (see e.g., Fig. 9 in [89])

[TABLE]

whose $n$th derivative at $R=0$ is of order $|a|^{-n}$. Hence, a new coordinate $\sigma$ is introduced instead

[TABLE]

with a mesh-refinement parameter [90, 89]

[TABLE]

in which the $n$th derivative of the function (37) at $R=0$ merely goes as $(\ln|a|)^{n}$. We monitor behaviors of ${}_{s}\hat{\psi}^{[lm]}$ during evolutions and set the parameter $\kappa$ by substituting ${}_{s}\hat{\psi}^{[lm]}$ into eq. (39) once ${}_{s}\hat{\psi}^{[lm]}$ evolves into a configuration similar with eq. (37). 333${}_{-1}\hat{\psi}^{[lm]}(T,R)$ enters the tail phase typically earlier than ${}_{-2}\hat{\psi}^{[lm]}(T,R)$.

Using the AnMR (38), we map $N^{\prime}+1$ standard Chebyshev-Gauss-Lobatto collocation points $\{\sigma_{i}^{\text{Cheb}}\}$ within $\sigma\in[-1,1]$ into the interval $R\in[0,R_{\text{H}}]$. The AnMR collocation points $\{R_{i}^{\text{AnMR}}\}=\{R(\sigma_{i}^{\text{Cheb}})\}$ are more dense near $R=0$ than $\{R_{i}^{\text{Cheb}}\}$, 444A drawback of AnMR is that $\{R_{i}^{\text{AnMR}}\}$ are sparse near the other boundary $R=L^{2}/2M$, for which we do not employ AnMR from the beginning of evolutions. as shown in Fig. 1. When transforming from $\{R_{i}^{\text{Cheb}}\}$ to $\{R_{i}^{\text{AnMR}}\}$, we interpolate the field ${}_{s}\hat{\psi}^{[lm]}$ in terms of Chebyshev polynomials of the first kind via its grid values $\{{}_{s}\hat{\psi}^{[lm]}(R_{i}^{\text{Cheb}})\}$ to obtain its values at the collocation points $\{R_{i}^{\text{AnMR}}\}$. 555Altering $\{R_{i}^{\text{AnMR}}\}$ by readjusting $\kappa$ also involves such a spectral interpolation, which will significantly increase computing costs if one readjusts $\kappa$ at each time step. Thus, updating $\kappa$ before the spectral convergence is spoiled will suffice. As shown in Fig. 2, spectral coefficients $c_{n}$ of the numerical solution $\{{}_{s}\hat{\psi}^{[lm]}(T,R_{i}^{\text{AnMR}})\}$ converge faster than those of $\{{}_{s}\hat{\psi}^{[lm]}(T,R_{i}^{\text{Cheb}})\}$, which reduces computing costs to get accurate tail behaviors.

III.2 First-order time reduction

With the above spatial discretization, the spatial derivative operator $\partial_{R}$ is replaced by a differentiation matrix $\bm{D}$ [91]. Then, the original PDEs of fields ${}_{s}\hat{\psi}^{[lm]}(T,R)$ turn into a system of coupled second-order ordinary differential equations (ODEs) of the value of the fields at the $N^{(\prime)}+1$ collocation points $\{{}_{s}\hat{\psi}^{[lm]}(T,R_{i})\}$, where $N^{(\prime)},R_{i}$ refer to $N,R_{i}^{\text{Cheb}}$ and $N^{\prime},R_{i}^{\text{AnMR}}$ respectively.

To reduce the second-order ODEs, we empirically find an auxiliary variable666Introducing $\partial_{T}\hat{\psi}^{[lm]}$ as an auxiliary variable will give almost the same numerical results when the DoubleFloat64 floating-point numbers are employed. However, the time evolution with this choice ($\partial_{T}\hat{\psi}^{[lm]}$) of the auxiliary variable might be unstable at late time due to increased roundoff error, if one works with double-precision floating-point numbers.

[TABLE]

can result in stable time evolutions at late times. Then, the final ODEs to be solved read

[TABLE]

where

[TABLE]

and the $(2N^{(\prime)}+2)\times(2N^{(\prime)}+2)$ matrix ${}_{s}\bm{L}^{[l]}$ can be derived from the master equation (9)

[TABLE]

with

[TABLE]

and

[TABLE]

Here, $\textsf{diag}[C]$ denotes an $(N^{(\prime)}+1)\times(N^{(\prime)}+1)$ diagonal matrix $\bm{C}$ with diagonal elements $\bm{C}_{ii}:=C(R_{i})$.

III.3 Time-symmetric integration

We solve the ODEs system (41) via a fourth-order time-symmetric integration method based on the Hermite rule [92, 93]. The time-symmetric integration method, compared with Runge-Kutta methods of the same order, is free of Courant limit, introduces smaller truncation error, and preserves Noether charges over long time periods [92, 93], which makes the time-symmetric integration method ideal for long time numerical evolution in BH perturbation theory.

We will specify the time-symmetric evolution method to solve an ODEs system of the following form

[TABLE]

where $\bm{L}$ is a matrix independent of $T$ and $\vec{u}$, and $\vec{\mathcal{S}}(T)$ is a known vector function. The problem can be converted to computing a series of integrals

[TABLE]

which are approximated by the Hermite rule

[TABLE]

with $\delta T=T_{n+1}-T_{n}$.777Numerical results shown in next section IV are obtained with $\delta T=2^{-10}$. Here, the time derivative is denoted by an over-dot, and the value of a vector function $\vec{f}(T)$ at a certain time $T_{i}$ is denoted by a subscript $\vec{f}(T_{i}):=\vec{f}_{i}$. Note that eq. (47) with the Hermite rule (48) is generally an implicit scheme to solve $\vec{u}_{n+1}$ for given $\vec{u}_{n}$. However, we will show that one can construct an explicit evolution scheme for ODEs of the form (46).

The time derivatives in (48) can be determined by

[TABLE]

Then, (48) becomes

[TABLE]

Thus, we obtain the time-symmetric evolution scheme for the ODEs (46)

[TABLE]

Finally, we rewrite the scheme (51) as a form of matrix-vector multiplication and addition to reduce round-off error at each time step

[TABLE]

For the ODEs of $\vec{v}^{[lm]}$ in (41), eq. (52) is obviously an explicit time evolution scheme. Then, for the ODE of $\vec{u}^{[lm]}$, the source term

[TABLE]

is not an explicit vector function of $T$, but $\vec{\mathcal{S}}_{n+1}(\vec{v}^{[l_{1}m_{1}]}_{n+1};\vec{v}^{[l_{2}m_{2}]}_{n+1})$ can be determined by $\vec{v}_{n}^{[lm]}$ explicitly. In addition, the time derivative of the source term

[TABLE]

can also be determined by $\vec{v}_{n}^{[lm]}$ explicitly. Hence, we obtain an explicit time evolution scheme for the whole ODEs system (41).

IV Numerical results

In this section, we discuss choices of initial data and mode coupling first. Then we present our main result, i.e., the power law of second-order tails, followed by its implications for multi-messenger observations. Finally, we describe numerical checks of our numerics.

IV.1 Initial data and mode coupling

It is known that the decay rates of linear tails depend on initial data, e.g., the angular profile, the time derivative of the field (stationary or not) or the radial profile (compact support or not) [94, 95, 96, 32].

For the first-order electromagnetic perturbation $\phi_{2}^{(1)}$, we employ the initial data that are commonly used in the literature.

Specifically, the angular profile is described by spin-weighted spherical harmonics with a single multipole $l$ for simplicity, i.e.,

[TABLE]

where we assume ${}_{-1}\hat{\psi}^{[lm]}\equiv{}_{-1}\hat{\psi}^{[l]}$ for different azimuthal numbers $m$.

For the time derivative of the field, we only consider the so-called stationary initial data, i.e.,

[TABLE]

Note that there is no regular solution $\psi(R)$ to an ODE $C_{TR}\partial_{R}\psi+{}_{s}C_{T}\psi=0$, for which the auxiliary variable ${}_{-1}P^{[l]}$ (40) does not vanish initially. Consequently, the power law of ${}_{-1}\hat{\psi}^{[l]}$ is determined by the late-time behavior of the retarded Green’s function $\mathcal{G}(T-T^{\prime},R,R^{\prime})$, following from the homogeneous solution888We also verified that the power law, with the initial data such that ${}_{s}P^{[l]}(T=0,R)\equiv 0$, is one less than that obtained with the so-called stationary initial data (56). This indicates that it is governed by the late-time behavior of $\mathcal{G}_{,T^{\prime}}$. In this sense, it is more appropriate to refer to ${}_{s}P^{[l]}(T=0,R)\equiv 0$ as “stationary initial data”.

[TABLE]

to the (1+1)-dimensional PDE (9) with ${}_{s}\not{\Delta}$ being replaced by the eigenvalue $-(l-s)(l+s+1)$, where $\mathcal{G}(T-T^{\prime},R,R^{\prime})$ is a solution to the PDE with a Dirac-delta source term $\delta(T-T^{\prime})\delta(R-R^{\prime})$ and is subject to the condition that $\mathcal{G}(T-T^{\prime},R,R^{\prime})=0$ for $T<T^{\prime}$.

The radial profile is given by

[TABLE]

with

[TABLE]

The results shown in this work are obtained with $R_{c}=R_{\text{H}}/2$ and $w=R_{\text{H}}/10$.

In this paper, we refer to the initial data (58) with

• •

$(k,b)=(1,0)$ as compact support, 999Even if a Gaussian wave packet is not strictly compact support and the boundary value $G(R=0;1,0)=\mathrm{e}^{-(R_{0}/w)^{2}}=\mathrm{e}^{-25}\simeq 1.4\times 10^{-11}$ is much larger than the round-off error of the DoubleFloat64 floating-point numbers $2^{-104}\simeq 4.9\times 10^{-32}$, the decay rates of linear tails for compact support initial data are successfully reproduced (see Fig. 8).

• •

$(k,b)=(0,1)$ as non-compact support.

These types of initial data are usually employed to study linear tails. Besides, we also take initial data employed in the literature concerned with nonlinear tails [62, 63], i.e., ingoing data $(\Delta{}_{-1}\hat{\psi}^{[l]}=0$),

[TABLE]

and outgoing data $(D{}_{-1}\hat{\psi}^{[l]}=0)$,

[TABLE]

The above types of initial data, plus zero initial data

[TABLE]

are employed for the second-order gravitational perturbation $\Psi_{4}^{(2)}$.

The pure multipole initial data (55) simplify the mode coupling, as the coupling of $l$ and $m$ in the source terms eqs.˜31, 32 and 33 can be separated. Different mode coupling channels $(m_{1},m_{2})$ such that $m_{1}-m_{2}=m_{3}$ just differ by a constant 101010Thus, we do not find exceptions, e.g., the coupling of $l=m=2$ to $l=-m=2$ reported in [62], to the power law (69) described in Sec. IV.2.

[TABLE]

in the source terms. Thus, we only need to consider mode coupling of $(l_{1}=l,l_{2}=l)$ such that $l_{3}=2,3,\cdots,2l$ for the pure multipole case.

We checked that the power law of second-order tails does not depend on the mode coupling channels $(l_{1},l_{2})\to l_{3}$ and that it also does not depend on the initial data of ${}_{-1}\hat{\psi}^{[l]}$ and ${}_{-2}\hat{\psi}^{[l_{3}]}$ both when $l_{3}\geq 4$.

IV.2 The power law of second-order tails

The power law decay of a field $\psi$ is monitored by local power index (LPI)

[TABLE]

such that the LPI is a constant $p$ if the field $\psi\sim T^{p}$.

Our main results are shown in the Fig. 3. The LPIs of ${}_{-2}\hat{\psi}^{[l]}$ for $l=4,5,6$ suggest a power law of the form

[TABLE]

and this power law, as mentioned above, is independent of the initial data, indicating that in the case $l\geq 4$ the homogeneous part (57) of ${}_{-2}\hat{\psi}^{[l]}$ cannot dominate over its inhomogeneous part

[TABLE]

This result supports a dominant position of nonlinear tails over their linear counterparts, in line with [62, 63, 64]. To see this more clearly, we solve both the inhomogeneous and homogeneous (namely ${}_{-2}\hat{\mathcal{S}}=0$) BPT equations (9) under the same initial data. As shown in Fig. 4, the source-driven perturbation overshadows the source-free one at late times, thereby determining how a perturbed BH approaches a stationary state and playing a critical role in fundamental theoretical issues, e.g., strong cosmic censorship. [97, 98, 99, 100, 101, 102, 103, 104].

When $l=2$ and $3$, however, the power law of second-order perturbation ${}_{-2}\hat{\psi}^{[l]}$ depends on the initial data of ${}_{-1}\hat{\psi}$ and ${}_{-2}\hat{\psi}^{[l]}$ both. Interested in the effect of source term, we focus on results with the zero initial data for ${}_{-2}\hat{\psi}^{[l]}$, as shown in Fig. 5. We find that for ${}_{-1}\hat{\psi}$ with the compact support initial data, ${}_{-2}\hat{\psi}^{[l]}$ decays as $T^{-2l-3}$ at fixed spatial position and $T^{-l-4}$ at future null infinity $\mathcal{I}^{+}$, which are the same as the linear Price’s law. However, the non-compact support initial data of ${}_{-1}\hat{\psi}$, forming an extended source at the beginning, lead to a slower tail that conforms to the power law (67).

Recent analytical results [62, 63, 64], formulated in the Regge-Wheeler-Zerilli (RWZ) formalism, indicate that the behavior of the source term $Q$ at large radii $Q\sim 1/r^{\beta}$ plays a key role in determining the power law of nonlinear tails. It was shown that the inhomogeneous part of Regge-Wheeler/Zerilli master function $\Psi$, induced by a compact outgoing source term $Q\sim 1/r^{\beta}$, decays as

[TABLE]

at fixed spatial position [63]. 111111The homogeneous part of $\Psi$, concerned with the initial data, was ignored in [63]. Note that the curvature perturbation $\Psi_{4}$ is related to the metric perturbation $\Psi$ by the Chandrasekhar transformation [79, 105]. The decay rates of $\Psi_{4}$ and $\Psi$ are the same near the black hole, but the former is two less than the later at null infinity [30]. Our numerical result, i.e., the power law $T^{-2l-2}$ at fixed spatial position when $l\geq 4$ for compact support ${}_{-1}\hat{\psi}(T=0,R)$, coincides with the analytical prediction (69) for the case $\beta=6$. 121212See [58] for a derivation of the electromagnetic source in RWZ formalism. In addition, to our knowledge, the dependence of the second-order power law on the initial data of the first-order quantities had not been discovered before. However, the power law $T^{-l-3}$ at future null infinity $\mathcal{I}^{+}$, or $T^{-l-1}$ in the metric perturbation picture, is inconsistent with the analytical prediction

[TABLE]

presented in [62], 131313But our finding is consistent with their numerical result of second-order tails for the self-coupling of $l=m=2$. which means that their analytical approach omits some features [63].

IV.3 Implications for multi-messenger observations

Despite still being strongly suppressed by the QNMs and neglected by any near-future observatories, the nonlinearly source-driven tails can be significantly amplified for binary BHs with high eccentricities, prospectively reaching the detection threshold of upcoming detectors [106, 107, 108, 109]. Hence, it would be valuable to discuss the implications of our results for multi-messenger observations in this subsection.

Foremost, the electromagnetically sourced nonlinear tails differ from gravitationally sourced nonlinear tails when $l=2,3$, because the latter, with a source $Q\sim 1/r^{2}$ [110, 111], decay as $T^{-2l-2}$ at fixed spatial position when $l\geq 2$ even if the first-order quantities are initially compact support [62]141414We conjecture that the gravitationally sourced nonlinear tails decay as $u^{-l-1}$ in the metric perturbation picture, albeit this was not explicitly reported in [62].. This difference, combined with the analysis of quadratic QNM [76], could help to identify astrophysical origins of GWs in the multi-messenger observations and offers a novel and complementary mechanism to detect black holes in the Milky Way.

As depicted in Fig. 6, the LPIs of second-order tails show a similar distance-dependence with linear tails [112, 113, 114, 29], i.e., the decay rates $p(T,R)$ at finite distance vary monotonously between that at null infinity $\mathcal{I}^{+}$ and that at event horizon $\mathcal{H}^{+}$. However, it is the decay rate at null infinity $\mathcal{I}^{+}$ that is relevant for astronomical observations, due to extremely distant astronomical distances [29].

As shown in Fig. 7, both the evolutions of ${}_{-1}\hat{\psi}$ and ${}_{-2}\hat{\psi}$ consist of three distinct phases: initial transient, quasi-normal ringing, and polynomial tail decay. During the quasi-normal ringing phase, we extract the quasi-normal frequencies with the matrix pencil method [115]. The spectrum of the first-order quantities ${}_{-1}\hat{\psi}$ matches the prediction of linear BH perturbation theory. Besides the linear QNMs, the quadratic QNMs with frequencies twice that of ${}_{-1}\hat{\psi}$’s linear mode are also found, which come from the free propagation part $\mathcal{G}_{F}$ of the Green’s function $\mathcal{G}$ in (68) [116]. Moreover, the peaks of GW waveform typically arrive at $\mathcal{I}^{+}$ later than those of their parent EM waveform in our simulations, indicating that EM events could be a forecast of their offspring GW events. These timing and spectral signatures provide practical search priors for multi-messenger campaigns: the EM peak time can trigger time-gated and stacked searches for long-lived ringdown tails in GW data, thereby lowering detection thresholds [117, 118, 119].

IV.4 Numerical checks

Our numerics are checked by successful reproduction of the power law for the first-order Maxwell scalars $\phi_{2}^{(1)},\phi_{1}^{(1)}$ and $\phi_{0}^{(1)}$ (see Fig. 8), as well as spectral accuracy of evolutionary variables (see Fig. 10).

An accurate reconstruction of $\phi_{0}^{(1)}$ requires quite a few collocation points, due to a LPI splitting151515The LPI splitting is that a field with positive spin weight decays with three different rates at $\mathcal{H}^{+},\mathcal{I}^{+}$ and the bulk [27, 25]. and a greater difference of LPIs between $\mathcal{I}^{+}$ and finite radii. However, we find that even if the number of collocation points is not large enough to reproduce the power law of $\phi_{0}^{(1)}$ at late times, the power law of ${}_{-2}\hat{\psi}$ is not affected.

The rescaled source term ${}_{-2}\hat{\mathcal{S}}^{[lm]}$ (30), shown in Fig. 9, resembles a Dirac-delta distribution in the compactified $R$ coordinate, whose peak approaches $R=0$ at a coordinate speed much less than the outgoing characteristic $c_{+}(R=0)=-1/8$ (see Table 2). This Dirac-delta-like profile can severely spoil the convergence of spectral methods. Nevertheless, with the AnMR collocation points (38), spectral coefficients of the source term ${}_{-2}\hat{\mathcal{S}}$, as well as those of the evolutionary variables ${}_{-1}\hat{\psi}$ and ${}_{-2}\hat{\psi}$, indicate an exponential convergence rate of the spectral expansions (see Fig. 10). In addition, no obvious aliasing is found in the nonlinear source term ${}_{-2}\hat{\mathcal{S}}$, for which we do not apply any filter in our evolutions.

V Conclusion and outlook

In this work, we investigate second-order gravitational perturbations under the Schwarzschild spacetime that is perturbed by a transient electromagnetic wave packet. Electromagnetically sourced nonlinear tails are found by numerically solving the inhomogeneous BPT equation. Accurate numerical results are obtained by efficient numerical methods, including AnMR and time-symmetric integration. When $l\geq 4$, the nonlinear tails of curvature perturbations decay as $t^{-2l-2}$ at fixed spatial position, slower than the linear Price’s law $t^{-2l-3}$, and $u^{-l-3}$ at null infinity, which are independent of mode coupling and initial data. While $l=2$ and $3$, given that real matter is compact, the nonlinear tails decay at the same rate as the linear tails under the compact support initial data. These tail behaviors are a bit different from the nonlinear tails sourced by first-order gravitational perturbations when $l=2$ and $3$, because the latter still decay as $t^{-2l-2}$ [62]. This difference could be used to identify astrophysical origins of GWs in the multi-messenger observations, providing a novel and complementary mechanism to detect black holes in the Milky Way.

One direction for future efforts is to extend the current analysis to Kerr spacetime, where nonlinear tails may reveal richer phenomena, e.g., mode projections or intermediate “splitting” of decay rates [121]. Given the rise in computing costs, it is essential to benchmark the efficiency of various numerical methods. In addition to the AnMR and time-symmetric integration methods used in this article, other approaches, such as the multi-domain spectral method, fully spectral method [89], and spectral decomposition method [122], are also suitable for studying polynomial tails.

Masked by the stronger contribution of QNMs to ringdown waveform modeling, tails have generally been regarded as negligible for near-future GW detections. However, the tail effect could be significantly amplified in inspirals with high eccentricities, potentially detectable by upcoming detectors [106, 107, 108, 109]. The power law of these eccentricity-induced tails seems to exhibit competition between source-free tails and source-driven ones [109]. Thus, another direction for future efforts is to investigate the electromagnetically sourced nonlinear tails in such an extreme mass ratio inspiral system comprised of a neutron star and a supermassive BH.

Acknowledgements.

This work is partly supported by the National Key Research and Development Program of China (Grant No. 2021YFC2203001). This work is supported in part by the National Natural Science Foundation of China under Grants No. 12035016, No. 12075026, No. 12275350, No. 12375048, No. 12375058, No. 12361141825, No. 12447182 and No. 12575047.

Appendix A Tetrad and spin coefficients

The tetrad employed in this work reads

[TABLE]

[TABLE]

[TABLE]

The unperturbed tetrad vectors (71) and (72) are chosen along the repeated principal null directions of the Weyl tensor, so that one can obtain

[TABLE]

and

[TABLE]

following from the Goldberg-Sachs theorem [123].

The only non-zero Weyl scalar is

[TABLE]

and the non-zero spin coefficients are

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

Finally, the vanishing spin coefficients are

[TABLE]