Search

We present {\sc Venice}, an operator splitting algorithm to integrate a numerical model on a hierarchy of timescales. {\sc Venice} allows a wide variety of different physical processes operating a different scales to be coupled on individual and adaptive time-steps. It therewith mediates the development of complex multi-scale and multi-physics simulation environments with a wide variety of independent components. The coupling between various physical models and scales is dynamic, and realized through (Strang) operators splitting using adaptive time steps. We demonstrate the functionality and performance of this algorithm using astrophysical models of a stellar cluster, first coupling gravitational dynamics and stellar evolution, then coupling internal gravitational dynamics with dynamics within a galactic background potential, and finally combining these models while also introducing dwarf galaxy-like perturbers. These tests show numerical convergence for decreasing coupling timescales, demonstrate how {\sc Venice} can improve the performance of a simulation by shortening coupling timescales when appropriate, and provide a case study of how {\sc Venice} can be used to gradually build up and tune a complex multi-physics model. Although the examples couple complete numerical models, {\sc Venice} can also be used to efficiently solve systems of stiff differential equations., Comment: Accepted for publication by A&A

Located at the core of the Galactic Centre, the S-star cluster serves as a remarkable illustration of chaos in dynamical systems. The long-term chaotic behaviour of this system can be studied with gravitational $N$-body simulations. By applying a small perturbation to the initial position of star S5, we can compare the evolution of this system to its unperturbed evolution. This results in the two solutions diverging exponentially, defined by the separation in position space $\delta_{r}$, with an average Lyapunov timescale of $\sim$420 yr, corresponding to the largest positive Lyapunov exponent. Even though the general trend of the chaotic evolution is governed in part by the supermassive black hole Sagittarius $\rm A^{*}$ (Sgr $\rm A^{*}$), individual differences between the stars can be noted in the behaviour of their phase-space curves. We present an analysis of the individual behaviour of the stars in this Newtonian chaotic dynamical system. The individuality of their behaviour is evident from offsets in the position space separation curves of the S-stars and the black hole. We propose that the offsets originate from the initial orbital elements of the S-stars, where Sgr $\rm A^{*}$ is considered in one of the focal points of the Keplerian orbits. Methods were considered to find a relationship between these elements and the separation in position space. Symbolic regression turns out to provide the clearest diagnostics for finding an interpretable expression for the problem. Our symbolic regression model indicates that $\left\langle\delta_r\right\rangle \propto e^{2.3}$, implying that the time-averaged individual separation in position space increases rapidly with the initial eccentricity of the S-stars., Comment: Accepted for publication in A&A

Recently we estimated that about 5 percent of supermassive black hole triple systems are fundamentally unpredictable. These gargantuan chaotic systems are able to exponentially magnify Planck length perturbations to astronomical scales within their interaction timescale. These results were obtained in the zero angular momentum limit, which we naively expected to be the most chaotic. Here, we generalise to triple systems with arbitrary angular momenta by systematically varying the initial virial ratio. We find the surprising result that increasing the angular momentum enhances the chaotic properties of triples. This is not only explained by the longer life times, allowing for a prolonged exponential growth, but also the maximum Lyapunov exponent itself increases. For the ensemble of initially virialised triple systems, we conclude that the percentage of unpredictable supermassive black hole triples increases to about 30 percent. A further increase up to about 50 percent is reached when considering triples on smaller astrophysical scales. Fundamental unpredictability is thus a generic feature of chaotic, self-gravitating triple populations., Comment: 13 pages, 11 figures, 2 tables. Submitted to MNRAS

Simulating the evolution of the gravitational N-body problem becomes extremely computationally expensive as N increases since the problem complexity scales quadratically with the number of bodies. We study the use of Artificial Neural Networks (ANNs) to replace expensive parts of the integration of planetary systems. Neural networks that include physical knowledge have grown in popularity in the last few years, although few attempts have been made to use them to speed up the simulation of the motion of celestial bodies. We study the advantages and limitations of using Hamiltonian Neural Networks to replace computationally expensive parts of the numerical simulation. We compare the results of the numerical integration of a planetary system with asteroids with those obtained by a Hamiltonian Neural Network and a conventional Deep Neural Network, with special attention to understanding the challenges of this problem. Due to the non-linear nature of the gravitational equations of motion, errors in the integration propagate. To increase the robustness of a method that uses neural networks, we propose a hybrid integrator that evaluates the prediction of the network and replaces it with the numerical solution if considered inaccurate. Hamiltonian Neural Networks can make predictions that resemble the behavior of symplectic integrators but are challenging to train and in our case fail when the inputs differ ~7 orders of magnitude. In contrast, Deep Neural Networks are easy to train but fail to conserve energy, leading to fast divergence from the reference solution. The hybrid integrator designed to include the neural networks increases the reliability of the method and prevents large energy errors without increasing the computing cost significantly. For this problem, the use of neural networks results in faster simulations when the number of asteroids is >70., Comment: Accepted for publication in the Journal of Computational Physics

We investigate the chaotic behavior of the S-star cluster in the Galactic center using precise $N$-body calculations, free from round-off or discretization errors. Our findings reveal that chaos among the Galactic center S-stars arises from close encounters, particularly among pairs and near the massive central body. These encounters induce perturbations, causing sudden changes in the orbital energies of the interacting stars. Consequently, neighboring solutions experience roughly exponential growth in separation. We propose a theory of "punctuated chaos" that describes the S-star cluster's chaotic behavior. This phenomenon results from nearly linear growth in the separation between neighboring orbits after repeated finite perturbations. Each participating star's orbit experiences discrete, abrupt changes in energy due to the perturbations. The cumulative effect of these events is further amplified by the steady drift in orbital phase. In the Galactic center, perturbations originate from coincidental encounters occurring within a distance of $\aplt 100$\,au between at least two stars (in some cases, three stars). Our model satisfactorily explains the observed exponential growth in the 27 S-star cluster. We determine that the S-star system has a Lyapunov time scale of approximately 462 +/-74 years. For the coming millennium, chaos in the S-star cluster will be driven mainly by a few of the closest orbiting stars: S2, S5, S6, S8, S9, S14, S18, S31, S21, S24, S27, S29, and S38., Comment: Accepted for publication in MNRAS

Dynamical chaos is a fundamental manifestation of gravity in astrophysical, many-body systems. The spectrum of Lyapunov exponents quantifies the associated exponential response to small perturbations. Analytical derivations of these exponents are critical for understanding the stability and predictability of observed systems. This essay presents a new model for chaos in systems with eccentric and crossing orbits. Here, exponential divergence is not a continuous process but rather the cumulative effect of an ever-increasing linear response driven by discrete events at regular intervals, i.e., punctuated chaos. We show that long-lived systems with punctuated chaos can magnify Planck length perturbations to astronomical scales within their lifetime, rendering them fundamentally indeterministic., Comment: Essay written for the Gravity Research Foundation 2023 Awards for Essays on Gravitation (Honorable Mention). Accepted by International Journal of Modern Physics D

The fundamental laws of physics are time-symmetric, but our macroscopic experience contradicts this. The time reversibility paradox is partly a consequence of the unpredictability of Newton's equations of motion. We measure the dependence of the fraction of irreversible, gravitational N-body systems on numerical precision and find that it scales as a power law. The stochastic wave packet reduction postulate then introduces fundamental uncertainties in the Cartesian phase space coordinates that propagate through classical three-body dynamics to macroscopic scales within the triple's lifetime. The spontaneous collapse of the wave function then drives the global chaotic behavior of the Universe through the superposition of triple systems (and probably multi-body systems). The paradox of infinitesimal granularity then arises from the superposition principle, which states that any multi-body system is composed of an ensemble of three-body problems., Comment: Accepted for publication in IC M2 Physical Sciences AIP Conference Proceedings

The Laser Interferometer Space Antenna (LISA) will be a transformative experiment for gravitational wave astronomy, and, as such, it will offer unique opportunities to address many key astrophysical questions in a completely novel way. The synergy with ground-based and space-born instruments in the electromagnetic domain, by enabling multi-messenger observations, will add further to the discovery potential of LISA. The next decade is crucial to prepare the astrophysical community for LISA's first observations. This review outlines the extensive landscape of astrophysical theory, numerical simulations, and astronomical observations that are instrumental for modeling and interpreting the upcoming LISA datastream. To this aim, the current knowledge in three main source classes for LISA is reviewed; ultracompact stellar-mass binaries, massive black hole binaries, and extreme or intermediate mass ratio inspirals. The relevant astrophysical processes and the established modeling techniques are summarized. Likewise, open issues and gaps in our understanding of these sources are highlighted, along with an indication of how LISA could help making progress in the different areas. New research avenues that LISA itself, or its joint exploitation with upcoming studies in the electromagnetic domain, will enable, are also illustrated. Improvements in modeling and analysis approaches, such as the combination of numerical simulations and modern data science techniques, are discussed. This review is intended to be a starting point for using LISA as a new discovery tool for understanding our Universe.

The gravitational $N$-body problem, which is fundamentally important in astrophysics to predict the motion of $N$ celestial bodies under the mutual gravity of each other, is usually solved numerically because there is no known general analytical solution for $N>2$. Can an $N$-body problem be solved accurately by a neural network (NN)? Can a NN observe long-term conservation of energy and orbital angular momentum? Inspired by Wistom & Holman (1991)'s symplectic map, we present a neural $N$-body integrator for splitting the Hamiltonian into a two-body part, solvable analytically, and an interaction part that we approximate with a NN. Our neural symplectic $N$-body code integrates a general three-body system for $10^{5}$ steps without diverting from the ground truth dynamics obtained from a traditional $N$-body integrator. Moreover, it exhibits good inductive bias by successfully predicting the evolution of $N$-body systems that are no part of the training set., Comment: 7 pages, 2 figures, accepted for publication at the NeurIPS 2021 workshop "Machine Learning and the Physical Sciences"

We study the influence of relativity on the chaotic properties and dynamical outcomes of an unstable triple system; the Pythagorean three-body problem. To this end, we extend the Brutus N-body code to include Post-Newtonian pairwise terms up to 2.5 order, and the first order Taylor expansion to the Einstein-Infeld-Hoffmann equations of motion. The degree to which our system is relativistic depends on the scaling of the total mass (the unit size was 1 parsec). Using the Brutus method of convergence, we test for time-reversibility in the conservative regime, and demonstrate that we are able to obtain definitive solutions to the relativistic three-body problem. It is also confirmed that the minimal required numerical accuracy for a successful time-reversibility test correlates with the amplification factor of an initial perturbation. When we take into account dissipative effects through gravitational wave emission, we find that the duration of the resonance, and the amount of exponential growth of small perturbations depend on the mass scaling. For a unit mass <= 10 MSun, the system behavior is indistinguishable from the Newtonian case, and the resonance always ends in a binary and one escaping body. For a mass scaling up to 1e7 MSun, relativity gradually becomes more prominent, but the majority of the systems still dissolve. The first mergers start to appear for a mass of ~1e5 MSun, and between 1e7 MSun and 1e9 MSun all systems end prematurely in a merger. These mergers are preceded by a gravitational wave driven in-spiral. For a mass scaling >= 1e9 MSun, all systems result in a gravitational wave merger upon the first close encounter. Relativistic three-body encounters thus provide an efficient pathway for resolving the final parsec problem. The onset of mergers at the characteristic mass scale of 1e7 MSun potentially leaves an imprint in the mass function of supermassive black holes., Comment: Accepted by Phys Rev D. 10 pages, 6 figures

Context. As the importance of Gravitational Wave (GW) Astrophysics increases rapidly, astronomers in different fields and with different backgrounds can have the need to get a quick idea of which GW source populations can be detected by which detectors and with what measurement uncertainties. Aims. The GW-Toolbox is an easy-to-use, flexible tool to simulate observations on the GW universe with different detectors, including ground-based interferometers (advanced LIGO, advanced VIRGO, KAGRA, Einstein Telescope, and also customised designs), space-borne interferometers (LISA and a customised design), pulsar timing arrays mimicking the current working ones (EPTA, PPTA, NANOGrav, IPTA) and future ones. We include a broad range of sources such as mergers of stellar mass compact objects, namely black holes, neutron stars and black hole-neutron stars; and supermassive black hole binaries mergers and inspirals, Galactic double white dwarfs in ultra-compact orbit, extreme mass ratio inspirals and Stochastic GW backgrounds. Methods. We collect methods to simulate source populations and determine their detectability with the various detectors. The paper aims at giving a comprehensive description on the algorithm and functionality of the GW-Toolbox. Results. The GW-Toolbox produces results that are consistent with more detailed calculations of the different source classes and can be accessed with a website interface (gw-universe.org) or as a python package (https://bitbucket.org/radboudradiolab/gwtoolbox). In the future, it will be upgraded with more functionality., Comment: 21 pages, 26 figures, 4 tables. Submitted to Astronomy & Astrophysics; Website url: gw-universe.org; comments on v2: made some revisions according to referee's comments; added Cosmic Explorer; comments on v3: make some more clarification in the caveats section, on the waveform and uncertainty estimation

Context. The star formation rates and stellar densities found in young massive clusters suggest that these stellar systems originate from gas at densities n(H$_2$) $>10^7$ cm$^{-3}$. Until today, however, the physical characterization of this ultra high density material remains largely unconstrained in observations. Aims. We investigated the density properties of the star-forming gas in the OMC-1 region located in the vicinity of the Orion Nebula Cluster (ONC). Methods. We mapped the molecular emission at 652 GHz in OMC-1 as part of the APEX-SEPIA660 Early Science. Results. We detect bright and extended N$_2$H$^+$ (J=7-6) line emission along the entire OMC-1 region. Comparisons with previous ALMA data of the (J=1-0) transition and radiative transfer models indicate that the line intensities observed in this N$_2$H$^+$ (7-6) line are produced by large mass reservoirs of gas at densities n(H$_2$) $>10^7$ cm$^{-3}$. Conclusions. The first detection of this N$_2$H$^+$ (7-6) line at parsec-scales demonstrates the extreme density conditions of the star-forming gas in young massive clusters such as the ONC. Our results highlight the unique combination of sensitivity and mapping capabilities of the new SEPIA660 receiver for the study of the ISM properties at high frequencies., Comment: 8 pages, 5 Figures. Accepted for publication by A&A

The importance of computing in astronomy continues to increase, and so is its impact on the environment. When analyzing data or performing simulations, most researchers raise concerns about the time to reach a solution rather than its impact on the environment. Luckily, a reduced time-to-solution due to faster hardware or optimizations in the software generally also leads to a smaller carbon footprint. This is not the case when the reduced wall-clock time is achieved by overclocking the processor, or when using supercomputers. The increase in the popularity of interpreted scripting languages, and the general availability of high-performance workstations form a considerable threat to the environment. A similar concern can be raised about the trend of running single-core instead of adopting efficient many-core programming paradigms. In astronomy, computing is among the top producers of green-house gasses, surpassing telescope operations. Here I hope to raise the awareness of the environmental impact of running non-optimized code on overpowered computer hardware., Comment: Originated at EAS 2020 conference, sustainability session by https://astronomersforplanet.earth - published in Nature Astronomy, September 2020

Hierarchical code coupling strategies make it possible to combine the results of individual numerical solvers into a self-consistent symplectic solution. We explore the possibility of allowing such a coupling strategy to be non-intrusive. In that case, the underlying numerical implementation is not affected by the coupling itself, but its functionality is carried over in the interface. This method is efficient for solving the equations of motion for a self-gravitating system over a wide range of scales. We adopt a dedicated integrator for solving each particular part of the problem and combine the results to a self-consistent solution. In particular, we explore the possibilities of combining the evolution of one or more microscopic systems that are embedded in a macroscopic system. The here presented generalizations of Bridge include higher-order coupling strategies (from the classic 2nd order up to 10th-order), but we also demonstrate how multiple bridges can be nested and how additional processes can be introduced at the bridge time-step to enrich the physics, for example by incorporating dissipative processes. Such augmentation allows for including additional processes in a classic Newtonian N-body integrator without alterations to the underlying code. These additional processes include for example the Yarkovsky effect, dynamical friction or relativistic dynamics. Some of these processes operate on all particles whereas others apply only to a subset. The presented method is non-intrusive in the sense that the underlying methods remain operational without changes to the code (apart from adding the get- and set-functions to enable the bridge operator). As a result, the fundamental integrators continue to operate with their internal time step and preserve their local optimizations and parallelism. ... abridged ..., Comment: Accepted for publication in Communications in Nonlinear Science and Numerical Simulation (CNSNS) The associated software is part of the AMUSE framework and can be downloaded from http:www.amusecode.org

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., Comment: Accepted for publication in MNRAS. 7 pages, 4 figures

We present the smoothed-particle hydrodynamics simulation code, Bonsai-SPH, which is a continuation of our previously developed gravity-only hierarchical $N$-body code (called Bonsai). The code is optimized for Graphics Processing Unit (GPU) accelerators which enables researchers to take advantage of these powerful computational resources. Bonsa-SPH produces simulation results comparable with state-of-the-art, CPU based, codes, but using an order of magnitude less computation time. The code is freely available online and the details are described in this work., Comment: Updated intro and multi-GPU sections. 23 pages, 9 figures. Submission to SciPost

We present stellar_wind.py, a module that provides multiple methods of simulating stellar winds using smoothed particle hydrodynamics codes (SPH) within the astrophysical multipurpose software environment (AMUSE) framework. With the simple wind mode, we create SPH wind particles in a spherically symmetric shell. We inject the wind particles with a velocity equal to their terminal velocity. The accelerating wind mode is similar, but with this method particles can be injected with a lower initial velocity than the terminal velocity and they are accelerated away from the star according to an acceleration function. With the heating wind mode, SPH particles are created with zero initial velocity with respect to the star, but instead wind particles are given an internal energy based on the integrated mechanical luminosity of the star. This mode is designed to be used on longer timescales and larger spatial scales compared to the other two modes and assumes that the star is embedded in a gas cloud. For fast winds, we find that both the simple and accelerating mode can reproduce the desired velocity, density and temperature profiles. For slow winds, the simple wind mode is insufficient due to dominant hydrodynamical effects that change the wind velocities. The accelerating mode, with additional options to account for these hydrodynamical effects, can still reproduce the desired wind profiles. We test the heating mode by simulating both a normal wind and a supernova explosion of a single star in a uniform density medium. The stellar wind simulation results matches the analytical solution for an expanding wind bubble. The supernova simulation gives qualitatively correct results, but the simulated bubble expands faster than the analytical solution predicts. We conclude with an example of a triple star system which includes the colliding winds of all three stars., Comment: Accepted for publication in A&A

Scientific discovery is mediated by ideas that, after being formulated in hypotheses, can be tested, validated, and quantified before they eventually lead to accepted concepts. Computer-mediated discovery in astrophysics is no exception, but antiquated code that is only intelligible to scientists who were involved in writing it is holding up scientific discovery in the field. A bold initiative is needed to modernize astrophysics code and make it transparent and useful beyond a small group of scientists. (abridged), Comment: Published in Science

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., Comment: Accepted for publication in Communications in Nonlinear Science and Numerical Simulation. Calculations are performed with Brutus as part of the AMUSE framework. Simultion data will become available online

Scientific discovery via numerical simulations is important in modern astrophysics. This relatively new branch of astrophysics has become possible due to the development of reliable numerical algorithms and the high performance of modern computing technologies. These enable the analysis of large collections of observational data and the acquisition of new data via simulations at unprecedented accuracy and resolution. Ideally, simulations run until they reach some pre-determined termination condition, but often other factors cause extensive numerical approaches to break down at an earlier stage. In those cases, processes tend to be interrupted due to unexpected events in the software or the hardware. In those cases, the scientist handles the interrupt manually, which is time-consuming and prone to errors. We present the Simulation Monitor (SiMon) to automatize the farming of large and extensive simulation processes. Our method is light-weight, it fully automates the entire workflow management, operates concurrently across multiple platforms and can be installed in user space. Inspired by the process of crop farming, we perceive each simulation as a crop in the field and running simulation becomes analogous to growing crops. With the development of SiMon we relax the technical aspects of simulation management. The initial package was developed for extensive parameter searchers in numerical simulations, but it turns out to work equally well for automating the computational processing and reduction of observational data reduction., Comment: 13 pages, 6 figures, Accepted for publication in Publications of the Astronomical Society of the Pacific (PASP)

Astrophysical direct $N$-body methods have been one of the first production algorithms to be implemented using NVIDIA's CUDA architecture. Now, almost seven years later, the GPU is the most used accelerator device in astronomy for simulating stellar systems. In this paper we present the implementation of the Sapporo2 $N$-body library, which allows researchers to use the GPU for $N$-body simulations with little to no effort. The first version, released five years ago, is actively used, but lacks advanced features and versatility in numerical precision and support for higher order integrators. In this updated version we have rebuilt the code from scratch and added support for OpenCL, multi-precision and higher order integrators. We show how to tune these codes for different GPU architectures and present how to continue utilizing the GPU optimal even when only a small number of particles ($N < 100$) is integrated. This careful tuning allows Sapporo2 to be faster than Sapporo1 even with the added options and double precision data loads. The code runs on a range of NVIDIA and AMD GPUs in single and double precision accuracy. With the addition of OpenCL support the library is also able to run on CPUs and other accelerators that support OpenCL., Comment: 15 pages, 7 figures. Accepted for publication in Computational Astrophysics and Cosmology

We describe the political and technical complications encountered during the astronomical CosmoGrid project. CosmoGrid is a numerical study on the formation of large scale structure in the universe. The simulations are challenging due to the enormous dynamic range in spatial and temporal coordinates, as well as the enormous computer resources required. In CosmoGrid we dealt with the computational requirements by connecting up to four supercomputers via an optical network and make them operate as a single machine. This was challenging, if only for the fact that the supercomputers of our choice are separated by half the planet, as three of them are located scattered across Europe and fourth one is in Tokyo. The co-scheduling of multiple computers and the 'gridification' of the code enabled us to achieve an efficiency of up to $93\%$ for this distributed intercontinental supercomputer. In this work, we find that high-performance computing on a grid can be done much more effectively if the sites involved are willing to be flexible about their user policies, and that having facilities to provide such flexibility could be key to strengthening the position of the HPC community in an increasingly Cloud-dominated computing landscape. Given that smaller computer clusters owned by research groups or university departments usually have flexible user policies, we argue that it could be easier to instead realize distributed supercomputing by combining tens, hundreds or even thousands of these resources., Comment: Accepted for publication in IEEE conference on ERRORs

We describe the implementation and performance of the ${\rm P^3T}$ (Particle-Particle Particle-Tree) scheme for simulating dense stellar systems. In ${\rm P^3T}$, the force experienced by a particle is split into short-range and long-range contributions. Short-range forces are evaluated by direct summation and integrated with the fourth order Hermite predictor-corrector method with the block timesteps. For long-range forces, we use a combination of the Barnes-Hut tree code and the leapfrog integrator. The tree part of our simulation environment is accelerated using graphical processing units (GPU), whereas the direct summation is carried out on the host CPU. Our code gives excellent performance and accuracy for star cluster simulations with a large number of particles even when the core size of the star cluster is small.

We have simulated, for the first time, the long term evolution of the Milky Way Galaxy using 51 billion particles on the Swiss Piz Daint supercomputer with our $N$-body gravitational tree-code Bonsai. Herein, we describe the scientific motivation and numerical algorithms. The Milky Way model was simulated for 6 billion years, during which the bar structure and spiral arms were fully formed. This improves upon previous simulations by using 1000 times more particles, and provides a wealth of new data that can be directly compared with observations. We also report the scalability on both the Swiss Piz Daint and the US ORNL Titan. On Piz Daint the parallel efficiency of Bonsai was above 95%. The highest performance was achieved with a 242 billion particle Milky Way model using 18600 GPUs on Titan, thereby reaching a sustained GPU and application performance of 33.49 Pflops and 24.77 Pflops respectively., Comment: 12 pages, 4 figures, Published in: 'Proceeding SC '14 Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis'. Gordon Bell Prize 2014 finalist

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., Comment: 22 pages, 15 digures

We review the recent optimizations of gravitational $N$-body kernels for running them on graphics processing units (GPUs), on single hosts and massive parallel platforms. For each of the two main $N$-body techniques, direct summation and tree-codes, we discuss the optimization strategy, which is different for each algorithm. Because both the accuracy as well as the performance characteristics differ, hybridizing the two algorithms is essential when simulating a large $N$-body system with high-density structures containing few particles, and with low-density structures containing many particles. We demonstrate how this can be realized by splitting the underlying Hamiltonian, and we subsequently demonstrate the efficiency and accuracy of the hybrid code by simulating a group of 11 merging galaxies with massive black holes in the nuclei., Comment: Accepted for publication in IEEE Computer

We present a novel method for efficient direct integration of gravitational N-body systems with a large variation in characteristic time scales. The method is based on a recursive and adaptive partitioning of the system based on the connected components of the graph generated by the particle distribution combined with an interaction-specific time step criterion. It uses an explicit and approximately time-symmetric time step criterion, and conserves linear and angular momentum to machine precision. In numerical tests on astrophysically relevant setups, the method compares favourably to both alternative Hamiltonian-splitting integrators as well as recently developed block time step-based GPU-accelerated Hermite codes. Our reference implementation is incorporated in the HUAYNO code, which is freely available as a part of the AMUSE framework., Comment: 12 pages, 8 figures, accepted for A&A

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., Comment: ApJ Letters (accepted for publication)

Astronomical phenomena are governed by processes on all spatial and temporal scales, ranging from days to the age of the Universe (13.8,Gyr) as well as from km size up to the size of the Universe. This enormous range in scales is contrived, but as long as there is a physical connection between the smallest and largest scales it is important to be able to resolve them all, and for the study of many astronomical phenomena this governance is present. Although covering all these scales is a challenge for numerical modelers, the most challenging aspect is the equally broad and complex range in physics, and the way in which these processes propagate through all scales. In our recent effort to cover all scales and all relevant physical processes on these scales we have designed the Astrophysics Multipurpose Software Environment (AMUSE). AMUSE is a Python-based framework with production quality community codes and provides a specialized environment to connect this plethora of solvers to a homogeneous problem solving environment., Comment: 10 pages, Submitted to: Philosophical Transactions A

We present the open source Astrophysical Multi-purpose Software Environment (AMUSE, www.amusecode.org), a component library for performing astrophysical simulations involving different physical domains and scales. It couples existing codes within a Python framework based on a communication layer using MPI. The interfaces are standardized for each domain and their implementation based on MPI guarantees that the whole framework is well-suited for distributed computation. It includes facilities for unit handling and data storage. Currently it includes codes for gravitational dynamics, stellar evolution, hydrodynamics and radiative transfer. Within each domain the interfaces to the codes are as similar as possible. We describe the design and implementation of AMUSE, as well as the main components and community codes currently supported and we discuss the code interactions facilitated by the framework. Additionally, we demonstrate how AMUSE can be used to resolve complex astrophysical problems by presenting example applications., Comment: 23 pages, 25 figures, accepted for A&A

We review the implementation of individual particle time-stepping for N-body dynamics. We present a class of integrators derived from second order Hamiltonian splitting. In contrast to the usual implementation of individual time-stepping, these integrators are momentum conserving and show excellent energy conservation in conjunction with a symmetrized time step criterion. We use an explicit but approximate formula for the time symmetrization that is compatible with the use of individual time steps. No iterative scheme is necessary. We implement these ideas in the HUAYNO (available online at www.amusecode.org) code and present tests of the integrators and show that the presented integration schemes shows good energy conservation, with little or no systematic drift, while conserving momentum and angular momentum to machine precision for long term integrations., Comment: 20 pages, 9 figures, accepted for publication in New Astronomy

We introduce a general-purpose framework for interconnecting scientific simulation programs using a homogeneous, unified interface. Our framework is intrinsically parallel, and conveniently separates all component numerical modules in memory. This strict separation allows automatic unit conversion, distributed execution of modules on different cores within a cluster or grid, and orderly recovery from errors. The framework can be efficiently implemented and incurs an acceptable overhead. In practice, we measure the time spent in the framework to be less than 1% of the wall-clock time. Due to the unified structure of the interface, incorporating multiple modules addressing the same physics in different ways is relatively straightforward. Different modules may be advanced serially or in parallel. Despite initial concerns, we have encountered relatively few problems with this strict separation between modules, and the results of our simulations are consistent with earlier results using more traditional monolithic approaches. This framework provides a platform to combine existing simulation codes or develop new physical solver codes within a rich "ecosystem" of interchangeable modules., Comment: 36 pages, 12 figures and source code (accepted for publication in Computer Physics Communications)

In this short review we present the developments over the last 5 decades that have led to the use of Graphics Processing Units (GPUs) for astrophysical simulations. Since the introduction of NVIDIA's Compute Unified Device Architecture (CUDA) in 2007 the GPU has become a valuable tool for N-body simulations and is so popular these days that almost all papers about high precision N-body simulations use methods that are accelerated by GPUs. With the GPU hardware becoming more advanced and being used for more advanced algorithms like gravitational tree-codes we see a bright future for GPU like hardware in computational astrophysics., Comment: To appear in: European Physical Journal "Special Topics" : "Computer Simulations on Graphics Processing Units" . 18 pages, 8 figures

We present a gravitational hierarchical N-body code that is designed to run efficiently on Graphics Processing Units (GPUs). All parts of the algorithm are executed on the GPU which eliminates the need for data transfer between the Central Processing Unit (CPU) and the GPU. Our tests indicate that the gravitational tree-code outperforms tuned CPU code for all parts of the algorithm and show an overall performance improvement of more than a factor 20, resulting in a processing rate of more than 2.8 million particles per second., Comment: 5 pages, 2 figures. Proceedings of "Advances in Computational Astrophysics: methods, tools and outcomes", June 13-17, 2011, Cefalu, Sicily, Italy, eds. Capuzzo Dolcetta, Limongi, Tornambe and Giobbi

Most recent progress in understanding the dynamical evolution of star clusters relies on direct N-body simulations. Owing to the computational demands, and the desire to model more complex and more massive star clusters, hardware calculational accelerators, such as GRAPE special-purpose hardware or, more recently, GPUs (i.e. graphics cards), are generally utilised. In addition, simulations can be accelerated by adjusting parameters determining the calculation accuracy (i.e. changing the internal simulation time step used for each star). We extend our previous thorough comparison (Anders et al. 2009) of basic quantities as derived from simulations performed either with STARLAB/KIRA or NBODY6. Here we focus on differences arising from using different hardware accelerations (including the increasingly popular graphic card accelerations/GPUs) and different calculation accuracy settings. We use the large number of star cluster models (for a fixed stellar mass function, without stellar/binary evolution, primordial binaries, external tidal fields etc) already used in the previous paper, evolve them with STARLAB/KIRA (and NBODY6, where required), analyse them in a consistent way and compare the averaged results quantitatively. For this quantitative comparison, we apply the bootstrap algorithm for functional dependencies developed in our previous study. In general we find very high comparability of the simulation results, independent of the used computer hardware (including the hardware accelerators) and the used N-body code. For the tested accuracy settings we find that for reduced accuracy (i.e. time step at least a factor 2.5 larger than the standard setting) most simulation results deviate significantly from the results using standard settings. The remaining deviations are comprehensible and explicable., Comment: 14 pages incl. 3 pages with figures and 4 pages of tables (analysis results), MNRAS in press

We present a data format for the output of general N-body simulations, allowing the presence of individual time steps. By specifying a standard, different N-body integrators and different visualization and analysis programs can all share the simulation data, independent of the type of programs used to produce the data. Our Particle Stream Data Format, PSDF, is specified in YAML, based on the same approach as XML but with a simpler syntax. Together with a specification of PSDF, we provide background and motivation, as well as specific examples in a variety of computer languages. We also offer a web site from which these examples can be retrieved, in order to make it easy to augment existing codes in order to give them the option to produce PSDF output., Comment: 5 pages; submitted to New Astronomy

We describe AMUSE, the Astrophysical Multipurpose Software Environment, a programming framework designed to manage multi-scale, multi-physics simulations in a hierarchical, extensible, and internally consistent way. Constructed as a collection of individual modules, AMUSE allows computational tools for different physical domains to be easily combined into a single task. It facilitates the coupling of modules written in different languages by providing inter-language tools and a standard programming interface that represents a balance between generality and computational efficiency. The framework currently incorporates the domains of stellar dynamics, stellar evolution, gas dynamics, and radiative transfer. We present some applications of the framework and outline plans for future development of the package., Comment: 8 pages, 4 figures; Invited talk; to appear in the Proceedings of "Advances in Computational Astrophysics: methods, tools and outcomes," Cefal\`u (Sicily, Italy), June 13-17, 2011

The young star clusters we observe today are the building blocks of a new generation of stars and planets in our Galaxy and beyond. Despite their fundamental role we still lack knowledge about the conditions under which star clusters form and the impact of these often harsh environments on the evolution of their stellar and substellar members. We demonstrate the vital role numerical simulations play to uncover both key issues. Using dynamical models of different star cluster environments we show the variety of effects stellar interactions potentially have. Moreover, our significantly improved measure of mass segregation reveals that it can occur rapidly even for star clusters without substructure. This finding is a critical step to resolve the controversial debate on mass segregation in young star clusters and provides strong constraints on their initial conditions., Comment: 8 pages, 4 figures; to appear in the proceedings of "Stellar Clusters and Associations - A RIA workshop on Gaia", 23-27 May 2011, Granada, Spain

We present parallel algorithms for constructing and traversing sparse octrees on graphics processing units (GPUs). The algorithms are based on parallel-scan and sort methods. To test the performance and feasibility, we implemented them in CUDA in the form of a gravitational tree-code which completely runs on the GPU.(The code is publicly available at: http://castle.strw.leidenuniv.nl/software.html) The tree construction and traverse algorithms are portable to many-core devices which have support for CUDA or OpenCL programming languages. The gravitational tree-code outperforms tuned CPU code during the tree-construction and shows a performance improvement of more than a factor 20 overall, resulting in a processing rate of more than 2.8 million particles per second., Comment: Accepted version. Published in Journal of Computational Physics. 35 pages, 12 figures, single column

We present a new very fast tree-code which runs on massively parallel Graphical Processing Units (GPU) with NVIDIA CUDA architecture. The tree-construction and calculation of multipole moments is carried out on the host CPU, while the force calculation which consists of tree walks and evaluation of interaction list is carried out on the GPU. In this way we achieve a sustained performance of about 100GFLOP/s and data transfer rates of about 50GB/s. It takes about a second to compute forces on a million particles with an opening angle of $\theta \approx 0.5$. The code has a convenient user interface and is freely available for use\footnote{{\tt http://castle.strw.leidenuniv.nl/software/octgrav.html}}., Comment: 9 pages, 8 figures. Accepted for publication at International Conference on Computational Science 2010

Understanding the universe is hampered by the elusiveness of its most common constituent, cold dark matter. Almost impossible to observe, dark matter can be studied effectively by means of simulation and there is probably no other research field where simulation has led to so much progress in the last decade. Cosmological N-body simulations are an essential tool for evolving density perturbations in the nonlinear regime. Simulating the formation of large-scale structures in the universe, however, is still a challenge due to the enormous dynamic range in spatial and temporal coordinates, and due to the enormous computer resources required. The dynamic range is generally dealt with by the hybridization of numerical techniques. We deal with the computational requirements by connecting two supercomputers via an optical network and make them operate as a single machine. This is challenging, if only for the fact that the supercomputers of our choice are separated by half the planet, as one is located in Amsterdam and the other is in Tokyo. The co-scheduling of the two computers and the 'gridification' of the code enables us to achieve a 90% efficiency for this distributed intercontinental supercomputer., Comment: Accepted for publication in IEEE Computer

The onset of runaway stellar collisions in young star clusters is more likely to initiate with an encounter between a binary and a third star than between two single stars. Using the initial conditions of such three-star encounters from direct $N$-body simulations, we model the resulting interaction by means of Smoothed Particle Hydrodynamics (SPH). We find that, in the majority of the cases considered, all three stars merge together, and in such three star mergers, the hydrodynamic simulations reveal that: (1) mass lost as ejecta can be a considerable fraction of the total mass in the system (up to $\sim25$%); (2) due to asymmetric mass loss, the collision product can sometimes receive a kick velocity that exceeds 10 km/s, large enough to allow the collision product to escape the core of the cluster; and (3) the energy of the ejected matter can be large enough (up to $\sim 3\times 10^{50}$ erg) to remove or disturb the inter cluster gas appreciably., Comment: 23 pages, 24 figures, Submitted to MNRAS

N-body simulations are widely used to simulate the dynamical evolution of a variety of systems, among them star clusters. Much of our understanding of their evolution rests on the results of such direct N-body simulations. They provide insight in the structural evolution of star clusters, as well as into the occurrence of stellar exotica. Although the major pure N-body codes STARLAB/KIRA and NBODY4 are widely used for a range of applications, there is no thorough comparison study yet. Here we thoroughly compare basic quantities as derived from simulations performed either with STARLAB/KIRA or NBODY4. We construct a large number of star cluster models for various stellar mass function settings (but without stellar/binary evolution, primordial binaries, external tidal fields etc), evolve them in parallel with STARLAB/KIRA and NBODY4, analyse them in a consistent way and compare the averaged results quantitatively. For this quantitative comparison we develop a bootstrap algorithm for functional dependencies. We find an overall excellent agreement between the codes, both for the clusters' structural and energy parameters as well as for the properties of the dynamically created binaries. However, we identify small differences, like in the energy conservation before core collapse and the energies of escaping stars, which deserve further studies. Our results reassure the comparability and the possibility to combine results from these two major N-body codes, at least for the purely dynamical models (i.e. without stellar/binary evolution) we performed. (abridged), Comment: accepted for publication in MNRAS, database can be found here: http://www.phys.uu.nl/~anders/data/NBODY_STARLAB_Comparison/ and http://members.galev.org/nbody/NBODY_STARLAB_Comparison/

We present Sapporo, a library for performing high-precision gravitational N-body simulations on NVIDIA Graphical Processing Units (GPUs). Our library mimics the GRAPE-6 library, and N-body codes currently running on GRAPE-6 can switch to Sapporo by a simple relinking of the library. The precision of our library is comparable to that of GRAPE-6, even though internally the GPU hardware is limited to single precision arithmetics. This limitation is effectively overcome by emulating double precision for calculating the distance between particles. The performance loss of this operation is small (< 20%) compared to the advantage of being able to run at high precision. We tested the library using several GRAPE-6-enabled N-body codes, in particular with Starlab and phiGRAPE. We measured peak performance of 800 Gflop/s for running with 10^6 particles on a PC with four commercial G92 architecture GPUs (two GeForce 9800GX2). As a production test, we simulated a 32k Plummer model with equal mass stars well beyond core collapse. The simulation took 41 days, during which the mean performance was 113 Gflop/s. The GPU did not show any problems from running in a production environment for such an extended period of time., Comment: 13 pages, 9 figures, accepted to New Astronomy