Analytic Kludge Waveforms for Extreme Mass Ratio Inspirals of Charged Object around Kerr-Newman Black Hole

The KN BH is a stationary, axisymmetrical, and asymptotically flat solution of the Einstein-Maxwell equation. In Boyer-Lindquist coordinates, the KN metric can be written as

	$\displaystyle ds^{2}$	$\displaystyle=$	$\displaystyle\frac{\Sigma}{\Delta}dr^{2}+\Sigma d\theta^{2}+\frac{\sin^{2}\theta}{\Sigma}[(r^{2}+\tilde{a}^{2})d\phi-\tilde{a}dt]^{2}$		(1)
		$\displaystyle-$	$\displaystyle\frac{\Delta}{\Sigma}[\tilde{a}\sin^{2}\theta d\phi-dt]^{2},$		(2)

where

	$\displaystyle\Sigma(r,\theta)=r^{2}+\tilde{a}^{2}\cos^{2}\theta,$		(3)
	$\displaystyle\Delta(r)=r^{2}-2Mr+\tilde{a}^{2}+\tilde{Q}^{2},$		(4)

and $M$ and $\tilde{Q}$ are the mass and electric charge of the BH, and $\tilde{a}$ is specific angular momentum. The electromagnetic potential is given by

$$A=\frac{\tilde{Q}r}{\Sigma}(dt-\tilde{a}\sin^{2}\theta d\phi).$$

(5)

The first-order equations of motion for a timelike charged particle with mass $m$ and electric charge $\tilde{q}$ are given by Hackmann:2013pva

	$\displaystyle\Sigma\frac{dt}{d\tau}$	$\displaystyle=$	$\displaystyle\frac{r^{2}+M^{2}a^{2}}{\Delta}P-a(aE\sin^{2}\theta-L_{z})M^{2},$		(6)
	$\displaystyle\Sigma\frac{dr}{d\tau}$	$\displaystyle=$	$\displaystyle\pm\sqrt{R},$		(7)
	$\displaystyle\Sigma\frac{d\theta}{d\tau}$	$\displaystyle=$	$\displaystyle\pm M\sqrt{\Theta},$		(8)
	$\displaystyle\Sigma\frac{d\phi}{d\tau}$	$\displaystyle=$	$\displaystyle\frac{aM}{\Delta}P-aME+\frac{ML_{z}}{\sin^{2}\theta},$		(9)

where

	$\displaystyle R$	$\displaystyle=$	$\displaystyle P^{2}-\Delta[r^{2}+M^{2}(L_{z}-aE)^{2}+M^{2}C],$
	$\displaystyle\Theta$	$\displaystyle=$	$\displaystyle C-\left[(1-E^{2})a^{2}+\frac{L_{z}^{2}}{\sin^{2}\theta}\right]\cos^{2}\theta,$		(10)

with

$$P=E(r^{2}+M^{2}a^{2})-M^{2}aL_{z}-MqQr.$$

(11)

In the equations of motion there are three constants of motion. $\tilde{E}$ and $\tilde{L}_{z}$ are the conserved total energy and component of angular momentum parallel to symmetry axis, respectively, and $\tilde{C}$ is the Carter constant Carter:1968rr . Note that in above expressions, we have used dimensionless quantities defined by

	$\displaystyle E=\frac{\tilde{E}}{m},\quad a=\frac{\tilde{a}}{M},\quad q=\frac{\tilde{q}}{m},\quad L_{z}=\frac{\tilde{L}_{z}}{mM},$
	$\displaystyle C=\frac{\tilde{C}}{m^{2}M^{2}},\quad Q=\frac{\tilde{Q}}{M}.$		(12)

One can see that the distinction between the above equations and the geodesic equations for a neutral particle moving in Kerr spacetime only reflects in the functions $\Delta$ and $P$.

In the AK model, if the spin of the CO is neglected, an EMRI event is completely specified by $14$ degrees of freedom. To obtain the inspiral orbits of the charged CO in the KN background, we append two charge parameter $Q$ and $q$ to the original AK model parameter space, such that

	$\displaystyle\lambda^{i}\equiv~{}$	$\displaystyle\Big{[}\lambda^{1},\cdots,\lambda^{16}\Big{]}$
	$\displaystyle\quad=\Big{[}$	$\displaystyle m,M,a,q,Q,e_{\rm LSO},\tilde{\gamma}_{0},\Phi,\cos\theta_{S},\phi_{S},\cos\lambda_{I},\alpha_{0},$
		$\displaystyle\cos\theta_{K},\phi_{K},D,t_{0}\Big{]},$		(13)

where the definition and meaning of each parameter can be found in Barack:2003fp .

The trajectories of charged CO are roughly treated as quasi-Keplerian ellipses, which are characterized by the eccentricity $e$, and the radial orbital frequency $\nu$. The instantaneous phase of the CO in the orbit is specified by the mean anomaly $\Phi$. The orientation of the orbit is described by three angles, $\lambda_{I}$, the inclination angle of the orbital plane with respect to the BH’s spin direction $\hat{S}$, $\tilde{\gamma}$, the angle from pericenter to $\hat{L}\times\hat{S}$ and $\alpha$ describing the direction of $\hat{L}$ around $\hat{S}$, where $\hat{S}$ is a unit vector of BH’s spin and $\hat{L}$ is a unit vector of the orbital angular momentum. The rate of change $\dot{\Phi}$, the orbital plane precession (also known as Lense-Thirring precession) frequency $\dot{\alpha}$ and the angular rate $\dot{\alpha}+\dot{\tilde{\gamma}}$ of periapsis precession are closely related to the fundamental frequencies of the CO’s orbit by

	$\displaystyle\dot{\Phi}$	$\displaystyle=$	$\displaystyle 2\pi\nu=\Omega_{r},$		(14)
	$\displaystyle\dot{\tilde{\gamma}}$	$\displaystyle=$	$\displaystyle\Omega_{\theta}-\Omega_{r},$		(15)
	$\displaystyle\dot{\alpha}$	$\displaystyle=$	$\displaystyle\Omega_{\phi}-\Omega_{\theta},$		(16)

where $\Omega_{r}$, $\Omega_{\theta}$ and $\Omega_{\phi}$ denote the fundamental frequencies of radial, polar and azimuthal motion, respectively. The closed form of these fundamental frequencies for Kerr BH orbits was first obtained by Schmidt Schmidt:2002qk by employing the elegant action-angle variable formalism of the Hamilton-Jacobi theory. Later on, combining Schmidt’s description and using the Mino time Mino:2003yg , Drasco and Hughes Drasco:2003ky derived the fundamental frequencies and showed the construction of the frequency domain representation of arbitrary functions of orbits. In this work we would like to follow the steps in Drasco:2003ky and also Fujita:2009bp ; Gair:2011ym to derive the analytical expressions of these three fundamental frequencies in the weak-field regime.

First of all, in terms of the dimensionless time variable, i.e. the so-called Mino’s time Mino:2003yg

$$\frac{d}{d\lambda}=\frac{\Sigma}{M}\frac{d}{d\tau},$$

(17)

the equations of motion for the charged CO now become

	$\displaystyle\frac{dr}{d\lambda}$	$\displaystyle=$	$\displaystyle\frac{\pm\sqrt{R}}{M},$		(18)
	$\displaystyle\frac{d\theta}{d\lambda}$	$\displaystyle=$	$\displaystyle\pm\sqrt{\Theta},$		(19)
	$\displaystyle\frac{dt}{d\lambda}$	$\displaystyle=$	$\displaystyle T_{r}+T_{\theta},$		(20)
	$\displaystyle\frac{d\phi}{d\lambda}$	$\displaystyle=$	$\displaystyle\Phi_{r}+\Phi_{\theta},$		(21)

where

$$T_{r}=\frac{r^{2}+M^{2}a^{2}}{M\Delta}P,\quad T_{\theta}=-a(aE\sin^{2}\theta-L_{z})M$$

(22)

$$\Phi_{r}=\frac{a}{\Delta}P-aE,\quad\Phi_{\theta}=\frac{L_{z}}{\sin^{2}\theta}.$$

(23)

One can see that the $r$ and $\theta$ motions are now apparently decoupled. Next, we parameterize the orbit in terms of the (dimensionless) semi-latus rectum $p$, the eccentricity $e$ and a phase angle $\psi$ via

$$r=\frac{Mp}{1+e\cos\psi},$$

(24)

where $\psi$ varies from $0$ to $2\pi$ as $r$ goes through a complete cycle. The two turning points of the radial motion,

$$r_{a}=\frac{Mp}{1+e},\quad r_{p}=\frac{Mp}{1-e},$$

(25)

are apoapsis and periapsis of the elliptic orbits, respectively. Moreover, the third orbital parameter is the turning point of the polar motion, $\theta_{tp}$, which is also called the inclination angle, since this is equivalent to $\lambda_{I}$ Gair:2011ym ¹¹1Another usual orbital inclination angle is defined in terms of the Carter constant by $\cos\lambda_{I}=\frac{L_{z}}{\sqrt{C+L_{z}^{2}}}$, which in the limit $a\to 0$ becomes $\cos\lambda_{I}=L_{z}/L$ and one has $\theta_{tp}+(\mathrm{sgn}\,L_{z})\lambda_{I}=\pi/2$, where $L$ is the total orbital angular momentum.. It is useful to express the constants of motion as functions of these three orbital parameters. This is implemented by solving the three equations

$$R(r_{a})=R(r_{p})=\Theta(\cos\theta_{tp})=0.$$

(26)

The asymptotic form of these constants of motion in the weak-field regime are given by

	$\displaystyle E$	$\displaystyle=$	$\displaystyle 1+\frac{1}{2p}(e^{2}-1)(1-qQ)+\mathcal{O}(p^{-2}),$		(27)
	$\displaystyle L_{z}$	$\displaystyle=$	$\displaystyle\sqrt{p}\sqrt{1-qQ}\sin\theta_{tp}+\mathcal{O}(p^{-1/2}),$		(28)
	$\displaystyle C$	$\displaystyle=$	$\displaystyle p(1-qQ)\cos^{2}\theta_{tp}+\mathcal{O}(p^{0}).$		(29)

Note that the terms at higher order in the $1/p$ expansion are not explicitly shown here, but they are important in the following calculations.

For bound orbits, $r(\lambda)$ and $\theta(\lambda)$ become periodic functions. Then from eq. (18), the fundamental period for the radial motion with respect to $\lambda$ is given by

$$\Lambda_{r}=\int_{0}^{\Lambda_{r}}d\lambda=2\int_{r_{a}}^{r_{p}}\frac{Mdr}{\sqrt{R}}=\int_{0}^{2\pi}\frac{d\psi}{\sqrt{V_{\psi}}},$$

(30)

where we have transformed the variable of the integral from $r$ to $\psi$, as the integral is easier to perform with the latter. The potential $V_{\psi}$ can be obtained through $R$ and (24).

Similarly, for the polar motion, the fundamental period is given by

$$\Lambda_{\theta}=4\int_{\theta_{tp}}^{\pi/2}\frac{d\theta}{\sqrt{\Theta}}=\int_{0}^{2\pi}\frac{d\chi}{\sqrt{V_{\chi}}},$$

(31)

where we have introduced the variable $\chi$ via $\cos^{2}\theta=\cos^{2}\theta_{tp}\cos^{2}\chi$, such that as $\chi$ varies from $0$ to $2\pi$, $\theta$ oscillates through its full range of motion, from $\theta_{tp}$ to $\pi-\theta_{tp}$ and back Drasco:2003ky . Thus, the angular frequencies of the radial and the polar motion with respect to $\lambda$ then become

$$\omega_{r}=\frac{2\pi}{\Lambda_{r}},\quad\omega_{\theta}=\frac{2\pi}{\Lambda_{\theta}}.$$

(32)

For the azimuthal motion, the equation (21) is the sum of a function of $r$ and a function of $\theta$, which allows us to define the frequencies of the coordinate $\phi$ with respect to $\lambda$ as

$$\omega_{\phi}=\left<\frac{d\phi}{d\lambda}\right>_{\lambda}=\left<\Phi_{r}\right>_{\lambda}+\left<\Phi_{\theta}\right>_{\lambda},$$

(33)

where

$$\left<\Phi_{r}\right>_{\lambda}=\frac{1}{\Lambda_{r}}\int\Phi_{r}d\lambda=\frac{1}{\Lambda_{r}}\int_{0}^{2\pi}\frac{\Phi_{r}}{\sqrt{V_{\psi}}}d\psi,$$

(34)

and

$$\left<\Phi_{\theta}\right>_{\lambda}=\frac{1}{\Lambda_{\theta}}\int_{0}^{2\pi}\frac{\Phi_{\theta}}{\sqrt{V_{\chi}}}d\chi.$$

(35)

Analogously, for the motion in $t$, the equation (20) is also the sum of a function of $r$ and a function of $\theta$, so we can define the frequencies of the coordinate $t$ with respect to $\lambda$ as

$$\omega_{t}=\frac{1}{\Lambda_{r}}\int_{0}^{2\pi}\frac{T_{r}}{\sqrt{V_{\psi}}}d\psi+\frac{1}{\Lambda_{\theta}}\int_{0}^{2\pi}\frac{T_{\theta}}{\sqrt{V_{\chi}}}d\chi.$$

(36)

Up to now, the above fundamental frequencies were written with respect to the Mino time $\lambda$, the frequencies with respect to the distant observer time, i.e. the coordinate time $t$, are obtained by Schmidt:2002qk ; Drasco:2003ky

$$\Omega_{r}=\frac{\omega_{r}}{\omega_{t}},\quad\Omega_{\theta}=\frac{\omega_{\theta}}{\omega_{t}},\quad\Omega_{\phi}=\frac{\omega_{\phi}}{\omega_{t}}.$$

(37)

Explicitly, the asymptotic form of these fundamental frequencies in the weak-field regime are given by

	$\displaystyle\Omega_{r}$	$\displaystyle=$	$\displaystyle\frac{\sqrt{1-qQ}}{M}\left(\frac{1-e^{2}}{p}\right)^{3/2}+\frac{3\sqrt{1-qQ}(qQ-4)}{4M}\left(\frac{1-e^{2}}{p}\right)^{5/2}$		(38)
		$\displaystyle+$	$\displaystyle\mathcal{O}(p^{-3}),$
	$\displaystyle\Omega_{\theta}$	$\displaystyle=$	$\displaystyle\frac{\sqrt{1-qQ}}{M}\left(\frac{1-e^{2}}{p}\right)^{3/2}+\frac{a(1-qQ)(e^{2}-1)^{3}\sin\theta_{tp}}{Mp^{3}}$		(39)
		$\displaystyle+$	$\displaystyle\frac{1}{M}\frac{\left.Q\left(q^{2}Q-3q+2Q\right)-3e^{2}\left(q^{2}Q^{2}-5qQ+4\right)\right)}{4(1-e^{2})\sqrt{1-qQ}}\left(\frac{1-e^{2}}{p}\right)^{5/2}$
		$\displaystyle+$	$\displaystyle\mathcal{O}(p^{-7/2}),$
	$\displaystyle\Omega_{\phi}$	$\displaystyle=$	$\displaystyle\frac{\sqrt{1-qQ}}{M}\left(\frac{1-e^{2}}{p}\right)^{3/2}$		(40)
		$\displaystyle+$	$\displaystyle\frac{a(2-qQ)(1-e^{2})^{3/2}+a(1-qQ)(e^{2}-1)^{3}\sin\theta_{tp}}{Mp^{3}}$
		$\displaystyle+$	$\displaystyle\frac{-Q(-3q+2Q+q^{2}Q)+3e^{2}(4-5qQ+q^{2}Q^{2})}{4M(1-e^{2})\sqrt{1-qQ}}\left(\frac{1-e^{2}}{p}\right)^{5/2}$
		$\displaystyle+$	$\displaystyle\mathcal{O}(p^{-7/2}).$

Therefore, we obtain the perihelion precession frequency and the orbital plane precession frequency at the leading order of the $1/p$ expansion

	$\displaystyle\dot{\tilde{\gamma}}$	$\displaystyle=$	$\displaystyle\frac{1}{M}\frac{6-6Qq-Q^{2}+q^{2}Q^{2}}{2(1-e^{2})\sqrt{1-qQ}}\left(\frac{1-e^{2}}{p}\right)^{5/2},$		(41)
	$\displaystyle\dot{\alpha}$	$\displaystyle=$	$\displaystyle\frac{a(2-qQ)(1-e^{2})^{3/2}}{Mp^{3}}.$		(42)

Comparing with the uncharged equations, the modification due to the charges in the system can be summarized as the electric force between the CO and the KN BH in the form $qQ$ and the deformation of the metric from the charge of the KN BH in the form $Q^{2}$. As we will see in the next subsection, the other contribution of the charges to the waveforms stems from the dipole electromagnetic radiation in the form $(Q-q)^{2}$.

Besides the fundamental frequencies, the other piece of the AK waveform is the change rates of the eccentricity $e$ and the radial orbital frequency $\nu$ with respect to the coordinate time, which are related to the energy flux and the angular momentum flux due to the gravitational radiation and electromagnetic radiation.

For the gravitational radiation, the standard quadrupole formulas of the energy flux and the angular momentum flux were already derived by Peters and Mathews Peters:1963ux ; Peters:1964zz

$$\left<\frac{dE}{dt}\right>=\frac{1}{5\mu}\left<\frac{d^{3}Q_{ij}}{dt^{3}}\frac{d^{3}Q^{ij}}{dt^{3}}-\frac{1}{3}\frac{d^{3}Q^{i}_{i}}{dt^{3}}\frac{d^{3}Q^{j}_{j}}{dt^{3}}\right>,$$

(43)

and

$$\left<\frac{dL_{i}}{dt}\right>=\frac{2}{5\mu M}\epsilon_{ijk}\left<\frac{d^{2}Q_{jm}}{dt^{2}}\frac{d^{3}Q^{km}}{dt^{3}}\right>,$$

(44)

where $\mu$ is the reduced mass $\mu=mM/(m+M)\simeq m$, $E$ and $L_{i}$ are the dimensionless energy and angular momentum appeared in the equations of motion of the charged CO. $Q_{ij}$ is the familiar quadrupole moment tensor of mass (also called the inertia tensor)

$$Q_{ij}=\mu x^{i}x^{j},$$

(45)

where $x^{i}$ is the relative position vector between the charged CO and the charged central BH, and in the weak-field regime one has $x^{i}=(r\cos\phi\sin\theta,r\sin\phi\sin\theta,r\cos\phi)$. In addition, the angle-brackets mean the average over one cyclic motion in $r$, which via (24), can be turned into the integral for $\psi$,

$$\left<\frac{dE}{dt}\right>=\frac{1}{T}\int_{0}^{2\pi}\frac{dE}{dt}\frac{d\psi}{\dot{\psi}},\quad T=\int_{0}^{2\pi}\frac{d\psi}{\dot{\psi}}.$$

(46)

In the following, without ambiguity we will just use $dE/dt$ and $dL_{z}/dt$ to denote the averaged ones. Then from above formulas and the equations of motion of the CO, we obtain the energy flux and angular momentum flux loss of the charged particle due to the gravitational radiation

	$\displaystyle\frac{dE}{dt}$	$\displaystyle=$	$\displaystyle-\frac{32\eta}{5Mp^{5}}(1-qQ)^{3}(1-e^{2})^{3/2}$		(47)
			$\displaystyle\times\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right),$
	$\displaystyle\frac{dL_{z}}{dt}$	$\displaystyle=$	$\displaystyle-\frac{32\eta\sin\theta_{tp}}{5Mp^{7/2}}(1-qQ)^{5/2}(1-e^{2})^{3/2}$		(48)
			$\displaystyle\times\left(1+\frac{7}{8}e^{2}\right),$

where the assumption that we have made is same to Barack:2003fp in which $\theta_{tp}$ is constant at leading order in the $1/p$ expansion. To simplify the expressions we have introduced the symmetric mass ratio,

$$\eta\equiv\frac{mM}{(m+M)^{2}}\simeq\frac{m}{M}.$$

(49)

In terms of the relation between $p$ and the radial orbital frequency at leading order

$$p=\frac{(1-e^{2})(1-qQ)^{1/3}}{(2\pi M\nu)^{2/3}},$$

(50)

the above two equations can be written as

	$\displaystyle\frac{dE}{dt}$	$\displaystyle=$	$\displaystyle-\frac{32\eta}{5M}(2\pi M\nu)^{10/3}(1-qQ)^{4/3}(1-e^{2})^{-7/2}$		(51)
			$\displaystyle\times\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right),$
	$\displaystyle\frac{dL_{z}}{dt}$	$\displaystyle=$	$\displaystyle-\frac{32\eta\sin\theta_{tp}}{5M}(2\pi M\nu)^{7/3}(1-qQ)^{4/3}(1-e^{2})^{-2}$		(52)
			$\displaystyle\times\left(1+\frac{7}{8}e^{2}\right).$

On the other hand, from the asymptotic form of the energy and angular momentum in the weak-field regime (28) and (29), we can obtain their change rates with respect to distant observer time

	$\displaystyle\frac{dE}{dt}$	$\displaystyle=$	$\displaystyle-\frac{(2\pi M)^{2/3}(1-qQ)^{2/3}}{3\nu^{1/3}}\frac{d\nu}{dt},$		(53)
	$\displaystyle\frac{dL_{z}}{dt}$	$\displaystyle=$	$\displaystyle-\frac{(1-qQ)^{2/3}\sin\theta_{tp}}{3(2\pi M\nu)^{1/3}\sqrt{1-e^{2}}\nu}$		(54)
			$\displaystyle\times\left(3e\nu\frac{de}{dt}+(1-e^{2})\frac{d\nu}{dt}\right).$

We should also take the orbit decay due to the electromagnetic radiation into account. The formulas of the energy and angular emission rates for the dipole electromagnetic radiation²²2Here we use the convention in Landau:1971 , where the Lagrangian of the electromagnetic field is $\mathcal{L}=-\frac{1}{16\pi}F_{\mu\nu}F^{\mu\nu}+A_{\mu}J^{\mu}$ and the equation of motion is $\nabla_{\mu}F^{\mu\nu}=-4\pi J^{\mu}$. are given by Landau:1971 and also Liu:2020cds ; Christiansen:2020pnv ; Liu:2020bag ; Liu:2022wtq

$$m\frac{dE_{EM}}{dt}=\frac{2}{3}\mu^{2}(Q-q)^{2}\left<\frac{d^{2}x^{i}}{dt^{2}}\frac{d^{2}x_{i}}{dt^{2}}\right>,$$

(55)

and

$$mM\frac{dL_{z}^{EM}}{dt}=\frac{2}{3}\mu^{2}(Q-q)^{2}\epsilon^{\jmath kl}\left<\frac{dx_{k}}{dt}\frac{d^{2}x_{l}}{dt^{2}}\right>,$$

(56)

where as before the angle-brackets denote the average over one period in $r$ motion, and $\jmath=z$. It should be noted that the above mentioned $E^{EM}$ and $L_{z}^{EM}$ are set to be dimensionless, their relation with the original ones are the same as (II.1). Then for a charged particle orbiting a KN BH, the energy flux and angular momentum flux loss due to the electromagnetic radiation at the leading order in $1/p$ expansion are given by

	$\displaystyle\frac{dE_{EM}}{dt}$	$\displaystyle=$	$\displaystyle-\frac{(2+e^{2})\eta^{2}}{3(1-e^{2})^{5/2}}(2\pi mM\nu)^{8/3}(Q-q)^{2}$		(57)
			$\displaystyle\times(1-qQ)^{2/3},$
	$\displaystyle\frac{dL_{z}^{EM}}{dt}$	$\displaystyle=$	$\displaystyle-\frac{2m\sin\theta_{tp}}{3M^{3}(1-e^{2})}(2\pi M\nu)^{5/3}(Q-q)^{2}$		(58)
			$\displaystyle\times(1-qQ)^{2/3}.$

Combining the equations (51), (52), (54), (54), (57) and (58), we have the evolution equations of orbital eccentricity $e$ and $\nu$ due to the gravitational and electromagnetic radiation

	$\displaystyle\frac{d\nu}{dt}$	$\displaystyle=$	$\displaystyle\frac{96}{10\pi}\frac{\eta}{M^{2}}(2\pi M\nu)^{11/3}(1-qQ)^{2/3}(1-e^{2})^{-7/2}$		(59)
			$\displaystyle\left(1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\right)$
			$\displaystyle+\frac{\eta(2+e^{2})(Q-q)^{2}(2\pi M\nu)^{3}(1-e^{2})^{-5/2}}{2\pi M^{2}},$
	$\displaystyle\frac{de}{dt}$	$\displaystyle=$	$\displaystyle-\frac{e\eta}{15M}(2\pi M\nu)^{8/3}(1-qQ)^{2/3}(1-e^{2})^{-5/2}(304+121e^{2})$		(60)
			$\displaystyle-\frac{e\eta}{M}(Q-q)^{2}(2\pi M\nu)^{2}(1-e^{2})^{-3/2}$

From the right hand sides of these equations, we find that the contribution from the dipole electromagnetic radiation is lower than that from the gravitational radiation by a factor of $(2\pi M\nu)^{2/3}$, which for a Keplerian orbit corresponds to $v^{2}$, where $v$ is the orbital velocity of the CO. This means that the correction due to the dipole electromagnetic radiation appears at $-1$ PN order in the waveforms and becomes prominent at the early stage of the inspiral of the CO where $v$ is small. This is verified by numerical relativity simulations of the coalescence of the charged binary BHs with comparable masses Bozzola:2020mjx , where it was found that the greatest difference between charged and uncharged BHs arises in the earlier inspiral. In this case, despite the fact that the AK model is not accurate enough to produce EMRI template waveforms in the strong field region, the behavior of the electric charges may well be captured by this model.

In terms of the radial frequency $\nu$ the other two equations (41) and (42) are expressed as

	$\displaystyle\frac{d\tilde{\gamma}}{dt}$	$\displaystyle=$	$\displaystyle\frac{(6-6Qq-Q^{2}+q^{2}Q^{2})}{(1-qQ)^{4/3}(1-e^{2})}\pi\nu(2\pi M\nu)^{2/3},$		(61)
	$\displaystyle\frac{d\alpha}{dt}$	$\displaystyle=$	$\displaystyle\frac{4a(2-qQ)\pi^{2}M\nu^{2}}{(1-e^{2})^{3/2}(1-qQ)}.$		(62)

From these four evolution equations we can see that the equations keep invariant under the operation $q\to-q$ and $Q\to-Q$, which means we cannot simultaneously determine the sign of $q$ and $Q$. Hereafter, without loss of generality, we will always assume that the charge of the MBH is positive, and let the sign of the charge of the orbiting particle to be free. For gravitational waveform generated by two charged compact objects moving on a Keplerian orbit in a plane, the charges appear only in the form $(Q-q)^{2}$ and $qQ$, which hinders the unique determination of the two charges. However, the relativistic effects considered in the AK waveform breaks the degeneracy between these two terms. As a consequence, once the sign of $Q$ is provided, we can not only identify the value of $Q$ but also uniquely discern both the magnitude and the sign of $q$. Even the dipole electromagnetic radiation disappears when $Q=q$, this conclusion still works.

Up to now, we have obtained the leading order evolution equations for the relevant orbital parameters when both the CO and central BH are charged. In analogy with the construction of EMRI waveforms in alternative theories of gravity, e.g. Zi:2021pdp , we combine these leading order corrected equations with those higher-order PN equations in the original AK model. Then the complete orbital evolution equations are given by

	$\displaystyle\dot{\Phi}~{}=~{}$	$\displaystyle 2\pi\nu,$
	$\displaystyle\dot{\nu}~{}=~{}$	$\displaystyle\frac{48}{5\pi}\frac{\eta}{M^{2}}X^{11/3}Y^{-9/2}\Big{\{}(1-qQ)^{2/3}Y\Big{(}1+\frac{73}{24}e^{2}+\frac{37}{96}e^{4}\Big{)}$
		$\displaystyle+X^{2/3}\Big{[}\frac{1273}{336}-\frac{2561}{336}e^{2}-\frac{3885}{128}e^{4}-\frac{13147}{5376}e^{6}\Big{]}$
		$\displaystyle-Xa\cos\lambda Y^{-1/2}\Big{[}\frac{73}{12}+\frac{1211}{24}e^{2}+\frac{3143}{96}e^{4}+\frac{65}{64}e^{6}\Big{]}\Big{\}}$
		$\displaystyle+\frac{\eta(Q-q)^{2}}{2\pi M^{2}}X^{3}(2+e^{2})Y^{-5/2},$
	$\displaystyle\dot{e}~{}=~{}$	$\displaystyle-\frac{e\eta}{15M}Y^{-7/2}X^{8/3}\Big{[}((1-qQ)^{2/3}+12X^{2/3})(304+121e^{2})Y$
		$\displaystyle-\frac{1}{56}X^{2/3}(133640+108984e^{2}+25211e^{4})\Big{]}$
		$\displaystyle+e\frac{\eta}{M}a\cos\lambda X^{11/3}Y^{-4}\left(\frac{1364}{5}+\frac{5032}{15}e^{2}+\frac{263}{10}e^{4}\right)$
		$\displaystyle-e\frac{\eta}{M}(Q-q)^{2}X^{2}Y^{-3/2},$
	$\displaystyle\dot{\tilde{\gamma}}~{}=~{}$	$\displaystyle\pi\nu X^{2/3}Y^{-1}\Big{[}(6-6qQ-Q^{2}+q^{2}Q^{2})(1-qQ)^{-4/3}$
		$\displaystyle+\frac{3}{2}X^{2/3}Y^{-1}(26-15e^{2})\Big{]}-12\pi\nu a\cos\lambda XY^{-3/2},$
	$\displaystyle\dot{\alpha}~{}=~{}$	$\displaystyle 4aM(\pi\nu)^{2}Y^{-3/2}(2-qQ)/(1-qQ)$		(63)

where dot denotes the derivative with respect to time and to avoid redundant expression, we have defined $Y=1-e^{2}$, $X=2\pi M\nu$. The equation for $\dot{\nu}$ and $\dot{e}$ are given accurately through $3.5$ PN order, the equations for $\dot{\tilde{\gamma}}$ and $\dot{\alpha}$ are accurate through $2$ PN order.

We would like to give a brief review of how do above evolving orbital parameters enter the AK EMRI waveforms. The general GW strain field at the detector is written as

$$h_{ij}(t)=A^{+}(t)H_{ij}^{+}(t)+A^{\times}(t)H^{\times}_{ij}(t),$$

(64)

where $H_{ij}^{+}$ and $H_{ij}^{\times}$ are the two polarization basis tensors constructed with the unit vector pointing from the detector to the source $\hat{n}$ and the unit vector along the CO’s orbital angular momentum $\hat{L}$,

$$H_{ij}^{+}(t)=\hat{p}_{i}\hat{p}_{j}-\hat{q}_{i}\hat{q}_{j},\quad H_{ij}^{\times}(t)=\hat{p}_{i}\hat{q}_{j}+\hat{q}_{i}\hat{p}_{j},$$

(65)

with

$$\hat{p}=\frac{\hat{n}\times\hat{L}}{|\hat{n}\times\hat{L}|},\quad\hat{q}=\hat{p}\times\hat{n},$$

(66)

and $A^{+}$ and $A^{\times}$ are the amplitudes of the two polarizations. The amplitudes of the two polarisations can be further written in terms of the Peters-Mathews harmonic decomposition as

	$\displaystyle A^{+}$	$\displaystyle\equiv\sum_{n}A_{n}^{+}$
		$\displaystyle=\sum_{n}-\Big{[}1+(\hat{L}\cdot\hat{n})^{2}\Big{]}\Big{[}a_{n}\cos 2\gamma-b_{n}\sin 2\gamma\Big{]}$
		$\displaystyle+c_{n}\Big{[}1-(\hat{L}\cdot\hat{n})^{2}\Big{]},$		(67)
	$\displaystyle A^{\times}$	$\displaystyle\equiv\sum_{n}A_{n}^{\times}$
		$\displaystyle=\sum_{n}2(\hat{L}\cdot\hat{n})\Big{[}b_{n}\cos 2\gamma+a_{n}\sin 2\gamma\Big{]},$		(68)

where $(a_{n},b_{n},c_{n})$ come from decomposition of the second time derivative of the inertia tensor $Q^{ij}$ into $n-$harmonics of the radial orbital frequency $\nu$ and are functions of $\nu$ and $e$ Peters:1963ux . Moreover, $\gamma$ is an azimuthal angle measuring the direction of pericentre with respect to the orthogonal projection of $\hat{n}$ onto the orbital plane, which further depends on $\tilde{\gamma}$ and $\alpha$.

Since the equilateral triangle detectors such as TQ can be used to construct two independent Michelson interferometers, the signal responded by such two interferometers can be written as:

$$h_{I,II}=\frac{\sqrt{3}}{2}\left(F^{+}_{I,II}h^{+}+F^{\times}_{I,II}h^{\times}\right),$$

(69)

where the antenna pattern function $F^{+,\times}_{I,II}$ of the detector depend on the orbits of satellites Cutler:1994ys . For TQ, the detailed information of the respond function for EMRI signal can be found in Fan:2020zhy .

At the final stage of EMRI, when the CO passed the boundary of stable orbits, it will plunge into the MBH directly in a short time. So we need to introduce a cutoff frequency to the waveform. When the CO is moving in the equatorial plane of the central BH, the cutoff is usually taken to be the last stable orbit (LSO). Here for simplicity we take the innermost stable circular orbit (ISCO) as the cutoff in our waveform. The radius of the ISCO results from the equations $R(r)=R^{\prime}(r)=R^{\prime\prime}(r)=0$. Although these equations can be solved analytically with Mathematica, the expression of the solution is very lengthy so we shall not show it here. Instead, the effects of the MBH charge $Q$ and CO charge $q$ on the ISCO radius can be demonstrated clearly in a graphical manner. Note that in this work we only consider the prograde orbits of the CO, since most of the detected events have prograde orbits Fan:2020zhy .

As shown in Fig. 1, $r_{\rm ISCO}$ changes gradually with the MBH charge $Q$ (left panel) and CO charge $q$ (right panel). This indicates that the CO gets a chance to orbit more circles in the vicinity of the KN BH, and the charged EMRI system radiates higher frequency GW signal than the neutral system. The effect of the CO charge on $r_{\rm ISCO}$ is more complicated. In this case, $r_{\rm ISCO}$ no longer depends on $q$ monotonically. Nevertheless, the turning point at which the monotonicity of $q$ changes increase with $Q$.