Gravity is one of the most fundamental forces which mediates the interaction of all the mass in the universe. Gravity also plays a very important role in the formation and evolution of everything from planetary systems to galaxies. The life of a star starts in a molecular cloud where dense regions of gas and dust collapse to form young stars. Many of these stars are surrounded by protoplanetary disks, from which planetary systems can eventually emerge. However, star formation is rarely an isolated process! Stars are typically born in clusters, which often results in many ejections and multi-stellar systems. Understanding how these systems evolve and interact over time is a really interesting problem, and because the collective motion of multiple bodies gets too difficult to solve merely analytically, computational astrophysics can help solve and visualize it in a better manner.
Introduction
Introduction
To study these interactions, I have used N-body simulations, which numerically model the gravitational forces acting between many stars over time.
This project has a progressive modeling approach; I started with a simple N-body simulation before introducing more realistic conditions. A small stellar cluster is generated at first using randomly assigned masses, positions and velocities to observe the fundamental gravitational interaction of the system. Then, I replace the arbitrary values with realistic masses from the Kroupa Initial Mass Function, initial separations, and initial velocities.
My main objective is to analyze and determine if binary star systems can be identified automatically by simply analyzing pairwise simulation data as the cluster evolves.
I took a brute-force approach, enumerating through all possible pairs to check if any are binary according to the binding energy criterion which states that the sum of the kinetic and potential energy of the pair must be negative. Also, I constructed graphs containing edges between binary candidates which provides a potential route to identifying trinary and/or higher-order systems.
Throughout the project some more experiments were conducted in order to evaluate how the number of bodies and the simulation time impacts the computational cost.
This project has a progressive modeling approach; I started with a simple N-body simulation before introducing more realistic conditions. A small stellar cluster is generated at first using randomly assigned masses, positions and velocities to observe the fundamental gravitational interaction of the system. Then, I replace the arbitrary values with realistic masses from the Kroupa Initial Mass Function, initial separations, and initial velocities.
My main objective is to analyze and determine if binary star systems can be identified automatically by simply analyzing pairwise simulation data as the cluster evolves.
I took a brute-force approach, enumerating through all possible pairs to check if any are binary according to the binding energy criterion which states that the sum of the kinetic and potential energy of the pair must be negative. Also, I constructed graphs containing edges between binary candidates which provides a potential route to identifying trinary and/or higher-order systems.
Throughout the project some more experiments were conducted in order to evaluate how the number of bodies and the simulation time impacts the computational cost.
Basic N-Body Simulation
Basic N-Body Simulation
The simulation begins by generating a set of initial conditions for the particles. Each particle is represented as an association containing a mass of 1 unit, and its is placed randomly in a three-dimensional position, and it is randomly assigned a three-dimensional velocity, with both positions and velocities sampled uniformly from the range [−1,1]:
Initializes the simulation by creating a list of 10 equal-mass stars with randomly assigned 3D positions and velocities.
In[]:=
data=Table[<|"Mass"->1,"Position"->RandomReal[{-1,1},3],"Velocity"->RandomReal[{-1,1},3]|>,10];
Random samples of the generated particles are then displayed using Dataset to verify that the initial conditions have been created correctly. The NBodySimulation function is subsequently used to simulate the motion of the particles under an inverse-square force law (1/r^2), which models gravitational interactions, over a time intervals of 5 units.
Displays a random sample of three stars from the generated dataset to verify the assigned masses, positions, and velocities.
In[]:=
Dataset[RandomSample[data,3]]
Out[]=
I simulated the motion of the particles under an inverse-square force law (1/r^2), which models gravitational interactions, over a time interval of 5 units:
nBodySim=NBodySimulation["InverseSquare",data,5];
It visualizes the trajectories of all particles by plotting their three-dimensional positions as functions of time. This provides a clear representation of how the particles evolve and interact throughout the simulation.
Extracts the position of every star from the simulation and plots their trajectories in three dimensions over the chosen time interval.
ParametricPlot3D[Evaluate[nBodySim[All,"Position",t]],{t,0,5},Boxed->False,Axes->True]
Out[]=
Realistic N-Body Simulation
Realistic N-Body Simulation
This section replaces the inverse square law with Newtonian gravity, and uses realistic masses, separations and velocities that one might expect in a proto-stellar cluster. The inverse square law essentially was using the formula 1/r^2 which gives a vague approximation of how the bodies would interact, but the use of “Newtonian” takes into account the true value of the gravitational constant. I also simulated over a much longer timescale years), by which time one would expect to see binary/trinary systems forming. The resulting trajectories help us see how gravitational interactions become more complex as the stars exchange energy and interact with each other.
(-
5
10
6
10
The function below is designed for repeat N-body experiments in which the number of bodies and the simulation time are inputs. Each star is assigned a realistic mass, a realistic positions, and a realistic velocity using astronomical units, before the system simulates it under the Newtonian gravity. The computation time is recorded in order to gauge what can be feasibly in the time that I had.
Defines a reusable function for running N-body simulations with different numbers of stars and simulation times. The function generates the initial conditions it executes the simulation it measures the runtime, and stores the results for later analysis.
ClearAll[runNBodyExperiment]runNBodyExperiment[nBodies_Integer,simTime_Quantity]:=Module[{data,simulationTime,time,sim,result},data=Table[<|"Mass"->Quantity[1,"SolarMass"],"Position"->Quantity[RandomReal[{-0.1,0.1},3],"Parsecs"],"Velocity"->Quantity[RandomReal[{-1,1},3],"Kilometers"/"Seconds"]|>,nBodies];simulationTime=QuantityMagnitude[UnitConvert[simTime,"Seconds"]];{time,sim}=AbsoluteTiming[NBodySimulation["Newtonian",data,simulationTime]];Echo[Style["Simulation complete in "<>ToString[time]<>" s","Text"]];result=<|"Result"->sim,"Bodies"->nBodies,"SimTime"->simTime,"RunTime"->Quantity[time,"Seconds"]|>;result]
The following executes the runNBodyExperiment function using 20 bodies and a simulation time of 1 million years. The quantity function specifies the simulation time in physical units. The output of this is saved as a realisticSim variable for using it in future analysis.
This code runs the N-body simulation with 20 stars for a simulation time of 10^6 years and stores the results in realisticSim for later analysis and visualization.
realisticSim=runNBodyExperiment[20,Quantity[,"Years"]]
6
10
»
Simulation complete in 7.61009 s
Out[]=
ResultNBodySimulationData
,Bodies20,SimTime,RunTime
The code visualizes the trajectories of every single star which is generated by the simulation. The simulation is first converted from years to seconds so it can be used as the plotting variable and it extracts the position of each star from the NBodySimulation.
Extracts the trajectories of all stars from the simulation and plots their motion in three dimensions over the full simulation time, adding labels to identify each star more easily.
tMax=QuantityMagnitude[UnitConvert[realisticSim["SimTime"],"Seconds"]];paths=realisticSim["Result"][All,"Position",t];labelPoints=TablerealisticSim["Result"]i,"Position",RandomReal,tMax,{i,1,Length[paths]};ShowParametricPlot3DEvaluate[paths],{t,0,tMax},,Graphics3D[MapThread[Text[Style[#2,12],#1]&,{QuantityMagnitude/@labelPoints,Range[Length[paths]]}]],Boxed->False
tMax
2
Out[]=
Conduct a Series of Experiments
Conduct a Series of Experiments
I performed a systematic series of experiments to explore how my simulation changes as I varied the inputs; I can then explore how each input affects the computation time for the N-body simulation. The collected data is then analyzed using plots, that makes it possible to see how the runtime increases with both the size of the stars system and the duration of the simulation. This matrix summarizes the design-of-experiment (DOE):
Generates all combinations of the selected numbers of bodies and simulation times, creating a list of input parameters for the series of N-body experiments.
This creates all of the input combinations which would be used for the experimentation an evaluation.
This automatically runs simulation using the combination of bodies and simulation times generated in the last section. Each experiment in this execution is order and the results and the runtime are stored together in a single dataset. This gives a complete collection of simulations that can be compared and analyzed in the coming section.
This runs the N-body simulation for every combination of input parameters, stores the results in a table, and displays a random sample to verify that the simulations were completed successfully.
Here I grouped the simulation results by simulation time and plots the runtime as a function of the number of bodies. The curves here show the different simulation durations, and this makes it easier to compare how the computational performance changes.
Grouping of the simulation results by simulation time, extracts the number of bodies and corresponding runtimes, and plots each simulation duration as a separate curve to compare how runtime scales with system size.
The following code plots the runtime of the simulation as a function of the simulation time for each number of bodies.
The following code plots the runtime of the simulation as a function of the simulation time for each number of bodies.
This shows the relationship between the number of bodies, the simulation time and the running time using a 3D scatter plot. By representing all the three variables in a single figure it becomes easier to see the overall trend.
This shows the relationship between the number of bodies, the simulation time and the running time using a 3D scatter plot. By representing all the three variables in a single figure it becomes easier to see the overall trend.
In this plot, the size of the bubble represents the computation time:
The experiments in this section show that the runtime of an N-body simulation increases as both the number of bodies and the simulation time increases. Amongst the two variables the number of bodies has a greater influence on the execution time as every additional star brings more gravitational pull/interaction in the system which needs to be computed.
Stellar Mass Function
Stellar Mass Function
The real stellar systems do not have stars with exactly identical masses. Instead, the low mass stars are generally significantly a lot more common than high mass stars. In order to create a more physically realistic situation the simulation incorporates initial mass functions (IMFs) which describe the probability distribution of stellar masses at the time of formation. In this section I have compared the Salpeter IMF with the Kroupa IMF.
Salpeter IMF
Salpeter IMF
Defines the Salpeter Initial Mass Function as a normalized probability distribution and generates a random sample of stellar masses from it.
I took random samples from this distribution and displayed the results as a histogram, this illustrates that low-mass stars are more dominant in stellar population while the massive stars are comparatively rarer.
Kroupa IMF
Kroupa IMF
The Kroupa IMF is like an extension of the Salpeter model by using different power-law slopes over different mass ranges, this gives a closer match to the observations of real stellar population observations. A piecewise function and random stellar masses are generated to produce a histogram. The piecewise function is converted to a probability distribution and it is used to generate random stellar masses.
Defines the Kroupa IMF a piecewise function with different power-law slopes for different stellar mass ranges. The piecewise function is then converted into a normalized probability distribution, allowing realistic stellar masses to be randomly generated for the simulation.
The histogram shows that low-mass stars are far more common than high mass stars, which is consistent with what I observe with real stellar populations:
The code samples stellar masses from the Kroupa Initial Mass Function and displays their distribution as a logarithmic histogram.
Simulate Using Kroupa IMF
Simulate Using Kroupa IMF
The N-body simulation is paired with the stars of different mass rather than just assigning the same mass. The masses are randomly sampled from the Kroupa distribution, though the initial positions and velocities are generated beforehand. This way the simulation better represents real stellar systems, where low mass stars are relatively rare. The resulting trajectories can then be compared in the equal mass simulation to see how a realistic mass distribution influences the evolution of the system. Each star is randomly assigned a mass from the Kroupa IMF, along with random initial positions and velocities.
Initializes the simulation with 20 stars whose masses are sampled from the Kroupa IMF.
Here are the masses that were sampled:
The following code visualizes the trajectories of all the stars, it first determines the total simulation time and extracts the position of every star as a function of time, then a random point is selected along each trajectory to place a label, it makes the stars easier to identify. Finally, the trajectories are displayed as a 3D parametric plot with labels to distinguish the paths of the individual stars.
I extracted the trajectories of all stars from the Kroupa simulation and plots their motion in three dimensions, adding labels to identify each star.
2D Slices
2D Slices
Although the N-body simulation can be visualized in 3D, interpreting the motion of the individual stars can be a bit difficult because of the overlapping trajectories. To give a clearer view of the system, the simulation is projected onto the 3 plane. These 2D slices make it easier to examine the spatial distribution of the stars and identify patterns that might not be immediately visible in the 3D visualization.
Extracts the simulated trajectories and projects them onto a selected two-dimensional plane (XY, YZ or XZ).
Plotted the projected trajectories of the stars in the XY, YZ and XZ planes, allowing the structure of the system to be viewed from different perspectives.
Searching for Binaries
Searching for Binaries
The previous sections showed how stars move and how they interact under the gravitational influence. The next step is to determine if any of these interactions help in the formation of gravitationally bound binary systems. In this section, the binding energy criterion is used to test every pair of stars in the simulation. For all of the pairs the total energy is calculated by combining their kinetic and gravitational potential energies. If the total energy is negative then the pair is considered gravitationally bound. The bound pairs are then represented as a graph to identify binary systems, and their trajectories and binding energies are analyzed to determine whether these systems remain stable throughout the simulation.
Binding Energy Criterion
Binding Energy Criterion
The following code computes the total binding energy of each pair of stars and applies the binding energy criterion E < 0 to determine whether the pair is gravitationally bound. I defined a function for calculating the reduced mass of a pair of stars. This code calculates the reduced mass of each pair of stars for use in the binding energy calculations.
Definition of the function for calculating the reduced mass:
Calculated the relative kinetic energy of each pair of stars using their reduced mass and relative velocity:
Calculates the gravitational potential energy between each pair of stars using their masses and separation distance:
Here, I have calculate the total binding energy of a pair of stars by adding the relative kinetic energy and the gravitational potential energy of the two stars.
The binding energy criterion is applied to determine whether a pair of stars is gravitationally bound, if the total energy is negative, the stars are classified as gravitationally bound; otherwise, they are considered unbound.
Examination of every possible pair of stars and recording whether they are gravitationally bound over time, this is really important for data analysis with larger set of real data. It helps use determine binaries without any visualization of their evolution trajectories.
Uses Subsets to generate every possible pair of stars, evaluates each pair at the specified time steps, and stores the binding energy together with whether the pair is gravitationally bound.
Then I had created a function which creates a graph with edges between bodies that are found to be gravitationally bound. The criteria for being gravitationally bound is that at some instance, the pair energy is negative. I would like to make this criteria more robust in future.
Filters the results to keep only gravitationally bound pairs, then uses Graph and UndirectedEdge to connect the stars in each bound pair, producing a graph that highlights candidate binary and multiple-star systems.
Scans the simulation for gravitationally bound pairs using the selected time steps and then constructs a graph connecting the bound pairs, allowing candidate binary and multiple-star systems to be identified.
The graph identifies candidate gravitationally bound systems in the simulation. Binary systems appear as connected pairs of stars, while larger connected components indicate possible triple or higher-order systems. These candidates are analyzed further by examining their trajectories and orbital behavior to assess whether they remain gravitationally bound throughout the simulation.
In this section I had choose bodies 14 and 17 to test whether my scanning function is successful.
Extracts the connected components of the graph, with each component representing a candidate binary or multiple-star system identified from the binding-energy analysis.
The trajectories of stars 14 and 17 closely follow one another throughout the simulation, remaining at a similar separation, so it looks like a promising candidate binary system. Furthermore, I can plot the total energy of a this pair of stars over time, providing more evidence that I have a good binary candidate:
Evaluates the total binding energy of the selected pair at regular time intervals and plots its evolution throughout the simulation.
The plot shows the total binding energy of the selected pair as a function of time. Although the energy varies during the simulation, it remains negative throughout, indicating that the two stars stay gravitationally bound. This provides further evidence that the pair is a stable candidate binary system rather than a temporary close encounter. Hence, I concluded that the function works.
Periodic Model
Periodic Model
A periodic model is a mathematical function that describes quantities that repeat at regular intervals over time, making it useful for modelling orbital motion. I had fit a PeriodicModel to the relative position data of bodies 14 and 17, to perhaps derive an orbital period of this binary system and see if they overlap.
PeriodicModel is fitted to the relative separation of bodies 14 and 17 to estimate the orbital period of this candidate binary system and assess how well a periodic function describes its motion.
The contribution of a third body (trinary system) may contribute additional frequencies that a more sophisticated PeriodicModel could find.
The fitted periodic model closely follows the simulated separation between the two stars, indicating that the candidate binary exhibits approximately periodic orbital motion. Small differences between the curves are expected because the stars are also influenced by the gravitational pull of the other bodies in the simulation. I can also obtain a ModelFitReport describing the quality of the model:
Passes the relative position data and the fitted PeriodicModel to ModelFit with the “Report” option, which computes goodness-of-fit statistics.
Finally, I can find the orbital period (in years) of my candidate binary system:
The period on this binary is about 360439 years, this is a really long period and it might be indicating that it is a wide binary system, this information is not enough no conclude that it is a wide binary, but it is one of the indications that it can be a wide and further tests should be run on this to confirm the exact classification.
Conclusion
Conclusion
In this project, I demonstrated how computational modelling can be used to analyze the behavior of stellar systems as they interact gravitationally in ways that are not possible to visualize through analytical methods alone. Using both the visual simulation output as well as automated analysis, I identified promising binary candidates using the binding energy criterion. I estimated their orbital period and examined their orbital behavior by fitting a periodic model. The PeriodicModel function is a new feature in Wolfram Language, and I am glad that I could use it in this project. Therefore, computational methods proved to make the automated analysis much more efficient. Finally, I stress tested the function NBodySimulation during the project in order to understand the computational capacities of the function.
Future Works
Future Works
The follow up projects to this might include:
1) Extending the simulation to study a wider range of gravitationally bound systems, including wide binaries, compact binaries, and higher order systems.
2) Expanding the automated analysis to classify different systems will help in studying larger data sets more efficiently.
3) Increasing the number of simulated stars in order to get closer real observations.
4) Incorporating additional physical process in order improve the realism of the models.
1) Extending the simulation to study a wider range of gravitationally bound systems, including wide binaries, compact binaries, and higher order systems.
2) Expanding the automated analysis to classify different systems will help in studying larger data sets more efficiently.
3) Increasing the number of simulated stars in order to get closer real observations.
4) Incorporating additional physical process in order improve the realism of the models.
References
References
◼
Binary Formation Timescale: https://arxiv.org/pdf/2605.16657v1, https://astrobites.org/2016/04/26/the-formation-of-wide-binaries
Acknowledgements
Acknowledgements
First and foremost, I am deeply thankful to my mentor, Joseph Brennan. He has guided me throughout the project and he was always there to help me whenever I needed any help. I am also grateful to Stephen Wolfram for helping me decide the project and its directions. I would like to thank Divya Thota and Arnay Garhyan for reviewing my project and additionally I would also like to thank the all the teaching assistants and mentors for their assistance in my project. I am also very thankful to the program directors, Megan, Rory, Cyrus, and Eryn and everyone who made this program and project possible.
CITE THIS NOTEBOOK
CITE THIS NOTEBOOK
Formation of binary star systems
by Katyayani Singh
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3753850
by Katyayani Singh
Wolfram Community, STAFF PICKS, July 9, 2026
https://community.wolfram.com/groups/-/m/t/3753850