Numerical integration of motion

[felix00]

Hi,I'd like to build a simulation (realtime) of space ships near a black hole (neutral, still or rotating possibly). Key features would be:
1) the ships are test particles that do not affect the metric
a) possibly test rigid-bodies with GR consistent rotational DOF.
2) the ships can fire their rockets, changing their momentum, and actuate gyros to apply torque.
3) Eventually I'd like to have a toy universe where c is quite low and code a raytracer rendering system that includes lensing, aberrations and the lot. But this is if I manage with 1 and 2.

I did a lot of molecular dynamics in the past so I am familiar with discretisation of equation of motions and the algorithms, however I could not find any clear source for doing so in GR. I also have quite limited GR knowledge so I'm getting a bit lost in length derivations of very special cases.

Does anyone have more insight? Thanks for your help!

 

[pervect]

The basic equations you need for this are called the geodesic equations. The most direct form is when your spaceship has no thrust, then one can write (https://en.wikipedia.org/wiki/Solving_the_geodesic_equations)

$$\sum_{b,c=0..3} \frac{d^2x^a}{ds^2} + \Gamma^{a}{}_{bc}\frac{dx^b}{ds}\frac{dx^c}{ds} = 0$$

Here the ##x^a## are coordinates, for instance ##x^0 = t, x^1 = x, x^2=y, x^3=z##. The equations are in parametric form, so we are solving for a solution of the form
$$t(s), x(s), y(s), z(s)$$
or in the original notation
$$x^0(s), x^1(s), x^2(s), x^3(s)$$

There are 4 equations, corresponding to values for a of 0, 1,2, and 3. b and c are dummy indices, which one sums over, so there are 16 terms if you write out all the terms.

The ##\Gamma## coefficients are called "Christoffel symbols", and they can be calculated from the metric. The metric describes the space-time, you can use the plain Schwarzschild metric for starters, which would be a non-rotating black hole. If you want to use cartesian-like coordinates (x,y,z) rather than (r, ##\theta##, ##\phi##) you probably want to use not regular Schwarzschild coordinates, but isotropic coordinates. But it's much easier to stick with regular Schwarzschild coordinates and use ##r, \theta, \phi##.

In GR you can use any coordinates you like, which is a powerful feature of the theory that can be confusing - it makes no physical difference (just as it makes no difference whether you use Cartesian or polar coordinates in a freshman physics problem), but you have the responsibility of choosing coordinates that mean something to you, and you have the opportunity to confuse yourself over such issues. It doesn't sound like a major problem, but people often think about problems in terms which require specific features of the coordinates (like the ability to find distances by subtracting coordinates), which hinders their undertanding of the significance of a solution written in coordinates that do not have these special features. A classic example of this is to interpret "r" in Schwarzschild coordinates as being a radius- the issue is that subtracting r values does not give you the distance.

The Christoffel symbols are calculated from the metric coefficients, which describe the space-time. See for instance https://en.wikipedia.org/wiki/Christoffel_symbols

$$\Gamma_{abc} = \left(\frac{\partial g_{ca}}{\partial x^b} + \frac{\partial g_{cb}}{\partial x^a} - \frac{\partial g_{ab}}{\partial x^c} \right)$$

Unfortunately, you need to learn how to raise the index, as you actually need ##\Gamma^a{}_{bc}##, not ##\Gamma_{abc}##. Wiki has the forumla, but it uses tensor notation. And it'd be a bit long to explain.

For an overview of putting this all together to find an unpowered orbit, you might look at http://www.fourmilab.ch/gravitation/orbits/, "Orbits in Strongly Curved Space-time". The treatment follows MTW's treatment in "Gravitation" fairly closely. It has a java applet, which may be in somewhat of a state of disrepair, last time I tried to use it it raised a bunch of security warnings.

I haven't discussed in detail how to add thrust, basically rather than set the geodesic equation to zero, you'd set it to the thrust value. But I think a good first goal would be to write and solve the geodesic equations for an unpowered test particle orbit in a Schwarzschild space-time.

This is already too long, and it's barely started. But I hope it's of some help.

Sean Caroll's lecture notes on GR might be a good source to replace MTW as a formal reference - they're available online, but not an easy read. https://www.preposterousuniverse.com/grnotes/. But they might be rather technical.

 

[PAllen]

Just to add a small addition to pervect's write up, you can posit an arbitrary rocket path as any time like trajectory and then ask what what thrust would be needed to achieve it. You just take the absolute derivative by proper time of the 4-velocity along the path; this will give the acceleration measured by an accelerometer at each moment.

 

[felix00]

Thanks a lot! this is very helpful!
Just a few points are not clear, mostly because it has been ~10 years since my GR course (that I didn't fully understand anyway).
The whole high/low indexing confuses me a bit.

1) [ itex ] x^0 = t, x^1 = x, x^2=y, x^3=z [ /itex ]
I noticed the high indexes, meaning they are contravariant? I've been reading about it but still confuses me.
Are these the coordinates I should store in memory and use to calculate what an observer on a second ship will see?

2) what exactly are t and s?
I suppose one is the time measured by a clock on the ship, and the other is measured by an observer infinitely far from the black hole? Which one is which?

 

[vanhees71]

##t## is the coordinate time and ##s## the proper time of your rocket (i.e., the time an observer at rest all the time relative to the rocket).

You have to refresh your mind concerning GR only a bit. Just look for the equations of motion of a test particle subject only to gravity. It turns out to be the geodesics in spacetime, and their equation is given in the posting #2 by pervect. A very good source to learn GR is Landau&Lifshitz vol. 2.

 

[pervect]

Thanks a lot! this is very helpful!
Just a few points are not clear, mostly because it has been ~10 years since my GR course (that I didn't fully understand anyway).
The whole high/low indexing confuses me a bit.

1) ## x^0 = t, x^1 = x, x^2=y, x^3=z ##
I noticed the high indexes, meaning they are contravariant? I've been reading about it but still confuses me.
Are these the coordinates I should store in memory and use to calculate what an observer on a second ship will see?

The ##x^i## are coordinates. The metric is like a map or chart of space-time, and the coordinates tell you when and where you are on the map. But they don't directly tell you anything about what an observer on a second ship might see (or compute) -a second ship would take some different math which I haven't gone into at all. The usual approach to what I think you're after is called a "frame field", though this technique has the limitation that it only describes things very near the second ship. It's rather like making a small map that's "to scale" around the second ship. The large map usually doesn't even try to be "to scale", since it's impossible to do, anyway.

A useful ananlogy - if you think of the curved surface of the Earth, you can never draw a large map of the Earth on a flat sheet of paper (though you could represent it to scale on a globe). But you can draw a map of some small subsection of the Earth - say a city - which will be very accurate.

I would actually suggest using Schwarzschild coordinates for at least the first go-around, in which case ##x^0=t, x^1=r, x^2=\theta, x^3=\phi##.

2) what exactly are t and s?
I suppose one is the time measured by a clock on the ship, and the other is measured by an observer infinitely far from the black hole? Which one is which?

For the case of a spaceship, s is indeed proper time, and it can be thought of as what a clock on the ship would read. "t" is just a coordinate, rather than invoking some "observers at infinity", I'd just say that "t" is what you write down to tell when an event occurs. It's an agreed on and more or less arbitrary label, a standard one chooses to communicate to someone where and when an event is happening in an agreed on manner.

So when one write down t,r,theta, and phi, one is telling when and where on one's charts the space-ship is, at some proper time s that can be regarded as the time on the ships clock. Hence one has t(s), r(s), theta(s), and phi(s).

If you want to compute what a second ship linterally sees, you need to do some ray tracing. You can compute the path that light takes via the same formula (the geodesic equations), except that you can't interpret "s" as any sort of proper time in that case. The math remains unchanged though, light rays follow geodesics that satisfy the geodesic equations just light space-ships do. S is called an "affine parameter" when one is computing the path through space-time that light rays take.

 

[felix00]

Thanks again! I used the spherical coordinates, and cooked up a quick s-step integrator which is probably wrong/inaccurate.
I'll do some more research to improve it.

Nevertheless I got somewhat reasonable.
https://www.dropbox.com/s/x1bilwwsk861iyj/GR.png?dl=0
https://www.dropbox.com/s/x1bilwwsk861iyj/GR.png?dl=0
The blue trajectory was obtained with starting velocity {1,0,0,0.1} and the yellow with {1,0,0,0.3}. Both started in the same point {t=0, r=10, theta=0, phi=0}. For simplicity, 2GM = 1, M=1 and c=1, Schwarzschild metric.
Intuitively I understand what I got: the yellow ship had more tangential velocity so it ends up in a larger trajectory. The blue one probably goes too close and the numerics goes bad.

What is the meaning of the initial "time" velocity. I set it to 1 for both tests, but when set to 0 I get weird trajectories.
I also noticed that it changes over time by small amounts. Is it because of my bad integrator?

 

[vanhees71]

You must have a time-like trajectory, so you must start with initial conditions such that ##g_{\mu \nu} \dot{x}^{\mu} \dot{x}^{\nu}>0## (in the west-coast convention with the metric "mostly -").

 

[felix00]

I think my code has the -+++ convention, if by that you mean the sign of the diagonal matrix elements of g.
Regardless of my choice of initial dx^0/ds the result of g_mn v^m v^n (sorry I cannot figure out how to put nice tex code) is negative, provided that the initial radial velocity and theta velocity are 0.

 

[vanhees71]

Yes, then your initial four-velocity must be such that ##g_{\mu \nu} v^{\mu} v^{\nu}<0## (time-like geodesics).

Concerning typesetting in the forums, have a look at the LaTeX primer:

https://www.physicsforums.com/help/latexhelp/

 

[pervect]

The time component of the 4-velocity is commonly called time dilation. It tells you the relative rate at which your coordinate time advances per unit of proper time. Proper time is what you measure with your local onboard clock, and coordinate time is the number that tells you when something occurred according to your coordinate chart.

Something I didn't mention. In order for the parameter s to be able to be interpreted as proper time, you need to follow the convention for a -+++ metric signature that the magnitude of the four velocity ##v^\mu## is equal to -1, by which I mean that

$$\sum_{a,b=0..3} g_{ab} v^a v^b = -1$$

If you were using a +--- signature, the magnitude would be +1 rather than -1.

If you don't follow this convention, you will still compute geodesic paths correctly via your equations, but the parameter s will just be a parameter that is related to proper time by a simple linear relationship, but will not be proper time.

Example: If you have a stationary particle far away from any gravitating mass so that special relativity applies, the metric is diag(-1,1,1,1) with your sign convention, and we want the rate at which coordinate time advances to be equal to the rate at which proper time advances. This means that v=(1,0,0,0), and ##g_{ab} v^a v^b = -1##.

Things are a bit different for light-like (null) geodesics - if you have a light-like worldline, the magnitude of it's 4-velocity is always zero - there is no meaningful concept of proper time for a light-like geodesic.

 

[felix00]

Thanks a lot to everyone, your answers made it way clearer than my old course notes!
I managed to put together a simple Runge-Kutta integrator and got meaningful trajectories. Made my day!
However it seems quite slow. I am doubting it is possible to have real-time raytraced rendering in a curved spacetime.

 

[PAllen]

My understanding is, indeed, that those fancy GR Ray trace images you see require LOTS of cpu time per frame. Each such movie can be a substantial project to generate.