^†^†thanks: Corresponding author

More Nonlinearities?
II. A Short Guide to First- and Second-Order Electromagnetic Perturbations in the Schwarzschild background.

Fawzi Aly¹

mabbasal[AT]buffalo.edu Dejan Stojkovic¹

ds77[AT]buffalo.edu ¹HEPCOS, Physics Department, SUNY at Buffalo, Buffalo, New York, USA

Abstract

We study second-order electromagnetic perturbations in the Schwarzschild background and derive the effective source terms for Regge-Wheeler equation which are quadratic in first-order gravitational and electromagnetic perturbations. In addition to the induced mixed quadratic modes, we find that linear gravitational modes are also excited, with amplitudes dependent on the electromagnetic potential. A toy model involving a Dirac delta function potential demonstrates mixing of linear gravitational and electromagnetic perturbations with frequencies $\omega^{(1)}$ and $\Omega^{(1)}$ , resulting in the second-order QNM mixing in the electromagnetic field at $\Omega^{(2)}=\Omega^{(1)}+\omega^{(1)}$ . This complements prior work in [1] on the second-order gravitational perturbation mixing and highlights potential applications in multi-messenger astrophysics for systems observed by LIGO and upcoming LISA. We also study first-order perturbations due to a point charge and show it could be reduced to a one-dimensional path integral. Within the toy model, we investigate the first-order electromagnetic perturbation due to a radially free-falling single charge $q$ and radial dipole moment $p=q\eta$ , employing semi-analytical and numerical methods. For the dipole case, we show that the QNM perturbation is excited with a nearly constant amplitude. Future work will focus on incorporating mixing in more realistic potentials and exploring numerical approach in the context of rotating spacetimes.

Introduction

The study of quasinormal modes (QNMs) is fundamental to understanding how black holes (BHs) respond to perturbations. In general relativity (GR), QNMs play an even more critical role, as they are only parametrized by BH macroscopic properties: mass $M$ , angular momentum $J$ , and charge $Q$ —making them characteristic signatures of these enigmatic objects [2, 3, 4, 5]. In this context, the black hole spectroscopy program was established. While observing a single QNM provides a direct measurement of these parameters, detecting two or more QNMs offers a valuable consistency check. This enables testing the validity limits of black hole perturbation theory (BHPT) in modeling the late-stage ringdown phase [2, 6, 7].

While gravitational QNMs have garnered significant attention—particularly following the groundbreaking gravitational wave (GW) detections by LIGO and Virgo [8, 9]—electromagnetic QNMs remain a relatively under-explored area with intriguing theoretical implications, and potential astrophysical relevance. This is particularly true in studying neutron stars (NSs), and charged BH properties. Electromagnetic QNMs encode some information about how electromagnetic fields evolve and decay around vacuum BH spacetimes, for example like those considered as perturbations of Kerr type remnants after magnetized mergers [10]. Understanding these modes is essential for better interpretation of electromagnetic radiation spectrum emitted during astrophysical events, including black hole-neutron star (NSBH) and binary neutron star (BNS) mergers, gravitational collapse of unstable magnetars, magnetar flares, and accretion processes near BHs.

In the context of multi-messenger astronomy [11, 12, 13, 14], the simultaneous detection of gravitational and electromagnetic signals offers a unique opportunity to probe the strong-field regime and more thoroughly examine astrophysical events where electromagnetic effects are switched on. The gravitational and electromagnetic signals are expected to provide both complementary and overlapping information about the same astrophysical event, allowing for cross-checks. Although, in principle, understanding electromagnetic QNM spectra might lead to more accurate extraction of information in the aforementioned systems [15, 5] and enable further testing of GR and its coupling with electromagnetism [1], the state of the art in any realistic signal processing for BNS is to ignore the QNM-component. This is partially due to the complexity of magnetized mergers and lack of full description of how electromagnetic field behaves in the vicinity of the ejecta and in plasma environment.

In the context of numerical simulations, scenarios involving mergers of charged binary BHs or BNS, magnetar collapse, and NSBH systems without significant disruption offer a higher probability of resolving electromagnetic QNM behavior. Current and future simulations of charged or magnetized systems could shed light on the remnant’s ringdown for these scenarios, especially following recent attempts to identify electromagnetic QNMs, (i.e. in [16, 17]). Given that magnetic fields in astrophysical BNS mergers can reach up to $B\sim 10^{12}\,\text{T}$ [18], it is intriguing to explore the role of electromagnetic fields in these extreme conditions.

In [1], we explored whether electromagnetic fields can leave imprints on gravitational radiation, focusing more on extreme mass ratio inspiral (EMRI) scenarios [19, 20, 21, 22]. We reported that it was plausible for electromagnetic effects to compete with the second-order gravitational perturbations in a NSBH merger. With a lower limit estimate for astrophysical NS of mass $2M_{\odot}$ , and with a magnetic field of order $\sim 10^{7}\,\text{T}$ , in an inspiral around a supermassive black hole of mass $10^{6}M_{\odot}$ . However, in the context of stellar-mass mergers, the situation is more challenging in the absence of a proper way to define perturbation parameters, so perhaps further analytical and numerical investigations are needed.

However, since gravitational effects are expected to dominate in astrophysics, this complementary work aims to generically investigate the reverse question: Can we detect gravitational signatures in the electromagnetic radiation? Such phenomena should be important in analyzing and modeling electromagnetic counterparts of magnetized collapse or mergers in future.

The work in this paper addresses several questions. First, we provide a brief review of the Regge-Wheeler formalism for electromagnetic perturbations. We then study second-order electromagnetic perturbations, where the electromagnetic-gravitational interaction, with its first-order modes, acts as an effective source term. We solve this problem for a sample source term for the purpose of illustration. Following this, we go back to first order perturbation, and solve for the electromagnetic perturbations induced by a point charge and also an ideal dipole, after approximating the Regge-Wheeler potential with a Dirac delta function. The choice of an ideal dipole serves as an analogy for astrophysical systems with relatively small macroscopic charges. Finally, we report our numerical solutions for the perturbation due to the dipole.

The structure of this paper is as follows. In Section I, we outline the electromagnetic perturbation in Schwarzschild coordinates. Next, in Section II, we formulate the second-order perturbation and derive the source terms for the second-order master Regge-Wheeler formalism for spin $s=1$ . By the end of this section, we solve for a sample source using the Dirac delta function approximation to illustrate electromagnetic-gravitational mixing, similar to [1]. Then, we return to first-order perturbations in Section III and show that, in the time domain, the inhomogeneous solution due to a point charge can be reduced to a one-dimensional path integral, and we discuss several examples. In Sections IV and V, we investigate the first-order electromagnetic perturbation due to a radially free-falling point charge and then a radial dipole in the Dirac delta approximation using a semi-analytical approach. Next, we also solve for the perturbation due to the dipole numerically and comment on the solution in Section VI. Finally, we discuss the results of our analysis and potential directions for future work in Section VII. Throughout this work, we use natural units where $G=c=1$ .

I Preliminary

In this section, we briefly review the spherical decomposition of Maxwell’s equations on a spherically symmetric spacetime, such as Schwarzschild. We discuss how the Faraday tensor $F$ can be expressed in terms of perturbation quantities $\Phi_{lm}$ , as well as the current multipoles $j_{lm}$ , and then obtain the corresponding electromagnetic stress-energy tensor $T$ . These few steps will allow us to solve the toy model introduced in IV.

I.1 Maxwell equations

To derive the Regge-Wheeler equations for electromagnetic perturbations in a Schwarzschild spacetime [10, 23, 24, 25], the vector potential $A_{\mu}$ and the current vector $J_{\mu}$ are decomposed into even and odd components, each of which is expanded in terms of 4-vector spherical harmonics, as outlined in Appendix D. Consequently, the Faraday tensor $F_{\mu\nu}$ , and therefore Maxwell’s equations, are similarly divided into sectors of even and odd parity. The first order even and odd electromagnetic perturbations

		$\displaystyle{}_{e}\Phi^{(1)}_{lm}(t,x)=\frac{r^{2}}{\lambda}\left(\frac{\partial h_{lm}}{\partial t}-\frac{\partial g_{lm}}{\partial r}\right)$		(1)
		$\displaystyle{}_{o}\Phi^{(1)}_{lm}(t,x)=a_{lm}$		(1)

are governed by the same potential ${}_{s=1}V_{RW}=\frac{\lambda f(r)}{r^{2}}$ , where $\lambda=l(l+1)$ and $l$ represents the angular number, but differ in how the source term ${}_{J}\mathcal{S}^{(1)}_{o,e\,lm}(t,x)$ is derived from $J_{\mu}$ . The governing equation for these perturbations is expressed as:

		$\displaystyle\mathcal{L}_{EM}:=\partial_{xx}-\partial_{tt}-\,{}_{s=1}V_{RW}(x),$		(2)
		$\displaystyle\mathcal{L}_{EM}\,{}_{o,e}\Phi^{(1)}_{lm}(t,x)={}_{J}\mathcal{S}^{(1)}_{o,e\,lm}(t,x),$		(2)

where $x$ is the tortoise radial coordinate, related to the areal radial coordinate through $dr=f(r)dx$ , while $t$ is the time measured by a clock asymptotically far from the Schwarzschild black hole. In addition, ${}_{J}\mathcal{S}^{(1)}_{o,e\,lm}(t,x)$ is formulated in terms of the current vector multipoles $j_{(i)\,lm}$ ,where $i=1,2,3,4$ , with an example provided in section III.

This ODE governing the perturbations has two regular singularities (rank = 1) at $r=0$ and $r=2M$ , along with an irregular singularity (rank = 2) as $r\rightarrow\infty$ . This structure classifies it as a confluent Heun ODE [26, 19], which can be solved using confluent Heun functions $H_{C}$ or, in some cases, a narrower subset of confluent Heun polynomials [27]. The QNMs boundary conditions are imposed on the solution, and together with two initial conditions, the full solution can be specified, the reader can find all the details in Section II.A of [1].

I.2 Faraday tensor

After successfully separating and decoupling the system into its two degrees of freedom ${}_{o,e}\Phi^{(1)}_{lm}$ , the Faraday tensor $F_{\mu\nu}$ can now be expressed in terms ${}_{o,e}\Phi^{(1)}_{lm}$ and the current multipoles $j_{(i)\,lm}(t,x)$ for $l>0$ . However, before proceeding, let’s first address the monopole term $l=0$ of $F_{\mu\nu}$ .

I.2.1 Monopole Case

Given the multipole decomposition, the monopole case can be solved separately, as the potential vanishes. The coupled equations for the scalar field $\Phi_{00}$ are:

	$\displaystyle\frac{\partial\Phi_{00}}{\partial r}$	$\displaystyle=\frac{r^{2}\,j_{(1)\,00}}{f(r)}$		(3)
	$\displaystyle\frac{\partial\Phi_{00}}{\partial t}$	$\displaystyle=r^{2}f(r)\,j_{(2)\,00}.$		(3)

When the perturbation is due to a collection of discrete charges, the solution reduces to a summation over the individual charges $q_{i}$ , each multiplied by a step function dependent on their $r-r_{p_{i}}(t)$ argument

\Phi_{00}=\sum_{i}q_{i}\,\Theta[r-r_{p_{i}}(t)],

(4)

where $\Theta$ denotes the Heaviside step function and $r_{p_{i}}(t)$ is the radial coordinate of the $i$ -th charge parameterized by $t$ . This solution should resemble a static perturbation in the background spacetime at radii greater than the position of the charge distribution. Such a perturbation effectively charges the spacetime, transforming the region outside the charge distribution into a member of the Reissner–Nordström family of spacetimes[28, 29, 30]:

F_{\mu\nu 00}=\begin{pmatrix}0&\frac{Q_{\text{tot}}(t,r)}{r^{2}}&0&0\\ **&0&0&0\\ 0&0&0&0\\ 0&0&0&0\end{pmatrix}

(5)

where $Q_{\text{tot}}=\sum_{i}q_{i}\,\Theta[r-r_{p_{i}}(t)]$ , somewhat analogous to Gauss’s law in electrostatics in the exterior region. It is important to note that if the source includes no monopole contribution, this perturbation vanishes, as in the case of dipole moments when modeling astrophysical systems with relatively small macroscopic charge relative to the higher moments. Subsequently, any electromagnetic moments would be radiated away on a time scale comparable to that of the merger and post-merger phases ¹¹1We are assuming that the remnant spacetime exterior is vacuum. However, if there were accretion disks or plasma present, the situation might differ, and the same conclusions may not hold.. During the process when the BHPT is expected to hold, the ${}_{o,e}\Phi^{(1)}_{lm}$ should describe this radiation. Furthermore, if the electromagnetic moments are sufficiently large, this radiation could be substantial, albeit subdominant to gravitational radiation. In such a case, coupling effects between the two radiations must be considered, which is the purpose of this work and its companion [1].

I.2.2 Multipoles $l>0$

For higher multipoles $l>0$ , the odd components of the Maxwell tensor ${}_{o}F_{\mu\nu lm}$ in the multipole expansion take the following forms for a generic source term:

{}_{o}F_{lm}=\begin{pmatrix}0&0&-\csc\theta Y_{lm}^{\prime}\frac{\partial{}_{o}\Phi_{lm}}{\partial t}&-\csc\theta Y_{lm}^{\prime}\frac{\partial{}_{o}\Phi_{lm}}{\partial r}\\ 0&0&\sin\theta\dot{Y}_{lm}\frac{\partial{}_{o}\Phi_{lm}}{\partial t}&\sin\theta\dot{Y}_{lm}\frac{\partial{}_{o}\Phi_{lm}}{\partial r}\\ **&**&0&-\lambda\,\sin\theta Y_{lm}\,{}_{o}\Phi_{lm}\\ **&**&**&0\end{pmatrix}.

(6)

The even part of the Maxwell tensor could be divided into two contributions, which can be expressed as:

{}_{e}F_{\mu\nu lm}={}_{e1}F_{\mu\nu lm}+{}_{e2}F_{\mu\nu lm},

where the first contribution, ${}_{e1}F_{\mu\nu lm}$ , is written in terms of the even scalar perturbation ${}_{e}\Phi$ , while the second contribution, ${}_{e2}F_{\mu\nu lm}$ , is expressed in terms of the $j_{(1)\,lm}$ and $j_{(2)\,lm}$ multipole components of the even current vector.

{}_{e1}F_{lm}=\begin{pmatrix}0&\frac{\lambda}{r^{2}}Y_{lm}{}_{e}\Phi_{lm}&\dot{Y}_{lm}f\frac{\partial{}_{e}\Phi_{lm}}{\partial r}&Y^{\prime}_{lm}f\frac{\partial{}_{e}\Phi_{lm}}{\partial r}\\ **&0&\dot{Y}_{lm}\frac{\partial{}_{e}\Phi_{lm}}{\partial t}\frac{1}{f}&Y^{\prime}_{lm}\frac{\partial{}_{e}\Phi_{lm}}{\partial t}\frac{1}{f}\\ **&**&0&0\\ **&**&**&0\end{pmatrix},

(7)

and,

{}_{e2}F_{lm}=\begin{pmatrix}0&0&-\dot{Y}_{lm}\frac{r^{2}j_{(1)lm}}{\lambda}&-Y^{\prime}_{lm}\frac{r^{2}j_{(1)lm}}{\lambda}\\ **&0&-\dot{Y}_{lm}\frac{r^{2}j_{(2)lm}}{\lambda}&-Y^{\prime}_{lm}\frac{r^{2}j_{(2)lm}}{\lambda}\\ **&**&0&0\\ **&**&**&0\end{pmatrix}.

(8)

Finally we can express the Maxwell tensor $F_{\mu\nu}$ as

	$\displaystyle F_{\mu\nu}$	$\displaystyle=F_{\mu\nu\,00}+\sum_{l=1}^{\infty}\sum_{m=-l}^{l}{}_{e2}F_{\mu\nu\,lm}$		(9)
		$\displaystyle+\sum_{l=1}^{\infty}\sum_{m=-l}^{l}\left({}_{e1}F_{\mu\nu\,lm}+{}_{o}F_{\mu\nu\,lm}\right)$		(9)

So far, we haven’t made any specific assumptions about the perturbation source, yet we can make some general statements about the behavior of the Maxwell tensor in its current form.

•

While solving the system (3), it becomes clear that the monopole $\Phi_{00}$ neither exhibits radiative behavior nor satisfies the typical boundary conditions imposed on QNMs. As $t\rightarrow\infty$ on the clock of an asymptotic observer, the charge will freeze near the horizon, all radiative components (QNMs,tails,..etc) will dissipate. Ultimately, the system will asymptotically approach the Reissner-Nordström solution.
•

The term ${}_{e2}F_{\mu\nu\,lm}$ contains more information about the behavior of the electromagnetic field in the neighborhood of the charge, typically localized if the source $J_{\mu}$ behaves as expected in astrophysical cases. Consequently, when calculating the stress-energy tensor $T_{\mu\nu}$ , this component needs to be regularized for discrete charge distributions to avoid infinite energy density. We shall refer to ${}_{e2}F_{\mu\nu\,lm}$ as “internal” field from now on.
•

The components ${}_{e1}F_{\mu\nu\,lm}$ and ${}_{o}F_{\mu\nu\,lm}$ are radiative in nature and represent the QNM behavior of the electromagnetic field. These components are obtained by solving the inhomogeneous Regge-Wheeler equations (2), where $j_{(1)\,lm}$ and $j_{(2)\,lm}$ source the even component, and $j_{(4)\,lm}$ sources the odd component. The initial conditions on the Electromagnetic field will add non-trivial constraints on ${}_{o,e}\Phi^{(1)}_{lm}$ , this will be highly relevant if we are focusing more on stellar-mass mergers (i.e. NSBH or BNS mergers) ringdown phase.

I.3 Stress-energy tensor

Given the Faraday tensor $F_{\mu\nu}$ , we can proceed to calculate the electromagnetic stress-energy tensor $T_{EM}^{\mu\nu}$ , which describes the distribution and flow of electromagnetic energy and momentum in spacetime. In spherical coordinates, $T_{EM}^{\mu\nu}$ takes the following standard form:

T_{EM}^{\mu\nu}=\begin{pmatrix}\epsilon&S_{r}&S_{\theta}&S_{\phi}\\ S_{r}&P_{r}&\sigma_{r\theta}&\sigma_{r\phi}\\ S_{\theta}&\sigma_{r\theta}&P_{\theta}&\sigma_{\theta\phi}\\ S_{\phi}&\sigma_{r\phi}&\sigma_{\theta\phi}&P_{\phi},\\ \end{pmatrix}

(10)

Assuming that the components of the stress-energy tensor retain their familiar physical interpretations as they do in Minkowski spacetime, $\epsilon$ represents the energy density or electromagnetic energy per unit volume; $S_{r}$ , $S_{\theta}$ , and $S_{\phi}$ are the components of the Poynting vector, indicating the flow of electromagnetic energy in the radial, polar, and azimuthal directions, respectively; $P_{r}$ , $P_{\theta}$ , and $P_{\phi}$ correspond to the pressure in these directions; and the off-diagonal elements $\sigma_{r\theta}$ , $\sigma_{r\phi}$ , and $\sigma_{\theta\phi}$ describe the shear stresses in the electromagnetic field.

As discussed earlier, we will exclude the contribution from ${}_{e1}F_{\mu\nu\,lm}$ when calculating $T_{EM}^{\mu\nu}$ . Since the full expressions are lengthy, for completeness and illustration purposes, we will just write down the expressions for $\epsilon$ and $S_{r}$ while rest will be provided in a companion Mathematica Notebook:

$\displaystyle\epsilon$	$\displaystyle=\sum_{l,m}\sum_{l^{\prime},m^{\prime}}\frac{A_{lml^{\prime}m^{\prime}}}{2r^{2}f(r)}\bigg{[}{}_{e}\Phi_{lm}\frac{\partial{}_{o}\Phi_{l^{\prime}m^{\prime}}}{\partial r}-\frac{\partial{}_{e}\Phi_{lm}}{\partial r}{}_{o}\Phi_{l^{\prime}m^{\prime}}$	(11)
	$\displaystyle\quad+{}_{o}\Phi_{lm}\frac{\partial{}_{e}\Phi_{l^{\prime}m^{\prime}}}{\partial r}+\frac{\partial{}_{o}\Phi_{lm}}{\partial r}{}_{e}\Phi_{l^{\prime}m^{\prime}}\bigg{]}$
	$\displaystyle+\frac{B_{lml^{\prime}m^{\prime}}}{2r^{2}}\bigg{[}-\frac{\partial{}_{o}\Phi_{lm}}{\partial r}\frac{\partial{}_{o}\Phi_{l^{\prime}m^{\prime}}}{\partial r}-\frac{\partial{}_{e}\Phi_{lm}}{\partial r}\frac{\partial{}_{e}\Phi_{l^{\prime}m^{\prime}}}{\partial r}$
	$\displaystyle\quad-\frac{1}{f(r)^{2}}\frac{\partial{}_{o}\Phi_{lm}}{\partial t}\frac{\partial{}_{o}\Phi_{l^{\prime}m^{\prime}}}{\partial t}-\frac{1}{f(r)^{2}}\frac{\partial{}_{e}\Phi_{lm}}{\partial t}\frac{\partial{}_{e}\Phi_{l^{\prime}m^{\prime}}}{\partial t}\bigg{]}$
	$\displaystyle-\frac{\lambda^{2}Y_{lm}Y_{l^{\prime}m^{\prime}}\left({}_{e}\Phi_{lm}{}_{e}\Phi_{l^{\prime}m^{\prime}}+{}_{o}\Phi_{lm}{}_{o}\Phi_{l^{\prime}m^{\prime}}\right)}{2r^{4}f(r)}$
	$\displaystyle-\frac{Q_{tot}}{{r^{4}f(r)}}\sum_{l,m}\lambda Y_{lm}{}_{e}\Phi_{lm}-\frac{Q_{tot}^{2}}{2r^{4}f(r)}$

and

$\displaystyle S_{r}$	$\displaystyle=\frac{1}{r^{2}}\sum_{l,m}\sum_{l^{\prime},m^{\prime}}\bigg{\{}A_{lml^{\prime}m^{\prime}}\left[f(r)\frac{\partial{}_{e}\Phi_{lm}}{\partial r}\frac{\partial{}_{o}\Phi_{l^{\prime}m^{\prime}}}{\partial r}\right.$	(12)
	$\displaystyle\left.-\frac{1}{f(r)}\frac{\partial{}_{o}\Phi_{lm}}{\partial t}\frac{\partial{}_{e}\Phi_{l^{\prime}m^{\prime}}}{\partial t}\right]+B_{lml^{\prime}m^{\prime}}\left[\frac{\partial{}_{o}\Phi_{l^{\prime}m^{\prime}}}{\partial r}\frac{\partial{}_{o}\Phi_{lm}}{\partial t}\right.$
	$\displaystyle\left.+\frac{\partial{}_{e}\Phi_{l^{\prime}m^{\prime}}}{\partial t}\frac{\partial{}_{e}\Phi_{lm}}{\partial r}\right]\bigg{\}}$

where

\begin{gathered}A_{lml^{\prime}m^{\prime}}=\csc\theta\left[\dot{Y}_{lm}Y^{\prime}_{l^{\prime}m^{\prime}}-Y^{\prime}_{lm}\dot{Y}_{l^{\prime}m^{\prime}}\right],\\ B_{lml^{\prime}m^{\prime}}=\csc^{2}\theta\left[Y^{\prime}_{lm}Y^{\prime}_{l^{\prime}m^{\prime}}+\dot{Y}_{lm}\dot{Y}_{l^{\prime}m^{\prime}}\right]\end{gathered}

(13)

The main takeaways from this brief summary so far are as follows:

•

The stress-energy tensor $T_{EM}^{\mu\nu}$ involves both linear terms in ${}_{e}\Phi_{lm}$ and ${}_{o}\Phi_{lm}$ , as well as quadratic terms involving products of the form ${}_{e}\Phi_{lm}\times{}_{e}\Phi_{l^{\prime}m^{\prime}}$ , ${}_{o}\Phi_{lm}\times{}_{o}\Phi_{l^{\prime}m^{\prime}}$ , and ${}_{e}\Phi_{lm}\times{}_{o}\Phi_{l^{\prime}m^{\prime}}$ and their derivatives, summed over the indices $l$ , $l^{\prime}$ , $m$ , and $m^{\prime}$ . As a result, we can anticipate that the frequencies associated with the mixing modes will be either linear or quadratic as argued for in [1].
•

We can decompose equation (10) using, for example, the 4D tensor harmonics defined as in [1]. In the Laplace spectral space (or less formally in the Fourier space), this will allow us to study the spectral decomposition of the energy flux, shear, and density as functions of $r$ for each multipole $lm$ and frequency $\Omega_{n}$ similar to [31, 32]. Furthermore, we can examine how they couple to other radiative multipoles $l^{\prime}m^{\prime}$ of $\Phi_{l^{\prime}m^{\prime}}$ .
•

The spherical decomposition of $T_{EM}^{\mu\nu}$ in expression (9) is not fully compatible with gravitational perturbations. To resolve this, we propose projecting it onto the tensorial spherical harmonics provided in Appendix A of [1] or following ref [23] as demonstrated in Appendix C,D,E, and F in [1], ensuring compatibility with the Regge-Wheeler gravitational analysis[1].
•

As $t\rightarrow\infty$ , both $\Phi$ and its derivatives will vanish, and the total charge $Q_{\text{tot}}$ will appear frozen near the horizon from the perspective of an asymptotic observer. Consequently, only the monopole term, i.e., the Reissner–Nordström charge, will be measurable at “late times”. However, as different radiative modes (i.e. QNMs and tails) have different decaying rates, that suggest the we should be careful what time-scale we are interested in studying.

Refer to caption — Figure 1: Comparison of gravitational and electromagnetic QNMs in Schwarzschild spacetime. Gravitational frequencies $\omega^{(1)}_{ln}$ for $l=2,3$ and $n=0,1$ are compared with electromagnetic frequencies $\Omega^{(1)}_{ln}$ for $l=1,2$ and $n=0,1$ . For the quadratic QNMs, the real part of the frequencies is given by the sum $\Omega^{(2)}_{(i\times j)R}=\Omega^{(1)}_{(i)R}\pm\omega^{(1)}_{(j)R}$ of the linear modes, while the imaginary part, representing the reciprocal of the decay time, is given by $\Omega^{(2)}_{(i\times j)I}=\Omega^{(1)}_{(i)I}+\omega^{(1)}_{(j)I}$ as a sum of the corresponding linear terms. The blue squares and red dots represent the electromagnetic and gravitational QNMs, respectively, while cyan rhombuses indicate the quadratic QNMs.

II Second order electromagnetic perturbation

In [1], the coupling between electromagnetic and gravitational modes was investigated perturbatively in the framework of minimally coupling Einstein-Maxwell system. It was reported that, in the simplest case of mixing between the two fields, when the electromagnetic field sources gravitational perturbations, it excites gravitational modes that are linear and quadratic in the electromagnetic QNMs. Also amplitude of the linear gravitational QNMs get correction from electromagnetic field. As we go higher in mixing, more coupling between the modes is expected to take place, such that modes quadratic in both fields arise.

On the other hand, at first and second order in $\alpha$ , the Maxwell equation for the perturbed electromagnetic field is expressed as:

\begin{gathered}\alpha\nabla_{\mu}^{(0)}F^{(1)\mu\nu}=\alpha J^{(1)\nu}\\ \alpha^{2}\nabla_{\mu}^{(0)}F^{(2)\mu\nu}=\alpha^{2}J_{eff}^{(2)\,\nu},\end{gathered}

(14)

where $F^{(1,2)}_{\mu\nu}$ and $J^{(1,2)\nu}$ represent the first and second -order perturbation of the electromagnetic field and second-order current source term respectively. The covariant derivative $\nabla_{\mu}^{(0)}$ is taken with respect to the unperturbed Schwarzschild metric $g_{\mu\nu}^{(0)}$ , hence it’s straightforward to reduce the second order field perturbation to their corresponding Regge-Wheeler second odd ${}_{o}\Phi^{(2)}_{lm}$ and even ${}_{e}\Phi^{(2)}_{lm}$ perturbation:

\mathcal{L}_{EM}\,{}_{o,e}\Phi^{(2)}_{lm}(t,x)={}_{\text{eff}}\mathcal{S}^{(2)}_{o,e\,lm}(t,x).

At the second order, there is additional effective nonlinear contribution which $J_{\text{eff}}^{(2)\nu}$ represents coupling between the first-order electromagnetic field $F^{(1)}_{\mu\nu}$ and the first-order metric perturbation $g^{(1)}_{\mu\nu}$ . The effective second-order current, $J_{\text{eff}}^{(2)\nu}$ , is given by

$\displaystyle J_{\text{eff}}^{(2)\nu}$	$\displaystyle=-\frac{1}{2}\left(g^{(0)\gamma\mu}F^{(1)\alpha\nu}+g^{(0)\nu\gamma}F^{(1)\mu\alpha}\right)$	(15)
	$\displaystyle\times\left(\nabla_{\mu}^{(0)}g_{\alpha\gamma}^{(1)}+\nabla_{\alpha}^{(0)}g_{\mu\gamma}^{(1)}-\nabla_{\gamma}^{(0)}g_{\alpha\mu}^{(1)}\right),$
	$\displaystyle=\frac{1}{2}F^{(1)\nu\alpha}\nabla_{\alpha}^{(0)}g_{\mu}^{\mu\,(1)}$

It is straightforward to find the corresponding Regge-Wheeler sources corresponding to $J_{\text{eff}}^{(2)\nu}$ , through using vector spherical harmonics in Appendix D. We report the final result for even and odd sources respectively here

\begin{gathered}{}_{\text{eff}}S_{e\,lm}=\frac{1}{2\lambda}\left(\frac{\partial^{2}}{\partial t^{2}}-\frac{\partial^{2}}{\partial x^{2}}\right)\tilde{j}_{lm}\\ \tilde{j}_{lm}=\left[-2\left(\frac{\partial}{\partial x}+\bar{A}\right)\psi_{e\,l^{\prime}m^{\prime}}\right]\left[Q_{\text{tot}}\,\delta_{ll^{\prime}}\delta_{mm^{\prime}}+(-1)^{m}\,\lambda^{\prime\prime}\,\Phi_{e\,l^{\prime\prime}m^{\prime\prime}}\,B_{l,-m;l^{\prime},m^{\prime},l^{\prime\prime},m^{\prime\prime}}\right],\\ \bar{A}=\frac{2}{r(\lambda-2)+6M}\left[\left(\frac{1}{2}\lambda-1\right)\frac{1}{2}\lambda+\frac{3M}{r}\left(\frac{1}{2}\lambda-f\right)\right].\end{gathered}

(16)

\begin{gathered}{}_{\text{eff}}S_{o\,lm}=(-1)^{m}\frac{f}{2}\Bigg{\{}-C_{l,-m;l^{\prime},m^{\prime},l^{\prime\prime},m^{\prime\prime}}\frac{\lambda^{\prime\prime}\Phi_{o\,l^{\prime\prime}m^{\prime\prime}}}{r^{2}}\left[-2\left(\frac{\partial}{\partial x}+\bar{A}\right)\psi_{e\,l^{\prime}m^{\prime}}\right],\\ +g^{aa}\frac{\partial}{\partial x_{a}}\left[-2\left(\frac{\partial}{\partial x}+\bar{A}\right)\psi_{e\,l^{\prime}m^{\prime}}\right]\Bigg{[}C_{l,-m;l^{\prime\prime},m^{\prime\prime},l^{\prime},m^{\prime}}\frac{\partial\Phi_{o\,l^{\prime\prime}m^{\prime\prime}}}{\partial x_{a}}+D_{l,-m;l^{\prime\prime},m^{\prime\prime},l^{\prime},m^{\prime}},\left(\epsilon_{a}^{b}\frac{\partial\Phi_{e\,l^{\prime\prime}m^{\prime\prime}}}{\partial x_{b}}+\frac{r^{2}j_{(a)\,l^{\prime\prime}m^{\prime\prime}}}{\lambda^{\prime\prime}}\right)\Bigg{]}\Bigg{\}}\\ \epsilon_{a}^{b}=\begin{pmatrix}0&f\\ \frac{1}{f}&0\end{pmatrix}.\end{gathered}

(17)

Here, $B_{l,m,l^{\prime},m^{\prime},l^{\prime\prime},m^{\prime\prime}}$ , $C_{l,m,l^{\prime},m^{\prime},l^{\prime\prime},m^{\prime\prime}}$ , and $D_{l,m,l^{\prime},m^{\prime},l^{\prime\prime},m^{\prime\prime}}$ are the multipolar mixing coefficients provided in Appendix D of [1]. In this context, we have ignored the source term for the linear gravitational perturbation to simplify the expression. Including these terms would contribute a source component independent of $\Phi$ or $\psi$ , as well as the linear electromagnetic modes $\Phi^{(1)}$ , which depend on the specifics of the electromagnetic potential.

Generically, beyond linear electromagnetic perturbations, at higher orders $\Phi^{(j)},j\geq 2$ , gravitational perturbations will similarly source higher-order electromagnetic modes. A similar analysis, as outlined in Section II of [1], can be applied here. For instance, at second order in electromagnetic QNMs, we expect to observe phenomena that only arise with higher-level mixing between the two fields, as seen in the gravitational case. For example at second order in $\alpha$ , if a gravitational linear QNM has a frequency $\omega^{(1)}_{i}$ and a linear electromagnetic QNM has a frequency $\Omega^{(1)}_{j}$ , part of the quadratic electromagnetic spectrum will include frequencies quadratic in both gravitational-electromagnetic modes $\Omega^{(2)}_{i}=\Omega^{(1)}_{j}\pm\omega^{(1)}_{k}$ as shown in Fig. 1, however quadratic modes in electromagnetic will be absent.

In most astrophysical scenarios, such as gravitational collapse or mergers, gravitational effects are expected to dominate. In this context, searching for the mixing between electromagnetic and gravitational fields suggests that electromagnetic QNMs with gravitational imprints may be more interesting than the reverse, but the richness of typical BNS and other magnetized phenomena will make it extremely hard to resolve the aforementioned effect. Nevertheless, both approaches could serve as tests for minimal coupling within the framework of General Relativity or for more complex couplings in modified gravity theories. For example, if we consider a different coupling such as the Einstein-Maxwell-Dilaton theory [33, 34], the induced QNMs would exhibit different frequencies. Moreover, if we can imagine an astrophysical phenomenon where both effects are pronounced, it could present a valuable opportunity for multi-messenger astronomy[11].

Following [1], we briefly consider a sample source term within the Dirac delta potential approximation, detailed in section IV and section III of [1], and solve for the resulting second-order electromagnetic perturbation. If we assume a quadratic QNM source term, as in [35], we have

S=\frac{\Phi^{(1)}\Psi^{(1)}}{(1+\zeta|x|)^{2}}\propto\frac{\,e^{-\left(\frac{V_{\text{EM}}}{2}+\frac{V_{G}}{2}\right)(t-|x|)}}{(1+\zeta|x|)^{2}}.

(18)

Here we only focus on the QNM part. Typically, we should expect a correction to the amplitude of the electromagnetic linear modes and the excitation of new frequencies given by

\Omega^{(2)}=-i\frac{V_{G}+V_{EM}}{2}.

For $x>0$ and $x<0$ , respectively, we have the second-order perturbation $\Phi^{(2)}_{QNM}$ :

		$\displaystyle\Phi^{(2)}_{QNM}(t,x>0)=\frac{8}{\zeta^{2}V_{G}}e^{-\frac{uV_{\text{EM}}}{2}}\left\{V_{\text{EM}}e^{-\frac{V_{\text{EM}}}{\zeta}}\left[\operatorname{Ei}\left(\frac{uV_{\text{EM}}}{2}+\frac{V_{\text{EM}}}{\zeta}\right)-\operatorname{Ei}\left(\frac{V_{\text{EM}}}{\zeta}\right)\right]+\zeta\right\}$		(19)
		$\displaystyle\quad\quad\quad\quad+\frac{8}{\zeta^{2}V_{G}}e^{-u\frac{V_{\text{EM}}+V_{G}}{2}}\left\{(V_{\text{EM}}+V_{G})e^{-\frac{V_{\text{EM}}+V_{G}}{\zeta}}\left[\operatorname{Ei}\left(\frac{V_{\text{EM}}+V_{G}}{\zeta}\right)-\operatorname{Ei}\left((V_{\text{EM}}+V_{G})\left(\frac{u}{2}+\frac{1}{\zeta}\right)\right)\right]-\zeta\right\}$		(19)

		$\displaystyle\Phi^{(2)}_{QNM}(t,x<0)=\frac{8}{\zeta^{2}V_{G}}e^{-\frac{vV_{\text{EM}}}{2}}\left\{e^{-\frac{V_{\text{EM}}}{\zeta}}V_{\text{EM}}\left[\operatorname{Ei}\left(\frac{uV_{\text{EM}}}{2}+\frac{V_{\text{EM}}}{\zeta}\right)-\operatorname{Ei}\left(\frac{V_{\text{EM}}}{\zeta}\right)\right]+\zeta\right\}$		(20)
		$\displaystyle\quad\quad\quad\quad+\frac{8}{\zeta^{2}V_{G}}e^{-\frac{uV_{G}}{2}-\frac{vV_{\text{EM}}}{2}}\left\{(V_{\text{EM}}+V_{G})e^{-\frac{V_{\text{EM}}+V_{G}}{\zeta}}\left[\operatorname{Ei}\left(\frac{V_{\text{EM}}+V_{G}}{\zeta}\right)-\operatorname{Ei}\left((V_{\text{EM}}+V_{G})\left(\frac{u}{2}+\frac{1}{\zeta}\right)\right)\right]-\zeta\right\}$		(20)

III Perturbations due to point charges

In this section, we analyze the electromagnetic perturbation due to one or more point charges in Schwarzschild coordinates. Consider a particle moving along a worldline parameterized by an affine parameter $\tau$ . The four-velocity vector $U^{\mu}$ of the particle is given by

U^{\mu}=\left(\frac{dt_{p}(\tau)}{d\tau},\frac{d\bar{r}_{p}(\tau)}{d\tau},\frac{d\bar{\theta}_{p}(\tau)}{d\tau},\frac{d\bar{\phi}_{p}(\tau)}{d\tau}\right),

(21)

which can also be reparameterized in terms of the time coordinate $t_{p}$ as:

U^{\mu}=\frac{1}{\frac{dt_{p}}{d\tau}}\left(1,\frac{dr_{p}(t_{p})}{dt_{p}},\frac{d\theta_{p}(t_{p})}{dt_{p}},\frac{d\phi_{p}(t_{p})}{dt_{p}}\right).

The four-current vector $J_{\mu}$ associated with the charge $q$ can then be written as:

	$\displaystyle J_{\mu}$	$\displaystyle=\left\{f(r),-\frac{1}{f(r)}\frac{dr_{p}(t)}{dt},-r^{2}\frac{d\theta_{p}(t)}{dt},-r^{2}\sin^{2}\theta\frac{d\phi_{p}(t)}{dt}\right\}$		(22)
		$\displaystyle\quad\times q\,\frac{\delta(r-r_{p}(t))}{r^{2}}\sum_{lm}Y_{lm}(\theta,\phi)Y^{*}_{lm}(\theta_{p}(t),\phi_{p}(t)).$		(22)

This four-current vector can be further decomposed using vector spherical harmonics. The components of the current can be expressed as

	$\displaystyle J_{\mu}$	$\displaystyle=\left\{j_{(1)lm}Y_{lm},j_{(2)lm}Y_{lm},j_{(3)lm}\dot{Y}_{lm}+j_{(4)lm}\csc\theta\,Y^{\prime}_{lm},\right.$		(23)
		$\displaystyle\quad\left.j_{(3)lm}Y^{\prime}_{lm}-j_{(4)lm}\sin\theta\dot{Y}_{lm}\right\},$		(23)

where the coefficients $j_{(i)lm}$ are extracted by projecting over the sphere as outlined in D:

j_{(i)lm}=\int d\Omega\mathbf{Y}^{\mu*}_{(i)}J_{\mu}.

(24)

Consequently,

$\displaystyle j_{(1)lm}$	$\displaystyle=qf(r)\frac{\delta(r-r_{p})}{r^{2}}\,Y^{*}_{lm}(\theta_{p},\phi_{p}),$	(25)
$\displaystyle j_{(2)lm}$	$\displaystyle=-\frac{q}{f(r)}\frac{dr_{p}}{dt}\frac{\delta(r-r_{p})}{r^{2}}Y^{*}_{lm}(\theta_{p},\phi_{p}),$	(26)
$\displaystyle j_{(3)lm}$	$\displaystyle=\frac{q}{\lambda}\delta\left(r-r_{p}\right)\left\{im\frac{d\phi_{p}}{dt}-\frac{d\theta_{p}}{dt}\mathcal{A}_{lm}[Y_{lm}^{*}]\right\},$	(27)
$\displaystyle j_{(4)lm}$	$\displaystyle=\frac{q}{\lambda}\delta\left(r-r_{p}\right)\left\{im\frac{d\theta_{p}}{dt}\mathcal{C}_{lm}[Y_{lm}^{}]+\frac{d\phi_{p}}{dt}\mathcal{D}_{lm}[Y_{lm}^{}]\right\}.$	(28)

Here, $\mathcal{A}_{lm}$ , $\mathcal{C}_{lm}$ , and $\mathcal{D}_{lm}$ are operators that effectively act on the multipolar numbers of the spherical harmonics, inducing mulipolar mode-mixing defined as

$\displaystyle\mathcal{A}_{lm}[k_{lm}]$	$\displaystyle=\sum_{lm}\,k_{lm}\int d\Omega Y_{lm}\frac{\partial Y_{l^{\prime}m^{\prime}}^{*}}{\partial\theta}\,,$	(29)
$\displaystyle\mathcal{C}_{lm}[k_{lm}]$	$\displaystyle=\sum_{lm}\,k_{lm}\int d\Omega\csc\theta\,Y_{lm}Y_{l^{\prime}m^{\prime}}^{*}\,,$
$\displaystyle\mathcal{D}_{lm}[k_{lm}]$	$\displaystyle=\sum_{lm}\,k_{lm}\int d\Omega\sin\theta\,Y_{lm}\frac{\partial Y_{l^{\prime}m^{\prime}}}{\partial\theta}\,.$

The effective source term for the Regge-Wheeler equation for multipoles $l>0$ (even perturbations) ${}_{e}\Phi$ is given by:

{}_{J}\mathcal{S}_{\text{e}\,lm}=\frac{f(r)}{\lambda}\left(\frac{\partial}{\partial r}\left[r^{2}j_{(1)lm}\right]-r^{2}\frac{\partial j_{(2)lm}}{\partial t}\right).

(30)

For the point charge case, we obtain:

$\displaystyle{}_{J}\mathcal{S}_{\text{e}\,lm}=$	$\displaystyle\frac{q}{\lambda}\Bigg{\{}\Bigg{[}\,\dot{r}_{p}(t)\left(\dot{\theta}_{p}(t)\frac{\partial Y_{lm}^{}}{\partial\theta}-im\,\dot{\phi}_{p}(t)Y_{lm}^{}\right)$	(31)
	$\displaystyle+\left(f(r)\frac{f(r)}{dr}+\ddot{r}_{p}(t)\right)Y_{lm}^{*}\Bigg{]}\delta\left(r-r_{p}(t)\right)$
	$\displaystyle+\left[f(r)^{2}-\dot{r}_{p}^{2}(t)\right]\delta^{\prime}\left(r-r_{p}(t)\right)Y_{lm}^{*}\Bigg{\}}.$

The odd perturbation ${}_{o}\Phi$ is described by

{}_{J}\mathcal{S}_{\text{o}\,lm}=-f(r)j_{(4)lm}.

(32)

Subsequently, for the point particle case we have

	$\displaystyle{}_{J}\mathcal{S}_{\text{o}\,lm}=$	$\displaystyle-\frac{q}{\lambda}f(r)\delta\left(r-r_{p}(t)\right)\left\{im\frac{d\theta_{p}}{dt}\mathcal{C}_{lm}[Y_{lm}^{*}]\right.$		(33)
		$\displaystyle\left.+\frac{d\phi_{p}}{dt}\mathcal{D}_{lm}[Y_{lm}^{*}]\right\}.$		(33)

Finally, the general solution for ${}_{e,o}\Phi_{lm}$ is given by

{}_{e,o}\Phi_{lm}\left(r^{\prime},t^{\prime}\right)=\int_{-\infty}^{\infty}dt\int_{2M}^{\infty}dr\,G_{lm}\left(t,t^{\prime},r,r^{\prime}\right){}_{J}\mathcal{S}_{e,o\,lm},

(34)

where $G_{lm}(t,t^{\prime},r,r^{\prime})$ represents the Green’s function for the Regge-Wheeler master equation with spin $s=1$ , see Section II in [1] for further details. For completeness, we emphasize that for any signal observed at a specific event $(t,x)$ , the Green’s functions define a limited domain of support in the $(t^{\prime},x^{\prime})$ plane [36, 37]. There are 2D sub-regions in the $(t^{\prime},x^{\prime})$ space that contribute to the event at $(t,x)$ through convolution of the Green’s function with the source term $\mathcal{S}_{e,o\,lm}$ . In the case of a point particle, this contribution is restricted to a 1D curve in $(t^{\prime},x^{\prime})$ . Consequently, the problem reduces to a path integral along the worldline of the particle, which can be reparameterized over a segment of the path with respect to any variable that is one-to-one with the affine parameter on that segment.

For example, if we stick to $t_{p}$ as affine parameter, then for odd perturbation ${}_{o}\Phi_{lm}$ the solution becomes

	$\displaystyle{}_{o}\Phi_{lm}(r^{\prime},t^{\prime})=$	$\displaystyle-\frac{q}{\lambda}\int_{-\infty}^{\infty}dt_{p}\,f\left(r_{p}(t_{p})\right)\left\{im\frac{d\theta_{p}}{dt_{p}}\mathcal{C}_{lm}[Y_{lm}^{*}]\right.$		(35)
		$\displaystyle\left.+\frac{d\phi_{p}}{dt_{p}}\mathcal{D}_{lm}[Y_{lm}^{*}]\right\}\,G_{lm}\left(r^{\prime},t^{\prime},r_{p}(t_{p}),t_{p}\right).$		(35)

For the even perturbation ${}_{e}\Phi_{lm}$

_e	$\displaystyle\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=\frac{q}{\lambda}\int_{-\infty}^{\infty}dt_{p}\Bigg{\{}\left[\ddot{r}_{p}(t_{p})-f\left(r_{p}(t_{p})\right)\frac{df\left(r_{p}\right)}{dr_{p}}\right.$	(36)
	$\displaystyle\left.-im\,\frac{dr_{p}}{dt_{p}}\frac{d\phi_{p}}{dt_{p}}+m\cot\theta_{p}\,\frac{dr_{p}}{dt_{p}}\frac{d\theta_{p}}{dt_{p}}\right]Y_{lm}^{*}$
	$\displaystyle+\left[\sqrt{(l-m)(l+m+1)}e^{i\hat{\phi}}\frac{dr_{p}}{dt_{p}}\frac{d\theta_{p}}{dt_{p}}\right]Y_{l,m+1}^{*}\Bigg{\}}$
	$\displaystyle\times G_{lm}\left(r^{\prime},t^{\prime},r_{p}(t_{p}),t_{p}\right)+\frac{q}{\lambda}\int^{\infty}_{-\infty}dt_{p}\left.\frac{G_{lm}}{dr}\right\|_{r=r_{p}(t_{p})}$
	$\displaystyle\times\left[\dot{r}_{p}(t_{p})^{2}-f^{2}\left(r_{p}(t_{p})\right)\right]Y_{lm}^{*}.$

However, as we will see in the upcoming subsections, a well-chosen parameter will depend on the worldline being examined, and possibly on the $(t,x)$ coordinates in which we are aiming to examine the perturbations. A thoughtful choice of this parameter can significantly simplify the complexity of the computation.

III.1 Ideal dipole

The linearity of Maxwell’s equations extends naturally to the electromagnetic Regge-Wheeler equation. If the total current $J_{\mu}$ is the sum of the currents from two particles, i.e.,

J_{\mu}=J_{\mu}^{\text{\scriptsize particle(1)}}+J_{\mu}^{\text{\scriptsize particle(2)}},

(37)

then the resulting scalar field $\Phi$ is the superposition of the fields from each particle:

\Phi=\Phi^{\text{\scriptsize particle(1)}}+\Phi^{\text{\scriptsize particle(2)}}.

(38)

Here, $\Phi^{\text{\scriptsize particle(1)}}$ and $\Phi^{\text{\scriptsize particle(2)}}$ are the solutions to the Regge-Wheeler equation for each current. For a simple dipole structure, consider a particle with charge $q$ on worldline $\gamma_{1}$ parameterized by $r_{p1}$ , and another with charge $-q$ on worldline $\gamma_{2}$ parameterized by $r_{p2}$ , where $r_{p2}-r_{p1}\approx\eta+\mathcal{O}(\eta^{2})$ ²²2In case of radial free fall of a radial dipole $\mathcal{O}(\eta^{2})$ =0 and the expansion truncates at first order correction in $\eta$ .. Thus, the total field for the dipole could simply be written as a variation of a the perturbation due to a positive charge as

{}_{e,o}\Phi^{d}=\delta\,{}_{e,o}\Phi^{+}.

(39)

The complexity of solving this problem can be divided into three parts:

1.

Solving for the worldines of charges $q$ and $-q$ considering the electromagnetic force between the particles and self-interaction.
2.

Solving for the Green’s function $G_{lm}$ .
3.

Evaluating the one-dimensional integral to obtain ${}_{e,o}\Phi$ .

We will now consider some simple cases.

III.2 Examples

III.2.1 Point charge in radial free fall

Consider a charge $q$ freely falling radially from rest with $E=1$ and with angular coordinate $(\theta_{0},\phi_{0})$ , while $\frac{dr_{p}}{dt_{p}}=-\beta(r_{p})f(r_{p})$ and $j_{(3)lm}=j_{(4)lm}=0$ vanishes. We have

\displaystyle{}_{o}\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=0.

(40)

After re-paramterization of the integral in terms of $r_{p}$ , as outlined in Appendix C, the even part is

	$\displaystyle{}_{e}\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=\frac{-qY_{lm}^{*}}{\lambda}\Bigg{[}$	$\displaystyle\int dr_{p}\Bigg{\{}\left(-\frac{3f(r_{p})}{2\beta(r_{p})}\frac{df(r_{p})}{dr_{p}}\right)G_{lm}$		(41)
		$\displaystyle+\left.\frac{dG_{lm}}{dr}\right\|_{r=r_{p}}\left(-\frac{f(r_{p})^{2}}{\beta(r_{p})}\right)\Bigg{\}}\Bigg{]}.$		(41)

III.2.2 Point charge in a circular orbit

Consider a charge $q$ in a circular orbit at a constant radius $r_{p}$ , with angular coordinates $(\theta_{0},\phi_{p}(t_{p}))$ . The angular position $\phi_{p}(t_{p})$ is related to time $t_{p}$ by $t_{p}=\frac{\phi_{p}}{\bar{\Omega}}$ , where $\bar{\Omega}$ is the constant angular velocity of the particle, and $j_{(2)lm}=0$ . The odd-parity perturbation is given by

	$\displaystyle{}_{o}\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=\frac{q}{\lambda}imf(r_{0})$	$\displaystyle\int_{\phi_{pi}}^{\phi_{pf}}d\phi_{p}\,\mathcal{D}_{lm}\left[Y_{lm}^{*}\right]_{\theta=\theta_{0},\phi=\phi_{p}}$		(42)
		$\displaystyle\times G_{lm}\left(r^{\prime},t^{\prime},r_{0},\frac{\phi_{p}}{\Omega}\right).$		(42)

The even-parity perturbation is given by

		$\displaystyle{}_{e}\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=\frac{-q}{\lambda\Omega}\int_{\phi_{pi}}^{\phi_{pf}}d\phi_{p}\,Y^{*}_{lm}(\theta_{0},\phi_{p})$		(43)
		$\displaystyle\quad\times{\left[f\left.\frac{\partial f(r)}{\partial r}\right\|_{r=r_{0}}G_{lm}(t^{\prime},r^{\prime},\frac{\phi_{p}}{\Omega},r_{0})+f^{2}\left.\frac{\partial G_{lm}}{\partial r}\right\|_{r=r_{0}}\right]}$		(43)

III.2.3 Quasi-circular orbit

We can imagine a worldline where the particle’s motion is confined to a plane at $\theta=\frac{\pi}{2}$ while undergoing angular motion in the $\hat{\phi}$ direction, similar to that of a circular orbit. Simultaneously, the particle experiences radial motion with velocity $v$ , where $v\ll r\hat{\Omega}$ , starting from $r_{\text{inner}}$ and moving to $r_{\text{outer}}$ . We can utilize integrals (36,35) with $\frac{dr_{p}}{dt_{p}}$ and $\frac{d\phi_{p}}{dt_{p}}$ as constants, while $\frac{d\theta_{p}}{dt_{p}}=0$ . This approach approximates the quasi-circular worldline of the point particle.

IV Dirac delta function potential

Away from the light ring, the potential diminishes, and solutions become asymptotically outgoing (or ingoing) near the flat region (or horizon), which is approximately true for the Regge-Wheeler potential at large distances. This holds analytically for superlocalized potentials like the Dirac delta function $V_{EM}\delta(x)$ [35]. The Green’s function³³3See [38, 35, 3] for a detailed discussion of the flat Green’s function component $G_{F}$ , the QNM component $G_{QNM}$ , and the branch cut component $G_{B}$ . is expressed as:

\begin{gathered}G(t-t^{\prime},r,r^{\prime})=G_{F}(t-t^{\prime},r,r^{\prime})+G_{Q}(t-t^{\prime},r,r^{\prime}),\end{gathered}

(44)

The Green’s function is further written as

G(T,r,r^{\prime})=-\frac{1}{2}\left[\Theta(A)-\Theta(B)\right]-\frac{1}{2}e^{-\frac{V_{EM}}{2}B}\Theta(B),

(45)

where $\Theta$ is the Heaviside step function. The derivative of the Green’s function with respect to $r^{\prime}$ is:

	$\displaystyle f(r^{\prime})$	$\displaystyle\frac{\partial\,G(t-t^{\prime},r,r^{\prime\prime})}{\partial r^{\prime\prime}}\bigg{\|}_{r^{\prime\prime}=r^{\prime}}=\frac{1}{2}\left[s\delta(A)-\tilde{s}\delta(B)\right]$		(46)
		$\displaystyle\quad-\frac{\tilde{s}}{2}e^{-\frac{V_{EM}}{2}(B)}\left[\frac{V_{EM}}{2}\Theta(B)-\delta(B)\right],$		(46)

⁴⁴4Taking derivative of with Dirac delta function distribution should be done under the integration sign, that is why we are not canceling the two contributions of

\delta(B)

, in case higher derivatives are needed.

Here, $s(x-x^{\prime})$ and $\tilde{s}(x-x_{EM})$ are sign functions, taking $+$ for positive arguments and $-$ for negative arguments. The terms $A$ and $B$ encode the causal structure, defined as:

	$\displaystyle A$	$\displaystyle=t-t^{\prime}-\|x-x^{\prime}\|,$		(47)
	$\displaystyle B$	$\displaystyle=t-t^{\prime}-\|x\|-\|x^{\prime}\|,$		(47)

IV.1 Support regions

The causal structure becomes more evident and sharper with simpler potentials, such as the Dirac delta potential. In this case, the QNMs component of the Green’s function is supported by the past light cone of the point $(0,t-x)$ , which is defined by the region $B$ in the $(t^{\prime},x^{\prime})$ plane.

Focusing on point sources, a particle moving along a trajectory $\gamma$ contributes to the QNM signal observed at $(x,t)$ primarily due to the particle’s interaction with the potential along its worldline. Almost all of this contribution arises from the particle’s history in proximity of the potential. Specifically, if the particle starts from an asymptotically far region, $r\rightarrow\infty$ , its contribution begins at $t^{\prime}\rightarrow-\infty$ and ends when its trajectory exits the past light cone of the point $(0,t-x)$ in the $(t^{\prime},x^{\prime})$ plane.

In many cases, as seen in the literature [39, 35], the source’s contribution is considered to start at a finite time, particularly when gravitational QNMs source the quadratic QNMs, as in binary black hole merger scenarios. While starting the particle’s contribution at $t^{\prime}\rightarrow-\infty$ helps get rid of contribution from electromagnetic field I.C., it can be helpful to consider the particle (the source) starting its contribution at a finite time at some large radius $r_{0}$ , to simplify the analysis. By setting the beginning of time to $t=0$ , the source’s contribution is defined by its intersection with the region $B$ , bounded below by $t^{\prime}=0$ in the $(t^{\prime},x^{\prime})$ space. For further clarification, visual references are useful; see Figures 2 (a–f).

In the Figures 2, we imagine a point charge $q$ described by the solid green curve contributing to the signal observed at $x=3M$ . Initially, in (a), there is no QNM contribution at $t=10M$ , but as time progresses as shown in (b) and (c), the support region $B$ expands, capturing more of the particle’s worldline as it enters the region. In contrast, when the observation point $x$ gets further away from the potential (as shown in (f)) the support region shrinks down. In that regime, in the limit of $x\rightarrow\infty$ , the contribution only near the potential will excite QNM signal. If we assume that the particle started from $r\rightarrow\infty$ while $t\rightarrow-\infty$ , the worldline would always remain inside the region $B$ of any $(t,x)$ , and these past light cones $B$ would extend to infinity.

For $G_{Q}$ , the redness of the contours, increases on a light-like mesh, with inner past light cones appearing whiter than the outer ones. The outermost part represents the edge of the signal, while the innermost part corresponds to the signal’s end.

The integration limits are $(t,x)$ -dependent. As we have seen before, if we choose to parameterize the path integral, defined by the overlap of the green curve with the $B$ region (colored in red), using the parameter $r_{p}$ for certain geodesics⁵⁵5i.e. radial free fall the integration limits are $r_{pBmax}(t,x)$ and $r_{pBmin}(t,x)$ . For a particle starting at infinity, we generally have $r_{pBmin}(t,x)\rightarrow\infty$ . However, solving for $r_{pBmax}(t,x)$ requires solving for $B=0$ in terms of $r_{p}$ , given by

t-|x|=t^{\prime}(r_{p})+|x^{\prime}(r_{p})|,

ensuring the integration is properly constrained and accounts for the causal structure encoded in $B$ .

To study the flat contribution of the Green’s function, the situation is analogous except that support region is $A-B$ , where $A$ represents another past light cone centered at $(x,t)$ . The reader can observe that the flat part’s density plot, as in figures 2, is constant and colored in brown for $x>0$ (or cyan for $x<0$ ), as the value of $G_{F}$ inside $A-B$ remains a constant $=-1/2$ .

Another contribution we should account for is what, mathematically speaking, can be regarded as the evaluation of the derivative of a Dirac delta distribution under the integration sign. Physically, we can interpret this derivative as an additional contribution arising from the edges of the QNM region $B$ and the flat region $A-B$ , where the particle is about to penetrate those regions. This contribution should manifest at the edge of the signal. Hence, we will refer only on the integration with $G_{Q}$ as the QNM part.

Moreover, the problem is not fully symmetric around the potential peak $x=0$ . Although the Green’s function itself is symmetric, our particular source term, the worldline of the particle, freezes in proximity to the horizon $x\rightarrow-\infty$ . Therefore, only at late times, for equal $\Delta x$ around the potential near the horizon or near the asymptotic flat region, will a QNM be observed. The same holds true for flat contributions, see (d) and (f) for comparison.

V First-order electromagnetic perturbation due to point charges

We are now in a good position to proceed with solving for the electromagnetic perturbation for radial free fall case. Generically, the path integrals in III.2 may present challenges; therefore, the choice between a semi-analytical approximation or a numerical solution depends on which aspects of the solution are of primary interest.

V.1 Point charge in radial free fall

For simplicity, we assume homogeneous initial conditions on ${}_{e}\Phi^{(1)}_{lm}(r,t)$ :

	$\displaystyle{}_{e}\Phi^{(1)}_{lm}(r,t)$	$\displaystyle=\frac{-q}{\lambda}Y_{lm}^{*}(\theta_{0},\phi_{0})\hat{\Phi}^{(1)}(t,r),$		(48)
	$\displaystyle\hat{\Phi}^{(1)}(t,r)$	$\displaystyle=I_{1}(t,r)+I_{2}(t,r),$		(48)

it’s important to note that $\hat{\Phi}^{(1)}(t,r)$ in this particular case is independent of the multipole numbers $(l,m)$ . The integrals $I_{1}(t,r)$ and $I_{2}(t,r)$ are defined as follows⁶⁶6where $\Theta(t-|x|)$ could be factored out just for clarity, although it should still be embedded in the $B$ integrals regardless:

	$\displaystyle I_{1}$	$\displaystyle=\frac{3}{2}\int_{-\infty}^{\infty}dt^{\prime}\,G\,f(r_{p}(t^{\prime}))^{2}\frac{df(r_{p}(t^{\prime}))}{dr_{p}(t^{\prime})}$
		$\displaystyle\equiv\frac{3}{4}\bigg{[}I_{F1A}-\Theta I_{F1B}+\Theta I_{Q1B}e^{-\frac{V_{EM}}{2}(t-\|x\|)}\bigg{]},$		(49)

and

	$\displaystyle I_{2}$	$\displaystyle=\int_{-\infty}^{\infty}dt^{\prime}\,\frac{\partial G}{\partial r^{\prime}}\bigg{\|}_{r^{\prime}=r_{p}(t^{\prime})}f(r_{p}(t^{\prime}))^{3}$
		$\displaystyle\equiv-\frac{1}{2}\bigg{[}sI_{F2A}-\Theta\tilde{s}\frac{V_{EM}}{2}e^{-\frac{V_{EM}}{2}(t-\|x\|)}I_{Q2B}\bigg{]}.$		(50)

The “amplitude” $A(t,x)$ of the QNM part can be expressed as

A(t,x)=\frac{3}{4}I_{Q1B}+\frac{\tilde{s}}{4}V_{EM}I_{Q2B},

However, evaluating $I_{Q1B}$ and $I_{Q2B}$ may not be as straightforward as the other integrals, as the reader can see in Appendix B. While these integrals can be computed numerically, if we solve for $r_{pBmax}(t,x)$ numerically. For simplicity, we will proceed with I.C. on the trajectory that $t=0$ at $r_{pBmin}=r_{0}$ ⁷⁷7if we are interested in the case starting far from the black hole, i.e. $r_{0}\rightarrow\infty$ , this will require that $t$ runs $-\infty<t<\infty$ ., and also that $t$ is positive definite. Since the integration limits are often easier to handle in double null coordinates $u^{\prime}$ and $v^{\prime}$ , we can proceed with the following transformation of variables in the right half on the $(t^{\prime},x^{\prime})$ plane as

	$\displaystyle v$	$\displaystyle\mapsto v(r_{p})=t(r_{p})+x(r_{p}),$
	$\displaystyle-dv$	$\displaystyle\,f(1+\beta)=\frac{f}{\beta}dr_{p},$

while for the left half we have

	$\displaystyle u$	$\displaystyle\mapsto u(r_{p})=t(r_{p})-x(r_{p}),$
	$\displaystyle-du$	$\displaystyle\,f(1-\beta)=\frac{f}{\beta}dr_{p}.$

The integration limits are a bit tricky as we encounter three different cases:

1.

If we are interested in the case where the particle starts from a finite position $x^{\prime}(r_{0})>t-|x|$ , then particle never enters the past light cone of $(0,t-|x|)$ . In this case, it would have zero QNM contribution for $(t,x)$ .

The particle exits the past light cone before reaching the potential peak at $x=0$ . There is no contribution from left half of the light cone.

\begin{gathered}r_{pBmin}(t,r)=r_{0}\mapsto v=x(r_{0})\\ r_{pBmax}(t,r)\mapsto v=t-|x|.\end{gathered}

(51)

A_{2}(t,x)=\int_{t-|x|}^{x(r_{0})}dv^{\prime}\,(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{1}}{4}V_{EM}\right)e^{v^{\prime}\frac{V_{EM}}{2}},

(52)

The particle exits after passing $x=0$ the past light cone. Then right side contribution in $v^{\prime}$ :

\begin{gathered}r_{pBmin}(t,r)=r_{0}\mapsto v=x(r_{0})\\ \quad\quad v=t^{\prime\prime},\end{gathered}

(53)

where $t^{\prime\prime}$ represents the time at which the particle passes the peak of the potential at $x=0$ . The left side contribution in $u^{\prime}$ is

\begin{gathered}u=t^{\prime\prime}.\\ r_{pBmax(t,r)}\mapsto v=t-|x|.\end{gathered}

(54)

Hence, the integral should be divided into two halves before and after $x=0$

	$\displaystyle A_{3}(t,x)$	$\displaystyle=\int_{t-\|x\|}^{t^{\prime\prime}}du^{\prime}\,(\beta-1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{u^{\prime}\frac{V_{EM}}{2}}$		(55)
		$\displaystyle+\int_{t^{\prime\prime}}^{x(r_{0})}dv^{\prime}\,(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{v^{\prime}\frac{V_{EM}}{2}}.$		(55)

While the integration limits are now straightforward, this has resulted in a complex integrand, as was expected since the worldline is timelike. In other words, because it’s not easy to invert the function $u^{\prime}(r_{p})$ (or $v^{\prime}(r_{p})$ ) analytically. However, for some point $(t,x)$ such that $t\approx|x|$ , we can expand around $u=0$ ⁸⁸8When $x>0$ , we can use $u$ , for $x<0$ , we use $v$ . Physically we can find the initial behavior of the QNM at some position $x$ , however as $t$ increases, higher terms in the expansions need to be considered. A Taylor expansion of the segment of the function of $r_{p}(u)$ around $u=0$ up to the order $n$ can be integrated as follows

\int P_{n}(u)e^{u\frac{V_{EM}}{2}}=e^{u\frac{V_{EM}}{2}}\sum_{i=0}^{n}c_{i}\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{u^{i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}},

where $c_{i}$ are the coefficient of the expansion

(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)

around $u=0$ that could easily be obtained through applying the chain rule.

Back to case (2), evaluating the integrand above at the upper limit $x(r_{0})$ will contribute to the QNM part with amplitude that only depends on the I.C. of the worldline. However the lower limit $t-|x|$ will result in a flat contribution which is $(t,x)$ dependent.

		$\displaystyle e^{-\frac{V_{EM}}{2}(t-\|x\|)}A_{2}(t,x)=e^{-\frac{V_{EM}}{2}(t-\|x\|)}\mathbf{A}_{2}(r_{0})-\sum_{i=0}^{n}c_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{(t-\|x\|)^{i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right),$		(56)
		$\displaystyle\mathbf{A}_{2}(r_{0})=e^{x(r_{0})\frac{V_{EM}}{2}}\sum_{i=0}^{n}c_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{x(r_{0})^{i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right).$		(56)

The case (3) should be very similar to the case (2), and we will end up with QNM part with amplitude only dependent on the initial conditions for the particle and a flat part polynomial in $t-|x|$ . However, the coefficients of expansion will be different.

		$\displaystyle e^{-\frac{V_{EM}}{2}(t-\|x\|)}A_{3}(t,x)=e^{-\frac{V_{EM}}{2}(t-\|x\|)}\mathbf{A}_{3}(r_{0})-\sum_{i=0}^{n}f_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{(t-\|x\|)^{i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right),$		(57)
		$\displaystyle\mathbf{A}_{3}(r_{0})=e^{x(r_{0})\frac{V_{EM}}{2}}\sum_{i=0}^{n}d_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{x(r_{0})^{i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right)-e^{t^{\prime\prime}\frac{V_{EM}}{2}}\sum_{i=0}^{n}d_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{t^{\prime\prime\,i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right)$
		$\displaystyle+e^{t^{\prime\prime}\frac{V_{EM}}{2}}\sum_{i=0}^{n}f_{i}\left(\sum_{k=0}^{i}(-1)^{k}\frac{i!}{(i-k)!}\frac{t^{\prime\prime\,i-k}}{\left(\frac{V_{EM}}{2}\right)^{k}}\right).$

Finally, the Maxwell tensor $F$ describing this perturbation can be written as

F=F_{\text{monopole}}+{}_{e1}F_{\text{rad}}+{}_{e2}F_{\text{int}},

(58)

where the monopole contribution as well as the radiative and “internal” contribution are expressed in terms of $\Phi^{(1)}$ and the current. The monopole contribution is simply $F_{\text{monopole}\,tr}=q\frac{\Theta(r-r_{p}(t))}{r^{2}}$ . The radiative contribution is

\begin{aligned} {}_{e1}F_{\text{rad}}=q\begin{pmatrix}0&**&**&**\\ \left[\delta^{2}(\theta,\theta_{0},\phi,\phi_{0})-\frac{1}{4\pi}\right]\frac{\hat{\Phi}^{(1)}}{r^{2}}&0&**&**\\ -\frac{\partial G_{S^{2}}}{\partial\theta}f\frac{\partial\hat{\Phi}^{(1)}}{\partial r}&-\frac{\partial G_{S^{2}}}{\partial\theta}\frac{1}{f}\frac{\partial\hat{\Phi}^{(1)}}{\partial t}&0&0\\ -\frac{\partial G_{S^{2}}}{\partial\phi}f\frac{\partial\hat{\Phi}^{(1)}}{\partial r}&-\frac{\partial G_{S^{2}}}{\partial\phi}\frac{1}{f}\frac{\partial\hat{\Phi}^{(1)}}{\partial t}&0&0\end{pmatrix}\end{aligned}.

(59)

While internal contribution is given by

\displaystyle{}_{e2}F_{\text{int}}=q\begin{pmatrix}0&0&**&**\\ 0&0&**&**\\ -\frac{\partial G_{S^{2}}}{\partial\theta}f\delta(r-r_{p}(t))&-\frac{\partial G_{S^{2}}}{\partial\theta}\beta\,\delta(r-r_{p}(t))&0&0\\ -\frac{\partial G_{S^{2}}}{\partial\phi}f\delta(r-r_{p}(t))&-\frac{\partial G_{S^{2}}}{\partial\phi}\beta\,\delta(r-r_{p}(t))&0&0\end{pmatrix}.

(60)

V.2 Radial ideal dipole in radial free fall

For simplicity, we will focus on the scenario of radial free fall, and the radial dipole. Additionally, we will disregard the effects of electromagnetic interactions and self-interactions between the two particles along their worldlines. Instead, we consider two charges, $q$ and $-q$ , with a separation $\eta$ along their radial worldlines into the black hole, starting from rest at large distance positions $r_{0}$ and $r_{0}+\eta$ respectively with $E=1$ .

We can derive expressions for an ideal radial dipole in the limit as $\eta\rightarrow 0$ as a variation of this expression with respect to $r_{0}$ or equivalently with respect to $t^{\prime\prime}$ , denoted as $\delta A(t,x)$ . Focusing on the QNM contribution, we will proceed with similar cases as we had before. In the dipole case (2) only segments of wordline $x_{p}>0$ contribute to the trajectory

\displaystyle\frac{1}{\eta}\frac{\delta A_{2}}{\delta r_{0}}=\left[(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{v^{\prime}\frac{V_{EM}}{2}}\right]_{r=r_{0}}+\int_{t-|x|}^{x(r_{0})}dv^{\prime}\,\frac{d}{dr_{p}}\left[(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)\right]e^{v^{\prime}\frac{V_{EM}}{2}}.

(61)

$\displaystyle\frac{1}{\eta}\frac{\delta A_{3}}{\delta r_{0}}$	$\displaystyle=\frac{r_{0}^{3/2}}{2(r_{0}-2)}\left[(\beta-1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{u^{\prime}\frac{V_{EM}}{2}}-(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{v^{\prime}\frac{V_{EM}}{2}}\right]_{r=r_{0}}$	(62)
	$\displaystyle-f(r_{0})\left[(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)e^{v^{\prime}\frac{V_{EM}}{2}}\right]_{r=r_{0}}+\int_{t-\|x\|}^{t^{\prime\prime}}du^{\prime}\,\frac{d}{dr_{p}}\left[(\beta-1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)\right]e^{u^{\prime}\frac{V_{EM}}{2}}$
	$\displaystyle+\int_{t^{\prime\prime}}^{x(r_{0})}dv^{\prime}\,\frac{d}{dr_{p}}\left[(\beta+1)f\left(\frac{3}{4}\frac{df}{dr_{p}}+\frac{\tilde{s}}{4}V_{EM}\right)\right]e^{v^{\prime}\frac{V_{EM}}{2}},$

where $\frac{\delta t^{\prime\prime}}{\delta r_{0}}=\frac{r_{0}^{3/2}}{2(r_{0}-2)}\eta$ . Hence, the EM perturbation is

\displaystyle{}_{e}\Phi^{dip(1)}_{lm}(r,t)

\displaystyle=-\frac{q}{\lambda}Y_{lm}^{*}\eta\,\hat{\Phi}^{dip(1)}(r,t),

(63)

where $+$ and $-$ refer to positive and negative charges respectively, and $p=q\eta$ , assuming that the scale of change of $r_{pzero}$ is much larger than $p$ . In a region outside the dipole, there is no monopole contribution. The radiative contribution is represented by

\begin{aligned} {}_{e1}F_{\text{rad}}&=-p\begin{pmatrix}0&**&**&**\\ \left[\delta^{2}(\theta,\theta_{0},\phi,\phi_{0})-\frac{1}{4\pi}\right]\frac{\hat{\Phi}^{dip(1)}(t,r)}{r^{2}}&0&**&**\\ -\frac{\partial G_{S^{2}}}{\partial\theta}f\frac{\partial\hat{\Phi}^{dip(1)}(t,r)}{\partial r}&-\frac{\partial G_{S^{2}}}{\partial\theta}\frac{1}{f}\frac{\partial\hat{\Phi}^{dip(1)}(t,r)}{\partial t}&0&0\\ -\frac{\partial G_{S^{2}}}{\partial\phi}f\frac{\partial\hat{\Phi}^{dip(1)}(t,r)}{\partial r}&-\frac{\partial G_{S^{2}}}{\partial\phi}\frac{1}{f}\frac{\partial\hat{\Phi}^{dip(1)}(t,r)}{\partial t}&0&0\end{pmatrix}\end{aligned},

(64)

where $\delta^{2}(\theta,\theta_{0},\phi,\phi_{0})$ is the Dirac delta function at the radial line of free fall, while $G_{S^{2}}$ is the angular Green’s function on a $S^{2}$ sphere (the expressions can be found in [40]). There is also another contribution to the electromagnetic field, localized at the position of the points charges, which is referred as ${}_{e2}F_{\text{int}}$ , given by

\displaystyle{}_{e2}F_{\text{int}}

\displaystyle=-p\begin{pmatrix}0&0&**&**\\ 0&0&**&**\\ -\frac{\partial G_{S^{2}}}{\partial\theta}f\left[\delta(r-r_{p}(t))-\delta(r-r_{p}(t)-\eta)\right]&-\frac{\partial G_{S^{2}}}{\partial\theta}\beta\,\left[\delta(r-r_{p}(t))-\delta(r-r_{p}(t)-\eta)\right]&0&0\\ -\frac{\partial G_{S^{2}}}{\partial\phi}f\left[\delta(r-r_{p}(t))-\delta(r-r_{p}(t)-\eta)\right]&-\frac{\partial G_{S^{2}}}{\partial\phi}\beta\,\left[\delta(r-r_{p}(t))-\delta(r-r_{p}(t)-\eta)\right]&0&0\end{pmatrix}.

(65)

Similar to $F_{\text{monopole}}$ , the term ${}_{e2}F_{\text{int}}$ vanishes outside the dipole. Consequently, its contribution is ignored when calculating the electromagnetic stress-energy tensor $T_{EM}^{\mu\nu}$ , which will source the gravitational perturbation.

VI Numerical solution

In this section, we discuss the results of numerically integrating the contributions to the QNM and flat electromagnetic perturbations arising from a radially free-falling radial dipole. In the figures, the particle’s position is marked by a red dashed vertical line, while the peak of the potential, $V_{\text{EM}}=\frac{M}{10}$ , is also indicated. The dipole initially begins around $x_{0}=200M$ with a separation $\eta=1M$ .

For the QNM, initially at $t=0M$ , in Figure 3(a), there is no contribution to $\hat{\Phi}^{(1)\text{dip}}_{\text{QNM}}$ . As time progresses, a constant pulse emerges (see Figure 3(b)), propagating symmetrically from the peak of the potential and moving to the left and right. This event occurs when the particle reaches $x_{p}=175.5M$ at $t=231.25M$ . This pulse might be related to the fact that charge $q$ enters the light cone of $(0,t-x)$ , while the other charge, for a short period, remains outside. This highlights two artifacts in this problem: the modeling of the dipole as two discrete charges, and the sharp causality boundaries within the QNM region due to the Dirac delta potential approximation. Additionally, reducing $\eta$ would control the width of this pulse, while a more realistic, extended potential might reduce its amplitude.

Later, as the particle crosses the potential peak at $x=0$ , as shown in (e), around $t=1282M$ , part of $\Phi^{(1)\text{dip}}_{\text{QNM}}$ reflects, while another part transmits with nearly equal amplitudes. The transmitted component continues toward the horizon, following the particle’s trajectory, as shown in (f) and (g). When analyzed on a logarithmic scale, the slope of $\log\left[\Phi^{(1)\text{dip}}_{\text{QNM}}\right]$ can be approximately fitted with a line of slope $-\frac{V_{\text{EM}}}{2}$ , suggesting a QNM excitation with an almost constant amplitude, consistent with our expectations.

Around $t=1596M$ , the particle is nearly frozen at $r=2M$ , and all QNM contributions to the perturbation vanish completely, as also shown in (h). After some time, around $t=1630M$ and $t=1670M$ , two additional pulses and a shock wave are excited at $x=0$ and begin propagating leftward and rightward, as shown in (i). The second shock wave might suggest that another QNM is excited at late times.

For the flat part, as shown in Figure 4(a), (b), and (c), there is a residual field contribution that does not decay and remains constant across the entire spacetime. The field flips sign across the particle, which is expected for a dipole moment. However, this nonlocal behavior is unrealistic, attributable to the fact that we are solving the wave equation in a 1D variable $x$ , while replacing the entire potential with a Dirac delta. This excludes the centrifugal effects, which represent the impact of the other two coordinates. Consequently, if a better approximation were made that preserved the Minkowski limit, we would expect this field to be local and to decay as $r^{-3}$ away from the potential, as a dipole would in flat spacetime. This suggests that a more accurate approximation could be adopted in future studies if we aim to explore the perturbation more thoroughly.

Furthermore, as the particle approaches the potential, the field builds residual strength near the potential, suggesting that some part of the field will be reflected due to the barrier, as shown in (d) and (e). In (f), we observe a portion of this field being reflected. We anticipate that in a more realistic scenario, as the source of the electromagnetic field moves closer to the potential peak, some of its local (non-radiative) fields would begin to decouple and radiate away, as shown in (g) in our simplified model. Once this shockwave has traveled far enough, all flat contributions vanish, as shown in (i).

VII Discussion & Conclusion

In this work, we have studied the electromagnetic perturbations at first $\Phi^{(1)}$ and second order $\Phi^{(2)}$ , using the Regge-Wheeler formalism in the Schwarzschild spacetime. We began by briefly summarizing relevant literature on the first-order electromagnetic perturbations in Schwarzschild coordinates, then reformulated the Faraday tensor $F_{\mu\nu}$ and electromagnetic stress-energy tensor $T_{\mu\nu}$ in terms of the even and odd perturbation scalars ${}_{o,e}\Phi^{(1)}_{lm}$ , as well as the current multipoles $j_{(i)\,lm}$ .

For the second-order perturbations, we derived the effective source terms for the second-order Regge-Wheeler master equation (with $s=1$ ). These sources are quadratic in the first-order gravitational perturbations (both even and odd), given by the Zerilli and Regge-Wheeler master equations (with $s=2$ ), and the electromagnetic perturbations (both odd and even). We accounted only for the electromagnetic current $J$ in these sources, which contributes to induced modes that are linear in gravitational perturbations with an amplitude dependent on the details of the electromagnetic potential. However, for brevity, we did not consider the perturber’s effective stress-energy source term $T^{\mu\nu}_{G}$ in the gravitational perturbations, which could be of some relevance for EMRI studies. We note that including these terms would induce a part of the source in the second-order electromagnetic perturbation that does not depend on first-order perturbations but rather on the physical parameters of the perturber, such as its mass $M$ and charge $q$ , along with a linear electromagnetic perturbation whose amplitude depends on both electromagnetic and gravitational details.

Furthermore, by replacing the Regge-Wheeler potential with a Dirac delta function as well as considering a sample source similar to those used in [1, 35], and representing a product of linear gravitational and linear electromagnetic perturbations of frequencies $\omega^{(1)}=-i\frac{V_{G}}{2}$ and $\Omega^{(1)}=-i\frac{V_{EM}}{2}$ respectively, we illustrated the mixing that takes place at the second-order electromagnetic QNM $\Phi^{(2)}_{QNM}$ with a frequency of $\Omega^{(2)}=-i\frac{V_{EM}+V_{G}}{2}$ .

Studying the mixing of the second-order electromagnetic perturbations complements our study of the mixing of second-order gravitational perturbations in [1]. We showed that, when examining the minimally coupled Einstein-Maxwell field equations perturbatively, second-order corrections are induced that are quadratic in the electromagnetic perturbations gravity sector and quadratic in the electromagnetic and gravitational perturbations for electromagnetic sector. Besides its theoretical significance, this mixing could be relevant for multi-messenger astrophysics, particularly for stellar mass and EMRI mergers observed by Advanced LIGO and the upcoming LISA. For more discussion on that we refer the reader to [1].

Moreover, we noticed that changing the coupling details, for example by considering the Maxwell-Dilaton-Einstein equations, induces a different mixing spectrum. In this context, this mixing phenomenon is relevant not only for testing the minimal coupling principle but also for probing modified theories of gravity.

For the first-order perturbation, we reduced the calculations of the perturbation from a point charge source to a one-dimensional path integral. A suitable choice of parameterization can simplify some of these integrals. By approximating the Regge-Wheeler potential as a Dirac delta function and focusing on two cases — a single charge $q$ and a radial dipole moment $p=q\eta$ — we found that the flat part of the perturbation is easy to derive analytically, though the QNM part requires a semi-analytical approximation or numerical implementation. We carried out both methods in this work. For the dipole case, based on the numerical solution, we concluded that the QNM perturbation is excited with an almost constant amplitude.

Nevertheless, the Dirac delta function approximation presents some limitations, as discussed in the numerical integration section. By approximating the full potential with a Dirac delta, we neglect the centrifugal correction, which encapsulates the Minkowski limit of the field’s behavior and effectively reduces the problem to one dimension. In this 1D context, the field does not decay as a function of $r$ , which is unrealistic for a 3D problem. This same issue leads to singular multipole coefficients in the stress-energy tensor reported in Appendix A of [1]. On the other hand, a less trivial approximation would be harder to handle analytically and would not yield simple modes like those with the Dirac delta function. Thus, the choice of approximation should depend on the specific aspects intended for investigation. In our case, as we are interested in mixing, a more realistic potential would be beneficial, and in future work we intend to revisit the mixing problem in Teukolsky equations in the Kerr background, using numerical approach.

Acknowledgements.

We acknowledge the partial support of D.S. through the U.S. National Science Foundation under Grant No. PHY-2310363. Additionally, we express our gratitude to Macarena Lagos, Ariadna Ribes Metidieri, Ahmed Elshahawy, A.K. Gorbatsievich, and Stanislav Komarov for their insightful discussions. Special thanks to Jacob Fields and Ish Gupta for their invaluable guidance and feedback on this manuscript. We also sincerely appreciate the valuable comments provided by David Radice and Elias Most. Finally, we extend special thanks to Mahmoud A. Mansour for his invaluable time, support, and assistance throughout this project.

Appendix A Green’s function on $S^{2}$

The Green’s function $G_{S^{2}}$ on $S^{2}$ with a unit radius is

$\displaystyle\nabla_{S^{2}}G_{S^{2}}$	$\displaystyle=-Y_{00}(\theta,\phi)Y_{00}^{*}(\theta_{0},\phi_{0})+\frac{\delta(\theta-\theta_{0})\delta(\phi-\phi_{0})}{\sin\theta}$	(66)
	$\displaystyle=\sum_{l>0,m}Y_{lm}(\theta,\phi)Y_{lm}^{*}(\theta_{0},\phi_{0}),$
$\displaystyle G_{S^{2}}$	$\displaystyle=\sum_{l>0,m}\frac{-1}{\lambda}Y_{lm}(\theta,\phi)Y_{lm}^{*}(\theta_{0},\phi_{0})$
	$\displaystyle=\frac{1}{4\pi}\ln\left[1-\cos\gamma\right]$

where $\cos\gamma=\cos\theta\cos\theta_{0}+\sin\theta\sin\theta_{0}\cos(\phi-\phi_{0}),$ and $\gamma$ is the angle between the two points $(\theta,\phi)$ and $(\theta_{0},\phi_{0})$ .

\frac{\partial G_{S^{2}}}{\partial\phi}=-\frac{\sin\theta\sin\theta_{0}\sin(\phi-\phi_{0})}{4\pi(1-\cos\gamma)}

\frac{G_{S^{2}}}{\partial\theta}=\frac{-\cos\theta_{0}\sin\theta+\cos\theta\cos(\phi-\phi_{0})\sin\theta_{0}}{4\pi(1-\cos\gamma)}

Appendix B Integrals

Here, $\epsilon_{r_{pBzero}(t,r)}$ is a binary indicator that equals 1 if $B$ has a zero at a given $(t,r)$ , and 0 otherwise. The quantities $r_{pA_{\min}}(t,r)$ and $r_{pB_{\min}}(t,r)$ denote the minimal positive values of $A$ and $B$ , respectively, while $r_{pA_{\max}}(t,r)$ and $r_{pB_{\max}}(t,r)$ represent their corresponding maximal positive values.

For a given $(t,r)$ , if $A=0$ (or $B=0$ ), it follows that $r_{pA_{\min}}=0$ (or $r_{pB_{\min}}=0$ ). Both $A$ and $B$ are monotonic functions with respect to the parameter $r_{p}$ along the particle’s trajectory for a fixed $(t,r)$ point. The partial derivatives of $A$ and $B$ with respect to $r_{p}$ are positive definite

\frac{\partial A}{\partial r_{p}}=s\beta(r_{p})^{2}+\frac{1}{\beta(r_{p})}-\frac{s}{\beta(r_{p})+1}>0

(67)

\frac{\partial B}{\partial r_{p}}=-\tilde{s}f(r_{p})+\frac{1}{\beta(r_{p})}+\frac{\beta(r_{p})}{f(r_{p})}\left[1-\tilde{s}\beta(r_{p})\right]>0

(68)

Consequently, the causal structure of the Green’s function confines the particle’s support to a $(t,r)$ cone along its trajectory. In summary, the problem reduces to finding max, min and zero (if it exists) for $A=0$ and $B=0$ as a function of $(t,r)$ , as well as evaluating the two integrals $I_{Q1B}$ and $I_{Q2B}$ which resemble part of the contribution of the QNM mode.

\displaystyle I_{F1A}=\int_{2M}^{\infty}dr_{p}\frac{f\left(r_{p}\right)}{\beta\left(r_{p}\right)}\frac{df\left(r_{p}\right)}{dr_{p}}\Theta(A)=\left[\frac{2}{3}\beta^{3}-2\beta\right]\bigg{|}_{r_{pA_{\min}}(t,r)}^{r_{pA_{\max}}(t,r)},

(69)

\displaystyle I_{F1B}=\int_{2M}^{\infty}dr_{p}\frac{f\left(r_{p}\right)}{\beta\left(r_{p}\right)}\frac{df\left(r_{p}\right)}{dr_{p}}\Theta(B)=\left[\frac{2}{3}\beta^{3}-2\beta\right]\bigg{|}_{r_{pB_{\min}}(t,r)}^{r_{pB_{\max}}(t,r)},

(70)

\displaystyle I_{Q1B}=\int_{2M}^{\infty}dr_{p}\,\frac{f(r_{p})}{\beta(r_{p})}\frac{df(r_{p})}{dr_{p}}\,e^{\frac{V_{EM}}{2}\left(t(r_{p})+\left|x(r_{p})\right|\right)}\Theta(B)=\int_{r_{pB_{\min}}(t,r)}^{r_{pB_{\max}}(t,r)}dr_{p}\,\frac{f(r_{p})}{\beta(r_{p})}\frac{df(r_{p})}{dr_{p}}\,e^{\frac{V_{EM}}{2}\left(t(r_{p})+\left|x(r_{p})\right|\right)},

(71)

I_{F2A}=\int_{2M}^{\infty}dr_{p}\frac{f\left(r_{p}\right)}{\beta\left(r_{p}\right)}\delta(A)=\epsilon_{r_{pAzero}(t,r)}\frac{f\left(r_{pAzero}(t,r)\right)}{\beta\left(r_{pAzero}(t,r)\right)}

(72)

I_{Q2B}=\int_{2M}^{\infty}dr_{p}\frac{f\left(r_{p}\right)}{\beta\left(r_{p}\right)}e^{\frac{V_{EM}}{2}\left(t\left(r_{p}\right)+\left|x\left(r_{p}\right)\right|\right)}\Theta(B)=\int_{r_{pB_{\min}}(t,r)}^{r_{pB_{\max}}(t,r)}dr_{p}\frac{f(r_{p})}{\beta(r_{p})}\,e^{\frac{V_{EM}}{2}\left(t(r_{p})+\left|x(r_{p})\right|\right)},

(73)

Appendix C Re-parameterization of perturbation integral equation

We can parameterize the path integral in terms of $r_{p}$ as well:

	$\displaystyle{}_{o}\Phi_{lm}(r^{\prime},t^{\prime})=$	$\displaystyle-\frac{q}{\lambda}\int_{2M}^{\infty}dr_{p}\,f\left(r_{p}\right)\left\{im\frac{d\hat{\theta}}{dr_{p}}\mathcal{C}_{lm}[Y_{lm}^{*}]\right.$		(74)
		$\displaystyle\left.+\frac{d\hat{\phi}}{dr_{p}}\mathcal{D}_{lm}[Y_{lm}^{*}]\right\}\,G_{lm}\left(r^{\prime},t^{\prime},r_{p},t\left(r_{p}\right)\right).$		(74)

_e	$\displaystyle\Phi_{lm}\left(t^{\prime},r^{\prime}\right)=\frac{q}{\lambda}\int^{\infty}_{2M}dr_{p}\Bigg{\{}\left[\frac{d\dot{r}_{p}(r_{p})}{dr_{p}}-\frac{f\left(r_{p}\right)}{\dot{r}_{p}\left(r_{p}\right)}\frac{df\left(r_{p}\right)}{dr_{p}}\right.$	(75)
	$\displaystyle\left.-im\,\dot{r}\left(r_{p}\right)\frac{d\hat{\phi}}{dr_{p}}+m\cot\hat{\theta}\,\dot{r}\left(r_{p}\right)\frac{d\hat{\theta}}{dr_{p}}\right]Y_{lm}^{*}$
	$\displaystyle+\left[\sqrt{(l-m)(l+m+1)}e^{i\hat{\phi}}\dot{r}\left(r_{p}\right)\frac{d\hat{\theta}}{dr_{p}}\right]Y_{l,m+1}^{*}\Bigg{\}}$
	$\displaystyle\times G_{lm}\left(r^{\prime},t^{\prime},r_{p},t\left(r_{p}\right)\right)$
	$\displaystyle+\frac{q}{\lambda}\int^{\infty}_{2M}dr_{p}\left.\frac{G_{lm}}{dr}\right\|_{r=r_{p}}\left[\dot{r}_{p}-\frac{f^{2}\left(r_{p}\right)}{\dot{r}_{p}}\right]Y_{lm}^{*}.$

where $\hat{\phi}(r_{p})=\phi(t(r_{p}))$ and $\hat{\theta}(r_{p})=\theta(t(r_{p}))$ . This re-parameterization in terms of $r_{p}$ is valid as long as the particle undergoes radial motion in Schwarzschild coordinates. Still, if the particle passes through the same radial point at different moments in time, the integral needs to be divided into the one-to-on domains and be handled separately. On the other hand, for a purely circular trajectory, we should re-parameterize using the angular coordinate instead.

Appendix D 4D vector harmonics

In this appendix, spherical decomposition in [23] is extended to 4D vector representations, following the same approach of Lousto and Barack [41], but for vectors instead of tensors. The harmonics are defined on the two-dimensional sphere $S^{2}$ , where scalar spherical harmonics $Y_{lm}(\theta,\phi)$ orthogonality are extended to the 4D vector harmonics $\mathbf{Y}_{\mu\,lm}(\theta,\phi)$ . The vector harmonics $Y^{\mu}_{(1)\,lm},Y^{\mu}_{(2)\,lm},Y^{\mu}_{(3)\,lm}$ are categorized as even-parity, while $Y^{\mu}_{(4)\,lm}$ belongs to the odd-parity class. These are expressed as:

\mathbf{Y}^{\mu}_{(1)\,lm}=\left(\begin{array}[]{c}Y_{lm}\\ 0\\ 0\\ 0\end{array}\right),\quad\mathbf{Y}^{\mu}_{(2)\,lm}=\left(\begin{array}[]{c}0\\ Y_{lm}\\ 0\\ 0\end{array}\right),

(76)

\mathbf{Y}^{\mu}_{(3)\,lm}=\frac{1}{\lambda}\left(\begin{array}[]{c}0\\ 0\\ \dot{Y}_{lm}\\ \csc^{2}\theta Y_{lm}^{\prime}\end{array}\right),\quad\mathbf{Y}^{\mu}_{(4)\,lm}=\frac{1}{\lambda}\left(\begin{array}[]{c}0\\ 0\\ \csc\theta Y_{lm}^{\prime}\\ -\csc\theta\dot{Y}_{lm}\end{array}\right).

(77)

where $\dot{Y}_{lm}$ and $Y_{lm}^{\prime}$ represent the derivatives with respect to $\theta$ and $\phi$ , respectively. Their is $l$ and $m$ sub-index label on the vector harmonics $\mathbf{Y}^{\mu}_{(i)\,lm}$ , but suppressed for simplicity. Using the following two identities:

\int_{S^{2}}d\Omega\csc\theta\,(Y^{*\prime}_{lm}\,\dot{Y}_{lm}-Y^{\prime}_{lm}\,\dot{Y}^{*}_{lm})=0,

(78)

and

\int_{S^{2}}d\Omega(\csc^{2}\theta\,Y^{\prime}_{lm}\,Y^{*\prime}_{lm}+\dot{Y}_{lm}\,\dot{Y}^{*}_{lm})=\delta_{ll^{\prime}}\delta_{mm^{\prime}}\lambda.

(79)

The vector harmonics orthogonality condition is given by:

\int_{S^{2}}\mathbf{Y}^{\mu}_{(i)lm}\mathbf{Y}_{\mu(j)l^{\prime}m^{\prime}}\,d\Omega=\delta_{ij}\delta_{ll^{\prime}}\delta_{mm^{\prime}}N_{i}(r),

(80)

where $\sqrt{N_{i}(r)}$ represents the norm of the $i$ -th vector harmonic. The components of the even-parity vectors ${}_{e}\mathbf{J}$ and ${}_{e}\mathbf{A}$ are expressed as:

{}_{e}\mathbf{J}_{\mu}=\begin{pmatrix}j_{1}Y\\ j_{2}Y\\ j_{3}\dot{Y}\\ j_{3}Y^{\prime}\end{pmatrix}\quad{}_{e}\mathbf{A}_{\mu}=\begin{pmatrix}gY\\ hY\\ k\dot{Y}\\ kY^{\prime}\end{pmatrix}

(81)

Similarly, the components of the odd-parity vectors ${}_{o}\mathbf{A}_{\mu}$ and ${}_{o}\mathbf{J}_{\mu}$ are given by:

{}_{o}\mathbf{J}_{\mu}=\begin{pmatrix}0\\ 0\\ j_{4}\csc\theta Y^{\prime}\\ -j_{4}\sin\theta\,\dot{Y}\end{pmatrix}\quad{}_{o}\mathbf{A}_{\mu}=\begin{pmatrix}0\\ 0\\ a\csc\theta Y^{\prime}\\ -a\sin\theta\,\dot{Y}.\end{pmatrix}

(82)

The summation over the angular $l$ and azimuthal $m$ numbers is implied but not explicitly shown. The functions $j_{i},a,g,k$ are dependent on the coordinates $(t,r)$ and the sub-index $lm$ .

References

[1] Fawzi Aly, Mahmoud A. Mansour, and Dejan Stojkovic. More nonlinearities? electromagnetic and gravitational mode mixing in nsbh mergers, 2024. URL: https://arxiv.org/abs/2410.12775, arXiv:2410.12775.
[2] Emanuele Berti, Vitor Cardoso, and Andrei O Starinets. Quasinormal modes of black holes and black branes. Classical and Quantum Gravity, 26(16):163001, July 2009. URL: http://dx.doi.org/10.1088/0264-9381/26/16/163001, doi:10.1088/0264-9381/26/16/163001.
[3] Kostas D. Kokkotas and Bernd G. Schmidt. Quasi-normal modes of stars and black holes. Living Reviews in Relativity, 2(1), September 1999. URL: http://dx.doi.org/10.12942/lrr-1999-2, doi:10.12942/lrr-1999-2.
[4] Roman A Konoplya and Alexander Zhidenko. Quasinormal modes of black holes: From astrophysics to string theory. Reviews of Modern Physics, 83(3):793–836, 2011.
[5] Valeria Ferrari and Leonardo Gualtieri. Quasinormal modes and gravitational wave astronomy. General Relativity and Gravitation, 40(6):945–970, 2008.
[6] Peter James Nee, Sebastian H. Völkel, and Harald P. Pfeiffer. Role of black hole quasinormal mode overtones for ringdown analysis. Physical Review D, 108(4), August 2023. URL: http://dx.doi.org/10.1103/PhysRevD.108.044032, doi:10.1103/physrevd.108.044032.
[7] Vishal Baibhav, Mark Ho-Yeuk Cheung, Emanuele Berti, Vitor Cardoso, Gregorio Carullo, Roberto Cotesta, Walter Del Pozzo, and Francisco Duque. Agnostic black hole spectroscopy: Quasinormal mode content of numerical relativity waveforms and limits of validity of linear perturbation theory. Physical Review D, 108(10), November 2023. URL: http://dx.doi.org/10.1103/PhysRevD.108.104020, doi:10.1103/physrevd.108.104020.
[8] Abbott et al. (LIGO Scientific Collaboration-Virgo Collaboration-KAGRA Collaboration). Observation of gravitational waves from a binary black hole merger. Phys. Rev. Lett., 116:061102, Feb 2016. URL: https://link.aps.org/doi/10.1103/PhysRevLett.116.061102, doi:10.1103/PhysRevLett.116.061102.
[9] Abbott et al. (LIGO Scientific Collaboration-Virgo Collaboration-KAGRA Collaboration). Observation of gravitational waves from two neutron star–black hole coalescences. The Astrophysical Journal Letters, 915(1):L5, June 2021. URL: http://dx.doi.org/10.3847/2041-8213/ac082e, doi:10.3847/2041-8213/ac082e.
[10] Vitor Cardoso and José P. S. Lemos. Quasinormal modes of schwarzschild–anti-de sitter black holes: Electromagnetic and gravitational perturbations. Physical Review D, 64(8), September 2001. URL: http://dx.doi.org/10.1103/PhysRevD.64.084017, doi:10.1103/physrevd.64.084017.
[11] Andrii Neronov. Introduction to multi-messenger astronomy. Journal of Physics: Conference Series, 1263(1):012001, jun 2019. URL: https://dx.doi.org/10.1088/1742-6596/1263/1/012001, doi:10.1088/1742-6596/1263/1/012001.
[12] R. Abbott et. al. Multi-messenger observations of a binary neutron star merger*. The Astrophysical Journal Letters, 848(2):L12, oct 2017. URL: https://dx.doi.org/10.3847/2041-8213/aa91c9, doi:10.3847/2041-8213/aa91c9.
[13] Ben Margalit and Brian D. Metzger. Constraining the maximum mass of neutron stars from multi-messenger observations of gw170817. The Astrophysical Journal Letters, 850(2):L19, nov 2017. URL: https://dx.doi.org/10.3847/2041-8213/aa991c, doi:10.3847/2041-8213/aa991c.
[14] Marica Branchesi. Multi-messenger astronomy: gravitational waves, neutrinos, photons, and cosmic rays. Journal of Physics: Conference Series, 718(2):022004, may 2016. URL: https://dx.doi.org/10.1088/1742-6596/718/2/022004, doi:10.1088/1742-6596/718/2/022004.
[15] Nils Andersson and Kostas D Kokkotas. Quasinormal modes of stars and black holes. Physical Review Letters, 77(20):4134–4137, 1996.
[16] Elias R. Most, Antonios Nathanail, and Luciano Rezzolla. Electromagnetic emission from blitzars and its impact on non-repeating fast radio bursts. The Astrophysical Journal, 864(2):117, sep 2018. URL: https://dx.doi.org/10.3847/1538-4357/aad6ef, doi:10.3847/1538-4357/aad6ef.
[17] Elias R. Most, Andrei M. Beloborodov, and Bart Ripperda. Monster shocks, gamma-ray bursts, and black hole quasi-normal modes from neutron-star collapse. The Astrophysical Journal Letters, 974(1):L12, October 2024. URL: http://dx.doi.org/10.3847/2041-8213/ad7e1f, doi:10.3847/2041-8213/ad7e1f.
[18] BP Abbott et al. Gw170817: Observation of gravitational waves from a binary neutron star inspiral. Physical Review Letters, 119(16):161101, 2017.
[19] Fawzi Aly and Dejan Stojkovic. In horizon penetrating coordinates: Kerr black hole metric perturbation, construction and completion. Classical and Quantum Gravity, 40(23):235010, nov 2023. URL: https://dx.doi.org/10.1088/1361-6382/ad0495, doi:10.1088/1361-6382/ad0495.
[20] Alejandro Cárdenas-Avendaño and Carlos F. Sopuerta. Testing Gravity with Extreme-Mass-Ratio Inspirals, page 275–359. Springer Nature Singapore, 2024. URL: http://dx.doi.org/10.1007/978-981-97-2871-8_8, doi:10.1007/978-981-97-2871-8_8.
[21] S. Amaro et al. Astrophysics with the laser interferometer space antenna. Living Reviews in Relativity, 26(1), March 2023. URL: http://dx.doi.org/10.1007/s41114-022-00041-y, doi:10.1007/s41114-022-00041-y.
[22] Emanuele Berti, Vitor Cardoso, and Clifford M. Will. Gravitational-wave spectroscopy of massive black holes with the space interferometer lisa. Physical Review D, 73(6), March 2006. URL: http://dx.doi.org/10.1103/PhysRevD.73.064030, doi:10.1103/physrevd.73.064030.
[23] Karl Martel and Eric Poisson. Gravitational perturbations of the schwarzschild spacetime: A practical covariant and gauge-invariant formalism. Physical Review D, 71(10), May 2005. URL: http://dx.doi.org/10.1103/PhysRevD.71.104003, doi:10.1103/physrevd.71.104003.
[24] karl Martel. Particles and black holes: Time-domain integration of the equations of black-hole perturbation theory(Doctoral Dissertation). PhD thesis, University of Guelph, CANADA, 2004. URL: https://ui.adsabs.harvard.edu/abs/2004PhDT.........4M/abstract.
[25] Ulrich H. Gerlach and Uday K. Sengupta. Gauge-invariant perturbations on most general spherically symmetric space-times. Phys. Rev. D, 19:2268–2272, Apr 1979. URL: https://link.aps.org/doi/10.1103/PhysRevD.19.2268, doi:10.1103/PhysRevD.19.2268.
[26] A. Ronveaux. Heun’s differential equations, chapter Confluent Heun Equation, pages 87–126. Oxford University press, 2007.
[27] Roumen S. Borissov and Plamen P. Fiziev. Exact solutions of teukolsky master equation with continuous spectrum. Bulgarian Journal of Physics, 37:65–68, 2010. arXiv:0903.3617.
[28] Subrahmanyan Chandrasekhar. The Mathematical Theory of Black Holes. Clarendon Press, Jan 1983.
[29] Jerry B. Griffiths and Jiří Podolský. Exact Space-Times in Einstein’s General Relativity. Cambridge University Press, Aug 2012.
[30] Fawzi Aly and Dejan Stojkovic. On the generalization of the kruskal–szekeres coordinates: a global conformal charting of the reissner–nordström spacetime. Classical and Quantum Gravity, 41(13):135005, may 2024. URL: https://dx.doi.org/10.1088/1361-6382/ad4dff, doi:10.1088/1361-6382/ad4dff.
[31] J. Tiomno. Maxwell equations in a spherically symmetric black-hole background and radiation by a radially moving charge. Lett. Nuovo Cim., 5S2:851–855, 1972. doi:10.1007/BF02832806.
[32] R. Ruffini. Fully relativistic treatment of the bremsstrahlung radiation from a charge falling in a strong gravitational field. Phys. Lett. B, 41:334–338, 1972. doi:10.1016/0370-2693(72)90591-6.
[33] Costantino Pacilio and Richard Brito. Quasinormal modes of weakly charged einstein-maxwell-dilaton black holes. Physical Review D, 98(10), November 2018. URL: http://dx.doi.org/10.1103/PhysRevD.98.104042, doi:10.1103/physrevd.98.104042.
[34] Eric W. Hirschmann, Luis Lehner, Steven L. Liebling, and Carlos Palenzuela. Black hole dynamics in einstein-maxwell-dilaton theory. Physical Review D, 97(6), March 2018. URL: http://dx.doi.org/10.1103/PhysRevD.97.064032, doi:10.1103/physrevd.97.064032.
[35] Macarena Lagos and Lam Hui. Generation and propagation of nonlinear quasinormal modes of a schwarzschild black hole. Physical Review D, 107(4), February 2023. URL: http://dx.doi.org/10.1103/PhysRevD.107.044040, doi:10.1103/physrevd.107.044040.
[36] Lam Hui, Daniel Kabat, and Sam S.C. Wong. Quasinormal modes, echoes and the causal structure of the green’s function. Journal of Cosmology and Astroparticle Physics, 2019(12):020–020, December 2019. URL: http://dx.doi.org/10.1088/1475-7516/2019/12/020, doi:10.1088/1475-7516/2019/12/020.
[37] Nikodem Szpak. Quasinormal mode expansion and the exact solution of the cauchy problem for wave equations, 2004. URL: https://arxiv.org/abs/gr-qc/0411050, arXiv:gr-qc/0411050.
[38] Edward W. Leaver. Spectral decomposition of the perturbation response of the schwarzschild geometry. Phys. Rev. D, 34:384–408, Jul 1986. URL: https://link.aps.org/doi/10.1103/PhysRevD.34.384, doi:10.1103/PhysRevD.34.384.
[39] Satoshi Okuzumi, Kunihito Ioka, and Masa-aki Sakagami. Possible discovery of a nonlinear tail and second-order quasinormal modes in black hole ringdown. Physical Review D, 77(12), June 2008. URL: http://dx.doi.org/10.1103/PhysRevD.77.124018, doi:10.1103/physrevd.77.124018.
[40] Fawzi Aly and Dejan Stojkovic. Electromagnetic radiation on schwarzschild background from radial free falling ideal dipole, 2024. in preparation.
[41] Leor Barack and Carlos O. Lousto. Perturbations of schwarzschild black holes in the lorenz gauge: Formulation and numerical implementation. Physical Review D, 72(10), November 2005. URL: http://dx.doi.org/10.1103/PhysRevD.72.104026, doi:10.1103/physrevd.72.104026.

More Nonlinearities? II. A Short Guide to First- and Second-Order Electromagnetic Perturbations in the Schwarzschild background.