Introduction to Astronomy - Kepler Problem

Kepler Problem A mass m orbits a much larger mass M in a non-circular orbit. We want to find r(t) and θ(t). mass m orbiting a much larger mass M in a noncircular orbit Start with F = ma in the two dimensions of the orbital plane. [ F_r = -\frac{GmM}{r^2} = -m(\frac{d\theta}{dt})^2r + m\frac{d^2r}{dt^2} ] [ F_\theta = 0 \Rightarrow \text{conservation of orbital angular momentum, L} ] but, [ L = mvr = mr\omega r = mr^2\omega = mr^2\left(\frac{d\theta}{dt}\right) ] Substitution 1: Eliminate the (dθ/dt) term from the equation for Fr, which in its present form, involves r, θ, and t. [ -\frac{GmM}{r^2} = -m\left(\frac{L}{mr^2}\right)^2r + m\frac{d^2r}{dt^2} ] or, [ \frac{d^2r}{dt^2} + \frac{GM}{r^2} - \frac{L^2}{m^2r^3} = 0 ] This equation turns out to be difficult to solve. We will therefore use the chain rule of differentiation to find an equation for r(θ) instead of r(t). This will yield the shape of the Keplerian orbit. From our conservation of angular momentum expression [L = m\left(\frac{d\theta}{dt}\right)r^{2}] we have [\frac{d\theta}{dt} = \frac{L}{mr}] Thus, [\frac{dr}{dt} = \left(\frac{dr}{d\theta}\right)\left(\frac{d\theta}{dt}\right) = \frac{L}{mr}\left(\frac{dr}{d\theta}\right)] Finally, [\frac{d^{2}r}{dt^{2}} = \frac{L^{2}}{m^{2}}\frac{d}{dt}\left[\frac{1}{r}\left(\frac{dr}{d\theta}\right)\right]] The equation for r(θ) can now be written as follows: With Substitution: #2 [\frac{L^{2}}{m^{2}}\frac{d}{d\theta}\left[\frac{1}{r}\left(\frac{dr}{d\theta}\right)\right] + \frac{GM}{r^{2}} - \frac{L^{2}}{m^{2}r^{3}} = 0] or [\frac{d}{d\theta}\left[\frac{1}{r}\left(\frac{dr}{d\theta}\right)\right] + \frac{GMm^{2}}{L^{2}} - \frac{1}{r} = 0] It is not obvious that this equation for r(θ) is easier to solve than our original equation for r(t), but it is! The trick is to make the substitution of variables r = 1/u With Substitution: #3 [\frac{d}{d\theta}\left[u^{2}\left(\frac{du}{d\theta}\right)\right] + \frac{GMm^{2}}{L^{2}} - u = 0] or [\frac{d^{2}u}{d\theta^{2}} + u = \frac{GMm^{2}}{L^{2}}] First, set the right-hand side equal to zero to yield an equation [\frac{d^2u}{d\theta^2} + u = 0] This equation has solutions of the form: [u(\theta) = B \cos(\theta + \beta)] where B and β are constants that are set by the boundary conditions. (Verify this solution for yourself) To find the solution with the right-hand side not equal to zero, simply add GMm²/L² to the solution found above. Thus [u(\theta) = \frac{1}{r(\theta)} = B \cos(\theta + \beta) + \frac{GMm^2}{L^2}] If we orient our plot of r(θ) so that the minimum and maximum values of r occur along the X axis, this fixes β = 0. Finally, [r(\theta) = \frac{L^2/GMm^2}{1 + B \cos{\theta}}] The constant B has yet to be determined. Our expression for r(θ) is that of an ellipse. Therefore, let us explore some of the properties of ellipses. The drawing on the left shows an ellipse with its semimajor and semiminor axes of lengths a and b, respectively. The drawing on the right shows the familiar geometrical construction of the same ellipse, where a string of fixed length is pulled taut around two fixed pins (defining the two foci) and a pencil at point P. The separation of the two foci are, by definition, equal to 2ea, where e is the orbital eccentricity. We can now derive the equation of an ellipse in r and θ coordinates. It is easy to show that the length of the string is 2a + 2ea = 2a(1+e). We now employ the law of cosines on the drawing on the right: [r'^2 + (2ea)^2 - 4ear\cos{\theta}] [a^2(1-e^2) = a^2(1-e)\cos{\theta}] Finally, we find the equation of an ellipse in polar coordinates: [r(\theta) = \frac{a(1-e^2)}{1 - e \cos{\theta}}] This is exactly the same form that we derived from the Keplerian equation of motion! By equating the numerators of these two expressions, we find: [a(1-e^2) = \frac{L^2}{GMm^2}] Thus, we see that for a given set of masses, if we choose the semimajor axis, a, and the orbital angular momentum, L, this sets the orbital eccentricity. This concludes the proof of Kepler's first law. Kepler's second and third laws easily follow: Recall that [L = mr^2\left(\frac{d\theta}{dt}\right)^2] is constant but, [\frac{1}{2}r^2d\theta] = the differential area swept out by the orbit [\Rightarrow \frac{dA}{dt} = \frac{1}{2}r^2\frac{d\theta}{dt} = \frac{L}{2m}] Finally, we have [\frac{dA}{dt} = \frac{L}{2m} = \text{constant}] Kepler's 2nd law If we integrate the expression for swept-out area around the orbit, we should arrive at the area of the ellipse. [\int_{0}^{P} \frac{dA}{dt} dt = \text{Area of ellipse} = \frac{L}{2m}P = \pi ab] [Area of ellipse in terms of semimajor and semiminor axes] From the drawing, one can see that a, b, and e are related by [\frac{b^2}{a^2} = 1 - e^2] Therefore, we can arrive at the following relation between orbital period, P, and semimajor axis, a: [\left(\frac{LP^2}{2m^2}\right)^2 = \pi^2 a^2 b^2 = \pi^2 a^4 (1-e^2) = \pi^2 a^3 \frac{L^2}{GMm^2}] where we have utilized the expression a(1-e²) = L²/GMm² from page 5. Finally, we are left with [\left(\frac{LP^2}{2m^2}\right)^2 = \frac{GM}{a^3}] or [\frac{P^2}{a^3} = \frac{4\pi^2}{GM}] This is exactly the same form as for a circular orbit, except here a is the semimajor axis. Time Dependence of the Orbital Motion Return to the original second-order differential equation in r(t) that we started with: [\frac{d^2r}{dt^2} + \frac{GM}{r^2} - \frac{L^2}{m^2r^3} = 0] Force Equation but, [\frac{d^2r}{dt^2}] can be written as [\frac{1}{2}\frac{d}{dt}(v_r^2)] where vr is the radial velocity. With this, we can integrate the above equation (with forces and accelerations) to yield an energy equation: [\frac{1}{2}v_r^2 - \frac{GM}{r} + \frac{L^2}{2m^2r^2} = \text{const}] At periastron, vr = 0, r = a(1-e), and L² = GMm²a(1-e²) (energy equation). Thus, const = -GM / 2a [\frac{1}{2}v_r^2 - \frac{GM}{r} + \frac{L^2}{2m^2r^2} = -\frac{GM}{2a}] [\frac{d\psi}{dt} = \sqrt{-\frac{GM}{a} + \frac{2GM}{r} - \frac{L^2}{m^2r^2}}] [dt = \frac{dr}{\sqrt{-\frac{GM}{a}r^2 + 2GMr - \frac{L^2}{m^2}}}] Now substitute L² = GMm²a(1-e²) [dt = \frac{rdr}{\sqrt{\frac{GM}{a}}\sqrt{-r^2 + 2ar - a^2(1-e^2)}}] [dt = \sqrt{\frac{a}{GM}} \frac{rdr}{\sqrt{1-e^2 - (r/a)^2}}] The equation in this form can be integrated by making the following (quite standard) trig substitution: Let (r-a) = -ae cos(u) [dr = ae \sin(u) du] [dt = \sqrt{\frac{a}{GM}} \frac{a(1-e\cos{u})(ae\sin{u})du}{\sqrt{a^2e^2 - a^2e^2\cos^2{u}}}] [dt = \frac{a^2e(1-e\cos{u})\sin{u}du}{\sqrt{\frac{GM}{a}}ae\sin{u}}] [\frac{2\pi dt}{P} = (1-e\cos{u})du] [M = \frac{2\pi(t-T_0)}{P} = u - e\sin{u}] (mean anomaly) (eccentric anomaly) And recall: [r = a(1-e\cos{u})] (from our trig substitution) The above two equations represent parametric solutions for r(t), namely r(u) and u(t). When fitting Keplerian orbits to observational data, one usually writes the parametric equations as x(u) and y(u) with the same u(t) equation. We can use our elliptical orbit solutions [r = \frac{a(1-e^2)}{1+e\cos{\theta}}] to find x(u) and y(u). (Note that there was a + sign in the denominator, which is opposite to the sign derived earlier. This simply flips the ellipse 180° in orientation.) From the equation for the ellipse, we find [r \cos{\theta} = x = a(1-e^2)/(1+e)] Substituting [r = a(1-e\cos{u})] from above, we find [x(u) = [a(1-e^2) - a(1-e\cos{u})]/e = ae(\cos{u}-e)/e] [x(u) = a(\cos{u}-e)] Similarly, [y(u) = a\sqrt{1-e^2}\sin{u}] and we had [\frac{2\pi(t-T)}{P} = u - e\sin{u}] These non-linear equations are used to fit orbits when x(t) and y(t) have been measured. Basically, for any observation time t, one solves the bottom expression by Newton's method for u. Then, x(u) and y(u) are evaluated from the top two equations and compared with the measurements. Acceleration in a Plane [\vec{a} = \frac{d\vec{v}}{dt}] So, let [\vec{F} = F\hat{r}], where [\hat{r}] is the unit vector in the radial direction Therefore, the velocity [\vec{v} = \frac{d\vec{r}}{dt} = \frac{dr}{dt}\hat{r} + r\frac{d\hat{\theta}}{dt}], by the chain rule [\frac{d\hat{\theta}}{dt} = \omega\hat{\theta}] So, [\vec{v} = \frac{dr}{dt}\hat{r} + \omega r\hat{\theta}] Finally, [m\vec{a} = m\frac{d\vec{v}}{dt} = m\left[\frac{d^2r}{dt^2}\hat{r} + \left(\frac{dr}{dt}\right)\left(\frac{d\hat{r}}{dt}\right) + \left(\frac{d\omega}{dt}\right)r\hat{\theta} + \omega\left(\frac{dr}{dt}\right)\hat{\theta} + \omega r\frac{d\hat{\theta}}{dt}\right]] From a derivation similar to the one above for [\frac{d\hat{r}}{dt}], we find [\frac{d\hat{\theta}}{dt} = -\omega\hat{r}] Collecting terms, we have [\vec{a} = \left[\left(\frac{d^2r}{dt^2}\right) - \omega^2r\right]\hat{r} + \left[\frac{d\omega}{dt} + 2\left(\frac{dr}{dt}\right)\omega\right]\hat{\theta}] Centripetal acceleration Coriolis acceleration