Excitation factors and exceptional points of Kerr–Newman quasinormal modes

Introduction

According to the black hole no-hair theorem, stationary, asymptotically flat black holes in four-dimensional Einstein–Maxwell theory belong to the Kerr–Newman (KN) family, characterized by their mass $M$, angular momentum $J = Ma$, and electric charge $Q$ [1, 2, 3, 4, 5]. As an exact solution of the Einstein–Maxwell equations, the KN spacetime generalizes both the Kerr and Reissner–Nordström black holes, providing a framework for investigating the effects of rotation and electric charge on black hole geometry, thermodynamics, and dynamics. From an observational perspective, gravitational-wave observations and black hole imaging offer opportunities to test strong-field gravity and constrain possible deviations from the Kerr geometry [6, 7]. Although astrophysical black holes are expected to carry negligible net electric charge, the KN solution remains an important theoretical reference for studying the influence of electric charge on stationary spacetimes. Linear perturbation theory describes the dynamical response of black holes to small disturbances and offers a powerful probe of their geometries [8, 9]. Their resonant oscillations are characterized by quasinormal modes (QNMs), whose complex frequencies determine the oscillation rates and damping times of black hole ringdown. For given mode labels, these frequencies depend on the black hole background, while the amplitudes of the excited modes also depend on external sources or initial conditions [8, 10, 11, 12, 13].

Within the Newman–Penrose formalism [14], gravitational and electromagnetic perturbations of Kerr black holes satisfy decoupled, separable Teukolsky equations [15, 16]. For KN black holes, however, the background electromagnetic field couples gravitational and electromagnetic perturbations, making their separation considerably more challenging. Previous studies employed the Dudley–Finley approximation [17, 18], slow-rotation and small-charge expansions [19, 20, 21], analytical approximations and WKB methods [22, 23, 24, 25], and direct numerical solutions of the coupled equations [26, 27, 28]. These numerical studies of the coupled perturbation equations have provided evidence for linear mode stability and revealed eigenvalue repulsion near extremality.

Recently, Hintz made an important breakthrough by establishing the exact separation of the radiative gravitoelectromagnetic perturbation equations on KN backgrounds [29]. This resolves a longstanding difficulty in KN perturbation theory by separating the coupled gravitational and electromagnetic perturbations into radial and angular equations linked by a common separation constant $\lambda$. Using these equations, Destounis, Cardoso, and Hintz computed KN QNM frequencies in excellent agreement with previous independent numerical results [30]. This framework has also been applied to the tidal response of KN black holes [31] and extended to Kerr–Newman–(anti-)de Sitter spacetimes [32, 33].

While QNM frequencies determine the oscillation rates and damping times, they do not characterize the strengths of the corresponding perturbations. Excitation factors (EFs), which are related to the residues of QNM poles, provide complementary information about the source-independent response of black holes. Another important feature of the QNM spectrum is the existence of exceptional points (EPs), where both eigenvalues and their associated eigenfunctions coalesce in such non-Hermitian eigenvalue problems [34, 35, 36]. Near such degeneracies, QNM frequencies can become highly sensitive to variations in black hole parameters, and the corresponding EFs may exhibit strong enhancements [37]. During the preparation of this manuscript, we became aware of two related studies of KN EPs. One identified EPs and eigenstate hysteresis in the coupled gravitoelectromagnetic QNM spectrum of KN black holes [38]. The other investigated spectral instability near KN EPs through pseudospectra and condition numbers for neutral massless scalar fields and coupled gravitoelectromagnetic perturbations [39]. Studying EFs and EPs together therefore provides further insight into the spectral structure and resonant response of black holes.

In this work, we investigate the EFs and EPs of KN QNMs using the exact separated equations for coupled gravitoelectromagnetic perturbations. We first examine how rotation and electric charge affect the EFs of selected gravitational-led and electromagnetic-led modes, comparing their magnitudes and complex-plane trajectories. We then identify 42 EPs in the gravitational-led and electromagnetic-led sectors, revealing ordered overtone ladders and additional off-ladder coalescences. These results provide further insight into the excitation properties and non-Hermitian spectral structure of coupled gravitoelectromagnetic perturbations of KN black holes.

The paper is organized as follows. Section II introduces the separated perturbation equations, the corresponding boundary conditions, and the methods for computing QNMs, EFs, and EPs. Section III presents the numerical results, including the dependence of EFs on black hole spin and charge and the behavior of QNMs near EPs. Finally, Section IV summarizes our main findings and discusses possible future directions.

Method

Gravitoelectromagnetic perturbations of Kerr–Newman black holes

We consider an asymptotically flat Kerr–Newman black hole in Einstein–Maxwell theory. In geometrized units $G = c = 1$ and with metric signature $(+,-,-,-)$, the background metric in Boyer–Lindquist coordinates reads

where

Here $M$, $a=J/M$, and $Q>0$ denote the black hole mass, spin parameter, and electric charge, respectively. We restrict to subextremal configurations, $a^2+Q^2<M^2$, with the outer and inner horizons located at

For convenience, we use the dimensionless parameters $\hat{a}=a/M$ and $\hat{Q}=Q/M$. Hintz obtained an exact separation of the radiative gravitoelectromagnetic perturbation equations on the KN background introduced above [29]. QNM frequencies computed from these equations were found to agree with existing numerical results for the coupled system [30]. We adopt the separated variables and conventions of these works below.

Angular eigenvalue problem. We take the time and azimuthal dependence of the perturbations to be $e^{-i\omega t+im\phi}$, with $m\in\mathbb{Z}$. For each spin label $s=\pm1$, the angular amplitudes $\chi_1(\theta)$ and $\chi_2(\theta)$ satisfy

where $\lambda$ is the separation constant and $s=\pm1$ distinguishes the positive- and negative-spin descriptions. The functions $W(\theta)$ and $T(\theta)$ are given by

We allow both $\omega$ and $\lambda$ to be complex. Regularity at the poles $\theta=0$ and $\theta=\pi$ requires the angular functions to satisfy

where

Here $s=1,2$ labels the two angular components. For a given frequency $\omega$, the angular equation subject to these boundary conditions determines the allowed eigenvalues $\lambda$.

At $a=0$, the angular eigenvalues are known explicitly. For $\ell\geq2$ and $|m|\leq\ell$, they are

For $\ell=1$, when only the $b=-1$ angular eigenvalue can be allowed, $\lambda_{1,-1}^{(0)}=-3M/(2Q)$, with $m=-1,0,1$. As $a$ is varied from zero, we follow each angular eigenvalue with labels $(\ell,b)$. The modes with $b=+1$ and $b=-1$ reduce to the gravitational and electromagnetic perturbations of Kerr black hole, respectively, as black hole charge $Q\to0$. We therefore refer to them as gravitational-led and electromagnetic-led modes. The equations written above contain terms proportional to $1/Q$, so this limit cannot be obtained by directly setting $Q=0$. The required variable rescalings are given in [30] and summarized in Appendix B.

Radial equation. To write down the radial equations, we first define

Following Ref. [30], we introduce the radial functions

where $z_{1,2}^{(\varsigma)}$ are the radial amplitudes in the separated first-order equations and

For each $\varsigma=\pm1$, the radial function $w=w^{(\varsigma)}$ satisfies the following second-order equation wherever $h(r)\ne0$

where primes denote derivatives with respect to $r$. Expanding the derivatives and multiplying by $h(r)$, we obtain the equivalent equation

with

This form contains no explicit $1/h$ terms and is convenient for numerical calculations.

Boundary conditions and quasinormal modes

We now impose boundary conditions on the separated perturbation equations. The radial solutions are required to be purely ingoing at the event horizon and purely outgoing at spatial infinity. These conditions determine the complex QNM frequencies. We denote two solutions of the radial equation (14) for $w$ by $R_{\ell m}^{\mathrm{in}}(r)$ and $R_{\ell m}^{\mathrm{up}}(r)$, with the usual in/up notation [40]. The first is characterized by purely ingoing behavior at the event horizon,

while the second is characterized by purely outgoing behavior at infinity. For generic frequencies, the two solutions are linearly independent and have the following asymptotic forms

The functions $g_H^{\mathrm{in/out}}$ and $f_\infty^{\mathrm{in/out}}$ describe the ingoing and outgoing solutions near the horizon and at infinity, respectively. The coefficients $A_{\ell m}^{(\cdot)}$ and $C_{\ell m}^{(\cdot)}$ depend on $\omega$ and $\lambda$.

Near the event horizon, we introduce $x=1-r_+/r$. For $h(r_+)\ne0$ and $\rho_H\notin\mathbb{Z}$, the two independent solutions have the following forms

where

Here $\Omega_H=a/(r_+^2+a^2)$ and $\kappa_H=(r_+-r_-)/[2(r_+^2+a^2)]$ are the angular velocity and surface gravity of the event horizon, respectively. The transformation defining $w$ removes the ingoing phase and the associated powers of $\Delta$, so the ingoing solution corresponds to the regular branch $g_H^{\mathrm{in}}$. At spatial infinity, for $\omega\ne0$, the ingoing and outgoing solutions behave as

The coefficients $a_1^{\mathrm{in/out}}$ and the higher-order terms are determined by substituting these expansions into the radial equation.

The QNM frequencies $\omega_{\ell mn}$ are obtained by requiring the incoming wave amplitude of $R_{\ell m}^{\mathrm{in}}(r)$ to vanish at infinity

Here $\lambda(\omega)$ is the angular eigenvalue satisfying the regularity conditions on the chosen $(\ell,m,b)$ branch, and $n$ denotes the radial overtone number. The labels $b$ and $\varsigma$ are suppressed for simplicity.

The QNMs $\omega_{\ell mn}$ and separation constant $\lambda$ are determined simultaneously using angular Galerkin projection and radial Chebyshev collocation (see, e.g. [30]). Similar techniques have been employed in studies of Kerr angular eigenvalue problems [41, 42], and coupled gravitoelectromagnetic perturbations of KN black holes [26, 28, 43]. The numerical procedure is described in Appendices A–D.

Excitation factors and exceptional points in QNM spectrum

In addition to the QNM frequencies, the EFs characterize the source-independent strength of the corresponding perturbations. For the separated radial variable $w$, they are determined by the outgoing amplitude $A_{\mathrm{ref}}$ and the frequency derivative of the incoming amplitude $A_{\mathrm{inc}}$, with the boundary conditions specified in Sec. II B. For a simple QNM with $\left.\mathrm{d}A_{\ell mn}^{(\mathrm{inc})}/\mathrm{d}\omega\right|_{\omega_q} \ne0$, the excitation factor is defined in terms of the asymptotic amplitudes in (16):

The outgoing amplitude in the numerator is evaluated at $\omega= \omega_{\ell mn}$, and the frequency derivative includes the dependence of the angular eigenvalue $\lambda$ on $\omega$, with $a$ and $Q$ held fixed. The amplitude of a QNM in a physical ringdown signal depends on the external source or initial conditions, and its calculation requires reconstructing the corresponding radiation fields from the separated variables. In contrast, the EF $\mathcal{B}_{\ell mn}^{(w)}$ is determined by the homogeneous perturbation equations and their boundary conditions, independently of the source.

The ingoing and outgoing boundary conditions at the event horizon and infinity make the QNM eigenvalue problem non-self-adjoint, allowing for the existence of EPs. At an EP, two QNM frequencies and their associated eigenfunctions coalesce. For Kerr–Newman black holes, the rotation parameter $a$ and electric charge $Q$ provide two independent parameters for exploring such degeneracies.

For the fixed black hole parameters $(a,Q)$ and a chosen angular eigenvalue branch $\lambda(\omega;a,Q)$, the incoming-wave amplitude is defined as

The QNM frequencies correspond to the zeros of $\mathcal{A}$. A simple QNM at $\omega= \omega_q$ satisfies $\mathcal{A}(\omega_q) = 0$ and $D_\omega\mathcal{A}(\omega_q) \ne0$, where $D_\omega$ denotes differentiation along the angular eigenvalue branch. As the black hole parameters vary, two QNM frequencies may coalesce, producing a double zero of $\mathcal{A}$. A candidate second-order EP therefore satisfies

evaluated at $(\omega_{\mathrm{EP}},a_{\mathrm{EP}},Q_{\mathrm{EP}})$. The frequency derivative is taken along the angular eigenvalue branch $\lambda(\omega;a,Q)$, with $a$ and $Q$ held fixed:

The connection between EPs and EFs follows from (22), which contains the inverse frequency derivative $(D_\omega\mathcal{A})^{-1}$ for a simple QNM. As two QNM frequencies approach an EP, this derivative in Eq. (25) tends to zero. Similar discussions of EPs near avoided crossings have been observed in Kerr QNMs [37]. The accuracy of the computed QNM and EP frequencies is assessed by comparison with previous results in Appendix E.

Results

Spin and charge dependence of the excitation factors

We investigate how rotation and electric charge affect the EFs of KN QNMs. Figure 1 and Figure 2 present the gravitational-led ($b=+1$) and electromagnetic-led ($b=-1$) sectors, respectively. In the gravitational-led sector, we consider $(\ell,m,n)=(2,2,0)$, $(2,2,1)$, $(3,3,0)$, and $(4,4,0)$, whereas the electromagnetic-led sector includes $(1,1,0)$, $(1,1,1)$, $(2,2,0)$, and $(2,2,1)$. The curves are shown for four fixed charges, $\hat{Q}=0.05$, $0.30$, $0.85$, and $0.90$, with the black hole spin $\hat{a}$ varied from zero toward the extremal bound $\hat{a}=\sqrt{1-\hat{Q}^{2}}$, while remaining in the subextremal region $\hat{a}^{2}+\hat{Q}^{2}<1$.

a. Gravitational-led sector. As shown in Figure 1, the EF magnitudes of all four selected gravitational-led modes decrease overall as the black hole spin increases at fixed charge. The strength of this suppression, however, depends strongly on the mode. Table I quantifies this difference for the fundamental and first-overtone modes in the $\ell=m=2$ sector at $\hat{Q}=0.05$. As $\hat{a}$ increases from $0$ to $0.98$, $\left|\mathcal{B}_{220}^{(w)}\right|$ decreases by a factor of approximately $36$, whereas $\left|\mathcal{B}_{221}^{(w)}\right|$ decreases by a factor of only about $4.5$.

The stronger suppression of the fundamental mode reverses the relative ordering of the two EFs: $\left|\mathcal{B}_{220}^{(w)}\right|$ is larger at $\hat{a}=0$, while $\left|\mathcal{B}_{221}^{(w)}\right|$ is larger at $\hat{a}=0.98$. Thus, black hole rotation can reorganize the relative excitation strengths of different overtones even within the same angular sector. The $(3,3,0)$ and $(4,4,0)$ modes show the same overall decrease with spin, although their rates of suppression are different. Figure 1 also reveals a nontrivial charge dependence. Over the spin interval covered by all four curves, the EF magnitudes generally increase as $\hat{Q}$ rises from $0.05$ to $0.30$. At higher charge, however, the $\hat{Q}=0.85$ and $0.90$ curves may cross, reversing their relative ordering as the spin varies.

b. Electromagnetic-led sector. As shown in Figure 2, the spin dependence of the electromagnetic-led modes shares some features with, but also differs from, that of the gravitational-led sector. At fixed charge, the EF magnitudes of the $(1,1,0)$, $(2,2,0)$, and $(2,2,1)$ modes decrease overall as the black hole spin increases. The $(1,1,1)$ mode, however, exhibits a nonmonotonic spin dependence: its EF magnitude initially increases, reaches a maximum at relatively large spin, and then decreases as the corresponding extremal bound is approached. The charge dependence is also mode dependent. Over the range of black hole spin considered for all four charges, the EF magnitudes of the $(1,1,0)$, $(2,2,0)$, and $(2,2,1)$ modes generally decrease as the charge increases. The $(1,1,1)$ mode again behaves differently: the $\hat{Q}=0.85$ and $0.90$ curves cross, so their relative ordering reverses as the spin varies.

$(\ell,m,n)$	$\hat{a}=0$	$\hat{a}=0.98$	$S$
$(2,2,0)$	$4.69\times10^{-3}$	$1.29\times10^{-4}$	$36.4$
$(2,2,1)$	$3.57\times10^{-3}$	$7.92\times10^{-4}$	$4.51$

Figure 3 compares the fundamental mode ($n=0$) and the first overtone ($n=1$) in the gravitational-led $\ell=m=2$ and electromagnetic-led $\ell=m=1$ sectors. The left panels show the ratio $|\mathcal{B}^{(w)}_{1}/\mathcal{B}^{(w)}_{0}|$ as a function of $\hat{a}$ at fixed charges. The right panels display the corresponding trajectories of $\mathcal{B}^{(w)}_{0}$ and $\mathcal{B}^{(w)}_{1}$ in the complex plane at $\hat{Q}=0.30$, with the arrows indicating increasing spin.

spin. For the low-charge cases $\hat{Q}=0.05$ and $0.30$, the ratio increases with spin in both perturbation sectors and eventually exceeds 1. Thus, at sufficiently rapid rotation, the EF magnitude of the first overtone becomes larger than that of the fundamental mode. The gravitational-led sector shows a similar increase at the higher charges considered here. In the electromagnetic-led sector, however, the $\hat{Q}=0.85$ and $0.90$ curves turn over and decrease sharply as their respective extremal bounds are approached, so that the fundamental-mode EF again becomes larger than the first-overtone EF. At $\hat{Q}=0.30$, the complex-plane trajectories in Figs. 3(b) and 3(d) provide phase information that is absent from the magnitudes alone. In the gravitational-led sector, the fundamental mode and first overtone trace pronounced loops as the black hole spin increases. The electromagnetic-led trajectories exhibit less pronounced winding, but the two modes still follow distinct paths.

Exceptional points of Kerr–Newman gravitoelectromagnetic perturbations

EPs provide a direct window into the non-Hermitian spectral structure of QNM problems. In the Kerr QNM spectrum, varying the single dimensionless spin parameter $a/M$ can lead to avoided crossings and eigenvalue repulsion between QNM branches [37, 44]. A generic EP, however, typically requires more than one independent control parameter. The KN family naturally supplies a two-dimensional parameter space through $\hat{a}$ and $\hat{Q}$, thereby offering a natural platform for investigating the exceptional-point structure of its QNM spectrum. Using the conditions introduced in Sec. II C, we identify 42 EPs in the coupled gravitoelectromagnetic QNM spectrum, as summarized in Table II. Twenty-three belong to the gravitational-led sector and nineteen to the electromagnetic-led sector. The gravitational-led EPs occur in the $(\ell, m) = (2, 2)$, $(3, 2)$, $(3, 3)$, and $(4, 3)$ modes, whereas the electromagnetic-led EPs are found in the $(1, 1)$, $(2, 1)$, and $(2, 2)$ modes. Additional EPs may exist for other mode pairs or in regions of parameter space not covered by the present search. We denote each EP by

EP label	$\hat{a}_{\mathrm{EP}}$	$\hat{Q}_{\mathrm{EP}}$	$M\omega_{\mathrm{EP}}$	$d_{\mathrm{EP}}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{2,2;\{0,1\}}$	0.304116510	0.952362666	$0.564398227-0.028304554\,\mathrm{i}$	$-5.185\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G},(0)}_{2,2;\{1,2\}}$	0.309806527	0.950727198	$0.570215504-0.025135809\,\mathrm{i}$	$-1.377\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G},(1)}_{2,2;\{1,2\}}$	0.209344101	0.971405377	$0.449330846-0.191881258\,\mathrm{i}$	$-1.255\times10^{-2}$
$\bullet\mathrm{EP}^{\mathrm{G},(2)}_{2,2;\{1,2\}}$	0.220095829	0.972984994	$0.454477546-0.185196350\,\mathrm{i}$	$-4.858\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{G},(0)}_{2,2;\{2,3\}}$	0.312388626	0.949922517	$0.572850039-0.023703417\,\mathrm{i}$	$-6.056\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G},(1)}_{2,2;\{2,3\}}$	0.225990537	0.972819474	$0.457397203-0.181809301\,\mathrm{i}$	$-2.551\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{2,2;\{3,4\}}$	0.313876728	0.949446103	$0.574366343-0.022870044\,\mathrm{i}$	$-3.350\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{2,2;\{4,5\}}$	0.314850778	0.949130064	$0.575357657-0.022319723\,\mathrm{i}$	$-2.111\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G},(0)}_{2,2;\{5,6\}}$	0.315541023	0.948904312	$0.576059338-0.021926955\,\mathrm{i}$	$-1.447\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G},(1)}_{2,2;\{5,6\}}$	0.749057505	0.532397483	$0.555905299-0.791341725\,\mathrm{i}$	$-1.555\times10^{-1}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,2;\{1,2\}}$	0.743448001	0.667234185	$0.972691674-0.081879676\,\mathrm{i}$	$-2.084\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,2;\{2,3\}}$	0.761028839	0.647930594	$0.975212331-0.076436630\,\mathrm{i}$	$-1.021\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,2;\{3,4\}}$	0.773694655	0.633079754	$0.977194724-0.073501994\,\mathrm{i}$	$-6.066\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,2;\{4,5\}}$	0.783206434	0.621438856	$0.978681255-0.071596994\,\mathrm{i}$	$-4.014\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,2;\{5,6\}}$	0.790616814	0.612078464	$0.979817398-0.070240484\,\mathrm{i}$	$-2.850\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,3;\{0,1\}}$	0.377451425	0.925982662	$0.994079166-0.012592716\,\mathrm{i}$	$-8.653\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,3;\{1,2\}}$	0.380844649	0.924626143	$0.999684309-0.010783874\,\mathrm{i}$	$-2.385\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G},(0)}_{3,3;\{2,3\}}$	0.382475757	0.923959822	$1.002389863-0.009969328\,\mathrm{i}$	$-1.054\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{G},(1)}_{3,3;\{2,3\}}$	0.299535626	0.952983464	$0.849371999-0.135671548\,\mathrm{i}$	$-2.101\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{3,3;\{3,4\}}$	0.383445163	0.923560493	$1.004002063-0.009496059\,\mathrm{i}$	$-5.822\times10^{-6}$
$\bullet\mathrm{EP}^{\mathrm{G}}_{4,3;\{1,2\}}$	0.667095777	0.744442401	$1.394804674-0.058219191\,\mathrm{i}$	$-7.887\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G},(0)}_{4,3;\{2,3\}}$	0.674811449	0.737727892	$1.398473807-0.053856357\,\mathrm{i}$	$-3.871\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{G},(1)}_{4,3;\{2,3\}}$	0.504647151	0.860790519	$1.244825997-0.196167093\,\mathrm{i}$	$-4.371\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{1,1;\{0,1\}}$	0.687265135	0.726012494	$0.472653256-0.023562859\,\mathrm{i}$	$-5.725\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{1,1;\{1,2\}}$	0.706931140	0.707171968	$0.475069610-0.020331190\,\mathrm{i}$	$-1.562\times10^{-4}$
$\bullet\mathrm{EP}^{\mathrm{E},(0)}_{1,1;\{2,3\}}$	0.716368444	0.697672665	$0.476161256-0.018872304\,\mathrm{i}$	$-6.911\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E},(1)}_{1,1;\{2,3\}}$	0.424335534	0.901822747	$0.394681893-0.198430498\,\mathrm{i}$	$-6.655\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{1,1;\{3,4\}}$	0.721974494	0.691892016	$0.476790542-0.018025020\,\mathrm{i}$	$-3.827\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{1,1;\{4,5\}}$	0.725714231	0.687978740	$0.477202514-0.017466386\,\mathrm{i}$	$-2.411\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{2,1;\{2,3\}}$	0.845724565	0.495798164	$0.558020601-0.321553435\,\mathrm{i}$	$-3.893\times10^{-2}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{2,1;\{3,4\}}$	0.893217748	0.425148005	$0.547326088-0.310235376\,\mathrm{i}$	$-2.141\times10^{-2}$
$\bullet\mathrm{EP}^{\mathrm{E},(0)}_{2,1;\{4,5\}}$	0.928040496	0.353599541	$0.540220533-0.303287275\,\mathrm{i}$	$-1.371\times10^{-2}$
$\bullet\mathrm{EP}^{\mathrm{E},(1)}_{2,1;\{4,5\}}$	0.594581447	0.783258533	$0.529534752-0.547071342\,\mathrm{i}$	$-3.298\times10^{-2}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{2,1;\{5,6\}}$	0.954235797	0.282599155	$0.534993813-0.292849976\,\mathrm{i}$	$-9.572\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{2,2;\{0,1\}}$	0.653891728	0.756531816	$0.918373632-0.010310056\,\mathrm{i}$	$-8.522\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E},(0)}_{2,2;\{1,2\}}$	0.662380552	0.749151818	$0.922230468-0.008685365\,\mathrm{i}$	$-2.356\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E},(1)}_{2,2;\{1,2\}}$	0.461271757	0.883991733	$0.793226261-0.144859314\,\mathrm{i}$	$-5.787\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E},(0)}_{2,2;\{2,3\}}$	0.666543573	0.745459097	$0.924092978-0.007959284\,\mathrm{i}$	$-1.040\times10^{-5}$
$\bullet\mathrm{EP}^{\mathrm{E},(1)}_{2,2;\{2,3\}}$	0.471332566	0.880489277	$0.795394060-0.137484232\,\mathrm{i}$	$-2.584\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E},(0)}_{2,2;\{3,4\}}$	0.669046534	0.743216659	$0.925204529-0.007539172\,\mathrm{i}$	$-5.733\times10^{-6}$
$\bullet\mathrm{EP}^{\mathrm{E},(1)}_{2,2;\{3,4\}}$	0.477535991	0.877784636	$0.797033638-0.133556710\,\mathrm{i}$	$-1.454\times10^{-3}$
$\bullet\mathrm{EP}^{\mathrm{E}}_{2,2;\{4,5\}}$	0.670728499	0.741700537	$0.925948100-0.007262762\,\mathrm{i}$	$-3.594\times10^{-6}$

where $X = G$ or $E$ denotes the gravitational-led or electromagnetic-led sector, respectively, $(\ell, m)$ gives the angular mode labels, and $(n_{1}, n_{2})$ identifies the two overtone branches whose QNM frequencies and eigenfunctions coalesce. If multiple EPs occur for the same sector and mode pair, the additional index $j$ distinguishes them. Figure 4 shows the distribution of the 42 EPs listed in Table II. Panel (a) gives their locations in the $(\hat{a}, \hat{Q})$ plane, where the solid black curve marks the extremal boundary $\hat{a}^{2} + \hat{Q}^{2} = 1$. Panel (b) displays the same EPs in the coordinates $(\vartheta, d_{\mathrm{EP}})$, with $\vartheta= \tan^{-1}(\hat{a}/\hat{Q})$ and $d_{\mathrm{EP}} = \hat{a}_{\mathrm{EP}}^{2} + \hat{Q}_{\mathrm{EP}}^{2} - 1$. Two features of this distribution are particularly noteworthy. First, both the gravitational-led and electromagnetic-led EPs exhibit ordered overtone-ladder structures, in which successive neighboring branches $\{n, n + 1\}$ coalesce. These ladders appear in several angular sectors, indicating a systematic organization of the EP spectrum across the two perturbation branches. For example, in the electromagnetic-led $\ell= m = 1$ sector, the overtone ladder takes the form

where the arrows indicate the ordering of the EPs with increasing overtone number. Along these ladders, increasing the overtone number is accompanied by increasing $\hat{a}_{\mathrm{EP}}$, decreasing $\hat{Q}_{\mathrm{EP}}$, and $d_{\mathrm{EP}}$ approaching zero. The higher-overtone members are therefore progressively closer to the extremal KN boundary.

A second notable feature is that the base labels $(X,\ell,m,n_{1},n_{2})$ do not always specify a unique EP. For a fixed sector, angular mode, and overtone pair, the two-dimensional KN parameter space may contain several distinct coalescences. When this occurs, we assign $j = 0$ to the EP belonging to the ordered overtone ladder and use $j = 1, 2, \ldots$ to distinguish the additional off-ladder EPs, ordered by increasing $\hat{a}_{\mathrm{EP}}$. For the gravitational-led $(2, 2;\{1,2\})$ pair, three distinct EPs are found: $\mathrm{EP}^{G,(0)}_{2,2;\{1,2\}}$ belongs to the near-extremal ladder, whereas $\mathrm{EP}^{G,(1)}_{2,2;\{1,2\}}$ and $\mathrm{EP}^{G,(2)}_{2,2;\{1,2\}}$ lie away from this sequence. Similarly, the electromagnetic-led $(2, 2;\{1,2\})$ pair possesses two EPs: $\mathrm{EP}^{E,(0)}_{2,2;\{1,2\}}$ lies on the near-extremal ladder, while $\mathrm{EP}^{E,(1)}_{2,2;\{1,2\}}$ is located off it. Thus, the KN parameter space contains not only ordered overtone ladders but also additional off-ladder coalescences carrying the same inherited overtone labels. The EP spectrum therefore has a more intricate organization than can be described by a single overtone ladder.

A closely related cascade of EPs has previously been found for massive scalar perturbations of near-extremal Kerr black holes, where successive overtone pairs become degenerate and the higher-overtone EPs accumulate toward the extremal limit [45]. Related exceptional-line and EP sequences have also been identified in the scalar Kerr–Newman spectrum [46]. The present results differ in several important respects. First, the EP ladders arise in the fully coupled gravitoelectromagnetic perturbations of the KN background, with the two control parameters supplied intrinsically by the black hole spin and charge rather than by an additional scalar-field mass. Second, whereas the displayed scalar-field ladders in Kerr black hole begin with the fundamental–first-overtone pair $\{0,1\}$, the lowest identified member of a gravitoelectromagnetic ladder in KN black hole depends on the angular sector. For example, the gravitational-led $(\ell,m)=(3,2)$ ladder begins with $\{1,2\}$, while the electromagnetic-led $(2,1)$ ladder begins with $\{2,3\}$. Moreover, the occurrence of additional off-ladder EPs with the same mode labels shows that the coupled spectrum is not organized by a single cascade alone. To our knowledge, these results provide the first evidence for overtone ladders, together with repeated on- and off-ladder coalescences, in the fully coupled Einstein–Maxwell QNM spectrum.

Taking $\mathrm{EP}^{E}_{1,1;\{0,1\}}$ as a representative example, Figure 5 examines the spectral structure around the coalescence and tests whether it exhibits the characteristic local branch-point structure of an EP. Using $\epsilon=\sqrt{1-\hat{a}^{2}-\hat{Q}^{2}}$ and $\vartheta=\tan^{-1}(\hat{a}/\hat{Q})$, we parameterize a closed loop around the EP by $u=\ln(\epsilon/\epsilon_{\mathrm{EP}})/r_{\epsilon}=\cos\phi$ and $v=(\vartheta-\vartheta_{\mathrm{EP}})/r_{\vartheta}=\sin\phi$, with $0\leq\phi\leq4\pi$. After one circuit, the two tracked QNM frequencies exchange branches; only after a second circuit does each frequency return to its initial branch. The squared frequency difference $D=[M(\omega_{A}-\omega_{B})]^{2}$ is invariant under the branch exchange, and its unwrapped phase increases by $2\pi$ per circuit. As an independent check of the local scaling of the frequency difference, we consider radial cuts in parameter space through the two EPs at their respective fixed values of $\vartheta=\vartheta_{\mathrm{EP}}$. Defining $\rho=|\ln(\epsilon/\epsilon_{\mathrm{EP}})|$, we fit

The fits yield $p=0.5028$ and $0.5029$ for the electromagnetic-led $\ell=m=1$ and $\ell=m=2$ $\{0,1\}$ pairs, respectively. Both values are very close to $1/2$, consistent with the square-root scaling of the frequency difference near an EP. The branch exchange, phase winding, and square-root scaling therefore provide mutually consistent evidence for the local EP structure.

Figure 6 compares the charge-dependent QNM spectra in the vicinity of the lower-overtone $\mathrm{EP}^{\mathrm{E}}_{1,1;\{1,2\}}$ and the higher-overtone $\mathrm{EP}^{\mathrm{E}}_{2,2;\{4,5\}}$. For the lower-overtone pair, varying $\hat{Q}$ at fixed $\hat{a}=\hat{a}_{\mathrm{EP}}$ traces the two branches through the EP at $\Delta\hat{Q}=0$, where they coalesce. For the higher-overtone pair, changing the fixed spin by only $\delta\hat{a}=4.0192\times10^{-8}$ noticeably changes how the two QNM frequencies vary with charge. The two complex frequencies become equal only when $\hat{a}=\hat{a}_{\mathrm{EP}}$. When the spin is shifted slightly away from this value, their real or imaginary parts may cross separately, but the full complex frequencies remain different.

Conclusions

We have investigated the EFs and EPs of KN QNMs using the exact separated gravitoelectromagnetic perturbation equations. The EFs shown in Figure 1 and Figure 2 exhibit distinct dependencies on black hole spin and charge across perturbation sectors and overtones. For the four selected gravitational-led modes, the EF magnitudes decrease overall as the black hole spin increases at fixed charge, but the suppression rates differ substantially between overtones and reverse the relative ordering of the fundamental and first-overtone factors in the $\ell=m=2$ sector. In the electromagnetic-led sector, the EF magnitudes of the $(1,1,0)$, $(2,2,0)$, and $(2,2,1)$ modes also decrease with increasing spin, whereas the $(1,1,1)$ mode exhibits nonmonotonic behavior. Their charge dependence is likewise mode dependent. The relative magnitudes, phases, and complex-plane trajectories further show that rotation can reorganize the excitation properties of different overtones. These results demonstrate that the excitation properties of KN QNMs depend on both the black hole hole parameters and the mode labels, rather than following a universal trend.

We have also identified 42 EPs in the coupled gravitoelectromagnetic QNM spectrum, with 23 in the gravitational-led sector and 19 in the electromagnetic-led sector, as summarized in Table II and Figure 4. Their distribution reveals two notable forms of spectral organization. First, both sectors contain ordered overtone ladders formed by successive coalescences of neighboring branches $\{n,n+1\}$, with the higher-overtone members moving progressively closer to the extremal KN boundary. The lowest identified member of a ladder depends on the angular sector and need not involve the fundamental mode. Second, the same sector, angular mode, and overtone pair may correspond to several distinct EPs. Thus, the EP spectrum contains both systematic overtone sequences and repeated coalescences that cannot be described by a single ladder alone. These results extend the overtone-ladder structure previously observed for probe scalar fields [45] to the fully coupled Kerr-Newman QNM spectrum. Figure 5 and Figure 6 confirm the characteristic square-root branch structure of the selected EPs and reveal the strong parameter sensitivity of higher-overtone QNMs. Together, these results demonstrate the rich non-Hermitian structure of the KN QNM spectrum.

Several directions deserve further investigation. A broader search over additional angular sectors and higher overtones is needed to determine how far the overtone ladders extend and to clarify how the accompanying off-ladder EPs are distributed relative to these ladders. An analytic treatment of the near-extremal regime may clarify why the higher-overtone members approach the extremal boundary and whether these sequences are related to photon-sphere, damped, or near-horizon mode families [27, 28]. An important extension is to study QNM excitation by specified external sources, such as a compact object plunging into a KN black hole from the innermost stable circular orbit. Combining the source-independent EFs with source-dependent contributions and the reconstruction of physical gravitational and electromagnetic perturbations would provide a more complete description of black hole ringdown signals. It would also be valuable to formulate the resonant response directly near an EP, where the coalescence of QNM frequencies and eigenfunctions may require a generalized modal description beyond a sum of isolated simple modes [47]. Such investigations would deepen our understanding of how gravitational-electromagnetic coupling shapes the dynamical response and spectral properties of rotating, charged black holes.

Acknowledgments

The work is funded by the National Natural Science Foundation of China with Grants No. 12405059 and No. 12505087, the Startup Fund for Advanced Talents of Nanchang University with Grant No. 28170256, Key Program of the Natural Science Foundation of Jiangxi Province under Grant No. 20232ACB201008 and the Ganpo Talent Support Program with No. gpjc20260105.

Angular formulation

We compute KN QNMs using the spectral method of Ref. [30], based on the exact separated gravitoelectromagnetic equations derived in Ref. [29]. The method combines an angular Galerkin expansion with radial Chebyshev collocation to solve simultaneously for the complex frequency $\omega$ and the separation constant $\lambda$. We outline the formulation and numerical procedure below.

We use the mode convention and the angular functions $W$ and $T$ defined in Sec. II A and Eq. (5). The label $\zeta=\pm1$ distinguishes the two spin descriptions. In the two-component basis used in Ref. [30], the restricted angular operator is

The constant transformation

brings its eigenvalue equation into the form of (4). We denote the resulting matrix differential operator by $\mathcal{B}(\omega)$. This change of basis preserves both the eigenvalues and the regularity conditions at the poles. It remains valid for complex $\omega$, although the angular operator is then generally non-self-adjoint.

At $\omega=0$, the rotation-dependent terms in the angular equations vanish, and the eigenvalues for $\ell\geq2$ are given by Eq. (8). The label $b$ identifies an eigenvalue branch continued from this limit and is distinct from the spin label $\zeta$. The correspondence between $b$ and the gravitational-led and electromagnetic-led sectors is established through the Kerr limit discussed below.

Sector identification in the Kerr limit

The apparent divergence of $\lambda$ as $Q\to0$ can be separated from its finite part by writing

Following Ref. [30], introduce $q=Q/M$ and rescale the angular and separated radial amplitudes according to

The branches of interest have finite $\widetilde{\lambda}$ as $q\to0$. These rescalings yield angular and radial equations with finite coefficients in this limit. At nonzero charge, the two rescaled angular problems are reparameterizations of the same spectrum; their distinct Kerr limits identify the two sectors.

The Kerr limit is taken after the sector-dependent field rescalings and the corresponding redefinition of the radial master variable described in [30]. It cannot be obtained by directly setting $Q = 0$ in Eq. (13), whose coefficients are singular in this limit. In particular, the limiting radial equations reduce to the Teukolsky equations with spin $s = 2\varsigma$ for $b = +1$ and $s = \varsigma$ for $b = -1$. With the convention $\Lambda_s = sA_{\ell m} + a^2\omega^2 - 2am\omega$, their separation constants satisfy

Thus $b = +1$ and $b = -1$ label the gravitational-led and electromagnetic-led modes, respectively. Away from $Q = 0$, both sectors retain gravitational–electromagnetic coupling.

Spectral discretization

Angular projection

For each angular component $k = 1, 2$, set $z = \cos\theta$ and $j_{\min,k} = \max(|m|, k)$, and use the pole exponents $\alpha_k$ and $\beta_k$ from (7), with the component index $s$ replaced by $k$. We approximate the regular solution by

where $L$ is the angular cutoff and

The Jacobi polynomials provide the regular part of each component, while the factors impose the required behavior at both poles. The power expansions with half-integer arise from the rescaled angular variables.

The basis is orthonormal with respect to

Writing $d = j-j_{\min,k}$, the normalization is

Here $k$ denotes the component and $d$ the polynomial degree; neither is the radial overtone index $n$.

Radial compactification

We start from the radial equation (14), with the coefficients in (15) and the background and radial functions defined in Section II A. The horizon and infinity conditions are those of Section II B. For $Q\omega\ne0$, we use the coordinate $x$ introduced in Section II B and extract the outgoing asymptotic factor with horizon normalization:

The horizon and infinity correspond to $x = 0$ and $x = 1$, respectively, and $F(r_+) = 1$.

To transform the derivatives, define

Then

where $\tau' = -2r_+/r^3$ and $g' = -4iM\omega/r^2$. After multiplying the transformed equation by $1-x$, we obtain

with

All radial functions are evaluated at $r = r_+/(1-x)$. The large-$r$ cancellations are performed analytically before evaluating the coefficient limits at $x = 1$.

Let $\Psi_{jk} = \Phi_{jk}e_k$, with $e_k$ the two-component coordinate vectors. Since $\mathcal{B}(\omega) = \mathcal{B}_0+\omega\mathcal{B}_1$, Galerkin projection gives

where

The vector inner product includes a sum over components. Matrix elements are evaluated by Gauss–Legendre quadrature with $4L+16$ nodes in $z$, retaining the measure in (C3); basis derivatives are evaluated analytically. The resulting matrix has dimension

All basis functions allowed by $j_{\min,k}$ are retained, including the additional lower-order function when $|m| = 1$. Although $G$ is affine in $\omega$, the selected eigenvalue $\lambda(\omega)$ generally depends nonlinearly on frequency.

Radial collocation and boundary conditions

We discretize Eq. (41) on $N+1$ Chebyshev–Lobatto points,

If $D$ denotes the differentiation matrix with respect to $x$, the radial matrix is

The endpoint rows use the finite limiting coefficients. Generically, $c_2(0)=c_2(1)=0$, and these rows reduce to first-order compatibility relations.

These relations implement the limiting differential equation, but smoothness alone does not establish the radiative boundary conditions for complex frequencies. In particular, after removing $F$, the competing asymptotic solution contains an $e^{-2i\omega r}$ factor and can be smooth and exponentially small at infinity for damped modes. As emphasized in Ref. [30], identifying physical QNMs therefore requires the prescribed asymptotic branch, together with resolution studies and independent frequency comparisons.

Simultaneous solution and mode tracking

A QNM corresponds to a common pair $(\omega,\lambda)$ for which both $G(\omega,\lambda)$ and $R(\omega,\lambda)$ admit nonzero null vectors. We locate these pairs using the normalized residual construction of Ref. [30]. For either matrix $\mathcal{M}=G$ or $R$, fixed row scaling is applied, and one equation is replaced by a normalization condition $v_k=1$. If row $j$ is omitted, this construction is

where $e_j$ and $e_k$ are coordinate vectors of the appropriate dimension. Provided $\widehat{\mathcal{M}}$ is invertible, $F_{\mathcal{M}}=0$ restores the omitted equation and yields a nontrivial null vector of $\mathcal{M}$. The normalization fixes an arbitrary eigenvector amplitude and adds no physical boundary condition.

The replacement row and normalization component are selected from the left and right singular vectors associated with the smallest singular value at the initial estimate. These choices, the row scalings, and the residual scalings remain fixed during each local search. We solve

as four real equations for the real and imaginary parts of $\omega$ and $\lambda$, using bounded nonlinear least squares followed by local root refinement.

For $\ell\ge2$, the angular branch is initialized at $\omega=0$ using (8). The $\ell=1$ electromagnetic-led branch is identified within the regular $\omega=0$ angular problem. The selected eigenpair is continued to the trial frequency along $\omega(\eta)=\eta\omega_{\mathrm{trial}}$, $0\le\eta\le1$, choosing successive eigenvectors by their largest normalized overlap. The angular branch assignment is checked again at the converged frequency.

During a scan in $(a,Q)$, solutions at neighboring parameter values provide initial estimates for the next calculation. The radial overtone label is inherited by continuation from an identified reference mode; it is not determined by the angular labels alone. Independent changes of $N$ and $L$ test the radial and angular truncations, respectively. Comparison with published QNM frequencies provides a further validation of the computed roots, as discussed in Appendix § 1.9.

Numerical check of our method

Table III compares our QNM and EP frequencies with published results. The QNMs are compared at the same background parameters and with the same mode labels. For the EPs, we use our results from Table II and compare the frequencies at the EP locations found in each study. The relative frequency difference is

$(\hat{a},\hat{Q};b)$	$M\omega$	$\epsilon_{\mathrm{QNM}}$
QNM benchmarks: Ref. [30]
$(0.2,0.4;+1)$	$0.4099175-0.088621\,i$	$1.854\times10^{-5}$
$(0.2,0.4;-1)$	$0.5150483-0.0955546\,i$	$3.209\times10^{-5}$
EP benchmarks: Ref. [38]
$(\hat{a},\hat{Q};b)$	$M\omega_{\mathrm{EP}}$	$\epsilon_{\mathrm{EP}}$
$(0.304116510,0.952362666;+1)$	$0.564398227-0.028304554\,i$	$9.537\times10^{-5}$
$(0.461271757,0.883991733;-1)$	$0.793226261-0.144859314\,i$	$1.986\times10^{-4}$

where $p$ and $r$ denote our result and the reference result, respectively. For the EPs, we use the published central values.

References