We combine analytical understanding of resonant dynamics in two-planet systems with machine learning techniques to train a model capable of robustly classifying stability in compact multi-planet systems over long timescales of $10^9$ orbits. Our Stability of Planetary Orbital Configurations Klassifier (SPOCK) predicts stability using physically motivated summary statistics measured in integrations of the first $10^4$ orbits, thus achieving speed-ups of up to $10^5$ over full simulations. This computationally opens up the stability constrained characterization of multi-planet systems. Our model, trained on $\approx 100,000$ three-planet systems sampled at discrete resonances, generalizes both to a sample spanning a continuous period-ratio range, as well as to a large five-planet sample with qualitatively different configurations to our training dataset. Our approach significantly outperforms previous methods based on systems' angular momentum deficit, chaos indicators, and parametrized fits to numerical integrations. We use SPOCK to constrain the free eccentricities between the inner and outer pairs of planets in the Kepler-431 system of three approximately Earth-sized planets to both be below 0.05. Our stability analysis provides significantly stronger eccentricity constraints than currently achievable through either radial velocity or transit duration measurements for small planets, and within a factor of a few of systems that exhibit transit timing variations (TTVs). Given that current exoplanet detection strategies now rarely allow for strong TTV constraints (Hadden et al., 2019), SPOCK enables a powerful complementary method for precisely characterizing compact multi-planet systems. We publicly release SPOCK for community use.

In this work we explore the thermodynamic aspects of dark energy in different scenarios: as a perfect fluid with constant and variable equation of state parameter; and as dissipative fluid described by a barotropic equation of state with bulk viscosity in the framework of the Eckart theory and the full Israel-Stewart theory. We explore cosmological solutions for a flat, homogeneous and isotropic universe. When modeled as a perfect fluid with a dynamical equation of state, $p=w(a)\rho$, the dark energy has an energy density, temperature and entropy well defined and an interesting result is that there is no entropy production even though been dynamical. For dissipative dark energy, in the Eckart theory two cases are studied: $\xi=const.$ and $\xi =(\beta/\sqrt{3}) \rho^{1/2}$; it is found that the entropy grows exponentially for the first case and as a power-law for the second. In the Israel-Stewart theory we consider a $\xi =\xi_0 \rho^{1/2}$ and a relaxation time $\tau = \xi/\rho$; an analytical big rip solution is obtained with a power-law entropy. In all cases is obtained a power-law relation between temperature and energy density. In order to maintain the second law of thermodynamics theoretical constraints for the equation of state are found in the different dark energy models studied. A barotropic dark fluid with $w<-1$ is thermodynamically difficult to support, but the overall effect of bulk viscosity in certain cases allows a phantom regime without thermodynamic anomalies.

Direct imaging surveys have found that long-period super-Jupiters are rare. By contrast, recent modeling of the widespread gaps in protoplanetary disks revealed by ALMA suggests an abundant population of smaller Neptune to Jupiter-mass planets at large separations. The thermal emission from such lower-mass planets is negligible at optical and near-infrared wavelengths, leaving only their weak signals in reflected light. Planets do not scatter enough light at these large orbital distances, but there is a natural way to enhance their reflecting area. Each of the four giant planets in our solar system hosts swarms of dozens of irregular satellites, gravitationally captured planetesimals that fill their host planets' spheres of gravitational influence. What we see of them today are the leftovers of an intense collisional evolution. At early times, they would have generated bright circumplanetary debris disks. We investigate the properties and detectability of such irregular satellite disks (ISDs) following models for their collisional evolution from Kennedy and Wyatt 2011. We find that ISD brightnesses peak in the $10-100$ AU range probed by ALMA, and can render planets detectable over a wide range of parameters with upcoming high-contrast instrumentation. We argue that future instruments with wide fields of view could simultaneously characterize the atmospheres of known close-in planets, and reveal the population of long-period $\sim$ Neptune-Jupiter mass exoplanets inaccessible to other detection methods. This provides a complementary and compelling science case that would elucidate the early lives of planetary systems.

Instabilities in compact planetary systems are generically driven by chaotic dynamics. This implies that an instability time measured through direct N-body integration is not exact, but rather represents a single draw from a distribution of equally valid chaotic trajectories. In order to characterize the "errors" on reported instability times from direct N-body integrations, we investigate the shape and parameters of the instability time distributions (ITDs) for ensembles of shadow trajectories that are initially perturbed from one another near machine precision. We find that in the limit where instability times are long compared to the Lyapunov (chaotic) timescale, ITDs approach remarkably similar lognormal distributions with standard deviations ~0.43 $\pm$ 0.16 dex, despite the instability times varying across our sample from $10^4-10^8$ orbits. We find excellent agreement between these predictions, derived from ~450 closely packed configurations of three planets, and a much wider validation set of ~10,000 integrations, as well as on ~20,000 previously published integrations of tightly packed five-planet systems, and a seven-planet resonant chain based on TRAPPIST-1, despite their instability timescales extending beyond our analyzed timescale. We also test the boundary of applicability of our results on dynamically excited versions of our Solar System. These distributions define the fundamental limit imposed by chaos on the predictability of instability times in such planetary systems. It provides a quantitative estimate of the intrinsic error on an N-body instability time imprinted by chaos, approximately a factor of 3 in either direction.

Symplectic integrators have made it possible to study the long-term evolution of planetary systems with direct N-body simulations. In this paper we reassess the accuracy of such simulations by running a convergence test on 20Myr integrations of the Solar System using various symplectic integrators. We find that the specific choice of metric for determining a simulation's accuracy is important. Only looking at metrics related to integrals of motions such as the energy error can overestimate the accuracy of a method. As one specific example, we show that symplectic correctors do not improve the accuracy of secular frequencies compared to the standard Wisdom-Holman method without symplectic correctors, despite the fact that the energy error is three orders of magnitudes smaller. We present a framework to trace the origin of this apparent paradox to one term in the shadow Hamiltonian. Specifically, we find a term that leads to negligible contributions to the energy error but introduces non-oscillatory errors that result in artificial periastron precession. This term is the dominant error when determining secular frequencies of the system. We show that higher order symplectic methods such as the Wisdom-Holman method with a modified kernel or the SABAC family of integrators perform significantly better in secularly evolving systems because they remove this specific term.

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.

Direct N-body simulations and symplectic integrators are effective tools to study the long-term evolution of planetary systems. The Wisdom-Holman (WH) integrator in particular has been used extensively in planetary dynamics as it allows for large timesteps at good accuracy. One can extend the WH method to achieve even higher accuracy using several different approaches. In this paper we survey integrators developed by Wisdom et al. (1996), Laskar & Robutel (2001, and Blanes et al. (2013). Since some of these methods are harder to implement and not as readily available to astronomers compared to the standard WH method, they are not used as often. This is somewhat unfortunate given that in typical simulations it is possible to improve the accuracy by up to six orders of magnitude (!) compared to the standard WH method without the need for any additional force evaluations. To change this, we implement a variety of high order symplectic methods in the freely available N-body integrator REBOUND. In this paper we catalogue these methods, discuss their differences, describe their error scalings, and benchmark their speed using our implementations.

The TRAPPIST-1 system has 7 known terrestrial planets arranged compactly in a mean motion resonant chain around an ultra-cool central star, some within the estimated habitable zone. Given their short orbital periods of just a few days, it is often presumed that the planets are tidally locked such that the spin rate is equal to that of the orbital mean motion. However, the compact, and resonant, nature of the system implies that there can be significant variations in the mean motion of these planets due to their mutual interactions. We show that such fluctuations can then have significant effects on the spin states of these planets. In this paper, we analyze, using detailed numerical simulations, the mean motion histories of the three planets that are thought to lie within or close to the habitable zone of the system: planets d, e, and f. We demonstrate that, depending on the strength of the mutual interactions within the system, these planets can be pushed into spin states which are effectively non-synchronous. We find that it can produce significant libration of the spin state, if not complete circulation in the frame co-rotating with the orbit. We also show that these spin states are likely to be unable to sustain long-term stability, with many of our simulations suggesting that the spin evolves, under the influence of tidal synchronization forces, into quasi-stable attractor states, which last on timescales of thousands of years.

The dark energy component of the universe still remains as a mystery, however, several papers based on observational data have shown that its equation of state may have an oscillatory behaviour. In this paper, we provide a general description for the dark-energy equation-of-state $w(z)$ in the form of Fourier series. This description generalises some previous dynamical dark energy models and is in agreement with the $w(z)$ reconstructions. We make use of a modified version of a simple and fast Markov Chain Monte Carlo code to constraint the model parameters. For the analysis we use data from supernovae type-Ia , baryon acoustic oscillations, $H(z)$ measurements and cosmic microwave background. We provide a comparison of the proposed model with $\Lambda$CDM, $w$CDM and the standard Taylor approximation. The Fourier series expansion of $w(z)$ is preferred from $\Lambda$CDM at more than $3\sigma$ significance level based on the improvement in the fit alone. We use the Akaike criteria to perform the model comparison and found that, even though there are extra parameters, there is a slight preference of the Fourier series compared with the $\Lambda$CDM model. The preferred shape of $w(z)$ found here puts in jeopardy the single scalar field models, as they as they cannot reproduce the crossing the phantom divide line $w=-1$.

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.

Crater counting on the Moon and other bodies is crucial to constrain the dynamical history of the Solar System. This has traditionally been done by visual inspection of images, thus limiting the scope, efficiency, and/or accuracy of retrieval. In this paper we demonstrate the viability of using convolutional neural networks (CNNs) to determine the positions and sizes of craters from Lunar digital elevation maps (DEMs). We recover 92% of craters from the human-generated test set and almost double the total number of crater detections. Of these new craters, 15% are smaller in diameter than the minimum crater size in the ground-truth dataset. Our median fractional longitude, latitude and radius errors are 11% or less, representing good agreement with the human-generated datasets. From a manual inspection of 361 new craters we estimate the false positive rate of new craters to be 11%. Moreover, our Moon-trained CNN performs well when tested on DEM images of Mercury, detecting a large fraction of craters in each map. Our results suggest that deep learning will be a useful tool for rapidly and automatically extracting craters on various Solar System bodies. We make our code and data publicly available at https://github.com/silburt/DeepMoon.git and https://doi.org/10.5281/zenodo.1133969 .

Many of the multi-planet systems discovered to date have been notable for their compactness, with neighbouring planets closer together than any in the Solar System. Interestingly, planet-hosting stars have a wide range of ages, suggesting that such compact systems can survive for extended periods of time. We have used numerical simulations to investigate how quickly systems go unstable in relation to the spacing between planets, focusing on hypothetical systems of Earth-mass planets on evenly-spaced orbits (in mutual Hill radii). In general, the further apart the planets are initially, the longer it takes for a pair of planets to undergo a close encounter. We recover the results of previous studies, showing a linear trend in the initial planet spacing between 3 and 8 mutual Hill radii and the logarithm of the stability time. Investigating thousands of simulations with spacings up to 13 mutual Hill radii reveals distinct modulations superimposed on this relationship in the vicinity of first and second-order mean motion resonances of adjacent and next-adjacent planets. We discuss the impact of this structure and the implications on the stability of compact multi-planet systems. Applying the outcomes of our simulations, we show that isolated systems of up to five Earth-mass planets can fit in the habitable zone of a Sun-like star without close encounters for at least $10^9$ orbits.

In this work we study a phenomenological non-gravitational interaction between dark matter and dark energy. The scenario studied in this work extends the usual interaction model proportional to the derivative of the dark component density adding to the coupling a non-linear term of the form $Q = \rho'/3(\alpha + \beta \rho)$. This dark sector interaction model could be interpreted as a particular case of a running vacuum model of the type $\Lambda(H) = n_0 + n_1 H^2 + n_2 H^4$ in which the vacuum decays into dark matter. For a flat FRW Universe filled with dark energy, dark matter and decoupled baryonic matter and radiation we calculate the energy density evolution equations of the dark sector and solve them. The different sign combinations of the two parameters of the model show clear qualitative different cosmological scenarios, from basic cosmological insights we discard some of them. The linear scalar perturbation equations of the dark matter were calculated. Using the CAMB code we calculate the CMB and matter power spectra for some values of the parameters $\alpha$ and $\beta$ and compare it with $\Lambda$CDM. The model modify mainly the lower multipoles of the CMB power spectrum remaining almost the same the high ones. The matter power spectrum for low wave numbers is not modified by the interaction but after the maximum it is clearly different. Using observational data from Planck, and various galaxy surveys we obtain the constraints of the parameters, the best fit values obtained are the combinations $\alpha = (3.7 \pm 7 )\times 10^{-4} $, $-(1.5\times10^{-5} {\rm eV}^{-1})^{4} \ll \beta < (0.07 {\rm eV}^{-1})^4$.

In single star systems like our own Solar system, comets dominate the mass budget of bodies that are ejected into interstellar space, since they form further away and are less tightly bound. However 1I/`Oumuamua, the first interstellar object detected, appears asteroidal in its spectra and in its lack of detectable activity. We argue that the galactic budget of interstellar objects like 1I/`Oumuamua should be dominated by planetesimal material ejected during planet formation in circumbinary systems, rather than in single star systems or widely separated binaries. We further show that in circumbinary systems, rocky bodies should be ejected in comparable numbers to icy ones. This suggests that a substantial fraction of additional interstellar objects discovered in the future should display an active coma. We find that the rocky population, of which 1I/`Oumuamua seems to be a member, should be predominantly sourced from A-type and late B-star binaries.

On February 6th, 2018 SpaceX launched a Tesla Roadster on a Mars-crossing orbit. We perform N-body simulations to determine the fate of the object over the next 15 Myr. The orbital evolution is initially dominated by close encounters with the Earth. While a precise orbit can not be predicted beyond the next several centuries due to these repeated chaotic scatterings, one can reliably predict the long-term outcomes by statistically analyzing a large suite of possible trajectories with slightly perturbed initial conditions. Repeated gravitational scatterings with Earth lead to a random walk. Collisions with the Earth, Venus and the Sun represent primary sinks for the Roadster's orbital evolution. Collisions with Mercury and Mars, or ejections from the Solar System by Jupiter, are highly unlikely. We calculate a dynamical half-life of the Tesla of approximately 15 Myr, with some 22%, 12% and 12% of Roadster orbit realizations impacting the Earth, Venus, and the Sun within one half-life, respectively. Because the eccentricities and inclinations in our ensemble increase over time due to mean-motion and secular resonances, the impact rates with the terrestrial planets decrease beyond a few million years, whereas the impact rate on the Sun remains roughly constant.

Hamiltonian systems such as the gravitational N-body problem have time-reversal symmetry. However, all numerical N-body integration schemes, including symplectic ones, respect this property only approximately. In this paper, we present the new N-body integrator JANUS, for which we achieve exact time-reversal symmetry by combining integer and floating point arithmetic. JANUS is explicit, formally symplectic and satisfies Liouville's theorem exactly. Its order is even and can be adjusted between two and ten. We discuss the implementation ofJANUS and present tests of its accuracy and speed by performing and analyzing long-term integrations of the Solar System. We show that JANUS is fast and accurate enough to tackle a broad class of dynamical problems. We also discuss the practical and philosophical implications of running exactly time-reversible simulations.

TRAPPIST-1 is a late M-dwarf orbited by seven Earth-sized planets with orbital period ratios near a chain of mean motion resonances. Due to uncertain system parameters, most orbital configurations drawn from the inferred posterior distribution are unstable on short timescales of $\sim$ 0.5 Myr, even when including the eccentricity damping effect of tides. By contrast, we show that most physically plausible resonant configurations generated through disk migration are stable even without tidal dissipation on timescales of at least 50 Myr ($10^{10}$ orbits), an increase of at least two orders of magnitude. This result, together with the remarkable chain of period ratios in the system, provide strong evidence for convergent migration naturally emplacing the system near an equilibrium configuration forced by the resonant chain. We provide an openly available database of physically plausible initial conditions for TRAPPIST-1 generated through parametrized planet-disk interactions, as well as bit-by-bit reproducible N-body integrations over $10^9-10^{10}$ orbits.

The Atacama Large Millimeter/submilimeter Array (ALMA) recently revealed a set of nearly concentric gaps in the protoplanetary disk surrounding the young star HL Tau. If these are carved by forming gas giants, this provides the first set of orbital initial conditions for planets as they emerge from their birth disks. Using N-body integrations, we have followed the evolution of the system for 5 Gyr to explore the possible outcomes. We find that HL Tau initial conditions scaled down to the size of typically observed exoplanet orbits naturally produce several populations in the observed exoplanet sample. First, for a plausible range of planetary masses, we can match the observed eccentricity distribution of dynamically excited radial velocity giant planets with eccentricities $>$ 0.2. Second, we roughly obtain the observed rate of hot Jupiters around FGK stars. Finally, we obtain a large efficiency of planetary ejections of $\approx 2$ per HL Tau-like system, but the small fraction of stars observed to host giant planets makes it hard to match the rate of free-floating planets inferred from microlensing observations. In view of upcoming GAIA results, we also provide predictions for the expected mutual inclination distribution, which is significantly broader than the absolute inclination distributions typically considered by previous studies.

The problem of cosmological production of gravitational waves is discussed in the framework of an expanding, spatially homogeneous and isotropic FRW type Universe with time-evolving vacuum energy density. The gravitational wave equation is established and its modified time-dependent part is analytically resolved for different epochs in the case of a flat geometry. Unlike the standard $\Lambda$CDM cosmology (no interacting vacuum), we show that gravitational waves are produced in the radiation era even in the context of general relativity. We also show that for all values of the free parameter, the high frequency modes are damped out even faster than in the standard cosmology both in the radiation and matter-vacuum dominated epoch. The formation of the stochastic background of gravitons and the remnant power spectrum generated at different cosmological eras are also explicitly evaluated. It is argued that measurements of the CMB polarization (B-modes) and its comparison with the rigid $\Lambda$CDM model plus the inflationary paradigm may become a crucial test for dynamical dark energy models in the near future.

The reproducibility of experiments is one of the main principles of the scientific method. However, numerical N-body experiments, especially those of planetary systems, are currently not reproducible. In the most optimistic scenario, they can only be replicated in an approximate or statistical sense. Even if authors share their full source code and initial conditions, differences in compilers, libraries, operating systems or hardware often lead to qualitatively different results. We provide a new set of easy-to-use, open-source tools that address the above issues, allowing for exact (bit-by-bit) reproducibility of N-body experiments. In addition to generating completely reproducible integrations, we show that our framework also offers novel and innovative ways to analyze these simulations. As an example, we present a high-accuracy integration of the Solar System spanning 10Gyrs, requiring several weeks to run on a modern CPU. In our framework we can not only easily access simulation data at predefined intervals for which we save snapshots, but at any time during the integration. We achieve this by integrating an on-demand reconstructed simulation forward in time from the nearest snapshot. This allows us to extract arbitrary quantities at any point in the saved simulation exactly (bit-by-bit), and within seconds rather than weeks. We believe that the tools we present in this paper offer a new paradigm for how N-body simulations are run, analyzed, and shared across the community.

We investigate the cosmological production of gravitational waves in a nonsingular flat cosmology powered by a "running vacuum" energy density described by $\rho_{\Lambda}\equiv\rho_{\Lambda}(H)$, a phenomenological expression potentially linked with the renormalization group approach in quantum field theory in curved spacetimes. The model can be interpreted as a particular case of the class recently discussed by Perico et al. (Phys. Rev. D {\bf 88}, 063531, 2013) which is termed complete in the sense that the cosmic evolution occurs between two extreme de Sitter stages (early and late time de Sitter phases). {The gravitational wave equation is derived and its time-dependent part numerically integrated since the primordial de Sitter stage. The generated spectrum of gravitons is also compared with the standard calculations where an abrupt transition, from the early de Sitter to the radiation phase, is usually assumed.} It is found that the stochastic background of gravitons is very similar to the one predicted by the cosmic concordance model plus inflation except at higher frequencies ($\nu \gtrsim 100$ kHz). This remarkable signature of a "running vacuum" cosmology combined with the proposed high frequency gravitational wave detectors and measurements of the CMB polarization (B-modes) may provide a new window to confront more conventional models of inflation.

The requirement that planetary systems be dynamically stable is often used to vet new discoveries or set limits on unconstrained masses or orbital elements. This is typically carried out via computationally expensive N-body simulations. We show that characterizing the complicated and multi-dimensional stability boundary of tightly packed systems is amenable to machine learning methods. We find that training an XGBoost machine learning algorithm on physically motivated features yields an accurate classifier of stability in packed systems. On the stability timescale investigated ($10^7$ orbits), it is 3 orders of magnitude faster than direct N-body simulations. Optimized machine learning classifiers for dynamical stability may thus prove useful across the discipline, e.g., to characterize the exoplanet sample discovered by the upcoming Transiting Exoplanet Survey Satellite (TESS). This proof of concept motivates investing computational resources to train algorithms capable of predicting stability over longer timescales and over broader regions of phase space.

Inspired by the recent Kepler discoveries of circumbinary planets orbiting nine close binary stars, we explore the fate of the former as the latter evolve off the main sequence. We combine binary star evolution models with dynamical simulations to study the orbital evolution of these planets as their hosts undergo common-envelope stages, losing in the process a tremendous amount of mass on dynamical timescales. Five of the systems experience at least one Roche-lobe overflow and common-envelope stages (Kepler-1647 experiences three), and the binary stars either shrink to very short orbits or coalesce; two systems trigger a double-degenerate supernova explosion. Kepler's circumbinary planets predominantly remain gravitationally bound at the end of the common-envelope phase, migrate to larger orbits, and may gain significant eccentricity; their orbital expansion can be more than an order of magnitude and can occur over the course of a single planetary orbit. The orbits these planets can reach are qualitatively consistent with those of the currently known post-common-envelope, eclipse-time variations circumbinary candidates. Our results also show that circumbinary planets can experience both modes of orbital expansion (adiabatic and non-adiabatic) if their host binaries undergo more than one common-envelope stage; multiplanet circumbinary systems like Kepler-47 can experience both modes during the same common-envelope stage. Additionally, unlike Mercury orbiting the Sun, a circumbinary planet with the same semi-major axis can survive the common envelope evolution of a close binary star with a total mass of 1 MSun.

In cosmology, phenomenologically motivated expressions for running vacuum are commonly parametrized as linear functions $\Lambda(H^2)$ or $\Lambda(R)$. Such kind of models assume an equation of state for vacuum given by $\,\overline P_\Lambda=-\,\overline\rho_\Lambda$, relating their background pressure $\,\overline P_\Lambda$ and mean energy density $\,\overline\rho_\Lambda\equiv\Lambda/8\pi G$. This equation of state requires that the dynamic for vacuum is due to the energy exchange with the material species. Most of the approaches to background level consider only the energy exchange between vacuum and the transient dominant material component of the universe. We extend such models assuming the running vacuum as the sum of independent contributions $\,\overline\rho_{\Lambda} =\sum_i\,\overline\rho_{\Lambda i}$, associated with (and interacting with) each of the $i$ material species. We derive the linear scalar perturbations for two running scenarios, modeling its cosmic evolution and identifying their different imprints on the cosmic microwave background anisotropies and the matter power spectrum. In the $\Lambda(H^2)$ scenario the running vacuum are coupled with all the material species in the universe, whereas the $\Lambda(R)$ description only leads to coupling between vacuum and the non-relativistic matter components; which produces different imprints of the two models on the matter power spectrum. A comparison with the Planck 2015 data was made in order to constrain the free parameters of the models. In the case of the $\Lambda(H^2)$ model, it was found that $\Omega_\Lambda=0.705\pm0.027$ and $H_0=69.6\pm2.9\, km\, Mpc^{-1}\, s^{-1}$, which diminish the tension with the low redshift expectations.

We present observations at optical wavelengths with the Cassini Spacecraft's Imaging Science System of the Phoebe ring, a vast debris disk around Saturn that seems to be collisionally generated by its irregular satellites. The analysis reveals a radial profile from 80-260 Saturn radii ($R_S$) that changes behavior interior to $\approx 110 R_S$, which we attribute to either the moon Iapetus sweeping up small particles, or to orbital instabilities that cause the ring to flare up vertically. Our study yields an integrated I/F at 0.635 $\mu$m along Saturn's shadow in the Phoebe ring's midplane from 80-250 $R_S$ of $2.7^{+0.9}_{-0.3} \times 10^{-9}$. We develop an analytical model for the size-dependent secular dynamics of retrograde Phoebe ring grains, and compare this model to the observations. This analysis implies that 1) the "Phoebe" ring is partially sourced by debris from irregular satellites beyond Phoebe's orbit and 2) the scattered light signal is dominated by small grains ($\lesssim 20\mu$m in size). If we assume that the Phoebe ring is generated through steady-state micrometeoroid bombardment, this implies a power-law size distribution with index $> 4$, which is unusually steep among solar system rings. This suggests either a steep size distribution of ejecta when material is initially released, or a subsequent process that preferentially breaks up large grains.

First-order variational equations are widely used in N-body simulations to study how nearby trajectories diverge from one another. These allow for efficient and reliable determinations of chaos indicators such as the Maximal Lyapunov characteristic Exponent (MLE) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). In this paper we lay out the theoretical framework to extend the idea of variational equations to higher order. We explicitly derive the differential equations that govern the evolution of second-order variations in the N-body problem. Going to second order opens the door to new applications, including optimization algorithms that require the first and second derivatives of the solution, like the classical Newton's method. Typically, these methods have faster convergence rates than derivative-free methods. Derivatives are also required for Riemann manifold Langevin and Hamiltonian Monte Carlo methods which provide significantly shorter correlation times than standard methods. Such improved optimization methods can be applied to anything from radial-velocity/transit-timing-variation fitting to spacecraft trajectory optimization to asteroid deflection. We provide an implementation of first and second-order variational equations for the publicly available REBOUND integrator package. Our implementation allows the simultaneous integration of any number of first and second-order variational equations with the high-accuracy IAS15 integrator. We also provide routines to generate consistent and accurate initial conditions without the need for finite differencing.

The James Webb Space Telescope (JWST) will provide unprecedented opportunities to observe the rings and small satellites in our solar system, accomplishing three primary objectives: 1) discovering new rings and moons, 2) unprecedented spectroscopy, and 3) time-domain observations. We give details on these science objectives and describe requirements that JWST must fulfill in order to accomplish the science objectives

Models of the dynamical evolution of the early solar system following the dispersal of the gaseous protoplanetary disk have been widely successful in reconstructing the current orbital configuration of the giant planets. Statistically, some of the most successful dynamical evolution simulations have initially included a hypothetical fifth giant planet, of ice giant mass, which gets ejected by a gas giant during the early solar system's proposed instability phase. We investigate the likelihood of an ice giant ejection event by either Jupiter or Saturn through constraints imposed by the current orbits of their wide-separation regular satellites Callisto and Iapetus respectively. We show that planetary encounters that are sufficient to eject an ice giant, often provide excessive perturbations to the orbits of Callisto and Iapetus making it difficult to reconcile a planet ejection event with the current orbit of either satellite. Quantitatively, we compute the likelihood of reconciling a regular Jovian satellite orbit with the current orbit of Callisto following an ice giant ejection by Jupiter of ~ 42% and conclude that such a large likelihood supports the hypothesis of a fifth giant planet's existence. A similar calculation for Iapetus reveals that it is much more difficult for Saturn to have ejected an ice giant and reconcile a Kronian satellite orbit with that of Iapetus (likelihood ~ 1%), although uncertainties regarding the formation of Iapetus, on its unusual orbit, complicates the interpretation of this result.

We present WHFast, a fast and accurate implementation of a Wisdom-Holman symplectic integrator for long-term orbit integrations of planetary systems. WHFast is significantly faster and conserves energy better than all other Wisdom-Holman integrators tested. We achieve this by significantly improving the Kepler-solver and ensuring numerical stability of coordinate transformations to and from Jacobi coordinates. These refinements allow us to remove the linear secular trend in the energy error that is present in other implementations. For small enough timesteps we achieve Brouwer's law, i.e. the energy error is dominated by an unbiased random walk due to floating-point round-off errors. We implement symplectic correctors up to order eleven that significantly reduce the energy error. We also implement a symplectic tangent map for the variational equations. This allows us to efficiently calculate two widely used chaos indicators the Lyapunov characteristic number (LCN) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO). WHFast is freely available as a flexible C package, as a shared library, and as an easy-to-use python module.

A recent ALMA image revealed several concentric gaps in the protoplanetary disk surrounding the young star HL Tau. We consider the hypothesis that these gaps are carved by planets, and present a general framework for understanding the dynamical stability of such systems over typical disk lifetimes, providing estimates for the maximum planetary masses. We collect these easily evaluated constraints into a workflow that can help guide the design and interpretation of new observational campaigns and numerical simulations of gap opening in such systems. We argue that the locations of resonances should be significantly shifted in massive disks like HL Tau, and that theoretical uncertainties in the exact offset, together with observational errors, imply a large uncertainty in the dynamical state and stability in such disks. This presents an important barrier to using systems like HL Tau as a proxy for the initial conditions following planet formation. An important observational avenue to breaking this degeneracy is to search for eccentric gaps, which could implicate resonantly interacting planets. Unfortunately, massive disks like HL Tau should induce swift pericenter precession that would smear out any such eccentric features of planetary origin. This motivates pushing toward more typical, less massive disks. For a nominal non-resonant model of the HL Tau system with five planets, we find a maximum mass for the outer three bodies of approximately 2 Neptune masses. In a resonant configuration, these planets can reach at least the mass of Saturn. The inner two planets' masses are unconstrained by dynamical stability arguments.

The Phoebe ring, Saturn's largest and faintest ring, lies far beyond the planet's well-known main rings. It is primarily sourced by collisions with Saturn's largest irregular satellite Phoebe, perhaps through stochastic macroscopic collisions, or through more steady micrometeoroid bombardment. The ring was discovered with the Spitzer Space Telescope at 24 $\mu$m and has a normal optical depth of $\sim 2 \times 10^{-8}$ (Verbiscer et al. 2009). We report the first observations of sunlight scattered off the Phoebe ring using the Cassini spacecraft's ISS camera at optical wavelengths. We find that material between $\approx 130-210$ Saturnian radii ($R_S$) from the planet produces an I/F of $1.7 \pm 0.1 \times 10^{-11}$ per $R_S$ of the line-of-sight distance through the disk. Combining our measurements with the Spitzer infrared data, we can place constraints on the ring-particles' light-scattering properties. Depending on the particles' assumed phase function, the derived single-scattering albedo can match either photometric models of Phoebe's dark regolith or brighter sub-surface material excavated by macroscopic impacts on Phoebe.

Fomalhaut b is currently the least massive, directly imaged exoplanet candidate. New observation epochs reveal this object to be moving on a highly eccentric orbit (Kalas et al. 2013), which sets important new constraints. I consider scenarios where Fomalhaut b is the only object interacting with the debris disk, and ones involving an additional unseen planet. I also investigate the possibility that Fomalhaut b is merely a transient dust cloud in light of the revised eccentric orbit. I argue that the scenario best able to match the observational constraints is a super-Earth Fomalhaut b surrounded by a vast cloud of dust that is generated by a population of irregular satellites, with an undetected $\sim$ Saturn-mass planet orbiting interior to the disk and driving the secular dynamics. Testable predictions are summarized that could differentiate between this scenario and other possibilities.

Dust grains generated by the Uranian irregular satellites will undergo chaotic large-amplitude eccentricity oscillations under the simultaneous action of radiation forces and the highly misaligned quadrupole potentials of the oblate planet and distant Sun. From a suite of orbital histories, we estimate collision proba- bilities of dust particles with the regular satellites and argue that this process may explain the observed hemispherical color asymmetries of the outermost four regular satellites of Uranus.

High-inclination circumplanetary orbits that are gravitationally perturbed by the central star can undergo Kozai oscillations---large-amplitude, coupled variations in the orbital eccentricity and inclination. We first study how this effect is modified by incorporating perturbations from the planetary oblateness. Tremaine et al. (2009) found that, for planets with obliquities > 68.875 degrees, orbits in the equilibrium local Laplace plane are unstable to eccentricity perturbations over a finite radial range, and execute large-amplitude chaotic oscillations in eccentricity and inclination. In the hope of making that treatment more easily understandable, we analyze the problem using orbital elements, confirming this threshold obliquity. Furthermore, we find that orbits inclined to the Laplace plane will be unstable over a broader radial range, and that such orbits can go unstable for obliquities less than 68.875 degrees. Finally, we analyze the added effects of radiation pressure, which are important for dust grains and provide a natural mechanism for particle semimajor axes to sweep via Poynting-Robertson drag through any unstable range. We find that generally the effect persists; however, the unstable radial range is shifted and small retrograde particles can avoid the instability altogether. We argue that this is occurs because radiation pressure modifies the equilibrium Laplace plane.

The leading face of Saturn's moon Iapetus, Cassini Regio, has an albedo only one tenth that on its trailing side. The origin of this enigmatic dichotomy has been debated for over forty years, but with new data, a clearer picture is emerging. Motivated by Cassini radar and imaging observations, we investigate Soter's model of dark exogenous dust striking an originally brighter Iapetus by modeling the dynamics of the dark dust from the ring of the exterior retrograde satellite Phoebe under the relevant perturbations. In particular, we study the particles' probabilities of striking Iapetus, as well as their expected spatial distribution on the Iapetian surface. We find that, of the long-lived particles (greater than about 5 microns), most particle sizes (greater than about 10 microns) are virtually certain to strike Iapetus, and their calculated distribution on the surface matches up well with Cassini Regio's extent in its longitudinal span. The satellite's polar regions are observed to be bright, presumably because ice is deposited there. Thus, in the latitudinal direction we estimate polar dust deposition rates to help constrain models of thermal migration invoked to explain the bright poles (Spencer & Denk 2010). We also analyze dust originating from other irregular outer moons, determining that a significant fraction of that material will eventually coat Iapetus--perhaps explaining why the spectrum of Iapetus' dark material differs somewhat from that of Phoebe. Finally we track the dust particles that do not strike Iapetus, and find that most land on Titan, with a smaller fraction hitting Hyperion. As has been previously conjectured, such exogenous dust, coupled with Hyperion's chaotic rotation, could produce Hyperion's roughly isotropic, moderate-albedo surface.