Analytical solutions of CPT-odd Maxwell equations in Schwarzschild spacetime

Lorentz symmetry (LS) is a fundamental symmetry in both general relativity (GR) [1] and the Standard Model [2] of particle physics. Some candidate theories of quantum gravity, such as certain formulations of loop quantum gravity [3, 4] and modifications of string theory [5, 6], allow for small deviations from exact Lorentz invariance at very high energies. In some scenarios, this could lead to extremely small Lorentz-violating effects at lower energies. However, Lorentz-violating signals at energy scales accessible in high-energy astrophysical observations ($ \sim 10^{11} \text{GeV} $ [7]) are expected to be extremely small and are generally suppressed by a small ratio involving the Planck scale $ 1.22 \times 10^{19} \text{GeV} $, as suggested by dimensional analysis and observational constraints.

Extremely small Lorentz violation (LV) effects may accumulate over long distances and at high energies in certain models, making them potentially detectable through terrestrial observations. Additionally, some experiments and astrophysical observations can test processes that are strictly forbidden in standard Lorentz-invariant (LI) physics but may occur in LV scenarios, such as vacuum birefringence [8, 9], photon decay [10], and photon splitting [11]. These observations have established stringent constraints on LV with high precision, particularly through high-energy observatories such as Pierre Auger and LHAASO [12, 13]. Moreover, most of these observatories primarily rely on multi-wavelength observations and long-distance photon propagation to study astrophysical events.

As a comprehensive framework of effective field theory, the Standard Model Extension is capable of describing both small deviations from LS in flat spacetime [6, 14] and violations of local LS in gravitational contexts [15], thus encompassing both high-energy phenomena in flat spacetime [6] and gravitational effects or particle motions in curved spacetime [16, 17]. In recent years, there has been increasing interest in probing LV in astrophysics through observations of CMB photons [18, 19], neutrinos [20], and gravitational waves [21]. A natural question is whether interesting LV effects manifest in intrinsically curved spacetime. Regarding the constraints on LV from studying the cosmological propagation of GRB or CMB photons, a key assumption is that spacetime is described by the Friedmann-Lemaître-Robertson-Walker (FLRW) metric, which is curved but conformally flat. Here, we investigate LV electrodynamics in the simplest non-conformally flat curved spacetime: the Schwarzschild metric.

There has been a growing number of studies on photon behavior in curved spacetime, particularly following the successful capture of black hole (BH) images by the Event Horizon Telescope (EHT) Collaboration [22]. The photon sphere and subrings at the edge of the black hole shadow may encode crucial information about potential new physics beyond GR [23]. This highlights the necessity of studying photon behavior in curved spacetime near BHs.

As a first attempt, we aim to investigate the asymptotic behavior of CPT-violating (CPTV) photons in the Schwarzschild geometry using the null formalism. Intuitively, employing the null formalism to study massless particles is, in some sense, analogous to describing the motion of massive particles — such as a gyroscope — using an orthonormal tetrad in its instantaneous rest frame, even though massless particles do not possess a rest frame. However, the underlying principle remains the same: null tetrads naturally accommodate massless particles that follow null trajectories. Moreover, the null formalism offers unique advantages in analyzing the asymptotic behavior of massless particles, particularly photons in this context. It significantly simplifies the description of both the tangent of the null geodesic for photons and gravitons, as well as their polarization states [24]. Consequently, it provides a coordinate-independent framework for studying photon dynamics with a clear geometric interpretation, especially in highly curved spacetimes. Notably, it also offers a physically transparent decomposition of the Faraday tensor into ingoing, outgoing, and Coulomb modes.

In the presence of CPT or Lorentz violation, contrary to conventional expectations, the asymptotic behavior of photons may be qualitatively altered [25]. For instance, vacuum birefringence can be understood in terms of the modified topology of the light cone structure, where different helicities of CPT-odd photons experience distinct causal cones [26].

Pioneering studies on exact solutions for LI photon fields as perturbations in given background geometries include vacuum solutions for photon fields in the Kerr spacetime [27] as well as solutions for a point charge near a Schwarzschild BH [28] and a Kerr BH [29]. Bičák et al. studied photon fields in curved spacetime within the Newman-Penrose (NP) framework [30], including the Schwarzschild [31], ReissnerNordström (R-N) [32], and Kerr [33] BH backgrounds. As a preliminary attempt, here we study the behavior of CPT-odd photons in the Schwarzschild geometry following a similar approach.

In Sec. II, we review the Newman-Penrose (NP) formalism and discuss some earlier studies on the LI Maxwell equations in curved spacetime using the NP framework. Then, we examine the CPT-violating Maxwell equations within the NP formalism in curved spacetime. Next, we present a method to solve the coupled Maxwell equations and provide special solutions in Sec. III. In the last section, we summarize our results and provide a short conclusion. In this paper, the signature of the metric tensor $ g_{\mu \nu} $ is chosen to be $ (+,-,-,-) $, and we use geometric units with $ \epsilon_0=\mu_0=c=G=\hbar=1 $. The notation conventions are as follows: spacetime indices are represented by Greek letters such as $ \mu, \nu, \rho $, while null tetrad indices are represented by Latin letters such as $ a, b, c $.

We study LV (more specifically CPTV) photon behavior within the photon sector of the minimal SME [6]. The action is given by

where $ \tilde{F}^{\alpha \beta}=\dfrac{1}{2}\epsilon ^{\mu \nu \alpha \beta } F_{\mu \nu } $, and $ (k_{AF})_\alpha $ is the CPTV coefficient [6] [34]. The coefficients $ (k_{AF})_\alpha $ are real and have mass dimension one. The equation of motion is

Using the null tetrad $ e_a^{\ \mu} =\left( l^{\mu},n^{\mu},m^{\mu},\bar{m}^{\mu}\right), a=1,2,3,4, $ correspond respectively to $ l,n,m, \bar{m} $, and the electromagnetic field tensor $ F_{\mu\nu} $ can be decomposed into three complex Newman-Penrose (NP) scalars,

Conversely, the electromagnetic field tensor $ F^{\mu\nu} $ can be expressed as

where $ a_{[\mu}b_{\nu]}:=\dfrac{1}{2}(a_{\mu}b_{\nu}-a_{\nu}b_{\mu}) $, and "$ \text{c.c.} $" denotes the complex conjugate.

This paper primarily focuses on photon fields propagating in the vicinity of a Schwarzschild BH, with the line element

where $ g(r)=1-\dfrac{2 M}{r} $. The corresponding null tetrad basis vectors are given by

with their corresponding covariant components given by

In terms of the null contravariant vectors, we project the derivatives into null directions:

following Ref. [30]. To simplify the calculations, we assume that the CPTV coefficient is spherically symmetric, consistent with the Schwarzschild background, and restrict our analysis to stationary electromagnetic fields. Consequently, the directional derivatives reduce to $ D=\partial_r $ and $ \Delta=-\dfrac{1}{2} g(r) \partial _r $. In other words, we focus on the behavior of static electric and magnetic fields in the presence of the CPT-odd term.

Similarly, we define the spin coefficients as $ \gamma_{abc} \equiv e_a^{\ \mu} e_{b\mu;\nu} e_c^{\ \nu} $, following the conventions in Ref. [35]. In the Schwarzschild metric, the nonzero spin coefficients are given by

It is important to note that the NP formalism employed here may not be fully applicable within a generic Lorentz-violating theory. In this paper, we implicitly adopt the test particle assumption, wherein the background metric is assumed to remain unaffected by the Lorentz-violating matter fields — specifically, the electromagnetic fields under consideration. Within this framework, the use of a complete and quasi-orthonormal null tetrad remains appropriate, given that it effectively captures the essential features of the quasi-null wavefront associated with CPT-violating electromagnetic fields. This is because the null tetrad can be regarded as a natural choice for describing massless particles.

However, for a more rigorous treatment — particularly when the back-reaction of matter fields on the spacetime metric is taken into account — the standard null tetrad may no longer suffice. In such cases, it may be necessary to generalize the framework, for example, by employing a quasi-null tetrad, as used in the analysis of gravitational wave polarization [36].

The LI Maxwell equations in the NP formalism have been derived in the appendix of Ref. [30] and in Chapter 1.8 of Chandrasekhar’s textbook [35]. For the CPTV contribution, the term $ (k_{AF})_{\mu} \widetilde{F}^{\mu\nu} $ in Eq. (2) can be projected onto the null tetrad basis as $ (k_{AF})_{a } \widetilde{F}^{a b} $, where $ \widetilde{F}^{a b}\equiv \dfrac{1}{2}\epsilon _{a b c d } F^{c d} $ and $ (k_{AF})_{a} \equiv(k_{AF})_{\mu}{\boldsymbol{e}}_a ^{\ \mu} \ (a=1, 2, 3, 4 $), which are tetrad components of the CPTV coefficient $ (k_{AF})_{\mu} $. For simplicity, we define $ (k_{AF})^{a}\equiv k^a $. As an example, for$ b=1 $, we obtain

Here, we use $ \epsilon_{1234} =i $, which follows from the definition of the complex null tetrad given in Eq. (6). The CPTV Maxwell equations in NP form are then given by

where $ J_l = l_\mu j^\mu $, $ J_n = n_{\mu} j^\mu $, etc.

Given that we assume that the CPTV coefficient $ k^\mu $ is spherically symmetric, a simple example is to consider only $ k^t \neq 0 $. Given that $ k^t = k^{1} l^{t} + k^{2} n^{t} \neq 0 $, this implies that $ k^{1} $ and $ k^{2} $ cannot be zero, while $ k^3 = k^4 = 0 $. By substituting the spin coefficients from Eq. (9) and the differential operators from Eq. (8), Eq. (11) simplifies to

where the differential operators $ \eth $ and $ \bar{\eth} $ are defined as

The eigenfunctions of $ \eth $ and $ \bar{\eth} $ are the spin-weighted spherical harmonics $ _{s}Y_{l m} $. Here, $ s=1, 0, -1 $ correspond to the spin weights of $ \Phi_0, \Phi_1, \Phi_2 $, respectively. For $ s=0 $, $ _{0}Y_{l m} = Y_{l m} $ are the standard spherical harmonics, and the indices $ s, l, m $ satisfy $ |m|\leq l $ and $ |s| \leq l $. For further details on the definition and properties of spin-weighted spherical harmonics, see Appendix A or Ref. [37].

To solve Eq. (12), we adopt the Teukolsky approach [38]. The key is to leverage the commutation relations [30] between differential operators (Eq. (8)) and spin coefficients (Eq. (9)) to decouple the four coupled Maxwell equations. Acting with a combination of the operators on Eq. (12), we obtain a set of partially decoupled equations:

Here,

For later convenience, we define $ \sum_{lm} \equiv \sum_{l=1}^{\infty} \sum_{m=-l}^{l}, $ and given that spherical symmetry is preserved in at least a special preferred reference frame, we expand the three complex scalars using spin-weighted spherical harmonics:

Inspection of Eq. (15) reveals that in the absence of CPTV coefficients, the three equations are decoupled. Given that the CPTV coefficients are experimentally constrained to be extremely small, say $ |k_{AF}|\le10^{-44} $GeV, [39–40], we may treat the CPTV terms on the right-hand side of Eq. (15) as perturbations.

Thus, the radial functions can be expanded in powers of CPTV coefficients $ k^1 $ and $ k^2 $ as

where the superscripts "(0)," "(1)," etc., indicate the corresponding order of $ k^1 $ and $ k^2 $. For example, $ R_{a|l m}^{(0)} $corresponds to the zeroth-order function without LV correction. The expansions of the NP scalars are thus given by

Keeping only the linear-order terms in the CPTV coefficients, Eq. (15) can be separated into two sets: the zeroth-order LI equations,

and the first-order equations with linear CPTV corrections,

The zeroth-order equations indicate that the external charge and current act as sources for the zeroth-order NP complex scalars, which are linear combinations of the components of the Faraday tensor. Similarly, the first-order equations show that the zeroth-order Faraday fields act as sources for the first-order CPTV corrections to the Faraday tensor.

We begin by solving the zeroth-order decoupled equations given in Eq. (20). By leveraging the orthogonality relations of spin-weighted spherical harmonics (Eq. (51)) and substituting Eqs. (8), (9) and (17) into Eq. (20), we obtain the following radial equations:

Here, $ R(r)'' $ and $ R(r)' $ denote the second- and first-order derivatives with respect to r, respectively. The source terms are given by

In the absence of source terms in Eq. (22), by introducing the variable transformation $ x \equiv r/(2M) $, the homogeneous equations can be rewritten in the standard form of hypergeometric equations:

The general solutions of the hypergeometric equations are

For $ l \neq 0 $, the linearly independent solutions for $ R_{a|l m} $ $ (a=0,1,2) $ are

For $ l=0 $, given that $ R_{0|l} $ and $ R_{2|l} $ correspond to spin-weight $ s=\pm1 $ and the spin-weighted functions are only defined for $ l\ge|s| $, $ R_{a|0} $ is not defined except for $ a=1 $, which corresponds to spin-weight $ s=0 $. The linearly independent solutions of $ R_{1|0} $ are

To fully characterize the solutions, we examine their asymptotic behaviors at both spatial infinity and near the event horizon. The asymptotic expressions of the solutions in the far-field regime are

Given that $ R^{(\text{I})}_{a|l} $ diverges for $ l > 0 $, only $ R^{(\text{II})}_{a|l} $ remains well-behaved at spatial infinity, ensuring an appropriate physical decay. The asymptotic expressions of the solutions near-horizon regime are

Thus, only $ R^{(\text{I})}_{a|l} $ solutions are regular near the event horizon. For $ l = 0 $, only $ R_{1|0} $ exists, and among its solutions, only $ R_{1|0}^{(\text{II})} $ remains well-behaved in both the far-field ($ x \to \infty $) and near-horizon ($ x \to 1 $) limits. In short, to have physically acceptable solutions, we have to ensure that the general solutions (Eq. (25)) exhibit proper and reasonable asymptotic behaviors both at infinity and near the horizon, namely $ R^{(\text{II})}_{a|l} $ and $ R^{(\text{I})}_{a|l} $ are chosen, respectively.

Next, we consider the non-homogeneous case of Eq. (22) in the presence of source terms. We assume that the source is localized within the finite region $ r_1 \leq r \leq r_2 $, where $ 2M < r_1 < r_2 < \infty $. In this range, we may let $ r_1 $ be sufficiently larger than the Schwarzschild radius of the compact object, say, $ r_1=10\,r_S=20M $, and $ r_2 $ far from infinity, i.e., $ r_2<<\infty $. This ensures that essential non-linear or curvature effects do not become dominant. More precisely, given that our analysis is primarily based on the test particle assumption, any back-reaction of electromagnetic perturbations on the background spacetime metric, as well as potential instability issues, are beyond the scope of this work. While these topics are indeed interesting and important, they are significantly more complex than the current study and are worth exploring in future research.

Based on the previously analyzed asymptotic behaviors, we use the fundamental solution set $ \{ R_{a|l}^{(\text{I})}, R_{a|l}^{(\text{II})} \} $ for $ a=0,1,2 $ to construct the general solution of Eq. (22). The radial solutions in different regions are provided next.

For $ l \neq 0 $, the solution takes different forms depending on the radial range:

• In the region $ 2M < r < r_1 $, where the solution is near the source but outside the event horizon, the general form is

where the coefficients $ u_{l m} $ correspond to $ a_{l m}, c_{l m}, e_{l m} $, respectively.

• In the region $ r > r_2 $, far from the source, the solution takes the form

where the coefficients $ v_{l m} $ correspond to $ b_{l m}, d_{l m}, f_{l m} $, respectively.

This piecewise formulation ensures that the solutions satisfy the appropriate boundary conditions at both spatial infinity and the event horizon while maintaining mathematical consistency across the defined radial domains.

There exists a special case for $ l = 0 $, where the radial function takes the form

and

Here, $ E_{a} $ and $ E_{b} $ are constants that, along with the previously introduced coefficients $ u_{l m} $ and $ v_{l m} $, will be determined using the method outlined below.

For the case of given sources in Eq. (22), we apply the method of variation of constants [41] (see also Appendix C). The corresponding particular solutions are obtained as follows:

where $ a=0,1,2 $. The function $ W\left( R_{a|l}^{(\text{I})}, R_{a|l}^{(\text{II})}, \xi\right) $ represents the Wronskian determinant of the two fundamental solutions $ R_{a|l}^{(\text{I})} $ and $ R_{a|l}^{(\text{II})} $ evaluated at ξ. By comparing Eq. (36) with Eq. (25), we obtain the following integral expressions for the expansion coefficients:

where ε is an infinitesimal positive constant and $ x_a\equiv r_a/(2M) $ , with $ a=1,2 $ , specify the external source region in radial direction. After performing explicit calculations, we find that the Wronskians are approximated as follows:

The detailed derivation of these expressions is provided in Appendix B. By substituting Eq. (38) into Eq. (37), we obtain the final expressions:

For coefficients $ E_{a} $ and $ E_{b} $, we refer to Eq. (12), from which we obtain $ E_{a} = 0 $ and $ E_{b} = \sqrt{\pi}\,e $, where e is the elementary charge obtained from integrating the charge aspect $ \Phi_1^0 $ of the Coulomb mode $ \Phi_1 $. Next, we consider the first-order correction, i.e., Eq. (21). By substituting the right-hand side of Eq. (19) into Eq. (21) and utilizing the properties of spin-weighted spherical harmonics (Eq. (51)), we obtain

Using the orthogonality of spin-weighted harmonics (Eq. (52)), we obtain the following equations for the radial components:

where we define a set of effective source terms $ J_{a|l m}^\text{LV} $ ($ a=0,1,2 $) induced by LV as follows:

From Eqs. (26)-(28), we observe that $ R_{a|l}^{(\text{II})} $ are functions of $ x^{-1} $. As $ x \rightarrow \infty $, we have $ x^{-1} \rightarrow 0 $, implying that $ R_{a|l}^{(\text{II})} \rightarrow (-x)^{-2-l} $ for $ a=0,1,2 $. Consequently, we approximate

By substituting these approximations and the spin coefficients into Eq. (44) and setting $ r=2Mx $, we derive the governing radial equations (Eq. (43)) up to the lowest order approximations of $ R_{a|l m}^{(0)} $ and their complex conjugate $ \bar{R}_{a|l m}^{(0)} $, where $ a=0,1,2 $.

Similar to the zeroth-order case, we express the general solution as

For the given LV sources $ J_{a|l m}^\text{LV} $ ($ a=0,1,2 $), following the procedure used in Eq. (36), we obtain the particular solution:

Comparing Eq. (47) with Eq. (46), we determine the following coefficients:

It can be seen that if the CPTV term is treated as an effective source term, the LV effect may be characterized by the effective currents $ J_{a \mid l m}^\text{LV}, a=0,1,2 $. In the case of point charges or other sources (see Ref. [31]), as $ x \rightarrow \infty $, the two quantities $ b_{l m}^{(1)} $ and $ f_{l m}^{(1)} $, corresponding to the spin-weight $ -1 $ and $ +1 $ modes, respectively, exhibit nearly identical asymptotic behavior. The primary difference between these spin-weight $ \mp 1 $ modes arises from the LV-induced effective currents $ J_{0 \mid l m}^\text{LV} $ and $ J_{2 \mid l m}^\text{LV} $, which is consistent with expectations.

The analytical solutions obtained in this section reveal that CPT-odd corrections induce significant mixing between different spin-weight modes of the electromagnetic field in the far-field regime ($ r \gg 2 M $), even though the electromagnetic field produced by static currents decays rapidly in this region. In the classical LI case, the zeroth-order solutions exhibit radial dependencies governed by the expansion of hypergeometric functions $ R_{a \mid l}^{(0)}(x) $, while the angular components are precisely described by spin-weighted spherical harmonics $ { }_s Y_{l m} $. The degeneracy in the behavior of different spin-weight modes is a direct consequence of the underlying spacetime symmetry. However, the introduction of the LV term breaks this symmetry, manifesting as perturbations to the original solutions through the equivalent source terms $ J_{a \mid l m}^\text{LV}, a=0,1,2 $. Specifically, the radial behavior $ R_{a \mid l m}^{(1)} $ still follows the power-law decay $ x^{-l-2} $ (where $ x=r / 2 M $), reflecting the suppression of electromagnetic multipole radiation. Notably, in the lowest-order case $ (l=0), R_{a \mid l m}^{(1)} \propto x^{-2} $, leading to $ \left|R_{a \mid l m}^{(1)}\right|^2 \sim r^{-4} $, which confirms the absence of net energy-momentum flux in the far-field region for electromagnetic fields generated by external static currents [42]. Furthermore, the linear relation $ R_{a \mid l m}^{(1)} \propto k^a \cdot R_{a \mid l m}^{(0)} $ (for $ a=1,2 $) in the far-field solutions indicates that LV effects could be extremely small, yet their cumulative impact might become significant in high-energy astrophysical environments, such as active galactic nucleus (AGN) jets.

Moreover, it is interesting to note that in time-dependent scenarios where radiation is present, the two helicity $ \pm $ modes correspond to different polarization states, making helicity-dependent effects particularly pronounced for higher multipole moments ($ l \geq 1 $). This suggests that LV effects may be more significant in radiation from higher multipole moments and could be constrained through cumulative effects in long-baseline photon propagation, such as polarization angle evolution in gamma-ray bursts (GRBs). Furthermore, the study of far-field behaviors in radiative cases reveals distinct deviations from LI electrodynamics [25]. For instance, logarithmic correction terms may arise owing to the absence of an additional derivative in the CPT-odd $ k_{A F} $ term compared to the LI Maxwell theory [25, 43]. Additionally, a form of energy flux cancellation between lower and higher frequency modes may occur to ensure the absence of net radiation for a charged particle moving at constant velocity [44]. For time-dependent situations, a nonzero net energy flux is expected, just as for conventional electrodynamics. A nontrivial example is provided by dipole radiation [45], though its polarization structure is non-perturbative in terms of the CPTV coefficients. In fact, a closer examination of Eqs. (12) reveals that the presence of the imaginary part of $ \Phi_1 $ in the first equation obstructs the separability of time and radial variables, u and r (for further details, see Eqs. (54c) in Ref. [25]). This issue may stem from the strong assumption of a constant radial component $ k^r $. A more physically plausible assumption — namely $ k^r\propto{\cal{O}}(1/r) $ — could resolve this difficulty. Under such a decay, the standard separability of u and r is recovered at sufficiently large radii, such as at null infinity. Then, the compatibility of asymptotic flatness with test particle assumption is ensured.

However, note that current constraints on dimension-3 CPTV coefficient $ k_\text{AF} $ are extremely stringent and have already reached $ |k_\text{AF}|\le 10^{-44} $ GeV. For a more detailed summary of the constraints, please see tables S3 and D16 in Ref. [46]. It is also important to distinguish between Lorentz violating constraints in curved spacetime from those in flat (or conformal flat) spacetime. The former may be subject to screen effect due to the minuscule nature of gravitational couplings [47]. While the existing constraints are so tight that further improvements through astrophysical observations may be impractical, this study still provides semi-analytical solutions to the CPTV Maxwell field equations in curved spacetime, which may be valuable from a theoretical standpoint.

In this paper, we investigate analytical solutions of the CPT-odd Maxwell equations in Schwarzschild spacetime using the NP formalism. By employing a perturbative approach, we treat the LI Maxwell equations as the zeroth-order approximation and incorporate the CPTV coefficients $ \left(k_{A F}\right)^\mu $ as first-order corrections. The electromagnetic field tensor is decomposed into three complex NP scalars $ \Phi_0, \Phi_1, \Phi_2 $, whose radial dependence is governed by hypergeometric functions, while the angular components are described by spin-weighted spherical harmonics. Specifically, the zeroth-order solutions preserve the spherical symmetry of the Schwarzschild metric, with multipole moments determined by hypergeometric radial functions. The introduction of CPTV coefficients $ \left(k_{A F}\right)^\mu $ induces anisotropic corrections by coupling different angular modes $ (l, m) $, even assuming that the CPTV coefficients are spherical symmetrically distributed. In the radiation case, this may also alter the polarization structure. If the radial behavior of radiation follows a similar scaling $ \left(R_{a \mid l m}^{(1)} \sim x^{-l-1}\right) $ with $ x=r / 2 M $, it may indicate the suppression of higher multipole radiation at large distances.

Although this is a very preliminary study on CPT-violating electrodynamics in curved spacetime, it represents the first attempt within the NP formalism. The interplay between geometric and LV effects may not only deepen our understanding of BH electrodynamics but also open new avenues for exploring physics beyond classical relativity. From a theoretical perspective, the results presented here provide a foundation for further investigations into CPT-violating and LV effects in curved spacetime. Future research could explore more complex spacetime backgrounds, such as Kerr spacetime, where frame-dragging and ergo-region dynamics may significantly influence LV effects. Extending this analysis to rotating black holes could reveal novel phenomena, such as LV-modified superradiance or potential imprints on BH shadow substructures observable by the Event Horizon Telescope. Moreover, the framework developed here is inherently adaptable to arbitrary orders of electric and magnetic multipole expansions. By extending the relation $ R_{a \mid l m}^{(1)} \propto x^{-l-2} $ to higher l, one could systematically analyze LV corrections to higher multipole moments, probing finer details of electromagnetic fields near compact objects.

We would like to acknowledge the valuable discussions with Zhanfeng Mai, and H. Wang would like to express his gratitude to Tengfei Li for his assistance during the calculation process.

The spin-weighted spherical harmonics could be defined in terms of the usual spherical harmonics as

It could also be represented as

From the definition above, the spin-weighted spherical harmonics have some useful properties:

The spin-weighted spherical harmonics obey the orthogonality condition

where $ \text{d} \Omega=\cos\theta\,\text{d}\theta\,\text{d}\phi $.

The Wronskian is a determinant constructed from n differentiable functions $ f_1, \ldots, f_n $ along with their first $ n-1 $ derivatives. Its explicit form is given by

where $ f_n^{(n-1)} $ is the $ (n-1)^\text{th} $ derivative of $ f_n $. For a general homogeneous differential equation

there exists a useful relation known as Abel's identity, which expresses the Wronskian of a set of solutions in terms of a known Wronskian at a reference point and the coefficients of the original differential equation:

This identity is particularly useful in computing the Wronskian of $ \Phi_0, \Phi_1, \Phi_2 $. For $ \Phi_0 $, using Eq. (24a), we obtain

By taking the limit $ x_0\rightarrow \infty $ and substituting $ R_{0|l}^{(\text{I})} $ and $ R_{0|l}^{(\text{II})} $ from Eq. (26), we obtain

For $ \Phi_1 $, using Eq. (24b), we obtain

By taking the limit $ x_0\rightarrow \infty $ and substituting $ R_{1|l}^{(\text{I})} $ and $ R_{1|l}^{(\text{II})} $ from Eq. (27), we obtain

For $ \Phi_2 $, using Eq. (24c), we obtain

Taking the limit $ x_0\rightarrow \infty $, and substituting $ R_{2|l}^{(\text{I})} $ and $ R_{2|l}^{(\text{II})} $ from Equation (28), we obtain

In this subsection, we introduce an approach [41] for solving non-homogeneous differential equations. Consider a general non-homogeneous ordinary differential equation of the form

where $ L(y) $ represents a linear differential operator. The general solution of Eq. (62) can be written as

where $ y_p $ is a particular solution of Eq. (62). When $ b(x) = 0 $, the equation reduces to a homogeneous form, whose general solution is given by $ c_1 y_1 + \cdots + c_n y_n $, where $ c_1, \cdots, c_n $ are constants.

To construct a particular solution $ y_p $ with a structure similar to the homogeneous solution, we assume the ansatz

where $ u_1, \cdots, u_n $ are functions rather than constants. These functions satisfy the following system of equations:

This system of equations can be solved using Cramer's rule, yielding

where $ W(y_1,\cdots,y_n)(t) $ denotes the Wronskian determinant of the fundamental solutions. Consequently, the particular solution $ y_p(x) $ can be expressed as