Circular orbits of charged particles

In the last couple of days I've been trying to add the Kerr-Newman metric to Gradus.jl. The Kerr-Newman metric describes a black hole spacetime with mass MM, spin aa, and charge QQ. The metric may be expressed in Boyer-Lindquist coordinates as

ds2=Δa2+sin2θΣdt2+ΣΔdr2+Σdθ2+sin2θΣ[(r2+a2)2a2Δsin2θ]dϕ22asin2θΣ(r2+a2Δ)dtdϕ, \text{d} s^2 = -\frac{\Delta - a^2 + \sin^2 \theta}{\Sigma} \text{d} t^2 + \frac{\Sigma}{\Delta} \text{d} r^2 + \Sigma \text{d} \theta^2 + \frac{\sin^2 \theta}{\Sigma} \left[\left(r^2 + a^2\right)^2 - a^2 \Delta \sin^2 \theta \right] \text{d} \phi^2 - 2 \frac{a \sin^2 \theta}{\Sigma} \left(r^2 + a^2 - \Delta \right) \text{d} t \text{d} \phi,


Σ=r2+a2cos2θ,andΔ=r22Mr+a2+Q2. \Sigma = r^2 + a^2 \cos^2 \theta, \quad \text{and} \quad \Delta = r^2 - 2 M r + a^2 + Q^2.

Implementing the metric itself is straight forward, allowing for a new class of geodesic to be traced, namely those with non-zero charge per unit mass qq. These are the particles which interact electromagnetically with the central singularity. We consider only time-like geodesics with charge.

These have been discussed in a number of papers, for example Schroven, Hackmann and Lämmerzahl (2017)[1], which studies in detail the innermost stable circular orbit (ISCO) for different charge configuration, and the effect this has on accreting matter of negative (electrons) and positive charge (protons/ions). There is also the earlier paper by Hackmann and Xu (2013)[2], which explores more generally the orbits of charged particles in the Kerr-Newman spacetime and classifies them to describe their behaviour. These works, and references therein, use the Hamilton-Jacobi formalism with the Carter constant to calculate an effective potential on the four-velocity components, and use this to infer stable (circular) orbits and the ISCO. I would instead like to study this from a step back, directly from the 2nd order geodesic equation, using modern numerical methods. This means adding an Ansatz for the Lorentz force

d2dλ2xμ=Γμαβμx˙αx˙β+qFμαμx˙α, \frac{\text{d}^2}{\text{d} \lambda^2} x^\mu = - \Gamma^{ \mu}_{\phantom{ \mu} \alpha\beta} \dot{x}^\alpha \dot{x}^\beta + q F^{ \mu}_{\phantom{ \mu} \alpha} \dot{x}^\alpha,

where Fμαμ F^{ \mu}_{\phantom{ \mu} \alpha} is the Faraday tensor, and the dot denotes differentiation with respect to affine time λ\lambda. The second-order ordinary differential equation may be solved with some initial position xμx^\mu and 3-velocity x˙\vec{\dot{x}}, as x˙t\dot{x}^t is constrained via

gαβx˙αx˙β=μ2 g_{\alpha \beta} \dot{x}^\alpha \dot{x}^\beta = -\mu^2

where μ\mu is the invariant mass.

In this post, I would like to study the class of circular charged orbits in the equatorial plane, beginning from eq. (3), to determine semi-analytic solutions for the angular momentum and energy of these orbits.

Faraday tensor

Starting from the four-potential[1]

Aμ=QrΣ(1dtasin2θdϕ), A_\mu = \frac{Q r}{\Sigma} \left(1 \text{d} t - a \sin^2 \theta \text{d} \phi \right),

the Faraday tensor is given by

Fμν=μAννAμ. F_{\mu\nu} = \partial_\mu A_\nu - \partial_\nu A_\mu.

This means there are only 4 unique non-zero components of the anti-symmetric tensor, namely Frt,Fθt,Frϕ,FrθF_{rt}, F_{\theta t}, F_{r\phi}, F_{r\theta}, with the other 4 being their anti-symmetric counterparts.

It is worth noting that with this potential, one must only consider the product of the particle and black hole charge qQqQ in the equations of motion, since there are no terms that depend solely on qq or QQ. There are therefore 2 classes of solution we can expect depending on the sign of qQqQ. Similarly in the metric, only terms with Q2Q^2 appear and consequently the parity is lost. We therefore restrict ourselves to Q0Q \geq 0 and allowing qq to take both negative and positive values without loss of generality.

Keplerian angular velocity

Following the method of my previous blog post, we will rework eq. (3) in the equatorial plane, subject to x˙r=x˙θ=0\dot{x}^r = \dot{x}^\theta = 0, and calculate the Keplerian angular velocity, defined

Ωϕ:=x˙ϕx˙t. \Omega_\phi := \frac{\dot{x}^\phi}{\dot{x}^t}.

The energy and angular momentum already have known solutions in terms of Ωϕ\Omega_\phi and metric terms. After expanding the Christoffel symbols and working through a little algebra (as in the previous post), one obtains

0=12(rgtt(x˙t)2+2rgtϕ x˙t x˙ϕ+rgϕϕ(x˙ϕ)2)+qgrr(Frtrx˙t+Frϕrx˙ϕ), 0 = \frac{1}{2}\left( \frac{\partial}{\partial r} g_{tt} (\dot{x}^t)^2 + 2 \frac{\partial}{\partial r} g_{t\phi} \, \dot{x}^t \dot{x}^\phi +\frac{\partial}{\partial r} g_{\phi \phi} (\dot{x}^\phi)^2 \right) + q g_{rr} \left( F^{ r}_{\phantom{ r} t} \dot{x}^t + F^{ r}_{\phantom{ r} \phi} \dot{x}^\phi \right),

and dividing by (x˙t)2(\dot{x}^t)^2, one obtains

0=12(rgtt+2rgtϕΩϕ+rgϕϕΩϕ2)+qgrrx˙t(Frtr+FrϕrΩϕ), 0 = \frac{1}{2}\left( \frac{\partial}{\partial r} g_{tt} + 2 \frac{\partial}{\partial r} g_{t\phi} \Omega_\phi +\frac{\partial}{\partial r} g_{\phi \phi} \Omega_\phi^2 \right) + q \frac{g_{rr}}{\dot{x}^t} \left( F^{ r}_{\phantom{ r} t} + F^{ r}_{\phantom{ r} \phi} \Omega_\phi \right),

which bears a factor 1/x˙t1/\dot{x}^t on the charge terms. The first term in this equation is a simple quadratic, but with these charge terms added, we require some way of eliminating or determining x˙t\dot{x}^t in order to solve for Ωϕ\Omega_\phi.

The only thing I can think to do here is to use eq. (4), and write

1x˙t=1μ(gtt+2gtϕΩϕ+gϕϕΩt2). \frac{1}{\dot{x}^t} = \frac{1}{\mu} \sqrt{- \left(g_{tt} + 2g_{t\phi} \Omega_\phi + g_{\phi \phi} \Omega_t^2 \right)}.

One could then balance eq. (9) to have one term either side of the equals sign, substitute the 1/x˙t1/\dot{x}^t, and square to obtain a quartic equation. Doing so also loses information about the sign of qq, but that is something we could add in later. Although this is in theory then analytically solvable, in reality the expression is horrendous, and I cannot find a nice simple way of reducing the quartic into e.g. the product of two quadratics, or depressing it to simplify the expression. There are no doubt other ways to tackle such an equation, but I am not versed in polynomial analysis and would have a hard time approaching the problem.

For now, I'd instead prefer a semi-analytic approach, and use a root finder to solve for the solutions. Since the overall form of the equation is leading order quadratic, and we can a priori select x˙t\dot{x}^t that goes forwards in time, we may therefore anticipate classes of solutions for Ωϕ\Omega_\phi, corresponding to product set of prograde and retrograde orbits with positively and negatively charged particles.

Implementing a short Julia function:

using Gradus

function functor_Ω(
    q = 0.0,
    μ = 1.0
    g, jacs = Gradus.metric_jacobian(m, rθ)
    # only want the derivatives w.r.t. r
    ∂rg = jacs[:, 1]

    x = SVector(0, rθ[1], rθ[2], 0)
    F = Gradus.maxwell_tensor(m, x)

    function f(ω)
        Δ = (ω^2 * ∂rg[4] + 2 * ω * ∂rg[5] + ∂rg[1])
        arg = -(ω^2 * g[4] + 2 * ω * g[5] + g[1]) / μ^2
        # check to ensure we can preserve forwards-in-time
        ut = if arg >= 0
            # return something distinct
            return NaN
        # combine with charge terms
        0.5 * Δ + (F[2, 4] * ω + F[2, 1]) * g[2] * q * ut

This function limits the domain of Ωϕ\Omega_\phi to solutions with 1ut0\frac{1}{u^t} \geq 0. We then set up some initial configuration:

m = KerrNewmanMetric(M = 1.0, a = 0.0, Q = 0.9)
rθ = SVector(3.0, π/2)

I'd like to note that this picture is unphysical. If a black hole were to have charge, it is likely that the charge is vanishingly small in the metric[1], and that the charge product qQqQ is instead non-zero. But there is nothing in the mathematics that prevents us from exploring this class of spacetimes, so we will press on.

Plotting the value of the constraint equation for different Ωϕ\Omega_\phi for two different spins:

Discontinuities are due to limited resolution (NΩϕ=300N_{\Omega_\phi}=300). The two solutions, prograde and retrograde, are clearly defined. When q>0.0|q| > 0.0, the possible range of Ωϕ\Omega_\phi is limited by the positivity constraint in eq. (10).

We can clearly identify the values of Ωϕ\Omega_\phi corresponding to each of the q={1,0,1}q = \{-1, 0, 1\} cases. Ambiguity is introduced if we move closer:

Examining the two cases again:

I am unsure what the interpretation here is: one could interpret that the direction of time has changed, as in eq. (10), there is no difference between qqq \rightarrow -q and ututu^t \rightarrow -u^t. Thus the transformation (q,ut)(q,ut)(q, u^t) \rightarrow (q, -u^t) is equivalent to the (q,ut)(-q, u^t) case, which puts negative charge orbits on the curve corresponding to positive charge orbits when the direction of time is changed.

We can analytically continue the curves to represent this with the modification to our function:

ut = sign(arg) * sqrt(abs(arg))

I have coloured the line corresponding to q=1q = 1 for a=0.4a = 0.4 orange to emphasize how the continuation is realized. I have chosen to adopt this in the Gradus.jl implementation. It also helps with the root finder for Ωϕ\Omega_\phi to have the function defined all along the domain (~;

Energy and ISCO

Here we will use the method outlined previously for determining the ISCO, using the same minima finding method in the energy expression. The energy and angular momentum are from the canonical four-momenta, which is modified under the influence of the electromagnetic potential, as described in Tursunov, Zdeněk, and Kološ (2016)[3] for an external magnetic field, and Carter (1968)[4] for the general class of Kerr metrics.

Quickly re-deriving the canonical momenta from the Lagrangian

L=12gμνx˙μx˙νqAσx˙σ, \mathcal{L} = \frac{1}{2} g_{\mu\nu} \dot{x}^\mu \dot{x}^\nu - q A_\sigma \dot{x}^\sigma,

with the momenta defined

πν:=Lx˙ν, \pi_\nu := \frac{\partial \mathcal{L}}{\partial \dot{x}^\nu},

and consequently

πν=gνμx˙μqAν=x˙νqAν. \pi_\nu = g_{\nu\mu}\dot{x}^\mu - q A_\nu = \dot{x}_\nu - q A_\nu.

There is here an ambiguity in the sign of AνA_\nu as it enters in the above equations. I am using the same convention as Hackmann et al. here, which is in contrast to e.g. Carter. I must admit I do not quite understand the origin of this discrepancy, but the differing results are dramatic.

We identify as usual E=πtE = -\pi_t and Lz=πϕL_z = \pi_\phi as constants of motion, and

E=(x˙tqAt),andLz=(x˙ϕqAϕ). \therefore \quad E = -(\dot{x}_t - q A_t), \quad \text{and} \quad L_z = (\dot{x}_\phi - q A_\phi).

We find x˙t\dot{x}_t and x˙ϕ\dot{x}_\phi using Ωϕ\Omega_\phi, as in the previous blog post, and can then plot the energy as a function of radius:

For the case of q=1q=1, we have E1E \approx 1 for the majority of orbits. The marginally bound orbit is defined as the circular orbit with E=1E = 1, so although these orbits are stable, they are not very tightly bound relative to the other setups. This has implications for a maximal charge product qQqQ for which circular orbits may be found. Numerical limitations furthermore prevent fully describing the energetic system at small radii, as can be seen by blatant discontinuities in the curves, where the root finder has struggled to find solutions for Ωϕ\Omega_\phi.

We can find the ISCO by determining the minima of the energy of the circular orbits, i.e.

Err=rISCO=0. \left.\frac{\partial E}{\partial r}\right\rvert_{r=r_\text{ISCO}} = 0.

Plotting the ISCO as a function of charge product qQqQ:

At qQ=1qQ = 1, both curves are divergent and rISCOr_\text{ISCO} \rightarrow \infty. For negative charge products the ISCO increases – this is equivalent to an additional attractive force, or an effective increase in gravitational potential, pulling on the particle. Orbits that were once stable are no longer, and hence rISCOr_\text{ISCO} increases.

For positive charge products, the ISCO initially decreases, reaching a minima, before increasing rapidly. The initial decrease may be thought of as the repulsion holding up the orbit and weakening the gravitational potential, allowing previously unstable orbits to become stable. The subsequent increase is interpreted as the EM repulsion overpowering gravity, suppressing the gravitational term in eq. (3).

Orbit finding approach

Using the method I previously outlined in another post, Gradus.jl can automatically find stable circular orbits in the equatorial plane by minimizing some stability measure. It requires only the modified geodesic equation eq. (3), and solves for some vϕv^\phi which minimizes the variations in Δxr\Delta x^r over the course of the integration. It is therefore a good numerical test of the semi-analytic approach to ensure the model is at least self consistent.

Plotting LzL_z against EE is a method for analyzing the stability and classifying different orbits. For example, the ISCO is a minimum of both EE and LzL_z, and therefore on the LzEL_z - E plane corresponds to a cusp in the curve representing circular orbits.

Plotting the analytic method (dashed) and the points determined by the circular orbit solver (crosses, connected by solid line):

Agreement is seen between the methods. This can also be used to more accurately determine the ISCO for the circular orbit solver, as sometimes the integrator is "too good" and classifies some unstable / marginally bound orbits as stable circular. These plots clearly identify the location of the ISCO energetically.

Line profiles

The effects of GR on the redshift is calculated using the ratio of energies

g:=kνuνobskμuμdisc, g := \frac{\left. k^\nu u_{\nu} \right\rvert_{\text{obs}}}{\left. k^\mu u_\mu \right\rvert_\text{disc}},

where we assume kνuνobs=1\left. k^\nu u_{\nu} \right\rvert_{\text{obs}} = -1 for simplicity, as the observer is assumed to be far enough from the singularity to be in flat space. Our calculation therefore simplifies to

using Tullio

function redshift_charged(q)
    function _f(m, gp, t)
        x = gp.u2
        g = Gradus.metric(m, x)
        u = CircularOrbits.fourvelocity(m, x[2]; q = q)
        @tullio o := -g[i, j] * gp.v2[i] * u[j]

Here it is already implemented in a form that returns a function that we can use to create a PointFunction to use with the integrator. We can then calculate moderately high resolution line profiles with:

# set charge and build point function
q = 1.0
redshift_pf = PointFunction(redshift_charged(q))
# position and disc configuration
u = SVector(0.0, 1000.0, deg2rad(40), 0.0)
d = GeometricThinDisc(0.0, 1000.0, π/2)
# set up binning and calculate over fine grid
bins = collect(range(0.1, 1.3, 200))
g, f = @time lineprofile(
    bins = bins,
    algorithm = BinnedLineProfile(),
    plane = PolarPlane(GeometricGrid(); Nr = 700, Nθ = 1300, r_max = 50.0),
    verbose = true,
    redshift_pf = predshift_pf,
    minrₑ = Gradus.isco(m, q = q),

Plotting these:

The effect of the energetics of the charged orbits is clear in the line profiles. For the case of q=1q=1, the circular orbit energy is pretty much constant for all stable orbits. This is particularly important around the ISCO radius, as that is where we can expect the maximal redshift. Consequently, the redshift of the line profile is much narrower. For the case of q=1q = -1 we see the opposite effect, where the energy of the stable orbits varies to a greater extent, and therefore the line profile is significantly broadened.

The effect is somewhat similar to changing the inclination of the observer:

The above is for the Schwarzschild case only.

A "more realistic" scenario anticipates QQ to be very small, but could still allow qQ1qQ \sim 1. We can approximate this by setting e.g. Q1016Q \approx 10^{-16} and q1016q \approx 10^{16}, so that QQ does not affect the spacetime, and enters the model only through the Lorentz force in the geodesic equation.

The Schwarzschild and qQ=0qQ = 0 case correspond exactly, but the effects of the charge are still very apparent.


[1] Schroven, Kris, Eva Hackmann, and Claus Lämmerzahl. ‘Relativistic Dust Accretion of Charged Particles in Kerr-Newman Spacetime’. Physical Review D 96, no. 6 (26 September 2017): 063015.
[2] Hackmann, Eva, and Hongxiao Xu. ‘Charged Particle Motion in Kerr-Newmann Space-Times’. Physical Review D 87, no. 12 (24 June 2013): 124030.
[3] Tursunov, Arman, Zdeněk Stuchlík, and Martin Kološ. ‘Circular Orbits and Related Quasiharmonic Oscillatory Motion of Charged Particles around Weakly Magnetized Rotating Black Holes’. Physical Review D 93, no. 8 (7 April 2016): 084012.
[4] B Carter (1968). Global structure of the Kerr family of gravitational fields. PHYS REV 174:5.