Contents
@@ -419,20 +421,233 @@Contents
4. Planet Packing#
-4.1. Scaling a planetary system#
-4.2. The Golden ratio#
-4.3. Adding planets#
+Gravity is a long-range force/effect that attracts all bodies with mass to each other. If one considers two masses (\(m_1\) and \(m_2\)) on nearby orbits around a central body \((m_o;\ m_1\ll m_o\ \&\ m_2\ll m_o)\), then it is expected that each mass will perturb the other mass’ orbit causing it to either speed up or slow down a little. However, if one body is substantially more massive or close enough, then the lower mass body can be scattered to a wide/eccentric orbit or escape the central body altogether. Through numerical experiments involving scattering, it leads to a seemingly simple question:
+-
+
How close can two (or more) planets be separated such that scattering does not occur?
+
This question requires a few things to be defined so that limited numerical experiments can be performed. For example, the above main question leads to four other sub-questions:
+-
+
How long is necessary to say that a scattering event will not occur?
+How do we measure the separation between orbits?
+Is there a scale-free way to define the previous two questions?
+Will our choice of initial orbital elements bias the potential outcomes? If so, how to overcome this bias?
+
4.1. Stellar lifetimes as a constraint on stability#
+Orbital stability often goes undefined, even though it can have multiple meanings. Typically, it implies that a mass will remain on a bound orbit indefinitely (i.e., Lagrange stability; see Hayashi et al. (2023) or Barnes & Greenberg (2006)). Using indefinitely as a timescale is not practical. Scattering events can also transport a mass to a wider orbit without fully expelling the mass from the system. To address the potential for scattering events, please review updating the standard stability formulae in the previous section. Otherwise, we look to natural constraints on time to arrive at a more practical definition for stability.
+One natural constraint on the timescale for stability is the lifetime of the system. To judge the lifetime of a system, we can look to astrophysics and the lifetimes of stars. Consider the stability of the Solar System, where we might define it to be stable if the planets remain on bound orbits for \({\sim}10\ {\rm Gyr}\), or the main-sequence lifetime of a typical \(\rm G2V\) star. Even this definition has problems because the planet Mercury has a some probability of undergoing an instability (e.g., Laskar & Gastineau (2009)) within the remainder of the Sun’s main-sequence lifetime.
+If Mercury is expelled from the Solar System, does that make the whole system unstable?
+The short answer is no because the remaining planets will adjust/exchange their angular momentum until a new equilibrium is found (see Laskar (1997)) and was worked out mathematically by Laplace in 1784. This may have occurred in the past, where the giant planets’ orbits were more compact (Quarles & Kaib (2019), Nesvorny et al. (2018), Nesvorny (2011)). The giant planets may have mutually scattered their orbits and resulted in the configuration we see today. An example simulation of this process can be found at www.billyquarles.com/latest-research.
+The main-sequence lifetime of stars can be determined numerically given that we have some accurate estimates for a given star’s composition (e.g., hydrogen, helium, and metal mass fraction). The details can be found in a chapter on stellar evolution in Modern Astrophysics. The general estimate of main-sequence stellar lifetime \(\tau_{\rm ms}\) is that it is proportional to the inverse-cube of the stellar mass, \(\tau_{\rm ms} \propto M^{-3}\). Through the proportionality, we can scale other stars relative to our Sun. For example, a \(10\ M_\odot\) star’s lifetime is \(10^{-3}\) times \(10\ {\rm Gyr}\), or \(1\ {\rm Myr}\). See the lecture from Jason Kendall below for a general overview.
+
+
+
4.2. Scaling a planetary system using the Hill Radius#
+In the three-body problem, the relative gravitational influence on a third body can be determined in which there are five special locations called the Lagrange points. These points define extrema in the gravitational potential, where the net gravitational force on the third body vanishes. Additionally, they correspond to high and low topological levels of pseudopotential \(U\). This pseudopotential exists within a rotated reference frame, where two massive bodies appear at fixed locations and a third body reacts to the presence of the massive primaries.
+Consider two massive primaries lying along the \(x\)-axis of a reference coordinate system with the center of mass (barycenter) at the origin. Since the barycenter lies at the origin, then the center of mass can be determined using:
+
+(4.1)#\[\begin{align}
+0 &= m_o x_o + m_1 x_1, \qquad {\rm where}\ m_o > m_1;\\
+a_{\rm bin} &= |x_o| + |x_1|, \\
+x_o &= -\frac{m_1}{m_o+m_1}a_{\rm bin}, \\
+x_o &= -\mu a_{\rm bin}, \quad {\rm and}\quad x_1 = (1-\mu)a_{\rm bin}.
+\end{align}\]
+Two intermediate constants: the binary semimajor axis \(a_{\rm bin}\) and mass ratio \(\mu\) are introduced to simplify the notation, but they also provide a convenient means to scale the problem. The above analysis can be extended into two dimensions using the pseudopotential \(U\):
+
+(4.2)#\[\begin{align}
+U(x,y) & = \frac{1-\mu}{r_o} + \frac{\mu}{r_1} + \frac{n^2}{2}\left(x^2+y^2\right), \\
+r_o &= \sqrt{(x+\mu)^2 + y^2}, \\
+r_1 &= \sqrt{[x-(1-\mu)]^2 + y^2}.
+\end{align}\]
+The first two terms are simply the gravitational potential relative to each of the masses (\(m_o\ \&\ m_1\)) measured by the respective radii (\(r_o\ \&\ r_1\)). But the last term arises from the rotated coordinate system as a centrifugal potential.
+
+
+Note
+The coordinates in the pseudopotential are also scaled by the binary semimajor axis \(a_{\rm bin}\). In these coordinates, the mean motion \(n = 1\). If you want to express everything in non-scaled coordinates, the \(\mu\) and \((1-\mu)\) need to be replaced with \(x_o\) and \(x_1\), respectively.
+Taking the partial derivatives of the pseudopotential \(U\) leads to deriving the forces within the rotated coordinate system (i.e., \(F = -\nabla U\)). Integrating the full equations of motion derived from the pseudopotential leads to the Jacobi Integral and the associated constant of motion, the Jacobi constant \(C_J\). The Jacobi constant is given mathematically as:
+
+(4.3)#\[\begin{align}
+C_J = 2U - v^2,
+\end{align}\]
+which is similar in appearance to an energy relation. It is not because we are considering the restricted problem (\(m_2\ll m_1 < m_o\)), where energy and angular momentum aren’t conserved.
+Why do we bother with the rotated coordinate system? Special boundaries arise that constrain the motion of the smaller mass. These boundaries are called the zero velocity contours (ZVC), which are defined when \(v=0\) as the name suggests. As a result, the Jacobi constant for the ZVC is \(C_J = 2U\). Using the ZVC, we can illustrate the extent of gravitational influence of each mass \(m_o\) and \(m_1\).
+Below is a python code that computes the ZVC using \((C_J = 3,\ 3.05,\ 3.1,\ 3.15)\) for the Earth-Moon system. Three of the Lagrange points (black dots) are collinear, while the remaining two are the triangular Lagrange points. Particles that start near \(L_4\) or \(L_5\) are bound within the lobe-shaped contours.
+Interior to the gray contours lie regions where either the Earth or Moon dominate in the gravitational tug-of-war and delineate each body’s Hill Sphere. The use of ‘sphere’ isn’t mathematically correct, where the generalization of the region to three dimensions is actually an ellipsoid. This can be seen in that the region around the Moon is more squashed in the vertical compared to the horizontal extent.
+
+
+
+
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+def potential(x,y):
+ ro = np.sqrt((x+mu)**2+y**2)
+ r1 = np.sqrt((x-(1.-mu))**2+y**2)
+ U = (1.-mu)/ro + mu/r1 + (x**2+y**2)/2.
+ return 2.*U
+
+mu = 1./81.
+xi = np.arange(-2.,2.,0.001)
+yi = np.arange(-2.,2.,0.001)
+xx,yy = np.meshgrid(xi,yi)
+zi = potential(xx,yy)
+L1 = (mu/3.)**(1./3)
+L3 = -(1.+(5*mu/12.))
+L4 = 0.5-mu
+col = 'k'
+mark = '.'
+ms = 15
+
+theta = np.arange(0,2.*np.pi+0.01,0.01)
+
+fig = plt.figure(figsize=(6,6),dpi=150)
+ax = fig.add_subplot(111,aspect='equal')
+
+#ax.plot((1.-mu)*np.cos(theta),(1.-mu)*np.sin(theta),'k-',lw=2)
+CS = ax.contour(xi,yi,zi,np.arange(3,3.2,0.05),cmap=plt.cm.Set1)
+ax.plot([-1.25,1.35],[0,0],'k--',ms=1.5)
+ax.plot(1.-L1,0,marker=mark,color=col,ms=ms)
+ax.text(1.-L1,0.065,'$L_1$',fontsize=20,horizontalalignment='center')
+ax.plot(1.+L1,0,marker=mark,color=col,ms=ms)
+ax.text(1.+L1+0.025,0.065,'$L_2$',fontsize=20,horizontalalignment='center')
+ax.plot(L3+mu,0,marker=mark,color=col,ms=ms)
+ax.text(L3-0.1+mu,0.065,'$L_3$',fontsize=20,horizontalalignment='center')
+ax.plot(L4,np.sqrt(3)/2.,marker=mark,color=col,ms=ms)
+ax.text(L4,np.sqrt(3)/2.+0.065,'$L_4$',fontsize=20,horizontalalignment='center')
+ax.plot(L4,-np.sqrt(3)/2.,marker=mark,color=col,ms=ms)
+ax.text(L4,-np.sqrt(3)/2.-0.135,'$L_5$',fontsize=20,horizontalalignment='center')
+
+ax.plot([L4,-mu],[np.sqrt(3.)/2.,0],'k--',lw=1.5)
+ax.plot([L4,1.-mu],[np.sqrt(3.)/2.,0],'k--',lw=1.5)
+ax.plot([L4,-mu],[-np.sqrt(3.)/2.,0],'k--',lw=1.5)
+ax.plot([L4,1.-mu],[-np.sqrt(3.)/2.,0],'k--',lw=1.5)
+
+ax.plot(-mu,0,'.',mfc='skyblue',mec='skyblue',ms=50,label='Earth')
+ax.plot(1.-mu,0,'.',mfc='darkblue',mec='darkblue',ms=20,label='Moon')
+#ax.plot(0,0,'kx',ms=5,label="CoM")
+
+ax.legend(bbox_to_anchor=(1.1, .5, .3, .5),numpoints=1,fontsize=20,frameon=False)
+ax.set_xlim(-1.3,1.3)
+ax.set_ylim(-1.3,1.3)
+ax.axis('off');
+
+
+
+
+The appearance of the contours change with the mass ratio, which can be seen below if we consider the Sun-Earth system instead. The inset panel zooms in on the region around the Earth and is scaled using the Hill radius \(R_H\). The Hill radius is defined using the mass ratio \(\mu\) and the binary separation \(a_{\rm bin}\):
+
+(4.4)#\[\begin{align}
+R_H = a_{\rm bin}\left(\frac{\mu}{3}\right)^{1/3}.
+\end{align}\]
+Note that this definition is the result of an approximation that assumes circular orbits, where an additional term \((1-e_{\rm bin})\) would be necessary if the binary orbit were elliptical. The definition effectively becomes
+
+(4.5)#\[\begin{align}
+R_H &\approx q_{\rm bin}\left(\frac{\mu}{3}\right)^{1/3}, \\
+q_{\rm bin} &= a_{\rm bin}(1-e_{\rm bin}).
+\end{align}\]
+
+
+
+
+
+
+import numpy as np
+import matplotlib.pyplot as plt
+from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+import matplotlib.colors as cm
+
+def Jacobi_const(x,y,xo,yo,x1,y1):
+ ro = np.sqrt((x-xo)**2+(y-yo)**2) #distance from particle to the Sun
+ r1 = np.sqrt((x-x1)**2+(y-y1)**2) #distance from particle to the Earth
+ phi = (1.-mu)*(ro**2/2. + 1./ro) + mu*(r1**2/2. + 1./r1)
+ return 2*phi
+
+G = 4*np.pi**2
+M_sun = 1 # 1 solar mass
+M_E = 3.0035e-6 #mass of Earth in M_sun
+mu = M_E/(M_E+M_sun)
+R_H = (mu/3)**(1./3.)
+x_E, y_E = (1-mu), 0
+x_S, y_S = -mu, 0
+
+r_levels = [0.3,0.6,0.9]
+for i in range(1,10):
+ r_levels.append(1+i*0.25*R_H)
+n_lev = len(r_levels)
+Z_levels = np.zeros(n_lev+1)
+for r in range(0,n_lev):
+ Z_levels[r] = Jacobi_const(r_levels[r],0,x_S,y_S,x_E,y_E)
+Z_levels[-1] = Jacobi_const(0.5-x_E,np.sqrt(3)/2,x_S,y_S,x_E,y_E)
+Z_levels.sort()
+print(Z_levels)
+
+fig = plt.figure(figsize=(5,5),dpi=150)
+ax = fig.add_subplot(111, aspect='equal')
+axins = inset_axes(ax, width="40%", height="40%",bbox_to_anchor=(0.25, 0.05, 1, 1), bbox_transform=ax.transAxes)
+axins.set_aspect('equal')
+
+ax.plot(x_S,y_S,'y*',ms=15)
+ax.plot(x_E,y_E,'b.',ms=10)
+axins.plot(0,0,'b.',ms=10)
+axins.plot(1,0,'rx',ms=5)
+axins.plot(-1,0,'rx',ms=5)
+
+xi = np.arange(x_E-1.5*R_H,x_E+1.5*R_H,0.01*R_H)
+yi = np.arange(y_E-1.1*R_H,y_E+1.1*R_H,0.01*R_H)
+xx,yy = np.meshgrid(xi,yi)
+Z = Jacobi_const(xx,yy,x_S,y_S,x_E,y_E)
+
+vmin, vmax = 3, 3.003
+my_cmap = plt.colormaps['Set1']
+norm = cm.Normalize(vmin=vmin, vmax=vmax)
+my_cmap.set_over('k')
+
+ax.contour(xx,yy,Z, levels=Z_levels,zorder=5,cmap=my_cmap,norm=norm)
+axins.contour((xx-x_E)/R_H,yy/R_H,Z, levels=Z_levels,zorder=5,cmap=my_cmap,norm=norm,linewidths=0.5)
+axins.set_xlabel("X ($R_H$)")
+axins.set_ylabel("Y ($R_H$)")
+
+xi = np.arange(-1.1,1.1,0.001)
+yi = np.arange(-1.1,1.1,0.001)
+xx,yy = np.meshgrid (xi,yi)
+Z = Jacobi_const(xx,yy,x_S,y_S,x_E,y_E)
+ax.contour(xx,yy,Z, levels=Z_levels,zorder=2,cmap=my_cmap,norm=norm,linewidths=1)
+ax.set_ylabel("Y (AU)")
+ax.set_xlabel("X (AU)")
+ax.set_yticks(np.arange(-1,1.5,0.5))
+ax.set_xticks(np.arange(-1,1.5,0.5));
+
+
+
+
+
+[3.00000347 3.0008897 3.00093672 3.00095947 3.00106035 3.00124328
+ 3.0012661 3.00147675 3.00175584 3.00240875 3.03227121 3.69332466
+ 6.75659148]
+
+The squashing in the vertical (see panel inset) is even more obvious in this regime.