astro-group-bristol / Gradus.jl

Extensible spacetime agnostic general relativistic ray-tracing (GRRT).
https://astro-group-bristol.github.io/Gradus.jl/dev/
GNU General Public License v3.0
18 stars 2 forks source link

Fix: refactor circular orbit calculations #96

Closed fjebaker closed 1 year ago

codecov-commenter commented 1 year ago

Codecov Report

Merging #96 (8055feb) into main (d763e67) will decrease coverage by 1.93%. The diff coverage is 63.25%.

:mega: This organization is not using Codecov’s GitHub App Integration. We recommend you install it so Codecov can continue to function properly for your repositories. Learn more

@@            Coverage Diff             @@
##             main      #96      +/-   ##
==========================================
- Coverage   72.10%   70.17%   -1.93%     
==========================================
  Files          45       49       +4     
  Lines        1606     1794     +188     
==========================================
+ Hits         1158     1259     +101     
- Misses        448      535      +87     
Impacted Files Coverage Δ
src/corona-to-disc/corona-models.jl 50.00% <0.00%> (-27.78%) :arrow_down:
src/line-profiles.jl 81.81% <ø> (+3.55%) :arrow_up:
src/orbits/orbit-interpolations.jl 90.90% <ø> (ø)
src/tracing/tracing.jl 90.19% <ø> (-3.05%) :arrow_down:
...ransfer-functions/cunningham-transfer-functions.jl 97.43% <ø> (ø)
src/corona-to-disc/disc-profiles.jl 58.22% <24.00%> (-7.24%) :arrow_down:
src/image-planes/grids.jl 60.86% <25.00%> (-39.14%) :arrow_down:
src/metrics/kerr-newman-ad.jl 46.25% <46.25%> (ø)
src/transfer-functions/types.jl 50.00% <50.00%> (ø)
src/accretion-geometry/geometry.jl 49.23% <60.00%> (+8.24%) :arrow_up:
... and 13 more

Help us with your feedback. Take ten seconds to tell us how you rate us. Have a feature suggestion? Share it here.

fjebaker commented 1 year ago

I'm trying to solve the circular orbits of charged particles in the Kerr-Newman spacetime analytically, and have made some progress, but am now concerned this method is not possible without using a non-linear solver.

In short, stemming from the method for calculating the Keplerian Frequency, writing the acceleration equation with the Lorentz force ansatz:

$$ \frac{\partial^2 x^\mu}{\partial \lambda^2} = - \Gamma^\mu{\phantom{\mu}\alpha \beta} \frac{\partial x^\alpha}{\partial \lambda} \frac{\partial x^\beta}{\partial \lambda} + e F^{\mu}{\phantom{\mu}\nu} \frac{\partial x^\nu}{\partial \lambda}, $$

which, by the same expansions in the blog post, with

$$ \frac{\partial x^r}{\partial \lambda} = \frac{\partial x^\theta}{\partial \lambda} = 0, $$

can be whittled down to

$$ 0 = \frac{1}{2} \frac{\partial x^\alpha}{\partial \lambda} \frac{\partial x^\beta}{\partial \lambda} \partialr g{\alpha \beta} + g{\sigma r} e F^{\sigma}{\phantom{\sigma}\mu} \frac{\partial x^\mu}{\partial \lambda}. $$

The method I'd been using is then to expand in terms of $t$ and $\phi$ components, and define the Keplerian velocity via

$$ \Omega := \frac{\dot{x}^\phi}{\dot{x}^t}, $$

where the dot denotes differentiation with respect to $\lambda$, and to solve for the Keplerian velocity purely in terms of metric factors -- and now also $F^r{\phantom{r}t}$ and $F^r{\phantom{r}\phi}$. The issue with this is that I made a mistake in thinking dividing by $(\dot{u}^t)^2$ would eliminate all velocity terms except for $\Omega$, however this is not the case, and one ends up with a function in the form

$$ 0 = \mathcal{P}(\Omega) + g{rr} \frac{e}{\dot{x}^t} \left(F^r{\phantom{r}t} + F^r_{\phantom{r}\phi}\Omega \right), $$

where $\mathcal{P}$ is a polynomial in $\Omega$ with coefficients derived from the metric.

One can use the constraint on velocity to find a formulation of $\dot{x}^t$ that depends on $\sqrt{\mathcal{\tilde{P}}(\Omega)}$, but the resulting expression I am unable to rearrange to solve for $\Omega$, once substituting the square root terms.

This may be accomplished numerically, with the full complex form written out and then using a root finder to determine $\Omega$, with some a priori condition to ensure either the prograde or retrograde solution is chosen.

The version implemented in the PR forgot that the $1/\dot{x}^t$ terms appear in the constraining equation, and solved the polynomial with coefficients modified by the Maxwell tensor terms. This, as one might expect, consistently calculates $\Omega$ too small compared to $\Omega$ derived through numerical solving methods (the setup here is arbitrary with $eQ \neq 0$):

orbit-diff-charge

fjebaker commented 1 year ago

One can make the constraint on $\Omega$ a quartic, in which case a general (but horrendously long) solution does indeed exist. I will investigate whether this quartic can be factored or depressed into a cubic to give a more tractable analytic solution.

fjebaker commented 1 year ago

Used a root finder to solve the angular velocity constraint equation for the Kerr-Newman metric. Am currently writing a short blog post about it to make sure it's all working as intended -- will link the post and add tests before merging.