Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Numerical Plot of Particle Orbits in Kerr Equatorial Plane

  1. Mar 23, 2017 #1
    In James Hartle's Gravity, pp. 318-319, Example 15.1, he considers the case of a particle that starts falling from infinity into a Kerr black hole, initially with no kinetic energy (e = 1) and initially moving radially (ℓ = 0). The particle's motion is constrained along the equatorial plane. He says we can integrate:
    [tex]\frac{d\phi}{dr} = \frac{d\phi/d\tau}{dr/d\tau} = -\frac{2Ma}{r\Delta} \Bigg[ \frac{2M}{r} \bigg( 1 - \frac{a^2}{r^2}\bigg)\Bigg] ^{-1/2} \tag{C-1}[/tex]...to get the orbit's shape ∅(r).

    (M is the mass of the Kerr Black Hole, a is its Angular Momentum per Mass M. e is the particle's energy-per-unit-rest-mass-at-infinity and ℓ is the particle's angular-momentum-per-unit-rest-mass-at-infinity. r is the coordinate radius of the Kerr metric and Δ = r^2 - 2Mr + a^2)

    1. The problem statement, all variables and given/known data
    I want to generalize this problem for orbits of a particle with e and ℓ values that I can freely set. Although I'm generalizing the problem for different values of e and ℓ, I still only want to work with particles constrained in the Kerr equatorial plane. I've worked out the equations for this generalization, and I've placed them under the Relevant equations of this template.

    The main thing I actually want to accomplish by generalizing this problem is to find out how I can obtain a particle with a spiral-like orbit by altering any of the Kerr Black hole parameters (M and a) or any of the particle parameters (e and ℓ).

    2. Relevant equations
    The following equations are derived from the geodesic equation for the Kerr metric, taken from Hartle's Gravity, page 318.
    \frac{d\phi}{d\tau} = \frac{1}{r^2 - 2Mr + a^2} \cdot \Bigg[ \bigg(1 - \frac{2M}{r} \bigg) ℓ + \bigg( \frac{2Ma}{r} \bigg) e \Bigg] \tag{1}
    [/tex] [tex]
    \frac{e^2 - 1}{2} = \frac{1}{2} \bigg( \frac{dr}{d\tau} \bigg)^2 - \frac{M}{r} + \frac{ℓ^2 - a^2 (e^2 - 1)}{2r^2} - \frac{M(ℓ-ae)^2}{r^3} \tag{2}
    [/tex]From these equations, I derived:
    \frac{dr}{d\tau} = -\sqrt{e^2 - 1 + \frac{2M}{r} - \frac{ℓ^2 - a^2 (e^2 - 1)}{r^2} + \frac{2M(ℓ-ae)^2}{r^3}} \tag{3}
    [/tex]Transforming the derivatives into finite differences, we get:
    \phi_{i+1} - \phi_i = \Delta \phi = \Delta \tau \cdot \Bigg[ \frac{d\phi}{d\tau} \Bigg] \tag{4}
    r_{i+1} - r_{i} = \Delta r = \Delta \tau \cdot \Bigg[ \frac{dr}{d\tau}\Bigg] \tag{5}[/tex]

    3. The attempt at a solution
    I gave several variables initial values:
    Code (Text):
    tau(1,1) = 0;
    rOFtau(1,1) = 30*M;
    phi(1,1) = 0;
    Dr = 1;
    I substituted Equations 1 and 3 into Equations 4 and 5 (respectively) and handled the integrations numerically as shown in the following code:
    Code (Text):
    while rOFtau(1,i) > 2*M && imag(rOFtau(1,i)) == 0 && abs(Dr) > 0

    % Iterate over Tau to get a set of Tau for plotting
        tau(1,i+1) = tau(1,i) + Dtau;

    % Iterate numerical integration over r to get a set of r for plotting
        Dr = Dtau * -(       e^2 - 1 ...
                            + 2*M/rOFtau(1,i) ...
                            - (l^2 - (a^2)*(e^2 - 1))/(rOFtau(1,i)^2) ...
                            + (2*M*(l-a*e)^2)/(rOFtau(1,i)^3) ...
        rOFtau(1,i+1) = rOFtau(1,i) + Dr;

    % Iterate numerical integration over Phi to get a set of Phi for plotting
        phi(1,i+1) = phi(1,i) + Dtau * ( (rOFtau(1,i)^2 - 2*M*rOFtau(1,i) + a^2)^(-1) * ...
                                               (    (1-(2*M/rOFtau(1,i)))*l + ...
                                                     (2*M*a/rOFtau(1,i))*e ...
        i = i + 1;
    I simply plot r against ∅ on a polar graph to get the shape of the orbit.

    Now, for the part where I start varying the parameters, I first tried:
    Code (Text):
    % Black Hole Parameters
    M = 20;
    a = 15;

    % Particle Parameters
    e = 1;
    l = 20;
    ...and I found that the orbit was too much of a plunge. (Image: https://scontent.fmnl4-2.fna.fbcdn....=086e26b46f326bb0509950e3f92842e6&oe=58D6F490)

    So at first, I thought acquiring a spiral-like orbit was a matter of increasing ℓ, so I increased it to ℓ = 100, and the program prompted me that it ignored the imaginary parts of the function. The orbit's angle ∅ varied more, but the loop terminated early because I let it stop if an imaginary part of r was found: (Image: https://scontent.fmnl4-2.fna.fbcdn....=82266e02168f425ad10ade24767114e5&oe=58D60994)

    I thought that this could be fixed by decreasing e (indicating an initial kinetic energy of the particle), but decreasing it to e = 0.9 gives r an imaginary part immediately after the first iteration. I've also tried altering a and M but the orbit barely changes.

    MATLAB is the program I use for coding.

    Any suggestions on how to proceed would be greatly appreciated. Or perhaps I misinterpreted a concept or two? Any insights or corrections would be appreciated as well. Thanks!
    Last edited: Mar 23, 2017
  2. jcsd
  3. Mar 24, 2017 #2
    Did you consider what relations between r and φ would give you a spiral and going from there?
  4. Mar 24, 2017 #3
    Thank you so much for the response, sir! I considered it when I was brainstorming for ways on how I could approach the problem. But I didn't know how I could mathematically relate r and φ to produce a general spiral orbit WHILE satisfying Equations 1 and 3.

    For example, off the top of my head, the first few spirals I could come up with were the Archimedes spiral (r = φ), and this other spiral from my calculus classes, r = sqrt(φ). But it was obvious that these relations would not satisfy Equations 1 and 3 for all r. So, I thought, even if I could come up with a general mathematical relation that pertains to a spiral, it would be useless if it doesn't satisfy the geodesic equation, thus, I abandoned the approach.

    Is there something I missed that could have made the approach fruitful?
  5. Mar 24, 2017 #4
    Not that's obvious to me. But it might be worth trying to narrow the range of possibilities first. Are you sure that spiral orbits are possible? Are there stable closed orbits in the Kerr metric? If so, it might be worth starting from one of those and trying to perturb it. ( Perhaps replace r by r.(1+ ) were k is small and see what happens.) Or consider generally r = f.φ, where for example f is a small-valued, slowly varying function. and see if you can put any constraints on f so as to satisfy the equations.
  6. Mar 25, 2017 #5
    Sir, I'm not actually sure. The only inspiration I had for imagining spiral orbits is the case of two slowly in-spiraling bodies. But now that I think about it, in that case, both bodies have a comparable gravitational effect on each other, while in this problem, the particle is assumed to have no curvature, and so it is only the black hole that affects the particle and not vice versa. This actually might explain why my first result was more of a plunge. Is it actually impossible?

    Okay, sir. Here goes: I start by assuming that f is a function of r and not φ because rotating the system a few degrees dφ mustn't yield a change in the solution because the spacetime is axisymmetric for φ. (This is justified because neither Equation 1 nor 3 depend on φ.) Then:
    [tex]r = \phi f(r) \tag{6}[/tex][tex]\partial r = \phi \partial f + f \partial\phi [/tex][tex]1 = \phi \frac{\partial f}{\partial r} + f \frac{\partial \phi}{\partial r}\tag{7}[/tex]But since both f and φ are only dependent on r:
    [tex]1 = \phi \frac{df}{dr} + f \frac{d\phi}{dr}\tag{8}[/tex][tex]\frac{d\phi/d\tau}{dr/d\tau} = \frac{1 - \phi \frac{df}{dr}}{f}\tag{9}[/tex] For f to be small enough and slowly varying, I suppose: [tex]\frac{df}{dr}<<1\tag{10a}[/tex][tex]f<<1\tag{10b}[/tex]I'm thinking I can substitute Equations 1 and 3 into Equation 9, and check if Equations 10 can be satisfied by selecting several parameters. But right now, I don't have much of an idea on how to execute this, exactly.

    Is this along the lines of what you were suggesting, sir? Should I alter Equation 10, or is it general enough? Please feel free to correct any mistakes I might have made. (This goes for anyone who might happen upon this thread as well. Please feel free to pitch in! :smile:) I'll try brainstorming on ways to continue this solution, but suggestions on how to proceed would be greatly appreciated as well. Thank you again for your response!
  7. Mar 25, 2017 #6
    I don't know.

    But it occurs to me that conservation of energy might be a problem. The only cases of spiral orbits I can think of immediately involve some change of orbital energy--drag from viscous forces, emission of gravitational waves or absorption of electromagnetic energy as in a cyclotron.
  8. Apr 1, 2017 #7
    Continuing the solution:

    I think the simplest way to solve it is to rearrange Equation 9 to get: [tex]
    \frac{f}{1 - \phi \frac{df}{dr}} = \frac{dr/d\tau}{d\phi/d\tau} \tag{11}
    [/tex] Since df/dr is very small,[tex]
    \frac{f}{1 - \phi \frac{df}{dr}} ≈ f \tag{12}
    [/tex]Then, substituting Equations 12 and 10b to Equation 11, we get:
    \frac{dr/d\tau}{d\phi/d\tau} ≈ f[/tex][tex]
    \frac{dr/d\tau}{d\phi/d\tau} << 1[/tex][tex]
    \frac{dr}{d\tau} << \frac{d\phi}{d\tau}\tag{13}[/tex]I suppose Equation 13 makes sense, because if I want a spiral orbit, the coordinate radius would have to change much more slowly than the angle.

    I checked whether or not Equation 13 can actually be satisfied by Equations 1 and 3:

    With careful selection of e and ℓ, I found that d∅/dτ can be greater than dr/dτ for some large values of r. (Example below)
    Code (Text):
    %% Sample values that satisfy dr/dτ < d∅/dτ
    M = 10;
    a = 5;

    r = 100000;
    l = 66666.67;
    e = 1.201730240307;
    But after letting r and ∅ evolve over time using Equations 1 and 3, r eventually decreased to a value for which Equation 3 became imaginary, and thus the entire solution became physically irrelevant, no matter how careful I was in selecting e and ℓ. After further analysis, it seemed to me that the only way that dr/dτ can remain real, while remaining less than d∅/dτ, is if ℓ and e are changing as well.

    And I believe this statement means that e and ℓ cannot remain constant in order for a spiral orbit to exist, which I think agrees with one of Sir John Park's previous statements:
    It would be greatly appreciated if someone could point out any mistakes I might have made in my solution, but right now, I'm almost convinced that a spiral-like particle orbit cannot exist in a Kerr spacetime's equatorial plane, and I hope to mark this thread as solved soon. Thanks!
  9. Jul 9, 2017 #8
    An infalling particle should spiral in like this:

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?