A new data structure for accelerating kinetic Monte Carlo method

The kinetic Monte Carlo simulation is a rigorous numerical approach to study the chemistry on dust grains in cold dense interstellar clouds. By tracking every single reaction in chemical networks step by step, this approach produces more precise results than other approaches but takes too much computing time. Here we present a method of a new data structure, which is applicable to any physical conditions and chemical networks, to save computing time for the Monte Carlo algorithm. Using the improved structure,the calculating time is reduced by 80 percent compared with the linear structure when applied to the osu-2008 chemical network at 10K. We investigate the effect of the encounter desorption in cold cores using the kinetic Monte Carlo model with an accelerating data structure. We found that the encounter desorption remarkably decreases the abundance of grain-surface H_2 but slightly influences the abundances of other species on the grain.     

High-order expansion of the solution of preliminary orbit determination problem

A method for high-order treatment of uncertainties in preliminary orbit determination is presented. The observations consist in three couples of topocentric right ascensions and declinations at three observation epochs. The goal of preliminary orbit determination is to compute a trajectory that fits with the observations in two-body dynamics. The uncertainties of the observations are usually mapped to the phase space only when additional observations are available and a least squares fitting problem is set up. A method based on Taylor differential algebra for the analytical treatment of observation uncertainties is implemented. Taylor differential algebra allows for the efficient computation of the arbitrary order Taylor expansion of a sufficiently continuous multivariate function. This enables the mapping of the uncertainties from the observation space to the phase space as high-order multivariate Taylor polynomials. These maps can then be propagated forward in time to predict the observable set at successive epochs. This method can be suitably used to recover newly discovered objects when a scarce number of measurements is available. Simulated topocentric observations of asteroids on realistic orbits are used to assess the performances of the method.     

Gravity assist space pruning based on differential algebra

In this paper a differential algebra version of the gravity assist space pruning algorithm is presented. The use of differential algebraic techniques is proposed to overcome the two main drawbacks of the existing algorithm, i.e., the steep increase of the number of function evaluations with the number of planets involved in the transfer, and the use of a bounding procedure that relies on Lipschitzian tolerances. Differential algebra allows us to process boxes in place of grid points, and to substitute pointwise evaluations of the constraint functions with their Taylor expansions. Thanks to the particular instance of multi-gravity assist problems dealt with, all the planet-to-planet legs can be treated independently, and forward and backward constraining can be applied. The proposed method is applied to preprocess the search space of sample interplanetary transfers and it also serves as a stepping stone towards a fully rigorous treatment of the pruning process based on Taylor models.     

How selection and weighting of astrometric observations influence the impact probability. The case of asteroid (99942) Apophis

The aim of this paper is to show that in the case of a low probability of asteroid collision with the Earth, the appropriate selection and weighting of the data are crucial for the impact investigation and for analysing the impact possibilities using extensive numerical simulations. By means of the Monte Carlo special method, a large number of 'clone' orbits have been generated. A full range of orbital elements in the six-dimensional parameter space, that is, in the entire confidence region allowed by the observational material, has been examined. On the basis of 1000 astrometric observations of (99942) Apophis, the best solutions for the geocentric encounter distance of  6.065 ± 0.081 R_⊕  (without perturbations by asteroids) or  6.064 ± 0.095 R_⊕  (including perturbations by the four largest asteroids) were derived for the close encounter with the Earth on 2029 April 13. The present uncertainties allow for special configurations ('keyholes') during this encounter that may lead to very close encounters in future approaches of Apophis. Two groups of keyholes are connected with the close encounter with the Earth in 2036 (within the minimal distance of  5.7736−5.7763 R_⊕  on 2029 April 13) and 2037 (within the minimal distance of  6.3359–6.3488 R_⊕  ). The nominal orbits for our most accurate models run almost exactly in the middle of these two impact keyhole groups. A very small keyhole for the impact in 2076 has been found between these groups at the minimal distance of 5.97347   R_⊕  . This keyhole is close to the nominal orbit. The present observations are not sufficiently accurate to eliminate definitely the possibility of impact with the Earth in 2036 and for many years after.     

Dealing with uncertainties in angles-only initial orbit determination

A method to deal with uncertainties in initial orbit determination (IOD) is presented. This is based on the use of Taylor differential algebra (DA) to nonlinearly map uncertainties from the observation space to the state space. When a minimum set of observations is available, DA is used to expand the solution of the IOD problem in Taylor series with respect to measurement errors. When more observations are available, high order inversion tools are exploited to obtain full state pseudo-observations at a common epoch. The mean and covariance of these pseudo-observations are nonlinearly computed by evaluating the expectation of high order Taylor polynomials. Finally, a linear scheme is employed to update the current knowledge of the orbit. Angles-only observations are considered and simplified Keplerian dynamics adopted to ease the explanation. Three test cases of orbit determination of artificial satellites in different orbital regimes are presented to discuss the feature and performances of the proposed methodology.     

Fast numerics for the spin orbit equation with realistic tidal dissipation and constant eccentricity

We present an algorithm for the rapid numerical integration of a time-periodic ODE with a small dissipation term that is \(C^1\) in the velocity. Such an ODE arises as a model of spin–orbit coupling in a star/planet system, and the motivation for devising a fast algorithm for its solution comes from the desire to estimate probability of capture in various solutions, via Monte Carlo simulation: the integration times are very long, since we are interested in phenomena occurring on timescales of the order of \(10^6\)–\(10^7\) years. The proposed algorithm is based on the high-order Euler method which was described in Bartuccelli et al. (Celest Mech Dyn Astron 121(3):233–260, 2015), and it requires computer algebra to set up the code for its implementation. The payoff is an overall increase in speed by a factor of about 7.5 compared to standard numerical methods. Means for accelerating the purely numerical computation are also discussed.     

Long-term evolution of Oort Cloud comets: methods and comparisons

Two long-term simulation methods for cometary orbits, a Monte Carlo method and a direct integration method, are compared with each other. The comparison is done in seven inclination and perihelion distance intervals, and shows differences in dynamical lifetime and capture probabilities for the following main reasons. We use a finite energy step approximation in the Monte Carlo method and the method considers only close approaches with the planets. The differences can be taken into account statistically and it is possible to calculate the correction factors for the capture probability and dynamical lifetime in the Monte Carlo method. Both corrections depend on the inclination and on the value of the minimum energy step. The capture probabilities of the short-period comets originating in the Oort Cloud are calculated by the corrected Monte Carlo method and compared with published results.     

Monte Carlo simulations of star clusters – II. Tidally limited, multimass systems with stellar evolution

A revision of Stodółkiewicz's Monte Carlo code is used to simulate evolution of large star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. A survey of the evolution of N -body systems influenced by the tidal field of a parent galaxy and by stellar evolution is presented. The process of energy generation is realized by means of appropriately modified versions of Spitzer's and Mikkola's formulae for the interaction cross-section between binaries and field stars and binaries themselves. The results presented are in good agreement with theoretical expectations and the results of other methods (Fokker–Planck, Monte Carlo and N -body). The initial rapid mass loss, resulting from stellar evolution of the most massive stars, causes expansion of the whole cluster and eventually leads to the disruption of less bound systems ( W ₀=3). Models with larger W ₀ survive this phase of evolution and then undergo core collapse and subsequent post-collapse expansion, like isolated models. The expansion phase is eventually reversed when tidal limitation becomes important. The results presented are the first major step in the direction of simulating evolution of real globular clusters by means of the Monte Carlo method.     

High order transfer maps for perturbed Keplerian motion

The paper presents a new semi-analytical technique for the propagation of near-Earth satellite motion. The approach uses differential algebra techniques to compute the high order expansion of the solution of the system’s ordinary differential equation for one orbital revolution, referred to as the transfer map. Once computed, a single high order transfer map (HOTM) can be reused to map an initial condition, or a set of initial conditions, forward in time for many revolutions. The only limiting factor is that the mapped objects must stay close to the reference orbit such that they remain within the region of validity of the HOTM. The performance of the method is assessed through a set of test cases in which both autonomous and non-autonomous perturbations are considered, including the case of continuously propelled trajectories.     

Roche tomography of cataclysmic variables – I. Artefacts and techniques

Roche tomography is a technique used for imaging the Roche-lobe-filling secondary stars in cataclysmic variables (CVs). In order to interpret Roche tomograms correctly, one must determine whether features in the reconstruction are real, or the result of statistical or systematic errors. We explore the effects of systematic errors using reconstructions of simulated data sets, and show that systematic errors result in characteristic distortions of the final reconstructions that can be identified and corrected. In addition, we present a new method of estimating statistical errors on tomographic reconstructions using a Monte Carlo bootstrapping algorithm, and show this method to be much more reliable than Monte Carlo methods which 'jiggle' the data points in accordance with the size of their error bars.     

Rassvet: Backward Monte Carlo radiative transfer in spherical-shell planetary atmospheres

Validation of global numerical models of planetary atmospheres requires simulating images and spectra from the IR to UV spectral regions in order to compare them with remote observations. This paper describes Rassvet, a 3-D spherical-shell backward Monte Carlo radiative transfer model developed for such simulations. It utilizes a new methodology for calculating atmospheric brightness in scattered sunlight by introducing the concept of an “effective emission source”. This allows for the accumulation of the scattered contribution along the entire path of a ray and the calculation of the atmospheric radiation when both scattered sunlight and thermal emission contribute to the remote measurement - which was not possible in previous models. A “polychromatic” algorithm is extended for applications with the backward Monte Carlo method and implemented in the model. It allows for the calculation of radiative intensity for several wavelengths simultaneously, resulting in improved efficiency. The capabilities of the model are demonstrated by simulating remote measurements from the atmosphere of Io.     

Modeling the spatial distribution of neutron stars in the Galaxy

In this paper we investigate the space and velocity distributions of old neutron stars (aged 10⁹ to 10¹⁰?yr) in our Galaxy. Galactic old Neutron Stars (NSs) population fills a torus-like area extending to a few tens kiloparsecs above the galactic plane. The initial velocity distribution of NSs is not well known, in this work we adopt a three component initial distribution, as given by the contribution of kick velocities, circular velocities and Maxwellian velocities. For the spatial initial distribution we use a Γ function. We then use Monte Carlo simulations to follow the evolution of the NSs under the influence of the Paczyński Galactic gravitational potential. Our calculations show that NS orbits have a very large Galactic radial expansion and that their radial distribution peak is quite close to their progenitors’ one. We also study the NS vertical distribution and find that it can well be described by a double exponential low. Finally, we investigate the correlation of the vertical and radial distribution and study the radial dependence of scale-heights.     

Orbits of short period comets captured from the Oort cloud

Oort cloud comets occasionally obtain orbits which take them through the planetary region. The perturbations by the planets are likely to change the orbit of the comet. We model this process by using a Monte Carlo method and cross sections for orbital changes, i.e. changes in energy, inclination and perihelion distance, in a single planet-comet encounter. The influence of all major planets is considered. We study the distributions of orbital parameters of observable comets, i.e. those which have perihelion distance smaller than a given value. We find that enough comets are captured from the Oort cloud in order to explain the present populations of short period comets. The median value of cos i for the Jupiter family is 0.985 while it is 0.27 for the Halley types. The results may explain the orbital features of short period comets, assuming that the active lifetime of a comet is not much greater than 400 orbital revolutions.     

Orbits of short period comets captured from the Oort cloud

Oort cloud comets occasionally obtain orbits which take them through the planetary region. The perturbations by the planets are likely to change the orbit of the comet. We model this process by using a Monte Carlo method and cross sections for orbital changes, i.e. changes in energy, inclination and perihelion distance, in a single planet-comet encounter. The influence of all major planets is considered. We study the distributions of orbital parameters of observable comets, i.e. those which have perihelion distance smaller than a given value. We find that enough comets are captured from the Oort cloud in order to explain the present populations of short period comets. The median value of cos i for the Jupiter family is 0.985 while it is 0.27 for the Halley types. The results may explain the orbital features of short period comets, assuming that the active lifetime of a comet is not much greater than 400 orbital revolutions.     

A steady-state model of NEA binaries formed by tidal disruption of gravitational aggregates

We present results of a simulation of a steady-state binary near-Earth asteroid (NEA) population. This study combines previous work on tidal disruption of gravitational aggregates [Walsh, K.J., Richardson, D.C., 2006. Icarus 180, 201-216] with a Monte Carlo simulation of NEA planetary encounters. Evolutionary effects include tidal evolution and binary disruption from close planetary encounters. The results show that with the best known progenitor (small Main Belt asteroids) shape and spin distributions, and current estimates of NEA lifetime and encounter probabilities, that tidal disruption should account for approximately 1-2% of NEAs being binaries. Given the best observed estimate of a ∼15% binary NEA fraction, we conclude that there are other formation mechanisms that contribute significantly to this population. We also present the expected distribution of binary orbital and physical properties for the steady-state binary NEAs formed by tidal disruption. We discuss the effects on binary fraction and properties due to changes in the least constrained parameters, and other possible effects on our model that could account for differences between the presented results and the observed binary population. Finally, we model possible effects of a significant population of binaries migrating to the near-Earth population from the Main Belt.     

Improving Bayesian analysis for LISA Pathfinder using an efficient Markov Chain Monte Carlo method

We present a parameter estimation procedure based on a Bayesian framework by applying a Markov Chain Monte Carlo algorithm to the calibration of the dynamical parameters of the LISA Pathfinder satellite. The method is based on the Metropolis-Hastings algorithm and a two-stage annealing treatment in order to ensure an effective exploration of the parameter space at the beginning of the chain. We compare two versions of the algorithm with an application to a LISA Pathfinder data analysis problem. The two algorithms share the same heating strategy but with one moving in coordinate directions using proposals from a multivariate Gaussian distribution, while the other uses the natural logarithm of some parameters and proposes jumps in the eigen-space of the Fisher Information matrix. The algorithm proposing jumps in the eigen-space of the Fisher Information matrix demonstrates a higher acceptance rate and a slightly better convergence towards the equilibrium parameter distributions in the application to LISA Pathfinder data. For this experiment, we return parameter values that are all within ～1σ of the injected values. When we analyse the accuracy of our parameter estimation in terms of the effect they have on the force-per-unit of mass noise, we find that the induced errors are three orders of magnitude less than the expected experimental uncertainty in the power spectral density.     

Application of high order expansions of two-point boundary value problems to astrodynamics

Two-point boundary value problems appear frequently in space trajectory design. A remarkable example is represented by the Lambert’s problem, where the conic arc linking two fixed positions in space in a given time is to be characterized in the frame of the two-body problem. Classical methods to numerically solve these problems rely on iterative procedures, which turn out to be computationally intensive in case of lack of good first guesses for the solution. An algorithm to obtain the high order expansion of the solution of a two-point boundary value problem is presented in this paper. The classical iterative procedures are applied to identify a reference solution. Then, differential algebra is used to expand the solution of the problem around the achieved one. Consequently, the computation of new solutions in a relatively large neighborhood of the reference one is reduced to the simple evaluation of polynomials. The performances of the method are assessed by addressing typical applications in the field of spacecraft dynamics, such as the identification of halo orbits and the design of aerocapture maneuvers.     

Element sets for high-order Poincaré mapping of perturbed Keplerian motion

The propagation and Poincaré mapping of perturbed Keplerian motion is a key topic in Celestial Mechanics and Astrodynamics, e.g., to study the stability of orbits or design bounded relative trajectories. The high-order transfer map (HOTM) method enables efficient mapping of perturbed Keplerian orbits using the high-order Taylor expansion of a Poincaré or stroboscopic map. The HOTM is only accurate close to the expansion point and therefore the number of revolutions for which the map is accurate tends to be limited. The proper selection of coordinates is of key importance for improving the performance of the HOTM method. In this paper, we investigate the use of different element sets for expressing the high-order map in order to find the coordinates that perform best in terms of accuracy. A new set of elements is introduced that enables extremely accurate mapping of the state, even for high eccentricities and higher-order zonal perturbations. Finally, the high-order map is shown to be very useful for the determination and study of fixed points and center manifolds of Poincaré maps.     

A systematic fitting scheme for caustic-crossing microlensing events

We outline a method for fitting binary-lens caustic-crossing microlensing events based on the alternative model parametrization proposed and detailed by Cassan. As an illustration of our methodology, we present an analysis of OGLE-2007-BLG-472, a double-peaked Galactic microlensing event with a source crossing the whole caustic structure in less than three days. In order to identify all possible models we conduct an extensive search of the parameter space, followed by a refinement of the parameters with a Markov Chain Monte Carlo algorithm. We find a number of low-  χ²  regions in the parameter space, which lead to several distinct competitive best models. We examine the parameters for each of them, and estimate their physical properties. We find that our fitting strategy locates several minima that are difficult to find with other modelling strategies and is therefore a more appropriate method to fit this type of event.     

On the periodicity hypothesis of the ages of large impact craters

An analysis is made of the periodicity hypothesis of the ages of large craters, based on the compilation by Grieve with the addition of recently identified craters. A method earlier proposed by Broadbent is used to derive a period, and the significance of the derived period is tested by a Monte Carlo experiment. In accordance with the result of Stothers, the ages of large craters  ( D >30 km)  are shown to exhibit a period close to 37.5 Myr. Monte Carlo experiments show, however, that the derived period is far from being statistically significant. A subset of crater data earlier adopted by Napier for the purpose of similar investigation is also tested, and it is shown that they also exhibit a similar period at an almost identical level of confidence. A brief discussion is made of the relation between the derived period and that associated with faunal mass extinctions.