# I Schwarzschild equation of motion: initial conditions

Tags:
1. May 28, 2016

### Yukterez

The Schwarzschild equation of motion, where coordinate length is differentiated by proper time is as far as I know

$$r''(t) = -\frac{G\cdot M}{r(t)^2} + r(t)\cdot{\theta}'(t)^2 - \frac{3\cdot G\cdot M\cdot{\theta}'(t)^2}{c^2}$$
$${\theta}''(t) = -2\cdot r'(t)/r(t)\cdot{\theta}'(t)$$

Now the question is, how do I choose my initial conditions for the angular and the radial initial velocities, for example, if I want to test the orbit of a test particle close to the speed of light $v_0=0.999c$ at the photon sphere $r_0=1.5r_s$?

I am not sure if my approach is correct, but when I try

$$\dot{\theta}_0\cdot r_0 = \frac{\xi_0}{ \color{red}{\sqrt{ 1-v_0^2/c^2}}\cdot \color{blue}{\sqrt{1-r_s/r_0}}}$$

with $\xi_0$ as the transversal component of $v_0$ and $\kappa_0$ as the radial component:

$$\dot{r}_0\ = \frac{\kappa_0\cdot \color{blue}{\sqrt{1-r_s/r_0}}}{ \color{red}{\sqrt{ 1-v_0^2/c^2}}}$$

so that $\xi = cos(\varphi_0), \kappa = sin(\varphi_0)$ with $\varphi$ as the launching angle from the shell, the orbits look like I would expect them, at least in the transversal direction:

But since radial and transverse directions have the factor $\color{blue}{\sqrt{1-r_s/r_0}}$ on different sides of the fraction, the inital angle is not conserved.

So are my initial conditions even correct, and if not, how would one transform the radial and the transversal components of the initial velocity $v_0$ to $\dot{r}(0)$ and $\dot{\varphi}(0)$?
_____________________________________________________

If it is of interest, the code for the plot is in Mathematica Syntax
Code (Text):
G = 1; M = 1; c = 1; rs = 2 G M/c^2; T = tMax G M/c^3;

r0 = 151/100 rs; vo = 999/1000 c; φ = Pi/2; θ0 = 0; tMax = 2;
vr0 = vo Cos[φ]/j*k; vθ0 = vo/r0 Sin[φ]/j/k;
d1 = T/10; d2 = d1; wp = 30; step = T/200;

j = Sqrt[1 - vo^2/c^2];
k = Sqrt[1 - rs/r0];

sol = NDSolve[{
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0] == vr0,
r[0] == r0,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0] == vθ0,
θ[0] == θ0,
τ'[t] == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0] == 0
}, {r, θ, τ}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

t[ξ_] :=
Quiet[χ /.
FindRoot[
Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ := Quiet[t[ι]];

x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]

s[text_] := Style[text, FontSize -> font];  font = 11;

(* Eigenzeit, Celerity *) Do[Print[
Rasterize[Grid[{{Show[Graphics[{
{Black, Circle[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 400, PlotRange -> 2 r0],
Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],
ParametricPlot[{x[η], y[η]}, {η, 0, т},
ColorFunction -> Function[{x, y, η},
Hue[0.85, 1, 0.5, Max[Min[(-т + (η + d1))/d1, 1], 0]]],
ColorFunctionScaling -> False],
ParametricPlot[{x[η], y[η]}, {η, 0, т},
ColorFunction -> Function[{x, y, η},
Hue[0, 1, 0.5, Max[Min[(-т + (η + d2))/d2, 1], 0]]],
ColorFunctionScaling -> False]]},
{Grid[{
{"   ", s["proper time"], " = ", s[N[т, 8]], s["GM/c³"]},
{"   ", s["coordinate time"], " = ", s[N[Evaluate[τ[т] /. sol][[1]], 8]], s["GM/c³"]},
{"   ", s["time dilation"], " = ", s[N[Evaluate[τ'[т] /. sol][[1]], 8]], s["dτ/dt"]},
{"   ", s["angle"], " = ", s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]], s["grad"]},
{"   ", s["radial distance"], " = ", s[N[Evaluate[r[т] /. sol][[1]], 8]], s["GM/c²"]},
{"   ", s["x-Axis"], " = ", s[N[x[т], 8]], s["GM/c²"]},
{"   ", s["y-Axis"], " = ", s[N[y[т], 8]], s["GM/c²"]}
}, Alignment -> Left]}}, Alignment -> Left]]
], {т, step, T, step}]

2. May 29, 2016

### m4r35n357

Finding initial conditions for these simulations is not trivial; I have devoted whole programs (some longer than the simulation itself) just to deriving useful sets. Coordinates and velocities are not particularly intuitive or useful for this purpose. For the Schwarzschild case (also Kerr, but that's for another day) you can use the constants of motion E and L (see equation 7.49 in Carroll for a circular orbit) and translate these into velocities when you initialize your solution.

I have provided links here many times to source code that does all this so I won't repeat it this time, but search for my posts regarding orbits if you are interested, it will save you a lot of time. I can probably help if you have a more specific question.

3. May 29, 2016

### m4r35n357

On further reflection it occurs that I might be able to help more with this, if you want to test out an orbit like the one you described. Let c = G = M = $\mu$ = 1. Then you should be able to get an orbit close to the photon sphere if you set:

a = 0.0
r = 3.001
E = 18.272630817366
L = 94.899952581663

You should get an unstable orbit where:

v = 0.999833 ($\gamma$ = 54.781)

Hope that helps (I set this up on the JavaScript simulator, you can use it as an independent check for your stuff - or let me know if I've got something wrong!).

Last edited: May 29, 2016
4. May 29, 2016

### Yukterez

I came to the conclusion that the right way to transform the v₀ into the initial conditions for the equation of state, r'(0) and θ'(0), must be

$$\dot{\theta}(0)\cdot r(0) = \frac{{v\perp}_0}{ \color{red}{\sqrt{ 1-v_0^2/c^2}}}$$
$$\dot{r}(0) = \frac{{v\parallel}_0\cdot \color{blue}{\sqrt{1-r_s/r_0}}}{ \color{red}{\sqrt{ 1-v_0^2/c^2}}}$$

When I convert the initial velocity components by that factors I get the same time dilation with the differential

$${\tau}'(t) = \frac{\sqrt{c^2\cdot r(t)-c^2\cdot r_s +r(t)\cdot r'(t)^2- r_s\cdot r(t)^2\cdot \theta '(t)^2+r(t)^3\cdot \theta '(t)^2}}{c\cdot \sqrt{r(t)- r_s } \cdot\sqrt{1-\frac{ r_s }{r(t)}}}$$

as I get it with

$$\frac{1}{\sqrt{1-v_0^2/c^2}} \cdot \frac{1}{\sqrt{1-r_s/r_0}}$$

see the test calculation where Out[7] matches Out[8]:

5. May 29, 2016

### Yukterez

How did you obtain that value? When I use the formula √(GM/r0)/√(1-rs/r0) I get v0 = 0.9995

6. May 30, 2016

### George Jones

Staff Emeritus
I just logged on to post that I do not get

When I logged on, I saw

This is what I get.

7. May 30, 2016

### Yukterez

Very fine, then the first problem is solved, and only one question remains.

Because the blue gravitational factor is only in the radial direction it seems that angles are not preserved, for example, if someone throws a ball with 45° locally, the initial angle seems flatter for a coordinate observer at infinity (that's what I suspect)? Or do I have to split the components in a way that the angle is preserved (I don't think so)?

My guess is that if I want to switch from proper circumference to proper radial distance (since the factor is not 2π anymore in Schwarzschild coordinates) I'd have to do something like multiply the old r coordinate with rs+∫(1/√(1-rs/R), R=rs..r) on a Plot?

Last edited: May 30, 2016
8. May 30, 2016

### m4r35n357

I don't see how the "potential" (blue) factor can produce a speed. I use the "kinetic" (red) factor:

$\gamma = \sqrt { \frac {1} {1 - v^2}}$, so that

$v = \sqrt {1 - \frac {1 } { \gamma^2} }$ (by $v$ I really mean $v / c$!)

where $\gamma = \frac {dt} {d\tau}$ from the simulation to give me that value of $v$.

Caveat:

We are using different methods to simulate and I must admit I don't fully understand your equations of motion. At first I thought they were the geodesic equations but I don't recognize the Christoffel symbols so now I'm not so sure.

Anyway, I think the $t$ in your equations is what I would refer to as $\tau$, and $\theta$ is what I would refer to as $\phi$. You would need to solve for $t$ (using another geodesic equation) and extract the intermediate $\frac {dt} {d\tau}$ in order to use my formula.

Are any of my assumptions mistaken?

Last edited: May 30, 2016
9. May 30, 2016

### George Jones

Staff Emeritus
This is what I think, too. Maybe Yukterez has a math background, Mathematicians and physicists tend to interchange the roles of $\phi$ and $\theta$ in spherical coordinates.

You're using $\gamma$ with respect to infinity, while Yukterez is using a $\gamma$ generated by a local physical velocity.

Suppose I stand on platform that uses a powerful rocket to keep it hovering constant $r_0$ and $\theta$ (i.e., $\phi$ for physics books). I then throw a ball that has initial speed $v_0$ (and that makes an initial angle $\varphi$ with respect to the platform). Then,

$$\gamma = \frac{1}{\sqrt{1 - v_0^2}} = \frac{dt}{d\tau} \sqrt{1 - \frac{2GM}{c^2 r_0}}.$$

When I get to work at about 9 Pacific time, I am going to have a look at Yukterez's equations of motion.

Last edited: May 30, 2016
10. May 30, 2016

### Yukterez

Here is the updated code (Mathematica Syntax), where θ is the longitude (the plane can always be rotated so that the lattitude is fixed), and φ the launching angle from the observers shell:

Code (Text):
G = 1; M = 1; c = 1; rs = 2 G M/c^2;

r0 = 151/100 rs; vo = 999/1000 c; φ = 0; θ0 = 0; tMax = 2; T = tMax G M/c^3;
vr0 = vo Sin[φ]/j*k; vθ0 = vo/r0 Cos[φ]/j;
d1 = T/10; d2 = d1; wp = 30; step = T/200;

j = Sqrt[1 - vo^2/c^2];
k = Sqrt[1 - rs/r0];

sol = NDSolve[{
r''[t] == -((G M)/r[t]^2) + r[t] θ'[t]^2 - (3 G M)/c^2 θ'[t]^2,
r'[0] == vr0,
r[0] == r0,
θ''[t] == -((2 r'[t] θ'[t])/r[t]),
θ'[0] == vθ0,
θ[0] == θ0,
τ'[t] == Sqrt[c^2 r[t] + r[t] r'[t]^2 - c^2 rs + r[t]^3 θ'[t]^2 - r[t]^2 rs θ'[t]^2]/(c Sqrt[r[t] - rs] Sqrt[1 - rs/r[t]]),
τ[0] == 0
}, {r, θ, τ}, {t, 0, T}, WorkingPrecision -> wp,
MaxSteps -> Infinity, Method -> Automatic,
InterpolationOrder -> All];

t[ξ_] :=
Quiet[χ /.
FindRoot[
Evaluate[τ[χ] /. sol][[1]] - ξ, {χ, 0},
WorkingPrecision -> wp, Method -> Automatic]];
Τ := Quiet[t[ι]];

x[t_] := (Sin[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]
y[t_] := (Cos[Evaluate[θ[t] /. sol]] Evaluate[r[t] /. sol])[[1]]

s[text_] := Style[text, FontSize -> font];  font = 11;

(* proper time *) Do[Print[
Rasterize[Grid[{{Show[Graphics[{
{Black, Circle[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 400, PlotRange -> 2 r0],
Graphics[{PointSize[0.01], Red, Point[{x[т], y[т]}]}],
ParametricPlot[{x[η], y[η]}, {η, 0, т},
PlotStyle->{LightGray}],
ParametricPlot[{x[η], y[η]}, {η, 0, т},
ColorFunction -> Function[{x, y, η},
Hue[0, 1, 0.5, Max[Min[(-т + (η + d2))/d2, 1], 0]]],
ColorFunctionScaling -> False]]},
{Grid[{
{" ", s["proper time"], " = ", s[N[т, 8]], s["GM/c³"]},
{" ", s["coordinate time"], " = ", s[N[Evaluate[τ[т] /. sol][[1]], 8]], s["GM/c³"]},
{" ", s["time dilation"], " = ", s[N[Evaluate[τ'[т] /. sol][[1]], 8]], s["dτ/dt"]},
{" ", s["angle"], " = ", s[N[Evaluate[(θ[т] /. sol) 180/Pi][[1]], 8]], s["degree"]},
{" ", s["radial distance"], " = ", s[N[Evaluate[r[т] /. sol][[1]], 8]], s["GM/c²"]},
{" ", s["x-Axis"], " = ", s[N[x[т], 8]], s["GM/c²"]},
{" ", s["y-Axis"], " = ", s[N[y[т], 8]], s["GM/c²"]}
}, Alignment -> Left]}}, Alignment -> Left]]
], {т, step, T, step}]

(* coordinate time *) Do[Print[
Rasterize[Grid[{{Show[Graphics[{
{Black, Circle[{0, 0}, rs]},
{Lighter[Gray], Dashed, Circle[{0, 0}, r0]}},
Frame -> True, ImageSize -> 400, PlotRange -> 2 r0],
Graphics[{PointSize[0.01], Red, Point[{x[Τ], y[Τ]}]}],
ParametricPlot[{x[η], y[η]}, {η, 0, Τ},
PlotStyle->{LightGray}],
ParametricPlot[{x[η], y[η]}, {η, 0, Τ},
ColorFunction ->
Function[{x, y, η},
Hue[0, 1, 0.5, Max[Min[(-Τ + (η + d2))/d2, 1], 0]]],
ColorFunctionScaling -> False]]},
{Grid[{
{" ", s["proper time"], " = ", s[N[Τ, 8]], s["GM/c³"]},
{" ", s["coordinate time"], " = ", s[N[ι, 8]], s["GM/c³"]},
{" ", s["time dilation"], " = ", s[N[Evaluate[τ'[Τ] /. sol][[1]], 8]], s["dτ/dt"]},
{" ", s["angle"], " = ", s[N[Evaluate[(θ[Τ] /. sol) 180/Pi][[1]], 8]], s["degree"]},
{" ", s["radial distance"], " = ", s[N[Evaluate[r[Τ] /. sol][[1]], 8]], s["GM/c²"]},
{" ", s["x-Axis"], " = ", s[N[x[Τ], 8]], s["GM/c²"]},
{" ", s["y-Axis"], " = ", s[N[y[Τ], 8]], s["GM/c²"]}
}, Alignment -> Left]}}, Alignment -> Left]]
], {ι, step, T, step}]

Last edited: May 30, 2016
11. May 30, 2016

### m4r35n357

I spent a good bit of time looking at this, and coming up with more complex approaches to calculating v (with and without being divided by the local "c") involving $g_{tt}$, but none seemed right (particularly wrt to the value being bounded between 0 and 1).

My pragmatic and empirical approach has been to look at the division of the speed with c (a light pulse instantaneously coincident with the particle and moving in the same direction), and assume the potential (position) terms cancel (so I am not really wrt to infinity any more?).

For example, all the cases I looked at of "knife-edge" orbits just outside the photon sphere for $a$ values of 1.0 ($r = 1.001$), 0.0 ($r = 3.001$) and -1.0 ($r = 4.001$) had $v/c$ tending to 1.0 (which seemed right). So pragmatically, if my answer is wrong I would expect the (adjusted for position-dependent potential, non-infinity) values of $v$ to be potentially either over or under 1.0. Neither seems right, so I kept the simpler solution, but I would appreciate comments on this assumption (I don't have the formal maths chops to make my reasoning any more rigorous).

If I am wrong and anyone knows of any references for calculating local $v / c$ in GR I would appreciate pointers, I couldn't find any.

NB my value of $v$ is purely due to post-processing of results, I would not be so cavalier about the simulation itself!

12. May 30, 2016

### George Jones

Staff Emeritus
I will try to expand on all this once I get to work. Trying to get my daughter and her visiting cousins out the door to school right now,

13. May 30, 2016

### Yukterez

For example, for a transversal initial motion right above the photon sphere at r=3GM/c²:
If you take the local velocity (proper distance differentiated by proper time) the v0 should be around c.
If you take the coordinate velocity (coordinate distances differentiated by coordiante time) it should be below c.
If you take the rapidity (coordinate distances differentiated by proper time) it should be above c.
In my equation of motion the latter is the case. That makes it very similar to Newtons equation, except the additional term - 3GM/c² θ'² in
r'' = -((G M)/r'²) + r θ'² - (3 G M)/c² θ'², θ''= -((2 r' θ')/r)
so I transform the initial velocity into the rapidity.
Then I differentiate coordinate time by proper time with
τ' = √[c² r + r r'² - c² rs + r³ θ'² - r² rs θ'²]/(c √[r - rs] √[1 - rs/r])
and finally get the positions at any given proper and coordinate time. From there I can easily transform into the system of any stationary shell observer at height R just by multiplying the coordinate time at infinity with √(1-rs/R).

Test:

When I place an observer at a shell at r0=10rs and let him throw a ball with an initial velocity of √(GM/r0) in an angle of 45° (measured locally) this ball approaches the horizon of the central mass after a revolution of 312.1 degrees, 434.323143 GM/c³ proper time and in coordinate time and shell time it only converges to r=rs, but never actually crosses it:

http://yukterez.net/org/newton.vs.einstein,45grad,v0=sqrt[GM,r0 [Broken]

On the left (Newtonian) and on the right (Einsteinian) the initial proper velocity and angle as measured by the observer at r0=10rs is the same, so the rapidity and energy is higher in the right animation. What results do you get for this configuration on your Java simulator?

Last edited by a moderator: May 7, 2017
14. May 30, 2016

### George Jones

Staff Emeritus
I am having trouble unpacking this paragraph. For example. what is $a$?

15. May 30, 2016

### m4r35n357

$a$ is the black hole spin parameter, normalized to the mass. The $r$ values I quoted are for particle orbits just outside the photon spheres at each of the $a$ values. Each example is supposed to be "close" to the actual photon sphere, where the value of $v/c = c/c$ has to be 1.0 by definition. My internal sense of physics justice says that there should not be a discontinuity in velocities at the photon sphere, but I can't express this in a more formal way. What I can do is try going closer, and verify that $v/c$ approaches 1.0 in practice. So far my prejudices have been confirmed. Does that make things clearer?

As I said, this is just post-processing and I will gladly chase up any analysis that I can find online. So far I have not found anything that calculates $v/c$ in curved spacetime.

16. May 30, 2016

### m4r35n357

Hmm, difficult to say as we are using different methods of setting up the problems! While I look at comparisons like this all the time, it is not trivial to set this up with the precise parameters you have given. I would need to set my radius to $10 r_s$ and experiment with the L factor until I measure (on the screen, and I have no protractor!) 45 degrees, so there would be no accurate figures to compare.

Can you give at least a brief outline how you derive your equations of motion, as I don't recognize them from the materials that I have read? Without understanding your equations, how you set the initial conditions (and you seem to have solved that problem already), and without Mathematica, I can't really follow what you are doing very well.

BTW did you try my simulator? It is HTML/JavaScript not Java ;)

17. May 30, 2016

### Yukterez

I don't think that the angle on the screen would be the same as the angle the observer on the shell would measure, since the Schwarzschild coordinates show the correct circumference, but the expandet length (the conversion factor between radius and circumference is not 2π, so the angle on the screen would look flatter than for the shell observer).

The equation uses the regular Schwarzschild coordinates as defined for an observer at infinity, but the proper time of the falling particle. If the velocity is c the rapidity is infinite. This allows an almost Newtonian form.

I tried, but I couldn't find the input parameters for the initial velocity and angle, but if you tell me how your input parameters are defined so I can convert my angle and velocity into them I'd love to! Do you work in cartesian, polar or spherical coordinates?

Last edited: May 30, 2016
18. May 30, 2016

### m4r35n357

Touche ;) We do seem to have diverged almost beyond the limits of compatibility. You seem to be invoking the local frame of the particle, whereas I am not using one (I don't expect anyone to be alive on it to measure angles!) so I don't really have anything to measure angles against.

As to the other matter, I just don't recognize those equations, in five years of looking at this I've never seen the Schwarzschild metric or its Christoffel symbols with a 3 in there. They look a bit like geodesic equations because of the second derivatives and velocity terms, but that 3 is a real spanner in the works. I really would like to understand what you are solving, not being funny here.

I might be able to compare with your results by varying the plummet until I match up the azimuthal angle of impact, I'll give that a go later..

19. May 30, 2016

### Yukterez

I played around with your simulator and at least I found out how to input a circular orbit at the photon sphere:

This gives an orbital velocity in coordinate distances by coordinate time of v0=0.57735, so v0/√[1-rs/r0] with rs=2 and r0=3 is c locally. So far so good, this is not so hard to convert between our observers, just a factor of √[1-rs/r0] here and there. Now I just have to find out how to zoom in (if I set G=M=c=1 for easier numerics the plot gets really small) and give some initial radial velocity additional to the angular velocity. What parameter of yours is responsible for the radial velocity component, or alternatively the launching angle of the test particle?

What I also don't understand is why your E and L show infinity in this example when L<100%: Screenshot

The 3 is all over the place, see
https://www.astro.umd.edu/~miller/teaching/astr498/lecture10.pdf
http://christian.vonschultz.se/forelant/gravitation_and_cosmology/2008-11-20.pdf

Last edited: May 30, 2016
20. May 30, 2016

### George Jones

Staff Emeritus
Okay, but this thread is about Shwarzschild, so let's keep $a=0$ everywhere.

For Yukterez's values of $r_0$, $v_0$, and $\varphi_0$

$$E = \frac{\sqrt{1 - \frac{2}{r_0}}}{\sqrt{1 - v_0^2}}$$

$$L = \frac{v_0 r_0 \sin \varphi_0}{\sqrt{1 - v_0^2}},$$

with $\varphi_0$ in radians.

Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted