The N-body problem refers to predicting the movement of N many bodies in space under gravity. We examine this problem with 4 bodies, namely: the Sun, Earth, Moon, and the New Moon which is our Moon’s own satellite. For this problem, functions like NBodySimulation, along with equations like the Hill Radius, Roche Limit and escape velocity, were used to determine the values of parameters: the mass, velocity and position of the New Moon. According to these values, it is demonstrated at the end what would occur when the New Moon is outside of these intervals with several parametric plots.
Introduction
Introduction
The N-Body problem investigates how multiple celestial bodies move due to gravitational forces. When dealing with two bodies, such as the Earth and Moon, we can accurately predict their trajectories using equations formulated by Kepler and Newton, such as F = G . However, when more than two bodies are involved—like planets or stars—the interactions become more complex due to chaotic and nonlinear behaviors. Early attempts by Euler and Lagrange involved using differential equations to predict these motions. While these equations provide valuable insights, they may not always include all details accurately. In our 4-Body Problem, we try to determine the parameters that might affect the 4-Body system and values of these parameters. Functions like the NBodySimulation and equations such as the Hill Radius are used here. These tools allow us to simulate the movement of the bodies and visualize how multiple bodies influence each other over time.
m
1
m
2
2
r
The 4-Body Problem and Defining the Variables
The 4-Body Problem and Defining the Variables
The function NBodySimulation takes the initial position and velocity vectors of N many bodies with their masses and generates a simulation of the motion of the bodies. Below, there is an example dataset and a simulation with 4 bodies.
In[]:=
exampleData=NBodySimulation["InverseSquare",{<|"Mass"1,"Position"{1,0,0},"Velocity"{0,0,0.5}|>,<|"Mass"2,"Position"{1,1,1},"Velocity"{0.5,0,1}|>,<|"Mass"3,"Position"{1,1,0},"Velocity"{0,0.5,0}|>,<|"Mass"4,"Position"{0,1,1},"Velocity"{0.5,0,0}|>},1]
Out[]=
NBodySimulationData
To plot the simulation in 3D, the function ParametricPlot3D is used. As I wanted to plot each body, “All” is used in the function. The simulation below shows the position of the bodies from 0 seconds to 1 second in 3D.
In[]:=
ParametricPlot3D[Evaluate[exampleData[All,"Position",t]],{t,0,1}]
Out[]=
As I will be simulating the Sun-Earth-Moon-NewMoon situation, we need to define the variables accordingly. For this problem, some parameters are determined: the New Moon’s distance to the other bodies, its mass, and its velocity.
The Masses of the Bodies
The Masses of the Bodies
The Sun' s, Earth' s and Moon' s masses are defined below. The New Moon' s mass is not defined here as it will be a variable, X parameter to see if the New Moon is in the orbit. Firstly, the bodies are defined to make it easier to use them in different code lines without using natural language input each time.
The shortcut of EntityProperty is used to get the mass of the bodies.
The Positions of the Bodies
The Positions of the Bodies
As it can be seen from the example data, the initial coordinates of the bodies need to be in vector form. Even though data for the net distance of the bodies to each other is accessible, EntityValue and EntityProperty are used to get the positions of the bodies in vector form. As the locations of the bodies change daily and EntityValue requires a date to get the locations, to get a constant value, instead of directly using , I specified the date and time inside the brackets.
DateObject[]
To get the coordinates, we need an origin. I determined the origin to be the Sun and got the HelioCoordinates of the bodies which means the positions of them relative to the Sun. As it is a Sun centered model, the Sun will be initially at {0,0,0}. For that, Table is applied to repeat 0-3 times in a list, and then Quantity function is used to give the coordinates a unit. For the Earth and the Moon, EntityValue and EntityProperty are used. This function gave the values in Astronomical Unit; therefore, UnitConvert is used to convert it to Meters.
In[]:=
sunPosition = Table[Quantity[0,"Meters"],3];earthPosition = UnitConvert[EntityValue[earth, EntityProperty["Planet","HelioCoordinates",{"Date"date}]],"Meters"];moonPosition = UnitConvert[EntityValue[moon, EntityProperty["PlanetaryMoon","HelioCoordinates",{"Date"date}]],"Meters"];
The Velocities of the Bodies
The Velocities of the Bodies
Along with the initial position coordinates, we also need the velocity of the bodies in vector form. Again, we can access the net velocity of the bodies easily; however, to get them in 3-D coordinates, the displacement of each body in one day is divided by one day. To get the displacement, the locations of the bodies on the next day are subtracted from the previous locations obtained with the date variable. A variable called “day” is also defined to get the time in one day. As the units become Meters/Day, UnitConvert is used to convert it to Meters/Second. The net velocity for the Earth is around its real value is 297,800 m/s, however, the net velocity obtained for the Moon was higher than it should be. Therefore, the values are scaled approximately to its real velocity which is around 1,022 m/s. As I got 28,299 m/s with my method, I divided it by 28.
In[]:=
day = ;sunVelocity = UnitConvert[sunPosition/day,"Meters"/"Seconds"];earthVelocity= UnitConvert[(EntityValue[earth, EntityProperty["Planet","HelioCoordinates",{"Date"date+day}]]-earthPosition)/(day),"Meters"/"Seconds"];moonVelocity = UnitConvert[(EntityValue[moon, EntityProperty["PlanetaryMoon","HelioCoordinates",{"Date"date+day}]]-moonPosition)/(day),"Meters"/"Seconds"]/28;
The Roche Limit
The Roche Limit
The Roche Limit is defined as the minimum distance a satellite can approach its primary body (the body it is in the orbit of) without the gravitational force between them overcoming the internal gravity of the satellite causing the satellite to break apart. The common formula is given by d ≈ R where
2×
ρ
1
ρ
2
◼
ρ
1
◼
ρ
2
◼
R: The Radius of The Primary Body
◼
d: Roche Limit
Another formula for the Roche Limit which I used is given by where is the mass of the primary body and is the mass of the satellite. For this problem, the primary body is the Moon and the satellite is the New Moon. Below, the radius of the Moon and a function for the Roche limit are defined.
d≈R
2×
m
1
m
2
m
1
m
2
In[]:=
rocheLimit[r_,m1_,m2_]:=r*(2*m1/m2)^(1/3);moonRadius=UnitConvert, "Meters";
The Hill Radius
The Hill Radius
The Hill Radius is the radius within which a body’s gravity dominates over another body’s. For instance, Earth’s Hill Radius is 1.5 million kilometers, meaning Earth’s gravity dominates the Sun’s within this range, allowing the Moon to orbit Earth. Beyond this distance, the Sun’s gravity would dominate. We need the New Moon to be within a range where the Moon’s gravity dominates over both the Sun’s and Earth’s. The Hill Radius is given by the formula r ≈ a×(1-e)× where:
3
m
1
3×
m
2
◼
a: The Semimajor Axis
◼
e: The Eccentricity of The Secondary Body
◼
m
1
◼
m
2
The formula is for a 3-Body problem. For a 4-body problem, the Moon’s Hill Radius first with Earth as the primary body is calculated. If this radius is entirely within Earth’s Hill Radius, the satellite will orbit the Moon. If not, we then calculate the Moon’s Hill Radius with the Sun as the primary body to determine the extent of the Moon’s gravitational dominance. Below is the function for the Hill Radius and the variables in the formula are defined . As the variables are already constant, we can get the Hill Radius value independent of the New Moon.
In[]:=
hillRadius[a_,e_,m1_,m2_]:=a*(1-e)*(m1/(3*m2))^(1/3);eccentricityMoon=;distanceMoonEarth=UnitConvert, "Meters";hillRadiusValue=hillRadius[distanceMoonEarth,eccentricityMoon,moonMass,earthMass]
Out[]=
According to Marvellous and Roy, the satellite needs to be well within the half of the Hill Radius to be stable in the orbit; therefore, the maximum distance becomes half of the Hill Radius value (Marvellous & Roy, 1981).
In[]:=
maxDistance=hillRadius[distanceMoonEarth,eccentricityMoon,moonMass,earthMass]/2
Out[]=
From these values, we get the maximum distance to be around 3.0510^7 meters which is well in the Earth’s Hill Radius so we don’t need the extra steps mentioned above.
Forces
Forces
We need the gravitational force exerted by the Moon on the New Moon to be greater than the force exerted by the Sun or the Earth on the New Moon to make sure it doesn’t get dragged to the other bodies. For this, the formula mentioned in the Introduction is used which is F = G . While is the mass of the New Moon, is the other body that is accounted for. From this,
m
1
m
2
2
r
m
2
m
1
◼
F
Moon-NewMoon
m
moon
m
newMoon
2
r
distanceMoon-NewMoon
◼
F
Earth-NewMoon
m
earth
m
newMoon
2
r
distanceEarth-NewMoon
◼
F
Sun-NewMoon
m
sun
m
newMoon
2
r
distanceSun-NewMoon
The force by the Moon needs to be greater than or equal to the other forces. Hence, from simplifying the equation:
The Velocity of the New Moon
The Velocity of the New Moon
The Mass of the New Moon
The Mass of the New Moon
Along with distance and velocity, the mass is also a factor that needs to be considered. According to Cruz, the satellite of a body should have at most 3 percent of the primary body (Cruz, 2020). The minimum mass it should have can be determined with maximum and minimum distances as the minimum distance depends on the mass of the New Moon and the maximum distance is a constant. The maximum distance should be larger than the minimum distance, therefore the value of mass can be found where the minimum distance would be larger and that value would be the minimum mass.
NSolveValues function is used to get the root of the equation. As the NSolveValues gives the value as a list, Part function is used to get the first element to take it as a Real number instead of a list.
From these equations, it can be seen that the mass of the New Moon should be approximately between 2.710^18 kg and 2.2010^21 kg.
Visualizations
Visualizations
A simulation, a slider, and several plots for the motion of the bodies are created with some functions such as NBodySimulation, Parametric3DPlot and so on.
Creating the 4-Body Simulation
Creating the 4-Body Simulation
A function called NBodyVisualization takes the number of objects, mass, position and velocity of each object, and the starting and finishing time of the simulation to represent the movements of the bodies.
The Module function is used to update the values in the slider and make an interface that joins the text, simulation and slider together. Just like the exampleData, a dataset for NBodySimulation needs to be defined. As the function will take the variables such as mass and position as a list (all N elements together), the Table function is used to associate each body to its properties. Quantity function is used to give proper units to the variables and the Norm function gives the net value such as net velocity. To display texts in a column, the Column function is used. To join the variable values with the text in one row, the Row function is used. After the simulation, we get the position plot of each body with the Plot and Norm functions. The Which functions below check if the parameters are within the right interval as explained above. However, the third Which function checks if any of the bodies collided by checking if there are any duplicates in the position vectors.
The Manipulate function is used to get the sliders for the simulation and NBodySimulation function defined above is used with the values of the bodies in lists.
Creating the Parametric Plots
Creating the Parametric Plots
By creating parametric plots, we can get the movement of each body throughout time. For these plots, I only included the Moon and the New Moon as, otherwise, the plane gets too large to include all the bodies and they become hard to spot. As each body needs to be defined in the dataset, the New Moon is defined as Mars, but it won’t have any effect as the variables of the New Moon are defined independent of the Mars definition (such as its mass is defined to be 0.02 of the Moon’s mass which is also in the required interval for it to be in the orbit).
The Parametric Plot for the New Moon at Long Distances
The Parametric Plot for the New Moon at Long Distances
A distance longer than the maximum distance is determined and added to the Moon’s position to get a vector position with the Quantity function as they can be added together only when they have the same unit; later a velocity value is also determined in the right interval. Then, I made it into a dataset and plotted it with the ParametricPlot3D and Graphics3D to get a representation of the bodies as spheres.
It can be seen from this simulation that when the bodies are too far away, they are not in the orbit.
The Parametric Plot for the New Moon at Short Distances
The Parametric Plot for the New Moon at Short Distances
This time, a distance shorter than minimum distance is determined and added to the Moon’s position vector.
The simulation time for this model was relatively much shorter as it can be seen from 9,407 seconds and 8.810^8 seconds. This is because when the bodies collide, the function gives an error and stops simulating the movements. Hence, we can see that a small distance causes bodies to collide.
The Parametric Plot for the New Moon at Enough Distances
The Parametric Plot for the New Moon at Enough Distances
A distance at the right interval is determined, added up, and made into a dataset again.
When the bodies are at the right interval, they don’t collide as can be seen from the simulation time, which is relatively large with 10^5 seconds. Also, the New Moon revolves around the Moon which can mean they are in the orbit.
The Parametric Plot for the New Moon at High Speeds
The Parametric Plot for the New Moon at High Speeds
For the velocity plot, I got the distance from the last plot as it is at enough distance but with a higher velocity than it should be. Below, we can see that, for this distance, the maximum velocity should be around 1,298 m/s. Therefore, I got a larger value from that.
As the velocity needs to be in vector form, the Moon’s velocity is multiplied with a number so that the net velocity is higher than the maximum velocity.
It can be seen from this plot that even though the bodies start at enough distance to each other, they don’t stay in the orbit and, as a result, escape because the velocity is higher than the escape velocity. They also don’t collide quickly like in the small distance plot.
Conclusion & Further Research
Conclusion & Further Research
For the Moon to have its own Moon in the 4-Body Problem, I examined some parameters: the mass, velocity and position. The maximum distance of the New Moon is obtained as around 3.0510^7 meters from the Hill Radius while the minimum distance depends on the mass of the New Moon from the Roche Limit. For the velocity, the maximum velocity it can have is determined to depend on its mass from the escape velocity formula.
According to studies, the maximum mass needs to be 3 percent of the mass of the Moon; we obtained the minimum mass from the equation which states that the maximum distance needs to be bigger than minimum distance as it only depends on the New Moon’s mass.
After obtaining the values, I made a simulation, plots and parametric plots. From the parametric plots, we saw that when the distance is too low, the two bodies collide shortly after while if they are too far away, they don’t stay in the orbit. Also, for the velocity, we can see that if the velocity is higher than the escape velocity, they don’t stay in orbit again.
For further research, instead of determining if the bodies are in orbit; eccentricity and conservation of momentum and energy can be used. As the eccentricity of a body in orbit is less than 1, it can be checked if the body is in the orbit. Additionally, energy and momentum are conserved when the body is in the orbit; hence, this can be determined from the conservation laws.
According to studies, the maximum mass needs to be 3 percent of the mass of the Moon; we obtained the minimum mass from the equation which states that the maximum distance needs to be bigger than minimum distance as it only depends on the New Moon’s mass.
After obtaining the values, I made a simulation, plots and parametric plots. From the parametric plots, we saw that when the distance is too low, the two bodies collide shortly after while if they are too far away, they don’t stay in the orbit. Also, for the velocity, we can see that if the velocity is higher than the escape velocity, they don’t stay in orbit again.
For further research, instead of determining if the bodies are in orbit; eccentricity and conservation of momentum and energy can be used. As the eccentricity of a body in orbit is less than 1, it can be checked if the body is in the orbit. Additionally, energy and momentum are conserved when the body is in the orbit; hence, this can be determined from the conservation laws.
Acknowledgments
Acknowledgments
◼
Macy Levin
References
References
◼
Burman, A. (2022). Stability in a Star-Planet-Moon three-body system. https://community.wolfram.com/groups/-/m/t/2575965.
◼
Marvellous, V. V., & Roy, A. E. (1981). Hill stability of satellite orbits. Celestial mechanics, 23(3), 269-275.
Cruz, J., & Dilão, R. (2020). Stability of circular orbits in the three body problem (Doctoral dissertation, Dissertation, Department of Physics, Tecnio Lisboa).
Cruz, J., & Dilão, R. (2020). Stability of circular orbits in the three body problem (Doctoral dissertation, Dissertation, Department of Physics, Tecnio Lisboa).
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Determining the conditions for the moon to have its own moon
by Damla Temur
Wolfram Community, STAFF PICKS, July 12, 2024
https://community.wolfram.com/groups/-/m/t/3215701
by Damla Temur
Wolfram Community, STAFF PICKS, July 12, 2024
https://community.wolfram.com/groups/-/m/t/3215701