With this report we provide users of an easy manual to facilitate the proper download and use of a sophisticated, high precision, few-body code originally developed by S. Mikkola, and later largely improved and implemented to treat a variety of cases. The code download can be done via the link https://drive.google.com/file/d/16FkVVR4Tk8eKhKMju2vQ9rlWI4Mpv01W/view The use of the code is free upon proper citation. The work is in progress and users are invited to help the authors to improve both the code and the user handbook.

Blazar OJ287 exhibits large thermal flares at least twice every 12 years. The times of these flares have been predicted successfully using the model of a quasi-Keplerian eccentric black hole binary where the secondary impacts the accretion disk of the primary, creating the thermal flares. New measurements of the historical light curve have been combined with the observations of the 2015 November/December flare to identify the impact record since year 1886, and to constrain the orbit of the binary. The orbit solution shows that the binary period, now 12.062 year, is decreasing at the rate of 36 days per century. This corresponds to an energy loss to gravitational waves that is 6.5 +- 4 % less than the rate predicted by the standard quadrupolar gravitational wave (GW) emission. We show that the difference is due to higher order gravitational radiation reaction terms that include the dominant order tail contributions.

The quasi-satellite (QS) phenomenon makes two celestial bodies to fly near each other (Mikkola et al. 2006) and that effect can be used also to make artificial satellites move in tandem. We consider formation flight of two or three satellites in low eccentricity near Earth orbits. With the help of weak ion thrusters it is possible to accomplish tandem flight. With ion thrusters it is also possible to mimic many kinds of mutual force laws between the satellites. We found that both a constant repulsive force or an attractive force that decreases with the distance are able to preserve the formation in which the eccentricities cause the actual relative motion and the weak thrusters keep the mean longitude difference small. Initial values are important for the formation flight but very exact adjustment of orbital elements is not important. Simplicity is one of our goals in this study and this result is achieved at least in the way that, when constant force thrusters are used, the satellites only need to detect the directions of the other ones to fly in tandem. A repulsive acceleration of the order of 10^-6 times the Earth attraction, is enough to effectively eliminate the disruptive effects of all the perturbations at least for a timescale of years.

We have carried out 200,000 N-body simulations of three identical stellar embryos with masses from a Chabrier IMF and embedded in a molecular core. The bodies are initially non-hierarchical and undergo chaotic motions, while accreting using Bondi-Hoyle accretion. The coupling of dynamics and accretion often leads to one or two dominant bodies controlling the center of the cloud core, while banishing the other(s) to the lower-density outskirts, leading to stunted growth. Eventually each system transforms either to a bound hierarchical configuration or breaks apart into separate single and binary components. The orbital motion is followed for 100 Myr. To illustrate the simulations we introduce the 'triple diagnostic diagram', which plots two dimensionless numbers against each other, representing the binary mass ratio and the mass ratio of the third body to the total system mass. Numerous freefloating BD binaries are formed in these simulations. The separation distribution function is in good correspondence with observations, showing a steep rise at close separations, peaking around 13 AU and then declining more gently. Unresolved BD triple systems may appear as wider BD binaries. Mass ratios are strongly peaked towards unity, as observed, but this is partially due to the initial assumptions. Eccentricities gradually increase towards higher values, due to the lack of viscous interactions in the simulations, which would both shrink the orbits and decrease their eccentricities. The main threat to newly born triple systems is internal instabilities, not external perturbations. Dynamical interactions in newborn triple systems of stellar embryos embedded in and accreting from a cloud core naturally form a population of freefloating BD binaries, and this mechanism may constitute a significant pathway for the formation of BD binaries. (Abstract abbreviated).

Coalescing black hole (BH) binaries forming in the dense core of globular clusters (GCs) are expected to be one the brightest sources of gravitational wave (GW) radiation for the next generation of ground-based laser interferometers. Favorable conditions for merger are initiated by the Kozai resonance in which the gravitational interaction with a third distant object, typically another BH, induces quasi-periodic variations of the inner BH binary eccentricity. In this paper we perform high precision N-body simulations of the long term evolution of hierarchical BH triples and investigate the conditions that lead to the merging of the BH binary and the way it might become an observable source of GW radiation. We find that the secular orbit average treatment, adopted in previous works, does not reliably describe the dynamics of these systems if the binary is orbited by the outer BH on a highly inclined orbit at a moderate distance. During the high eccentricity phase of a Kozai cycle the torque due to the outer BH can drive the binary to extremely large eccentricities in a fraction of the binary's orbital period. This occurs before relativistic terms become important to the evolution and allows the binary GW signal to reach large GW frequencies (>~10 Hz) at high eccentricities. We show that ~50 % of coalescing BH binaries driven by the Kozai mechanism in GCs will have eccentricities larger than 0.1, with 10 % of them being extremely eccentric, (1-e)<~10^-5, when they first chirp in the frequency band of ground based laser interferometers. This implies that a large fraction of such GW sources could be missed if conventional quasi-circular templates are used for analysis of GW detectors data. The efficient detection of all coalescing BH binaries in GCs will therefore require template banks of eccentric inspiral waveforms for matched-filtering and dedicated search strategies.

The formation of very wide binaries, such as the alpha Cen system with Proxima (also known as alpha Centauri C) separated from alpha Centauri (which itself is a close binary A/B) by 15000 AU, challenges current theories of star formation, because their separation can exceed the typical size of a collapsing cloud core. Various hypotheses have been proposed to overcome this problem, including the suggestion that ultra-wide binaries result from the dissolution of a star cluster - when a cluster star gravitationally captures another, distant, cluster star. Recent observations have shown that very wide binaries are frequently members of triple systems and that close binaries often have a distant third companion. Here we report Nbody simulations of the dynamical evolution of newborn triple systems still embedded in their nascent cloud cores that match observations of very wide systems. We find that although the triple systems are born very compact - and therefore initially are more protected against disruption by passing stars - they can develop extreme hierarchical architectures on timescales of millions of years as one component is dynamically scattered into a very distant orbit. The energy of ejection comes from shrinking the orbits of the other two stars, often making them look from a distance like a single star. Such loosely bound triple systems will therefore appear to be very wide binaries.

The blazar OJ287 is the most promising (and the only) case for an extragalactic binary black hole system inspiralling under the action of gravitational radiation reaction. At present, though it is not possible to directly observe the binary components, it is possible to observe the jet emanating form the primary black hole. We argue that the orbital motion of the secondary black hole is reflected in the wobble of the jet and demonstrate that the wobble is orbital position dependent. The erratic wobble of the jet, reported in Agudo et al. (2012), is analyzed by taking into account the binary nature of the system and we find that the erratic component of jet wobble is very small.

We examine the ability to test the black hole no-hair theorem at the 10% level in this decade using the binary black hole in OJ287. In the test we constrain the value of the dimensionless parameter q that relates the scaled quadrupole moment and spin of the primary black hole: q2 = -q 2 . At the present we can say that q = 1 \pm 0.3 (one), in agreement with General Relativity and the no-hair theorems. We demonstrate that this result can be improved if more observational data is found in historical plate archives for the 1959 and 1971 outbursts. We also show that the predicted 2015 and 2019 outbursts will be crucial in improving the accuracy of the test. Space-based photometry is required in 2019 July due the proximity of OJ287 to the Sun at the time of the outburst. The best situation would be to carry out the photometry far from the Earth, from quite a different vantage point, in order to avoid the influence of the nearby Sun. We have considered in particular the STEREO space mission which would be ideal if it has a continuation in 2019 or LORRI on board the New Horizons mission to Pluto.

Inspiral of compact stellar remnants into massive black holes (MBHs) is accompanied by the emission of gravitational waves at frequencies that are potentially detectable by space-based interferometers. Event rates computed from statistical (Fokker-Planck, Monte-Carlo) approaches span a wide range due to uncertaintities about the rate coefficients. Here we present results from direct integration of the post-Newtonian N-body equations of motion descrbing dense clusters of compact stars around Schwarzschild MBHs. These simulations embody an essentially exact (at the post-Newtonian level) treatment of the interplay between stellar dynamical relaxation, relativistic precession, and gravitational-wave energy loss. The rate of capture of stars by the MBH is found to be greatly reduced by relativistic precession, which limits the ability of torques from the stellar potential to change orbital angular momenta. Penetration of this "Schwarzschild barrier" does occasionally occur, resulting in capture of stars onto orbits that gradually inspiral due to gravitational wave emission; we discuss two mechanisms for barrier penetration and find evidence for both in the simulations. We derive an approximate formula for the capture rate, which predicts that captures would be strongly disfavored from orbits with semi-major axes below a certain value; this prediction, as well as the predicted rate, are verified in the N-body integrations. We discuss the implications of our results for the detection of extreme-mass-ratio inspirals from galactic nuclei with a range of physical properties.

We explore the origin of a population of distant companions (~1000 - 5000 AU) to Class I protostellar sources recently found by Connelley and co-workers, who noted that the companion fraction diminished as the sources evolved. Here we present N-body simulations of unstable triple systems embedded in dense cloud cores. Many companions are ejected into unbound orbits and quickly escape, but others are ejected with insufficient momentum to climb out of the potential well of the cloud core and associated binary. These loosely bound companions reach distances of many thousands of AU before falling back and eventually being ejected into escapes as the cloud cores gradually disappear. We use the term orphans to denote protostellar objects that are dynamically ejected from their placental cloud cores, either escaping or for a time being tenuously bound at large separations. Half of all triple systems are found to disintegrate during the protostellar stage, so if multiple systems are a frequent outcome of the collapse of a cloud core, then orphans should be common. Bound orphans are associated with embedded close protostellar binaries, but escaping orphans can travel as far as ~0.2 pc during the protostellar phase. The steep climb out of a potential well ensures that orphans are not kinematically distinct from young stars born with a less violent pre-history. The identification of orphans outside their heavily extincted cloud cores will allow the detailed study of protostars high up on their Hayashi tracks at near-infrared and in some cases even at optical wavelengths.

The long-term dynamics of Oort cloud comets are studied under the influence of both the radial and the vertical components of the Galactic tidal field. Sporadic dynamical perturbation processes are ignored, such as passing stars, since we aim to study the influence of just the axisymmetric Galactic tidal field on the cometary motion and how it changes in time. We use a model of the Galaxy with a disc, bulge and dark halo, and a local disc density, and disc scale length constrained to fit the best available observational constraints. By integrating a few million of cometary orbits over 1 Gyr, we calculate the time variable flux of Oort cloud comets that enter the inner Solar System, for the cases of a constant Galactic tidal field, and a realistically varying tidal field which is a function of the Sun's orbit. The applied method calculates the evolution of the comets by using first-order averaged mean elements. We find that the periodicity in the cometary flux is complicated and quasi-periodic. The amplitude of the variations in the flux are of order 30%. The radial motion of the Sun is the chief cause of this behaviour, and should be taken into account when the Galactic influence on the Oort cloud comets is studied.

We model the binary black hole system OJ287 as a spinning primary and a non-spinning secondary. It is assumed that the primary has an accretion disk which is impacted by the secondary at specific times. These times are identified as major outbursts in the light curve of OJ287. This identification allows an exact solution of the orbit, with very tight error limits. Nine outbursts from both the historical photographic records as well as from recent photometric measurements have been used as fixed points of the solution: 1913, 1947, 1957, 1973, 1983, 1984, 1995, 2005 and 2007 outbursts. This allows the determination of eight parameters of the orbit. Most interesting of these are the primary mass of $1.84\cdot 10^{10} M_\odot$, the secondary mass $1.46\cdot 10^{8} M_\odot$, major axis precession rate $39^\circ.1$ per period, and the eccentricity of the orbit 0.70. The dimensionless spin parameter is $0.28\:\pm\:0.01$ (1 sigma). The last parameter will be more tightly constrained in 2015 when the next outburst is due. The outburst should begin on 15 December 2015 if the spin value is in the middle of this range, on 3 January 2016 if the spin is 0.25, and on 26 November 2015 if the spin is 0.31. We have also tested the possibility that the quadrupole term in the Post Newtonian equations of motion does not exactly follow Einstein's theory: a parameter $q$ is introduced as one of the 8 parameters. Its value is within 30% (1 sigma) of the Einstein's value $q = 1$. This supports the $no-hair theorem$ of black holes within the achievable precision. We have also measured the loss of orbital energy due to gravitational waves. The loss rate is found to agree with Einstein's value with the accuracy of 2% (1 sigma).

The compact binary system in OJ287 is modelled to contain a spinning primary black hole with an accretion disk and a non-spinning secondary black hole. Using Post Newtonian (PN) accurate equations that include 2.5PN accurate non-spinning contributions, the leading order general relativistic and classical spin-orbit terms, the orbit of the binary black hole in OJ287 is calculated and as expected it depends on the spin of the primary black hole. Using the orbital solution, the specific times when the orbit of the secondary crosses the accretion disk of the primary are evaluated such that the record of observed outbursts from 1913 up to 2007 is reproduced. The timings of the outbursts are quite sensitive to the spin value. In order to reproduce all the known outbursts, including a newly discovered one in 1957, the Kerr parameter of the primary has to be $0.28 \pm 0.08$. The quadrupole-moment contributions to the equations of motion allow us to constrain the `no-hair' parameter to be $1.0\:\pm\:0.3$ where 0.3 is the one sigma error. This supports the `black hole no-hair theorem' within the achievable precision. It should be possible to test the present estimate in 2015 when the next outburst is due. The timing of the 2015 outburst is a strong function of the spin: if the spin is 0.36 of the maximal value allowed in general relativity, the outburst begins in early November 2015, while the same event starts in the end of January 2016 if the spin is 0.2

The spin and quadrupole moment of the supermassive black hole at the Galactic center can in principle be measured via astrometric monitoring of stars orbiting at milliparsec (mpc) distances, allowing tests of general relativistic "no-hair" theorems (Will 2008). One complicating factor is the presence of perturbations from other stars, which may induce orbital precession of the same order of magnitude as that due to general relativistic effects. The expected number of stars in this region is small enough that full N-body simulations can be carried out. We present the results of a comprehensive set of such simulations, which include a post-Newtonian treatment of spin-orbit effects. A number of possible models for the distribution of stars and stellar remnants are considered. We find that stellar perturbations are likely to obscure the signal due to frame-dragging for stars beyond ~0.5 mpc from the black hole, while measurement of the quadrupole moment is likely to require observation of stars inside ~0.2 mpc. A high fraction of stellar remnants, e.g. 10-Solar-mass black holes, in this region would make tests of GR problematic at all radii. We discuss the possibility of separating the effects of stellar perturbations from those due to GR.

Supermassive black holes are common in centers of galaxies. Among the active galaxies, quasars are the most extreme, and their black hole masses range as high as to $6\cdot 10^{10} M_\odot$. Binary black holes are of special interest but so far OJ287 is the only confirmed case with known orbital elements. In OJ287, the binary nature is confirmed by periodic radiation pulses. The period is twelve years with two pulses per period. The last four pulses have been correctly predicted with the accuracy of few weeks, the latest in 2007 with the accuracy of one day. This accuracy is high enough that one may test the higher order terms in the Post Newtonian approximation to General Relativity. The precession rate per period is $39^\circ.1 \pm 0^\circ.1$, by far the largest rate in any known binary, and the $(1.83\pm 0.01)\cdot 10^{10} M_\odot$ primary is among the dozen biggest black holes known. We will discuss the various Post Newtonian terms and their effect on the orbit solution. The over 100 year data base of optical variations in OJ287 puts limits on these terms and thus tests the ability of Einstein's General Relativity to describe, for the first time, dynamic binary black hole spacetime in the strong field regime. The quadrupole-moment contributions to the equations of motion allows us to constrain the `no-hair' parameter to be $1.0\:\pm\:0.3$ which supports the black hole no-hair theorem within the achievable precision.

We discuss the implementation of a new regular algorithm for simulation of the gravitational few-body problem. The algorithm uses components from earlier methods, including the chain structure, the logarithmic Hamiltonian, and the time-transformed leapfrog. This algorithmic regularization code, AR-CHAIN, can be used for the normal N-body problem, as well as for problems with softened potentials and/or with velocity-dependent external perturbations, including post-Newtonian terms, which we include up to order PN2.5. Arbitrarily extreme mass ratios are allowed. Only linear coordinate transformations are used and thus the algorithm is somewhat simpler than many earlier regularized schemes. We present the results of performance tests which suggest that the new code is either comparable in performance or superior to the existing regularization schemes based on the Kustaanheimo-Stiefel (KS) transformation. This is true even for the two-body problem, independent of eccentricity. An important advantage of the new method is that, contrary to the older KS-CHAIN code, zero masses are allowed. We use our algorithm to integrate the orbits of the S stars around the Milky Way supermassive black hole for one million years, including PN2.5 terms and an intermediate-mass black hole. The three S stars with shortest periods are observed to escape from the system after a few hundred thousand years.

The young stars near the supermassive black hole at the galactic center follow orbits that are nearly random in orientation and that have an approximately thermal distribution of eccentricities, N(e) ~ e. We show that both of these properties are a natural consequence of a few million years' interaction with an intermediate-mass black hole (IBH), if the latter's orbit is mildly eccentric and if its mass exceeds approximately 1500 solar masses. Producing the most tightly-bound S-stars requires an IBH orbit with periastron distance less than about 10 mpc. Our results provide support for a model in which the young stars are carried to the galactic center while bound to an IBH, and are consistent with the hypothesis that an IBH may still be orbiting within the nuclear star cluster.

We describe a novel N-body code designed for simulations of the central regions of galaxies containing massive black holes. The code incorporates Mikkola's 'algorithmic' chain regularization scheme including post-Newtonian terms up to PN2.5 order. Stars moving beyond the chain are advanced using a fourth-order integrator with forces computed on a GRAPE board. Performance tests confirm that the hybrid code achieves better energy conservation, in less elapsed time, than the standard scheme and that it reproduces the orbits of stars tightly bound to the black hole with high precision. The hybrid code is applied to two sample problems: the effect of finite-N gravitational fluctuations on the orbits of the S-stars; and inspiral of an intermediate-mass black hole into the galactic center.

In this paper we carry out a quantitative analysis of the three-body systems and map them as a function of decaying time and intial conguration, look at this problem as an example of a simple deterministic system, and ask to what extent the orbits are really predictable. We have investigated the behavior of about 200 000 general Newtonian three body systems using the simplest initial conditions. Within our resolution these cover all the possible states where the objects are initially at rest and have no angular momentum. We have determined the decay time-scales of the triple systems and show that the distribution of this parameter is fractal in appearance. Some areas that appear stable on large scales exhibit very narrow strips of instability and the overall pattern, dominated by resonances, reminds us of a traditional Maasai warrior shield. Also an attempt is made to recover the original starting conguration of the three bodies by backward integration. We find there are instances where the evolution to the future and to the past lead to different orbits, in spite of time symmetric initial conditions. This implies that even in simple deterministic systems there exists an Arrow of Time.

[Abridged] In galactic nuclei with sufficiently short relaxation times, binary supermassive black holes can evolve beyond their stalling radii via continued interaction with stars. We study this "collisional" evolutionary regime using both fully self-consistent N-body integrations and approximate Fokker-Planck models. The N-body integrations employ particle numbers up to 0.26M and a direct-summation potential solver; close interactions involving the binary are treated using a new implementation of the Mikkola-Aarseth chain regularization algorithm. Even at these large values of N, two-body scattering occurs at high enough rates in the simulations that they can not be simply scaled to the large-N regime of real galaxies. The Fokker-Planck model is used to bridge this gap; it includes, for the first time, binary-induced changes in the stellar density and potential. The Fokker-Planck model is shown to accurately reproduce the results of the N-body integrations, and is then extended to the much larger N regime of real galaxies. Analytic expressions are derived that accurately reproduce the time dependence of the binary semi-major axis as predicted by the Fokker-Planck model. Gravitational wave coalescence is shown to occur in <10 Gyr in nuclei with velocity dispersions below about 80 km/s. Formation of a core results from a competition between ejection of stars by the binary and re-supply of depleted orbits via two-body scattering. Mass deficits as large as ~4 times the binary mass are produced before coalescence. After the two black holes coalesce, a Bahcall-Wolf cusp appears around the single hole in one relaxation time, resulting in a nuclear density profile consisting of a flat core with an inner, compact cluster, similar to what is observed at the centers of low-luminosity spheroids.

We conducted extensive numerical experiments of equal mass three-body systems until they became disrupted. The system lifetimes, as a bound triple, and the Lyapunov times show a correlation similarto what has been earlier obtained for small bodies in the Solar System. Numerical integrations of several sets of differently randomised initial conditions produced the same relationship of the instability time and Lyapunov time. Marginal probability densities of the various times in the three-body experiments are also discussed. Our high accuracy numerical method for three-body orbit computations and Lyapunov time determinations is concisely described.

We consider secular perturbations of nearly Keplerian two-body motion under a perturbing potential that can be approximated to sufficient accuracy by expanding it to second order in the coordinates. After averaging over time to obtain the secular Hamiltonian, we use angular momentum and eccentricity vectors as elements. The method of variation of constants then leads to a set of equations of motion that are simple and regular, thus allowing efficient numerical integration. Some possible applications are briefly described.

Algorithmic regularization uses a transformation of the equations of motion such that the leapfrog algorithm produces exact trajectories for two-body motion as well as regular results in numerical integration of the motion of strongly interacting few-body systems. That algorithm alone is not sufficiently accurate and one must use the extrapolation method for improved precision. This requires that the basic leapfrog algorithm be time-symmetric, which is not directly possible in the case of velocity-dependent forces, but is usually obtained with the help of the implicit midpoint method. Here we suggest an alternative explicit algorithmic regularization algorithm which can handle velocity-dependent forces. This is done with the help of a generalized midpoint method to obtain the required time symmetry, thus eliminating the need for the implicit midpoint method and allowing the use of extrapolation.

A chain regularization method is combined with special purpose computer hardware to study the evolution of massive black hole binaries at the centers of galaxies. Preliminary results with up to N=260,000 particles are presented. The decay rate of the binary is shown to decrease with increasing N, as expected on the basis of theoretical arguments. The eccentricity of the binary remains small.

We investigate the dynamical decay of non-hierarchical accreting triple systems and its implications on the ejection model as Brown Dwarf formation scenario. A modified chain-regularization scheme is used to integrate the equations of motion, that also allows for mass changes over time as well as for momentum transfer from the accreted gas mass onto the bodies. We integrate an ensemble of triple systems within a certain volume with different accretion rates, assuming several prescriptions of how momentum is transferred onto the bodies. We follow their evolution until the systems have decayed. We analyze the end states and decay times of these systems and determine the fraction of Brown Dwarfs formed, their escape speeds as well as the semi-major axis distribution of the formed Brown Dwarf binaries. We find that the formation probability of Brown Dwarfs depends strongly on the assumed momentum transfer which is related to the motion of the gas. Due to ongoing accretion and consequent shrinkage of the systems, the median escape velocity is increased by a factor of 2 and the binary separations are decreased by a factor of 5 compared with non-accreting systems. Furthermore, the obtained semi-major axis distribution drops off sharply to either side of the median, which is also supported by observations. We conclude that accretion and momentum transfer of accreted gas during the dynamical decay of triple systems is able to produce the observed distribution of close binary Brown Dwarfs, making the ejection model a viable option as Brown Dwarf formation scenario.

Observational searches for asteroids orbiting near Earth's triangular Lagrange points face unique obstacles. A population of such asteroids would occupy a large projected area on the sky (possibly hundreds of square degrees) and is not favorably placed with respect to the Sun. Here we examine the properties of synthetic populations of Earth ``Trojans'' in order to aid in the optimization of observational searches for them. We find that the highest on-sky projected number densities are not located at the positions of the L4 and L5 points themselves, but rather a few degrees closer to the Sun. Also, asteroids on orbits about the L4 and L5 points typically brighten as the difference between their ecliptic longitude and that of the Sun increases owing to phase effects, but their number density on the sky concurrently falls rapidly.