$N$-body integrations are used to model a wide range of astrophysical dynamics, but they suffer from errors which make their orbits diverge exponentially in time from the correct orbits. Over long time-scales, their reliability needs to be established. We address this reliability by running a three-body planetary system over about $200$ e-folding times. Using nearby initial conditions, we can construct statistics of the long-term phase-space structure and compare to rough estimates of resonant widths of the system. We compared statistics for a wide range of numerical methods, including a Runge--Kutta method, Wisdom--Holman method, symplectic corrector methods, and a method by Laskar and Robutel. "Improving" an integrator did not increase the phase space accuracy, but simply increasing the number of initial conditions did. In fact, the statistics of a higher order symplectic corrector method were inconsistent with the other methods in one test.

Symplectic methods, in particular the Wisdom-Holman map, have revolutionized our ability to model the long-term, conservative dynamics of planetary systems. However, many astrophysically important effects are dissipative. The consequences of incorporating such forces into otherwise symplectic schemes is not always clear. We show that moving to a general framework of non-commutative operators (dissipative or not) clarifies many of these questions, and that several important properties of symplectic schemes carry over to the general case. In particular, we show that explicit splitting schemes generically exploit symmetries in the applied external forces which often strongly suppress integration errors. Furthermore, we demonstrate that so-called `symplectic correctors' (which reduce energy errors by orders of magnitude at fixed computational cost) apply equally well to weakly dissipative systems and can thus be more generally thought of as `weak splitting correctors.' Finally, we show that previously advocated approaches of incorporating additional forces into symplectic methods work well for dissipative forces, but give qualitatively wrong answers for conservative but velocity-dependent forces like post-Newtonian corrections. We release REBOUNDx, an open-source C library for incorporating additional effects into REBOUND N-body integrations, together with a convenient Python wrapper. All effects are machine-independent and we provide a binary format that interfaces with the SimulationArchive class in REBOUND to enable the sharing and reproducibility of results. Users can add effects from a list of pre-implemented astrophysical forces, or contribute new ones.

Symplectic integrators are a foundation to the study of dynamical $N$-body phenomena, at scales ranging from from planetary to cosmological. These integrators preserve the Poincar\'e invariants of Hamiltonian dynamics. The $N$-body Hamiltonian has another, perhaps overlooked, symmetry: it is smooth, or, in other words, it has infinite differentiability class order (DCO) for particle separations greater than $0$. Popular symplectic integrators, such as hybrid methods or block adaptive stepping methods do not come from smooth Hamiltonians and it is perhaps unclear whether they should. We investigate the importance of this symmetry by considering hybrid integrators, whose DCO can be tuned easily. Hybrid methods are smooth, except at a finite number of phase space points. We study chaotic planetary orbits in a test considered by Wisdom. We find that increasing smoothness, at negligible extra computational cost in particular tests, improves the Jacobi constant error of the orbits by about $5$ orders of magnitude in long-term simulations. The results from this work suggest that smoothness of the $N$-body Hamiltonian is a property worth preserving in simulations.

Symplectic integrators are the preferred method of solving conservative $N$-body problems in cosmological, stellar cluster, and planetary system simulations because of their superior error properties and ability to compute orbital stability. Newtonian gravity is scale free, and there is no preferred time or length scale: this is at odds with construction of traditional symplectic integrators, in which there is an explicit timescale in the time-step. Additional timescales have been incorporated into symplectic integration using various techniques, such as hybrid methods and potential decompositions in planetary astrophysics, integrator sub-cycling in cosmology, and block time-stepping in stellar astrophysics, at the cost of breaking or potentially breaking symplecticity at a few points in phase space. The justification provided, if any, for this procedure is that these trouble points where the symplectic structure is broken should be rarely or never encountered in practice. We consider the case of hybrid integrators, which are used ubiquitously in astrophysics and other fields, to show that symplecticity breaks at a few points are sufficient to destroy beneficial properties of symplectic integrators, which is at odds with some statements in the literature. We show how to solve this problem in the case of hybrid integrators by requiring Lipschitz continuity of the equations of motion. For other techniques, like time step subdivision, {consequences to this this problem are not explored here}, and the fact that symplectic structure is broken should be taken into account by $N$-body simulators, who may find an alternative non-symplectic integrator performs similarly.

Hybrid symplectic integrators such as MERCURY are widely used to simulate complex dynamical phenomena in planetary dynamics that could otherwise not be investigated. A hybrid integrator achieves high accuracy during close encounters by using a high order integration scheme for the duration of the encounter while otherwise using a standard 2nd order Wisdom-Holman scheme, thereby optimizing both speed and accuracy. In this paper we reassess the criteria for choosing the switching function that determines which parts of the Hamiltonian are integrated with the high order integrator. We show that the original motivation for choosing a polynomial switching function in MERCURY is not correct. We explain the nevertheless excellent performance of the MERCURY integrator and then explore a wide range of different switching functions including an infinitely differentiable function and a Heaviside function. We find that using a Heaviside function leads to a significantly simpler scheme compared to MERCURY, while maintaining the same accuracy in short term simulations.

Context. The TRAPPIST-1 system hosts seven Earth-sized, temperate exoplanets orbiting an ultra-cool dwarf star. As such, it represents a remarkable setting to study the formation and evolution of terrestrial planets that formed in the same protoplanetary disk. While the sizes of the TRAPPIST-1 planets are all known to better than 5% precision, their densities have significant uncertainties (between 28% and 95%) because of poor constraints on the planet's masses. Aims.The goal of this paper is to improve our knowledge of the TRAPPIST-1 planetary masses and densities using transit-timing variations (TTV). The complexity of the TTV inversion problem is known to be particularly acute in multi-planetary systems (convergence issues, degeneracies and size of the parameter space), especially for resonant chain systems such as TRAPPIST-1. Methods. To overcome these challenges, we have used a novel method that employs a genetic algorithm coupled to a full N-body integrator that we applied to a set of 284 individual transit timings. This approach enables us to efficiently explore the parameter space and to derive reliable masses and densities from TTVs for all seven planets. Results. Our new masses result in a five- to eight-fold improvement on the planetary density uncertainties, with precisions ranging from 5% to 12%. These updated values provide new insights into the bulk structure of the TRAPPIST-1 planets. We find that TRAPPIST-1\,c and e likely have largely rocky interiors, while planets b, d, f, g, and h require envelopes of volatiles in the form of thick atmospheres, oceans, or ice, in most cases with water mass fractions less than 5%.

Calculating the long term solution of ordinary differential equations, such as those of the $N$-body problem, is central to understanding a wide range of dynamics in astrophysics, from galaxy formation to planetary chaos. Because generally no analytic solution exists to these equations, researchers rely on numerical methods which are prone to various errors. In an effort to mitigate these errors, powerful symplectic integrators have been employed. But symplectic integrators can be severely limited because they are not compatible with adaptive stepping and thus they have difficulty accommodating changing time and length scales. A promising alternative is time-reversible integration, which can handle adaptive time stepping, but the errors due to time-reversible integration in astrophysics are less understood. The goal of this work is to study analytically and numerically the errors caused by time-reversible integration, with and without adaptive stepping. We derive the modified differential equations of these integrators to perform the error analysis. As an example, we consider the trapezoidal rule, a reversible non-symplectic integrator, and show it gives secular energy error increase for a pendulum problem and for a H\'{e}non---Heiles orbit. We conclude that using reversible integration does not guarantee good energy conservation and that, when possible, use of symplectic integrators is favored. We also show that time-symmetry and time-reversibility are properties that are distinct for an integrator.

The symplectic Wisdom-Holman map revolutionized long-term integrations of planetary systems. There is freedom in such methods of how to split the Hamiltonian and which coordinate system to employ, and several options have been proposed in the literature. These choices lead to different integration errors, which we study analytically and numerically. The Wisdom-Holman method in Jacobi coordinates and the method of Hernandez, H16, compare favorably and avoid problems of some of the other maps, such as incorrect center-of-mass position or truncation errors even in the one-planet case. We use H16 to compute the evolution of Pluto's orbital elements over 500 million years in a new calculation.

We study analytically and experimentally certain symplectic and time-reversible N-body integrators which employ a Kepler solver for each pair-wise interaction, including the method of Hernandez & Bertschinger (2015). Owing to the Kepler solver, these methods treat close two-body interactions correctly, while close three-body encounters contribute to the truncation error at second order and above. The second-order errors can be corrected to obtain a fourth-order scheme with little computational overhead. We generalise this map to an integrator which employs a Kepler solver only for selected interactions and yet retains fourth-order accuracy without backward steps. In this case, however, two-body encounters not treated via a Kepler solver contribute to the truncation error.

Formation of globular clusters (GCs), the Galactic bulge, or galaxy bulges in general, are important unsolved problems in Galactic astronomy. Homogeneous infrared observations of large samples of stars belonging to GCs and the Galactic bulge field are one of the best ways to study these problems. We report the discovery by APOGEE of a population of field stars in the inner Galaxy with abundances of N, C, and Al that are typically found in GC stars. The newly discovered stars have high [N/Fe], which is correlated with [Al/Fe] and anti-correlated with [C/Fe]. They are homogeneously distributed across, and kinematically indistinguishable from, other field stars in the same volume. Their metallicity distribution is seemingly unimodal, peaking at [Fe/H]~-1, thus being in disagreement with that of the Galactic GC system. Our results can be understood in terms of different scenarios. N-rich stars could be former members of dissolved GCs, in which case the mass in destroyed GCs exceeds that of the surviving GC system by a factor of ~8. In that scenario, the total mass contained in so-called "first-generation" stars cannot be larger than that in "second-generation" stars by more than a factor of ~9 and was certainly smaller. Conversely, our results may imply the absence of a mandatory genetic link between "second generation" stars and GCs. Last, but not least, N-rich stars could be the oldest stars in the Galaxy, the by-products of chemical enrichment by the first stellar generations formed in the heart of the Galaxy.

We apply one of the exactly symplectic integrators, that we call HB15, of \cite{HB15}, along with the Kepler problem solver of \cite{WH15}, to solve planetary system $N$-body problems. We compare the method to Wisdom-Holman methods (WH) in the \texttt{MERCURY} software package, the \texttt{MERCURY} switching integrator, and others and find HB15 to be the most efficient method or tied for the most efficient method in many cases. Unlike WH, HB15 solved $N$-body problems exhibiting close encounters with small, acceptable error, although frequent encounters slowed the code. Switching maps like \texttt{MERCURY} change between two methods and are not exactly symplectic. We carry out careful tests on their properties and suggest they must be used with caution. We then use different integrators to solve a 3-body problem consisting of a binary planet orbiting a star. For all tested tolerances and time steps, \texttt{MERCURY} unbinds the binary after 0 to 25 years. However, in the solutions of HB15, a time-symmetric Hermite code, and a symplectic Yoshida method, the binary remains bound for $>1000$ years. The methods' solutions are qualitatively different, despite small errors in the first integrals in most cases. Several checks suggest the qualitative binary behavior of HB15's solution is correct. The Bulirsch-Stoer and Radau methods in the \texttt{MERCURY} package also unbind the binary before a time of 50 years, suggesting this dynamical error is due to a \texttt{MERCURY} bug.

We derive and present a fast and accurate solution of the initial value problem for Keplerian motion in universal variables that does not use the Stumpff series. We find that it performs better than methods based on the Stumpff series.

We present a new symplectic integrator designed for collisional gravitational $N$-body problems which makes use of Kepler solvers. The integrator is also reversible and conserves 9 integrals of motion of the $N$-body problem to machine precision. The integrator is second order, but the order can easily be increased by the method of \citeauthor{yos90}. We use fixed time step in all tests studied in this paper to ensure preservation of symplecticity. We study small $N$ collisional problems and perform comparisons with typically used integrators. In particular, we find comparable or better performance when compared to the 4th order Hermite method and much better performance than adaptive time step symplectic integrators introduced previously. We find better performance compared to SAKURA, a non-symplectic, non-time-reversible integrator based on a different two-body decomposition of the $N$-body problem. The integrator is a promising tool in collisional gravitational dynamics.

The SDSS-III/APOGEE survey operated from 2011-2014 using the APOGEE spectrograph, which collects high-resolution (R~22,500), near-IR (1.51-1.70 microns) spectra with a multiplexing (300 fiber-fed objects) capability. We describe the survey data products that are publicly available, which include catalogs with radial velocity, stellar parameters, and 15 elemental abundances for over 150,000 stars, as well as the more than 500,000 spectra from which these quantities are derived. Calibration relations for the stellar parameters (Teff, log g, [M/H], [alpha/M]) and abundances (C, N, O, Na, Mg, Al, Si, S, K, Ca, Ti, V, Mn, Fe, Ni) are presented and discussed. The internal scatter of the abundances within clusters indicates that abundance precision is generally between 0.05 and 0.09 dex across a broad temperature range; within more limited ranges and at high S/N, it is smaller for some elemental abundances. We assess the accuracy of the abundances using comparison of mean cluster metallicities with literature values, APOGEE observations of the solar spectrum and of Arcturus, comparison of individual star abundances with other measurements, and consideration of the locus of derived parameters and abundances of the entire sample, and find that it is challenging to determine the absolute abundance scale; external accuracy may be good to 0.1-0.2 dex. Uncertainties may be larger at cooler temperatures (Teff<4000K). Access to the public data release and data products is described, and some guidance for using the data products is provided.

We map the distribution and properties of the Milky Way's interstellar medium as traced by diffuse interstellar bands (DIBs) detected in near-infrared stellar spectra from the SDSS-III/APOGEE survey. Focusing exclusively on the strongest DIB in the H-band, at ~1.527 microns, we present a projected map of the DIB absorption field in the Galactic plane, using a set of about 60,000 sightlines that reach up to 15 kpc from the Sun and probe up to 30 magnitudes of visual extinction. The strength of this DIB is linearly correlated with dust reddening over three orders of magnitude in both DIB equivalent width (W_DIB) and extinction, with a power law index of 1.01 +/- 0.01, a mean relationship of W_DIB/A_V = 0.1 Angstrom mag^-1, and a dispersion of ~0.05 Angstrom mag^-1 at extinctions characteristic of the Galactic midplane. These properties establish this DIB as a powerful, independent probe of dust extinction over a wide range of A_V values. The subset of about 14,000 robustly detected DIB features have an exponential W_DIB distribution. We empirically determine the intrinsic rest wavelength of this transition to be lambda_0 = 15,272.42 Angstrom, and then calculate absolute radial velocities of the carrier, which display the kinematical signature of the rotating Galactic disk. We probe the DIB carrier distribution in three dimensions and show that it can be characterized by an exponential disk model with a scaleheight of about 100 pc and a scalelength of about 5 kpc. Finally, we show that the DIB distribution also traces large-scale Galactic structures, including the central long bar and the warp of the outer disk.

Neutrino oscillations are only observable when the neutrino production, propagation and detection coherence conditions are satisfied. In this paper we consider in detail neutrino production coherence, taking \pi\to \mu \nu \ decay as an example. We compare the oscillation probabilities obtained in two different ways: (1) coherent summation of the amplitudes of neutrino production at different points along the trajectory of the parent pion; (2) averaging of the standard oscillation probability over the neutrino production coordinate in the source. We demonstrate that the results of these two different approaches exactly coincide, provided that the parent pion is considered as pointlike and the detection process is perfectly localized. In this case the standard averaging of the oscillation probability over the finite spatial extensions of the neutrino source (and detector) properly takes possible decoherence effects into account. We analyze the reason for this equivalence of the two approaches and demonstrate that for pion wave packets of finite width \sigma_{x\pi} the equivalence is broken. The leading order correction to the oscillation probability due to \sigma_{x\pi}\ne 0 is shown to be \sim [v_g/(v_g-v_\pi)]\sigma_{x\pi}/l_{osc}, where v_g and v_\pi \ are the group velocities of the neutrino and pion wave packets, and l_{osc} is the neutrino oscillation length.

We propose an orbifolded, warped, extra dimension scenario in which the visible brane is not parallel to the hidden brane. This leads automatically to Lorentz violation in the visible, four dimensional world. The background solution to the Einstein equations is a function of a parameter that can be identified with the amount of 'tilting' of the brane. The cosmological constant is found to coincide with the classic Randall-Sundrum value to first order in this tilt. Lorentz violating effects induced in the Standard Model are considered. We find that the strongest constraint on the tilt comes from determinations of the electron-proton mass ratio in six quasar spectra (four optical and two radio). Measurements of a third radio source could improve this by an order of magnitude.

The two dark sectors of the universe - dark matter and dark energy - may interact with each other. Background and linear density perturbation evolution equations are developed for a generic coupling. We then establish the general conditions necessary to obtain models free from early time non-adiabatic instabilities. As an application, we consider a viable universe in which the interaction strength is proportional to the dark energy density. The scenario does not exhibit "phantom crossing" and is free from instabilities, including early ones. A sizeable interaction strength is compatible with combined WMAP, HST, SN, LSS and H(z) data. Neutrino mass and/or cosmic curvature are allowed to be larger than in non-interacting models. Our analysis sheds light as well on unstable scenarios previously proposed.