Runaway collisions in dense clusters may lead to the formation of supermassive black hole (SMBH) seeds, and this process can be further enhanced by accretion, as recent models of SMBH seed formation in Population III star clusters have shown. This may explain the presence of supermassive black holes already at high redshift, $z>6$. However, in this context, mass loss during collisions was not considered and could play an important role for the formation of the SMBH seed. Here, we study the effect of mass loss, due to collisions of protostars, in the formation and evolution of a massive object in a dense primordial cluster. We consider both constant mass loss fractions as well as analytic models based on the stellar structure of the collision components. Our calculations indicate that mass loss can significantly affect the final mass of the possible SMBH seed. Considering a constant mass loss of 5% for every collision, we can lose between 60-80% of the total mass that is obtained if mass loss were not considered. Using instead analytical prescriptions for mass loss, the mass of the final object is reduced by 15-40%, depending on the accretion model for the cluster we study. Altogether, we obtain masses of the order of $10^4M_{\odot}$, which are still massive enough to be SMBH seeds.

Chaos is present in most stellar dynamical systems and manifests itself through the exponential growth of small perturbations. Exponential divergence drives time irreversibility and increases the entropy in the system. A numerical consequence is that integrations of the N-body problem unavoidably magnify truncation and rounding errors to macroscopic scales. Hitherto, a quantitative relation between chaos in stellar dynamical systems and the level of irreversibility remained undetermined. In this work we study chaotic three-body systems in free fall initially using the accurate and precise N-body code Brutus, which goes beyond standard double-precision arithmetic. We demonstrate that the fraction of irreversible solutions decreases as a power law with numerical accuracy. This can be derived from the distribution of amplification factors of small initial perturbations. Applying this result to systems consisting of three massive black holes with zero total angular momentum, we conclude that up to five percent of such triples would require an accuracy of smaller than the Planck length in order to produce a time-reversible solution, thus rendering them fundamentally unpredictable.

Since its formulation by Sir Isaac Newton, the problem of solving the equations of motion for three bodies under their own gravitational force has remained practically unsolved. Currently, the solution for a given initialization can only be found by performing laborious iterative calculations that have unpredictable and potentially infinite computational cost, due to the system's chaotic nature. We show that an ensemble of solutions obtained using an arbitrarily precise numerical integrator can be used to train a deep artificial neural network (ANN) that, over a bounded time interval, provides accurate solutions at fixed computational cost and up to 100 million times faster than a state-of-the-art solver. Our results provide evidence that, for computationally challenging regions of phase-space, a trained ANN can replace existing numerical solvers, enabling fast and scalable simulations of many-body systems to shed light on outstanding phenomena such as the formation of black-hole binary systems or the origin of the core collapse in dense star clusters.

To explain the observed population of supermassive black holes at z~7, very massive seed black holes or, alternatively, super-Eddington scenarios are needed to reach final masses of the order of 10^9 solar masses. A popular explanation for massive seeds has been the direct collapse model, which predicts the formation of a single massive object due to the direct collapse of a massive gas cloud. Simulations over the last years have however shown that such a scenario is very difficult to achieve. A realistic model of black hole formation should therefore take fragmentation into account, and consider the interaction between stellar-dynamical and gas-dynamical processes. We present here numerical simulations pursued with the AMUSE code, employing an approximate treatment of the gas. Based on these simulations, we show that very massive black holes of 10^4-10^5 solar masses may form depending on the gas supply and the accretion onto the protostars.

The search for extrasolar planets in the past decades has shown that planets abound in the Solar neighborhood. While we are still missing an Earth twin, the forthcoming space missions and ground-based instrumentation are already driven to achieve this goal. But, in order to fully understand the conditions for life appearing in the Solar System, we still miss some pieces of the planetary system jigsaw puzzle, namely a deeper understanding of the minor bodies. Trojans, moons, and comets are tracers of the formation and evolution processes of planetary systems. These missing pieces are also critical to understand the emergence and evolution of life over millions of years. With the large crop of planetary systems discovered so far and yet to be detected with the forthcoming missions, the hunt for minor bodies in extrasolar systems is a natural continuation of our search for real Solar System- and, in particular, Earth- analogs. This white paper is focused on detection of these minor components and their relevance in the emergence, evolution and survival of life.

Numerical solutions to Newton's equations of motion for chaotic self gravitating systems of more than 2 bodies are often regarded to be irreversible. This is due to the exponential growth of errors introduced by the integration scheme and the numerical round-off in the least significant figure. This secular growth of error is sometimes attributed to the increase in entropy of the system even though Newton's equations of motion are strictly time reversible. We demonstrate that when numerical errors are reduced to below the physical perturbation and its exponential growth during integration the microscopic reversibility is retrieved. Time reversibility itself is not a guarantee for a definitive solution to the chaotic N-body problem. However, time reversible algorithms may be used to find initial conditions for which perturbed trajectories converge rather than diverge. The ability to calculate such a converging pair of solutions is a striking illustration which shows that it is possible to compute a definitive solution to a highly unstable problem. This works as follows: If you (i) use a code which is capable of producing a definitive solution (and which will therefore handle converging pairs of solutions correctly), (ii) use it to study the statistical result of some other problem, and then (iii) find that some other code produces a solution S with statistical properties which are indistinguishable from those of the definitive solution, then solution S may be deemed veracious.

Collisions were suggested to potentially play a role in the formation of massive stars in present day clusters, and have likely been relevant during the formation of massive stars and intermediate mass black holes within the first star clusters. In the early Universe, the first stellar clusters were particularly dense, as fragmentation typically only occurred at densities above $10^9$cm$^{-3}$, and the radii of the protostars were enhanced due to the larger accretion rates, suggesting a potentially more relevant role of stellar collisions. We present here a detailed parameter study to assess how the number of collisions as well as the mass growth of the most massive object depends on the properties of the cluster, and we characterize the time evolution with three effective parameters, the time when most collisions occur, the duration of the collisions period, as well as the normalization required to obtain the total number of collisions. We apply our results to typical Population III (Pop.III) clusters of about $1000$M$_\odot$, finding that a moderate enhancement of the mass of the most massive star by a factor of a few can be expected. For more massive Pop.III clusters as expected in the first atomic cooling halos, we expect a more significant enhancement by a factor of $15-32$. We therefore conclude that collisions in massive Pop.III clusters were likely relevant to form the first intermediate mass black holes.

Models aiming to explain the formation of massive black hole seeds, and in particular the direct collapse scenario, face substantial difficulties. These are rooted in rather ad hoc and fine-tuned initial conditions, such as the simultaneous requirements of extremely low metallicities and strong radiation backgrounds. Here we explore a modification of such scenarios where a massive primordial star cluster is initially produced. Subsequent stellar collisions give rise to the formation of massive (10^4 - 10^5 solar mass) objects. Our calculations demonstrate that the interplay between stellar dynamics, gas accretion and protostellar evolution is particularly relevant. Gas accretion onto the protostars enhances their radii, resulting in an enhanced collisional cross section. We show that the fraction of collisions can increase from 0.1-1% of the initial population to about 10% when compared to gas-free models or models of protostellar clusters in the local Universe. We conclude that very massive objects can form in spite of initial fragmentation, making the first massive protostellar clusters viable candidate birth places for observed supermassive black holes.

Observations of the Orion-A integral shaped filament (ISF) have shown indications of an oscillatory motion of the gas filament. This evidence is based on both the wave-like morphology of the filament as well as the kinematics of the gas and stars, where the characteristic velocities of the stars require a dynamical heating mechanism. As proposed by Stutz and Gould (2016), such a heating mechanism (the "Slingshot") may be the result of an oscillating gas filament in a gas-dominated (as opposed to stellar-mass dominated) system. Here we test this hypothesis with the first stellar-dynamical simulations in which the stars are subjected to the influence of an oscillating cylindrical potential. The accelerating, cylindrical background potential is populated with a narrow distribution of stars. By coupling the potential to N-body dynamics, we are able to measure the influence of the potential on the stellar distribution. The simulations provide evidence that the slingshot mechanism can successfully reproduce several stringent observational constraints. These include the stellar spread (both in projected position and in velocity) around the filament, the symmetry in these distributions, and a bulk motion of the stars with respect to the filament. Using simple considerations we show that star-star interactions are incapable of reproducing these spreads on their own when properly accounting for the gas potential. Thus, properly accounting for the gas potential is essential for understanding the dynamical evolution of star forming filamentary systems in the era of Gaia.

White dwarfs (WDs) are powerful tools to study the evolutionary history of stars and binaries in the Galaxy. But do we understand their multiplicity from a theoretical point of view? This can be tested by a comparison with the sample of WDs within 20 pc, which is minimally affected by selection biases. From the literature, we compile the available information of the local WD sample with a particular emphasis on their multiplicity, and compare this to synthetic models of WD formation in single stars and binaries. As part of our population synthesis approach, we also study the effect of different assumptions concerning the star formation history, binary evolution, and the initial distributions of binary parameters. We find that the observed space densities of single and binary WDs are well reproduced by the models. The space densities of the most common WD systems (single WDs and unresolved WD-MS binaries) are consistent within a factor two with the observed value. We find a discrepancy only for the space density of resolved double WDs. We exclude that observational selection effects, fast stellar winds, or dynamical interactions with other objects in the Milky Way explain this discrepancy. We find that either the initial mass ratio distribution in the Solar neighbourhood is biased towards low mass-ratios, or more than ten resolved DWDs have been missed observationally in the 20 pc sample. Furthermore, we show that the low binary fraction of WD systems (~25%) compared to Solartype MS-MS binaries (~50%) is consistent with theory, and is mainly caused by mergers in binary systems, and to a lesser degree by WDs hiding in the glare of their companion stars. Lastly, Gaia will dramatically increase the size of the volume-limited WD sample, detecting the coolest and oldest WDs out to 50 pc. We provide a detailed estimate of the number of single and binary WDs in the Gaia sample.

According to Munoz-Gutierrez et al. (2015) the orbit of comet 1P/Halley is chaotic with a surprisingly small Lyapunov time scale of order its orbital period. In this work we analyse the origin of chaos in Halley's orbit and the growth of perturbations, in order to get a better understanding of this unusually short time scale. We perform N-body simulations to model Halley's orbit in the Solar System and measure the separation between neighbouring trajectories. To be able to interpret the numerical results, we use a semi-analytical map to demonstrate different growth modes, i.e. linear, oscillatory or exponential, and transitions between these modes. We find the Lyapunov time scale of Halley's orbit to be of order 300 years, which is significantly longer than previous estimates in the literature. This discrepancy could be due to the different methods used to measure the Lyapunov time scale. A surprising result is that next to Jupiter, also encounters with Venus contribute to the exponential growth in the next 3000 years. Finally, we note an interesting application of the sub-linear, oscillatory growth mode to an ensemble of bodies moving through the Solar System. Whereas in the absence of encounters with a third body the ensemble spreads out linearly in time, the accumulation of weak encounters can increase the lifetime of such systems due to the oscillatory behaviour.

We use self-consistent numerical simulations of the evolution and disruption of the Sun's birth cluster in the Milky Way potential to investigate the present-day phase space distribution of the Sun's siblings. The simulations include the gravitational N-body forces within the cluster and the effects of stellar evolution on the cluster population. In addition the gravitational forces due to the Milky Way potential are accounted for in a self-consistent manner. Our aim is to understand how the astrometric and radial velocity data from the Gaia mission can be used to pre-select solar sibling candidates. We vary the initial conditions of the Sun's birth cluster, as well as the parameters of the Galactic potential. We show that the disruption time-scales of the cluster are insensitive to the details of the non-axisymmetric components of the Milky Way model and we make predictions, averaged over the different simulated possibilities, about the number of solar siblings that should appear in surveys such as Gaia or GALAH. We find a large variety of present-day phase space distributions of solar siblings, which depend on the cluster initial conditions and the Milky Way model parameters. We show that nevertheless robust predictions can be made about the location of the solar siblings in the space of parallaxes ($\varpi$), proper motions ($\mu$) and radial velocities ($V_\mathrm{r}$). By calculating the ratio of the number of simulated solar siblings to that of the number of stars in a model Galactic disk, we find that this ratio is above 0.5 in the region given by: $\varpi \geq 5$mas, $4 \leq \mu \leq 6$masyr$^{-1}$, and $-2\leq V_\mathrm{r} \leq 0$kms$^{-1}$. Selecting stars from this region should increase the probability of success in identifying solar siblings through follow up observations [Abridged].

The general consensus in the N-body community is that statistical results of an ensemble of collisional N-body simulations are accurate, even though individual simulations are not. A way to test this hypothesis is to make a direct comparison of an ensemble of solutions obtained by conventional methods with an ensemble of true solutions. In order to make this possible, we wrote an N-body code called Brutus, that uses arbitrary-precision arithmetic. In combination with the Bulirsch--Stoer method, Brutus is able to obtain converged solutions, which are true up to a specified number of digits. We perform simulations of democratic 3-body systems, where after a sequence of resonances and ejections, a final configuration is reached consisting of a permanent binary and an escaping star. We do this with conventional double-precision methods, and with Brutus; both have the same set of initial conditions and initial realisations. The ensemble of solutions from the conventional simulations is compared directly to that of the converged simulations, both as an ensemble and on an individual basis to determine the distribution of the errors. We find that on average at least half of the conventional simulations diverge from the converged solution, such that the two solutions are microscopically incomparable. For the solutions which have not diverged significantly, we observe that if the integrator has a bias in energy and angular momentum, this propagates to a bias in the statistical properties of the binaries. In the case when the conventional solution has diverged onto an entirely different trajectory in phase-space, we find that the errors are centred around zero and symmetric; the error due to divergence is unbiased, as long as the time-step parameter, eta <= 2^(-5) and when simulations which violate energy conservation by more than 10% are excluded.

The conservation of energy, linear momentum and angular momentum are important drivers for our physical understanding of the evolution of the Universe. These quantities are also conserved in Newton's laws of motion under gravity \citep{Newton:1687}. Numerical integration of the associated equations of motion is extremely challenging, in particular due to the steady growth of numerical errors (by round-off and discrete time-stepping, \cite{1981PAZh....7..752B,1993ApJ...415..715G,1993ApJ...402L..85H,1994LNP...430..131M}) and the exponential divergence \citep{1964ApJ...140..250M,2009MNRAS.392.1051U} between two nearby solution. As a result, numerical solutions to the general N-body problem are intrinsically questionable \citep{2003gmbp.book.....H,1994JAM....61..226L}. Using brute force integrations to arbitrary numerical precision we demonstrate empirically that ensembles of different realizations of resonant 3-body interactions produce statistically indistinguishable results. Although individual solutions using common integration methods are notoriously unreliable, we conjecture that an ensemble of approximate 3-body solutions accurately represents an ensemble of true solutions, so long as the energy during integration is conserved to better than 1/10. We therefore provide an independent confirmation that previous work on self-gravitating systems can actually be trusted, irrespective of the intrinsic chaotic nature of the N-body problem.

We developed a Keplerian-based Hamiltonian splitting for solving the gravitational $N$-body problem. This splitting allows us to approximate the solution of a general $N$-body problem by a composition of multiple, independently evolved $2$-body problems. While the Hamiltonian splitting is exact, we show that the composition of independent $2$-body problems results in a non-symplectic non-time-symmetric first-order map. A time-symmetric second-order map is then constructed by composing this basic first-order map with its self-adjoint. The resulting method is precise for each individual $2$-body solution and produces quick and accurate results for near-Keplerian $N$-body systems, like planetary systems or a cluster of stars that orbit a supermassive black hole. The method is also suitable for integration of $N$-body systems with intrinsic hierarchies, like a star cluster with primordial binaries. The superposition of Kepler solutions for each pair of particles makes the method excellently suited for parallel computing; we achieve $\gtrsim 64\%$ efficiency for only $8$ particles per core, but close to perfect scaling for $16384$ particles on a $128$ core distributed-memory computer. We present several implementations in \texttt{Sakura}, one of which is publicly available via the AMUSE framework.