Gaussian processes (GPs) are commonly used as a model of stochastic variability in astrophysical time series. In particular, GPs are frequently employed to account for correlated stellar variability in planetary transit light curves. The efficient application of GPs to light curves containing thousands to tens of thousands of datapoints has been made possible by recent advances in GP methods, including the celerite method. Here we present an extension of the celerite method to two input dimensions, where, typically, the second dimension is small. This method scales linearly with the total number of datapoints when the noise in each large dimension is proportional to the same celerite kernel and only the amplitude of the correlated noise varies in the second dimension. We demonstrate the application of this method to the problem of measuring precise transit parameters from multiwavelength light curves and show that it has the potential to improve transit parameters measurements by orders of magnitude. Applications of this method include transit spectroscopy and exomoon detection, as well a broader set of astronomical problems.

We present an extensive study of spectroscopically-derived chemical abundances for M87 and its globular cluster (GC) system. Using observations from the Mitchell spectrograph at McDonald, LRIS at Keck, and Hectospec on the MMT, we derive new metallicity gradients from $\sim 2$ to $140$ kpc. We use a novel hierarchical statistical framework to simultaneously separate the GC system into subpopulations while measuring the metallicity gradients of those subpopulations. We create physically-motivated spectral stacks of the GC subpopulations by leveraging the output of this statistical framework to perform the first application of abundance tagging in a massive ETG to better constrain the origins of the GC subpopulations and, thus, the assembly history of M87. We find a metal-poor, $\alpha$-enhanced population of GCs in both in the inner and outer halo unanticipated by current cosmological simulations of galaxy evolution. We use the remarkably flat metallicity gradients we find for both the metal-rich and metal-poor GC subpopulations in the inner halo as tentative evidence that some amount of the metal-poor GCs formed directly in the halo of M87 at high redshift.

AU Microscopii (AU Mic) is the second closest pre main sequence star, at a distance of 9.79 parsecs and with an age of 22 million years. AU Mic possesses a relatively rare and spatially resolved3 edge-on debris disk extending from about 35 to 210 astronomical units from the star, and with clumps exhibiting non-Keplerian motion. Detection of newly formed planets around such a star is challenged by the presence of spots, plage, flares and other manifestations of magnetic activity on the star. Here we report observations of a planet transiting AU Mic. The transiting planet, AU Mic b, has an orbital period of 8.46 days, an orbital distance of 0.07 astronomical units, a radius of 0.4 Jupiter radii, and a mass of less than 0.18 Jupiter masses at 3 sigma confidence. Our observations of a planet co-existing with a debris disk offer the opportunity to test the predictions of current models of planet formation and evolution.

The rotational evolution of cool dwarfs is poorly constrained after around 1-2 Gyr due to a lack of precise ages and rotation periods for old main-sequence stars. In this work we use velocity dispersion as an age proxy to reveal the temperature-dependent rotational evolution of low-mass Kepler dwarfs, and demonstrate that kinematic ages could be a useful tool for calibrating gyrochronology in the future. We find that a linear gyrochronology model, calibrated to fit the period-Teff relationship of the Praesepe cluster, does not apply to stars older than around 1 Gyr. Although late-K dwarfs spin more slowly than early-K dwarfs when they are young, at old ages we find that late-K dwarfs rotate at the same rate or faster than early-K dwarfs of the same age. This result agrees qualitatively with semi-empirical models that vary the rate of surface-to-core angular momentum transport as a function of time and mass. It also aligns with recent observations of stars in the NGC 6811 cluster, which indicate that the surface rotation rates of K dwarfs go through an epoch of inhibited evolution. We find that the oldest Kepler stars with measured rotation periods are late-K and early-M dwarfs, indicating that these stars maintain spotted surfaces and stay magnetically active longer than more massive stars. Finally, based on their kinematics, we confirm that many rapidly rotating GKM dwarfs are likely to be synchronized binaries.

Mutual gravitation between a pulsating star and an orbital companion leads to a time-dependent variation in path length for starlight traveling to Earth. These variations can be used for coherently pulsating stars, such as the {\delta} Scuti variables, to constrain the masses and orbits of their companions. Observing these variations for {\delta} Scuti stars has previously relied on subdividing the light curve and measuring the average pulsation phase in equally sized subdivisions, which leads to under-sampling near periapsis. We introduce a new approach that simultaneously forward-models each sample in the light curve and show that this method improves upon current sensitivity limits - especially in the case of highly eccentric and short-period binaries. We find that this approach is sensitive enough to observe Jupiter mass planets around {\delta} Scuti stars under ideal conditions, and use gravity-mode pulsations in the subdwarf B star KIC 7668647 to detect its companion without radial velocity data. We further provide robust detection limits as a function of the SNR of the pulsation mode and determine that the minimum detectable light travel time amplitude for a typical Kepler {\delta} Scuti is around 2 s. This new method significantly enhances the application of light travel time variations to detecting short period binaries with pulsating components, and pulsating A-type exoplanet host stars, especially as a tool for eliminating false positives.

We determine rotation periods for 127 stars in the ~115 Myr old Blanco 1 open cluster using ~200 days of photometric monitoring with the Next Generation Transit Survey (NGTS). These stars span F5-M3 spectral types (1.2 $\gtrsim M \gtrsim$ 0.3 M$_{\odot}$) and increase the number of known rotation periods in Blanco 1 by a factor of four. We determine rotation periods using three methods: Gaussian process (GP) regression, generalised autocorrelation (G-ACF) and Lomb-Scargle (LS) periodograms, and find that GPs and G-ACF are more applicable to evolving spot modulation patterns. Between mid-F and mid-K spectral types, single stars follow a well-defined rotation sequence from ~2 to 10 days, whereas stars in photometric multiple systems typically rotate faster. This may suggest that the presence of a moderate-to-high mass ratio companion inhibits angular momentum loss mechanisms during the early pre-main sequence, and this signature has not been erased at ~100 Myr. The majority of mid-F to mid-K stars display evolving modulation patterns, whereas most M stars show stable modulation signals. This morphological change coincides with the shift from a well-defined rotation sequence (mid-F to mid-K stars) to a broad rotation period distribution (late-K and M stars). Finally, we compare our rotation results for Blanco 1 to the similarly-aged Pleiades: the single star populations in both clusters possess consistent rotation period distributions, which suggests that the angular momentum evolution of stars follows a well-defined pathway that is, at least for mid-F to mid-K stars, strongly imprinted by ~100 Myr.

emcee is a Python library implementing a class of affine-invariant ensemble samplers for Markov chain Monte Carlo (MCMC). This package has been widely applied to probabilistic modeling problems in astrophysics where it was originally published, with some applications in other fields. When it was first released in 2012, the interface implemented in emcee was fundamentally different from the MCMC libraries that were popular at the time, such as PyMC, because it was specifically designed to work with "black box" models instead of structured graphical models. This has been a popular interface for applications in astrophysics because it is often non-trivial to implement realistic physics within the modeling frameworks required by other libraries. Since emcee's release, other libraries have been developed with similar interfaces, such as dynesty (Speagle 2019). The version 3.0 release of emcee is the first major release of the library in about 6 years and it includes a full re-write of the computational backend, several commonly requested features, and a set of new "move" implementations.

We derive analytic, closed-form solutions for the light curve of a planet transiting a star with a limb darkening profile which is a polynomial function of the stellar elevation, up to arbitrary integer order. We provide improved analytic expressions for the uniform, linear, and quadratic limb-darkened cases, as well as novel expressions for higher order integer powers of limb darkening. The formulae are crafted to be numerically stable over the expected range of usage. We additionally present analytic formulae for the partial derivatives of instantaneous flux with respect to the radius ratio, impact parameter, and limb darkening coefficients. These expressions are rapid to evaluate, and compare quite favorably in speed and accuracy to existing transit light curve codes. We also use these expressions to numerically compute the first partial derivatives of exposure-time averaged transit light curves with respect to all model parameters. An additional application is modeling eclipsing binary or eclipsing multiple star systems in cases where the stars may be treated as spherically symmetric. We provide code which implements these formulae in C++, Python, IDL, and Julia, with tests and examples of usage.

Exoplanets orbiting pre-main sequence stars are laboratories for studying planet evolution processes, including atmospheric loss, orbital migration, and radiative cooling. V1298 Tau, a young solar analog with an age of 23 $\pm$ 4 Myr, is one such laboratory. The star is already known to host a Jupiter-sized planet on a 24 day orbit. Here, we report the discovery of three additional planets --- all between the size of Neptune and Saturn --- based on our analysis of K2 Campaign 4 photometry. Planets c and d have sizes of 5.6 and 6.4 $R_\oplus$, respectively and with orbital periods of 8.25 and 12.40 days reside 0.25% outside of the nominal 3:2 mean-motion resonance. Planet e is 8.7 $R_\oplus$ in size but only transited once in the K2 time series and thus has a period longer than 36 days, but likely shorter than 223 days. The V1298 Tau system may be a precursor to the compact multiplanet systems found to be common by the Kepler mission. However, the large planet sizes stand in sharp contrast to the vast majority of Kepler multis which have planets smaller than 3 $R_\oplus$. Simple dynamical arguments suggest total masses of $<$28 $M_\oplus$ and $<$120 $M_\oplus$ for the c-d and d-b planet pairs, respectively. The implied low masses suggest that the planets may still be radiatively cooling and contracting, and perhaps losing atmosphere. The V1298 Tau system offers rich prospects for further follow-up including atmospheric characterization by transmission or eclipse spectroscopy, dynamical characterization through transit-timing variations, and measurements of planet masses and obliquities by radial velocities.

In recent years, dedicated extreme-precision radial velocity (EPRV) spectrographs have produced vast quantities of high-resolution, high-signal-to-noise time-series spectra for bright stars. These data contain valuable information for the dual purposes of planet detection via the measured RVs and stellar characterization via the co-added spectra. However, considerable data analysis challenges exist in extracting these data products from the observed spectra at the highest possible precision, including the issue of poorly-characterized telluric absorption features and the common use of an assumed stellar spectral template. In both of these examples, precision-limiting reliance on external information can be sidestepped using the data directly. Here we propose a data-driven method to simultaneously extract precise RVs and infer the underlying stellar and telluric spectra using a linear model (in the log of flux). The model employs a convex objective and convex regularization to keep the optimization of the spectral components fast. We implement this method in wobble, an open-source python package which uses TensorFlow in one of its first non-neural-network applications to astronomical data. In this work, we demonstrate the performance of wobble on archival HARPS spectra. We recover the canonical exoplanet 51 Pegasi b, detect the secular RV evolution of the M dwarf Barnard's Star, and retrieve the Rossiter-McLaughlin effect for the Hot Jupiter HD 189733b. The method additionally produces extremely high-S/N composite stellar spectra and detailed time-variable telluric spectra, which we also present here.

We implement a search for exoplanets in campaigns zero through eight (C0-8) of the K2 extension of the Kepler spacecraft. We apply a modified version of the QATS planet search algorithm to K2 light curves produced by the EVEREST pipeline, carrying out the C0-8 search on $1.5 \times 10^5$ target stars with magnitudes in the range of Kp = 9-15. We detect 818 transiting planet candidates, of which 374 were undiscovered by prior searches, with {64,15,5,2,1} in {2,3,4,5,6}-planet multi-planet candidate systems, respectively. Of the new planets detected, 100 orbit M dwarfs, including one which is potentially rocky and in the habitable zone. 154 of our candidates reciprocally transit with our Solar System: they are geometrically aligned to see at least one Solar System planet transit. We find candidates which display transit timing variations and dozens of candidates on both period extremes with single transits or ultra-short periods. We point to evidence that our candidates display similar patterns in frequency and size-period relation as confirmed planets, such as tentative evidence for the radius gap. Confirmation of these planet candidates with follow-up studies will increase the number of K2 planets by up to 50%, and characterization of their host stars will improve statistical studies of planet properties. Our sample includes many planets orbiting bright stars amenable for radial velocity follow-up and future characterization with JWST. We also list the 579 eclipsing binary systems detected as part of this search.

We present a new age-dating technique that combines gyrochronology with isochrone fitting to infer ages for FGKM main-sequence and subgiant field stars. Gyrochronology and isochrone fitting are each capable of providing relatively precise ages for field stars in certain areas of the Hertzsprung-Russell diagram: gyrochronology works optimally for cool main-sequence stars, and isochrone fitting can provide precise ages for stars near the main-sequence turnoff. Combined, these two age-dating techniques can provide precise and accurate ages for a broader range of stellar masses and evolutionary stages than either method used in isolation. We demonstrate that the position of a star on the Hertzsprung- Russell or color-magnitude diagram can be combined with its rotation period to infer a precise age via both isochrone fitting and gyrochronology simultaneously. We show that incorporating rotation periods with 5% uncertainties into stellar evolution models improves age precision for FGK stars on the main sequence, and can, on average, provide age estimates up to three times more precise than isochrone fitting alone. In addition, we provide a new gyrochronology relation, calibrated to the Praesepe cluster and the Sun, that includes a variance model to capture the rotational behavior of stars whose rotation periods do not lengthen with the square-root of time, and parts of the Hertzsprung-Russell diagram where gyrochronology has not been calibrated. This publication is accompanied by an open source Python package, stardate, for inferring the ages of main-sequence and subgiant FGKM stars from rotation periods, spectroscopic parameters and/or apparent magnitudes and parallaxes.

During its two year prime mission the Transiting Exoplanet Survey Satellite (TESS) will perform a time-series photometric survey covering over 80% of the sky. This survey comprises observations of 26 24 x 96 degree sectors that are each monitored continuously for approximately 27 days. The main goal of TESS is to find transiting planets around 200,000 pre-selected stars for which fixed aperture photometry is recorded every two minutes. However, TESS is also recording and delivering Full-Frame Images (FFIs) of each detector at a 30 minute cadence. We have created an open-source tool, eleanor, to produce light curves for objects in the TESS FFIs. Here, we describe the methods used in eleanor to produce light curves that are optimized for planet searches. The tool performs background subtraction, aperture and PSF photometry, decorrelation of instrument systematics, and cotrending using principal component analysis. We recover known transiting exoplanets in the FFIs to validate the pipeline and perform a limited search for new planet candidates in Sector 1. Our tests indicate that eleanor produces light curves with significantly less scatter than other tools that have been used in the literature. Cadence-stacked images, and raw and detrended eleanor light curves for each analyzed star will be hosted on MAST, with planet candidates on ExoFOP-TESS as Community TESS Objects of Interest (CTOIs). This work confirms the promise that the TESS FFIs will enable the detection of thousands of new exoplanets and a broad range of time domain astrophysics.

We report the Transiting Exoplanet Survey Satellite (TESS) discovery of three terrestrial-sized planets transiting L 98-59 (TOI-175, TIC 307210830) -- a bright M dwarf at a distance of 10.6 pc. Using the Gaia-measured distance and broad-band photometry we find that the host star is an M3 dwarf. Combined with the TESS transits from three sectors, the corresponding stellar parameters yield planet radii ranging from 0.8REarth to 1.6REarth. All three planets have short orbital periods, ranging from 2.25 to 7.45 days with the outer pair just wide of a 2:1 period resonance. Diagnostic tests produced by the TESS Data Validation Report and the vetting package DAVE rule out common false positive sources. These analyses, along with dedicated follow-up and the multiplicity of the system, lend confidence that the observed signals are caused by planets transiting L 98-59 and are not associated with other sources in the field. The L 98-59 system is interesting for a number of reasons: the host star is bright (V = 11.7 mag, K = 7.1 mag) and the planets are prime targets for further follow-up observations including precision radial-velocity mass measurements and future transit spectroscopy with the James Webb Space Telescope; the near resonant configuration makes the system a laboratory to study planetary system dynamical evolution; and three planets of relatively similar size in the same system present an opportunity to study terrestrial planets where other variables (age, metallicity, etc.) can be held constant. L 98-59 will be observed in 4 more TESS sectors, which will provide a wealth of information on the three currently known planets and have the potential to reveal additional planets in the system.

Over the past century, major advances in astronomy and astrophysics have been largely driven by improvements in instrumentation and data collection. With the amassing of high quality data from new telescopes, and especially with the advent of deep and large astronomical surveys, it is becoming clear that future advances will also rely heavily on how those data are analyzed and interpreted. New methodologies derived from advances in statistics, computer science, and machine learning are beginning to be employed in sophisticated investigations that are not only bringing forth new discoveries, but are placing them on a solid footing. Progress in wide-field sky surveys, interferometric imaging, precision cosmology, exoplanet detection and characterization, and many subfields of stellar, Galactic and extragalactic astronomy, has resulted in complex data analysis challenges that must be solved to perform scientific inference. Research in astrostatistics and astroinformatics will be necessary to develop the state-of-the-art methodology needed in astronomy. Overcoming these challenges requires dedicated, interdisciplinary research. We recommend: (1) increasing funding for interdisciplinary projects in astrostatistics and astroinformatics; (2) dedicating space and time at conferences for interdisciplinary research and promotion; (3) developing sustainable funding for long-term astrostatisics appointments; and (4) funding infrastructure development for data archives and archive support, state-of-the-art algorithms, and efficient computing.

We determine the orbital eccentricities of individual small Kepler planets, through a combination of asteroseismology and transit light-curve analysis. We are able to constrain the eccentricities of 51 systems with a single transiting planet, which supplement our previous measurements of 66 planets in multi-planet systems. Through a Bayesian hierarchical analysis, we find evidence that systems with only one detected transiting planet have a different eccentricity distribution than systems with multiple detected transiting planets. The eccentricity distribution of the single-transiting systems is well described by the positive half of a zero-mean Gaussian distribution with a dispersion $\sigma_e = 0.32 \pm 0.06$, while the multiple-transit systems are consistent with $\sigma_e = 0.083^{+0.015}_{-0.020}$. A mixture model suggests a fraction of $0.76^{+0.21}_{-0.12}$ of single-transiting systems have a moderate eccentricity, represented by a Rayleigh distribution that peaks at $0.26^{+0.04}_{-0.06}$. This finding may reflect differences in the formation pathways of systems with different numbers of transiting planets. We investigate the possibility that eccentricities are "self-excited" in closely packed planetary systems, as well as the influence of long-period giant companion planets. We find that both mechanisms can qualitatively explain the observations. We do not find any evidence for a correlation between eccentricity and stellar metallicity, as has been seen for giant planets. Neither do we find any evidence that orbital eccentricity is linked to the detection of a companion star. Along with this paper we make available all of the parameters and uncertainties in the eccentricity distributions, as well as the properties of individual systems, for use in future studies.

In systems with detected planets, hot-Jupiters and compact systems of multiple planets are nearly mutually exclusive. We compare the relative occurrence of these two architectures as a fraction of detected planetary systems to determine the role that metallicity plays in planet formation. We show that compact multi-planet systems occur more frequently around stars of increasingly lower metallicities using spectroscopically derived abundances for more than 700 planet hosts. At higher metallicities, compact multi-planet systems comprise a nearly constant fraction of the planet hosts despite the steep rise in the fraction of hosts containing hot and cool-Jupiters. Since metal poor stars have been underrepresented in planet searches, this implies that the occurrence rate of compact multis is higher than previously reported. Due to observational limits, radial velocity planet searches have focused mainly on high-metallicity stars where they have a higher chance of finding giant planets. New extreme-precision radial velocity instruments coming online that can detect these compact multi-planet systems can target lower metallicity stars to find them.

We derive analytic, closed form, numerically stable solutions for the total flux received from a spherical planet, moon or star during an occultation if the specific intensity map of the body is expressed as a sum of spherical harmonics. Our expressions are valid to arbitrary degree and may be computed recursively for speed. The formalism we develop here applies to the computation of stellar transit light curves, planetary secondary eclipse light curves, and planet-planet/planet-moon occultation light curves, as well as thermal (rotational) phase curves. In this paper we also introduce STARRY, an open-source package written in C++ and wrapped in Python that computes these light curves. The algorithm in STARRY is six orders of magnitude faster than direct numerical integration and several orders of magnitude more precise. STARRY also computes analytic derivatives of the light curves with respect to all input parameters for use in gradient-based optimization and inference, such as Hamiltonian Monte Carlo (HMC), allowing users to quickly and efficiently fit observed light curves to infer properties of a celestial body's surface map.

Measuring the centroid of a spectral line is a common problem in astronomy. Many methods have been devised to overcome limitations due to either noise in the spectra or asymmetric profiles, the most common of which are the intensity weighted averages (first moment) or fits of analytical (typically Gaussian) profiles. If the spectral line can be considered a single component, we demonstrate that a simple quadratic fit to the pixel of maximum intensity and its two neighboring pixels provides a robust measure of the line centroid. This approach allows for a sub-velocity resolution precision on the line centroid, without be biases by noise or asymmetric features in the line profile and outperforming traditional methods in most situations.

We present the first kinematical detection of embedded protoplanets within a protoplanetary disk. Using archival ALMA observations of HD 163296, we demonstrate a new technique to measure the rotation curves of CO isotopologue emission to sub-percent precision relative to the Keplerian rotation. These rotation curves betray substantial deviations caused by local perturbations in the radial pressure gradient, likely driven by gaps carved in the gas surface density by Jupiter-mass planets. Comparison with hydrodynamic simulations shows excellent agreement with the gas rotation profile when the disk surface density is perturbed by two Jupiter mass planets at 83 au and 137 au. As the rotation of the gas is dependent on the pressure of the total gas component, this method provides a unique probe of the gas surface density profile without incurring significant uncertainties due to gas-to-dust ratios or local chemical abundances which plague other methods. Future analyses combining both methods promise to provide the most accurate and robust measures of embedded planetary mass. Furthermore, this method provides a unique opportunity to explore wide-separation planets beyond the mm continuum edge and to trace the gas pressure profile essential in modelling grain evolution in disks.

It has been demonstrated that the time variability of a star's brightness at different frequencies can be used to infer its surface gravity, radius, mass, and age. With large samples of light curves now available from Kepler and K2, and upcoming surveys like TESS, we wish to quantify the overall information content of this data and identify where the information resides. As a first look into this question we ask which stellar parameters we can predict from the brightness variations in red-giant stars data and to what precision, using a data-driven model. We demonstrate that the long-cadence (30-minute) Kepler light curves for 2000 red-giant stars can be used to predict their $T_{\rm eff}$ and $\log g$. Our inference makes use of a data-driven model of a part of the autocorrelation function (ACF) of the light curve, where we posit a polynomial relationship between stellar parameters and the ACF pixel values. We find that this model, trained using 1000 stars, can be used to recover the temperature $T_{\rm eff}$ to $<$100 K, the surface gravity to $<$ 0.1 dex, and the asteroseismic power-spectrum parameters $\rm \nu_{max}$ and $\rm \Delta{\nu}$ to $<11$ $\mu$Hz and $<0.9$ $\mu$Hz ($\lesssim$ 15\%). We recover $T_{\rm eff}$ from range of time-lags 0.045 $<$ $T_{\rm lag}$ $<$ 370 days and the $\log g$, $\rm \nu_{max}$ and $\rm \Delta{\nu}$ from the range 0.045 $<$ $T_{\rm lag}$ $<$ 35 days. We do not discover any information about stellar metallicity. The information content of the data about each parameter is empirically quantified using this method, enabling comparisons to theoretical expectations about convective granulation.

We analyze the evolution of the potentially habitable planet Proxima Centauri b to identify environmental factors that affect its long-term habitability. We consider physical processes acting on size scales ranging from the galactic to the stellar system to the planet's core. We find that there is a significant probability that Proxima Centauri has had encounters with its companion stars, Alpha Centauri A and B, that are close enough to destabilize an extended planetary system. If the system has an additional planet, as suggested by the discovery data, then it may perturb planet b's eccentricity and inclination, possibly driving those parameters to non-zero values, even in the presence of strong tidal damping. We also model the internal evolution of the planet, evaluating the roles of different radiogenic abundances and tidal heating and find that magnetic field generation is likely for billions of years. We find that if planet b formed in situ, then it experienced 169 +/- 13 million years in a runaway greenhouse as the star contracted during its formation. This early phase could remove up to 5 times as much water as in the modern Earth's oceans, possibly producing a large abiotic oxygen atmosphere. On the other hand, if Proxima Centauri b formed with a substantial hydrogen atmosphere (0.01 - 1% of the planet's mass), then this envelope could have shielded the water long enough for it to be retained before being blown off itself. After modeling this wide range of processes we conclude that water retention during the host star's pre-main sequence phase is the biggest obstacle for Proxima b's habitability. These results are all obtained with a new software package called VPLANET.

This research note presents a derivation and implementation of efficient and scalable gradient computations using the celerite algorithm for Gaussian Process (GP) modeling. The algorithms are derived in a "reverse accumulation" or "backpropagation" framework and they can be easily integrated into existing automatic differentiation frameworks to provide a scalable method for evaluating the gradients of the GP likelihood with respect to all input parameters. The algorithm derived in this note uses less memory and is more efficient than versions using automatic differentiation and the computational cost scales linearly with the number of data points.

Photometry from the Kepler mission is optimized to detect small, short duration signals like planet transits at the expense of long-term trends. This long-term variability can be recovered in photometry from the Full Frame Images (FFIs), a set of calibration data collected approximately monthly during the Kepler mission. Here, we present f3, an open-source package to perform photometry on the Kepler FFIs in order to detect changes in the brightness of stars in the Kepler field of view over long time baselines. We apply this package to a sample of 4,000 Sun-like stars with measured rotation periods. We find $\approx 10\%$ of these targets have long-term variability in their observed flux. For the majority of targets we identify the luminosity variations are either correlated or anticorrelated with the short-term variability due to starspots on the stellar surface. We find a transition between anticorrelated (starspot-dominated) variability and correlated (facula-dominated) variability between rotation periods of 15 and 25 days, suggesting the transition between the two modes is complete for stars at the age of the Sun. We also identify a sample of stars with apparently complete cycles, as well as a collection of short-period binaries with extreme photometric variation over the Kepler mission.

The target of many astronomical studies is the recovery of tiny astrophysical signals living in a sea of uninteresting (but usually dominant) noise. In many contexts (i.e., stellar time-series, or high-contrast imaging, or stellar spectroscopy), there are structured components in this noise caused by systematic effects in the astronomical source, the atmosphere, the telescope, or the detector. More often than not, evaluation of the true physical model for these nuisances is computationally intractable and dependent on too many (unknown) parameters to allow rigorous probabilistic inference. Sometimes, housekeeping data---and often the science data themselves---can be used as predictors of the systematic noise. Linear combinations of simple functions of these predictors are often used as computationally tractable models that can capture the nuisances. These models can be used to fit and subtract systematics prior to investigation of the signals of interest, or they can be used in a simultaneous fit of the systematics and the signals. In this Note, we show that if a Gaussian prior is placed on the weights of the linear components, the weights can be marginalized out with an operation in pure linear algebra, which can (often) be made fast. We illustrate this model by demonstrating the applicability of a linear model for the non-linear systematics in K2 time-series data, where the dominant noise source for many stars is spacecraft motion and variability.

Markov Chain Monte Carlo (MCMC) methods for sampling probability density functions (combined with abundant computational resources) have transformed the sciences, especially in performing probabilistic inferences, or fitting models to data. In this primarily pedagogical contribution, we give a brief overview of the most basic MCMC method and some practical advice for the use of MCMC in real inference problems. We give advice on method choice, tuning for performance, methods for initialization, tests of convergence, troubleshooting, and use of the chain output to produce or report parameter estimates with associated uncertainties. We argue that autocorrelation time is the most important test for convergence, as it directly connects to the uncertainty on the sampling estimate of any quantity of interest. We emphasize that sampling is a method for doing integrals; this guides our thinking about how MCMC output is best used.

Despite more than 20 years since the discovery of the first gas giant planet with an anomalously large radius, the mechanism for planet inflation remains unknown. Here, we report the discovery of EPIC228754001.01, an inflated gas giant planet found with the NASA K2 Mission, and a revised mass for another inflated planet, K2-97b. These planets reside on ~9 day orbits around host stars which recently evolved into red giants. We constrain the irradiation history of these planets using models constrained by asteroseismology and Keck/HIRES spectroscopy and radial velocity measurements. We measure planet radii of 1.31 +\- 0.11 Rjup and and 1.30 +\- 0.07 Rjup, respectively. These radii are typical for planets receiving the current irradiation, but not the former, zero age main sequence irradiation of these planets. This suggests that the current sizes of these planets are directly correlated to their current irradiation. Our precise constraints of the masses and radii of the stars and planets in these systems allow us to constrain the planetary heating efficiency of both systems as 0.03% +0.03%/-0.02%. These results are consistent with a planet re-inflation scenario, but suggest the efficiency of planet re-inflation may be lower than previously theorized. Finally, we discuss the agreement within 10% of stellar masses and radii, and planet masses, radii, and orbital periods of both systems and speculate that this may be due to selection bias in searching for planets around evolved stars.

Difference imaging or image subtraction is a method that measures differential photometry by matching the pointing and point-spread function (PSF) between image frames. It is used for the detection of time-variable phenomena. Here we present a new category of method---CPM Difference Imaging, in which differences are not measured between matched images but instead between image frames and a data-driven predictive model that has been designed only to predict the pointing, PSF, and detector effects but not astronomical variability. In CPM Difference Imaging each pixel is modelled by the Causal Pixel Model (CPM) originally built for modeling Kepler data, in which pixel values are predicted by a linear combination of other pixels at the same epoch but far enough away such that these pixels are causally disconnected, astrophysically. It does not require that the user have any explicit model or description of the pointing or point-spread function of any of the images. Its principal drawback is that---in its current form---it requires an imaging campaign with many epochs and fairly stable telescope pointing. The method is applied to simulated data and also the K2 Campaign 9 microlensing data. We show that CPM Difference Imaging can detect variable objects and produce precise differentiate photometry in a crowded field. CPM Difference Imaging is capable of producing image differences at nearly photon-noise precision.

The growing field of large-scale time domain astronomy requires methods for probabilistic data analysis that are computationally tractable, even with large datasets. Gaussian Processes are a popular class of models used for this purpose but, since the computational cost scales, in general, as the cube of the number of data points, their application has been limited to small datasets. In this paper, we present a novel method for Gaussian Process modeling in one-dimension where the computational requirements scale linearly with the size of the dataset. We demonstrate the method by applying it to simulated and real astronomical time series datasets. These demonstrations are examples of probabilistic inference of stellar rotation periods, asteroseismic oscillation spectra, and transiting planet parameters. The method exploits structure in the problem when the covariance function is expressed as a mixture of complex exponentials, without requiring evenly spaced observations or uniform noise. This form of covariance arises naturally when the process is a mixture of stochastically-driven damped harmonic oscillators -- providing a physical motivation for and interpretation of this choice -- but we also demonstrate that it can be a useful effective model in some other cases. We present a mathematical description of the method and compare it to existing scalable Gaussian Process methods. The method is fast and interpretable, with a range of potential applications within astronomical data analysis and beyond. We provide well-tested and documented open-source implementations of this method in C++, Python, and Julia.

Variability in the light curves of spotted, rotating stars is often non-sinusoidal and quasi-periodic --- spots move on the stellar surface and have finite lifetimes, causing stellar flux variations to slowly shift in phase. A strictly periodic sinusoid therefore cannot accurately model a rotationally modulated stellar light curve. Physical models of stellar surfaces have many drawbacks preventing effective inference, such as highly degenerate or high-dimensional parameter spaces. In this work, we test an appropriate effective model: a Gaussian Process with a quasi-periodic covariance kernel function. This highly flexible model allows sampling of the posterior probability density function of the periodic parameter, marginalising over the other kernel hyperparameters using a Markov Chain Monte Carlo approach. To test the effectiveness of this method, we infer rotation periods from 333 simulated stellar light curves, demonstrating that the Gaussian process method produces periods that are more accurate than both a sine-fitting periodogram and an autocorrelation function method. We also demonstrate that it works well on real data, by inferring rotation periods for 275 Kepler stars with previously measured periods. We provide a table of rotation periods for these 1132 Kepler objects of interest and their posterior probability density function samples. Because this method delivers posterior probability density functions, it will enable hierarchical studies involving stellar rotation, particularly those involving population modelling, such as inferring stellar ages, obliquities in exoplanet systems, or characterising star-planet interactions. The code used to implement this method is available online.

The TRAPPIST-1 system is the first transiting planet system found orbiting an ultra-cool dwarf star. At least seven planets similar to Earth in radius and in mass were previously found to transit this host star. Subsequently, TRAPPIST-1 was observed as part of the K2 mission and, with these new data, we report the measurement of an 18.77 d orbital period for the outermost planet, TRAPPIST-1h, which was unconstrained until now. This value matches our theoretical expectations based on Laplace relations and places TRAPPIST-1h as the seventh member of a complex chain, with three-body resonances linking every member. We find that TRAPPIST-1h has a radius of 0.727 Earth radii and an equilibrium temperature of 173 K. We have also measured the rotational period of the star at 3.3 d and detected a number of flares consistent with a low-activity, middle-aged, late M dwarf.

Given sparse or low-quality radial-velocity measurements of a star, there are often many qualitatively different stellar or exoplanet companion orbit models that are consistent with the data. The consequent multimodality of the likelihood function leads to extremely challenging search, optimization, and MCMC posterior sampling over the orbital parameters. Here we create a custom Monte Carlo sampler for sparse or noisy radial-velocity measurements of two-body systems that can produce posterior samples for orbital parameters even when the likelihood function is poorly behaved. The six standard orbital parameters for a binary system can be split into four non-linear parameters (period, eccentricity, argument of pericenter, phase) and two linear parameters (velocity amplitude, barycenter velocity). We capitalize on this by building a sampling method in which we densely sample the prior pdf in the non-linear parameters and perform rejection sampling using a likelihood function marginalized over the linear parameters. With sparse or uninformative data, the sampling obtained by this rejection sampling is generally multimodal and dense. With informative data, the sampling becomes effectively unimodal but too sparse: in these cases we follow the rejection sampling with standard MCMC. The method produces correct samplings in orbital parameters for data that include as few as three epochs. The Joker can therefore be used to produce proper samplings of multimodal pdfs, which are still informative and can be used in hierarchical (population) modeling. We give some examples that show how the posterior pdf depends sensitively on the number and time coverage of the observations and their uncertainties.

We present an update to the EVEREST K2 pipeline that addresses various limitations in the previous version and improves the photometric precision of the de-trended light curves. We develop a fast regularization scheme for third order pixel level decorrelation (PLD) and adapt the algorithm to include the PLD vectors of neighboring stars to enhance the predictive power of the model and minimize overfitting, particularly for faint stars. We also modify PLD to work for saturated stars and improve its performance on extremely variable stars. On average, EVEREST 2.0 light curves have 10-20% higher photometric precision than those in the previous version, yielding the highest precision light curves at all Kp magnitudes of any publicly available K2 catalog. For most K2 campaigns, we recover the original Kepler precision to at least Kp = 14, and to at least Kp = 15 for campaigns 1, 5, and 6. We also de-trend all short cadence targets observed by K2, obtaining even higher photometric precision for these stars. All light curves for campaigns 0-8 are available online in the EVEREST catalog, which will be continuously updated with future campaigns. EVEREST 2.0 is open source and is coded in a general framework that can be applied to other photometric surveys, including Kepler and the upcoming TESS mission.

In probabilistic (Bayesian) inferences, we typically want to compute properties of the posterior distribution, describing knowledge of unknown quantities in the context of a particular dataset and the assumed prior information. The marginal likelihood, also known as the "evidence", is a key quantity in Bayesian model selection. The Diffusive Nested Sampling algorithm, a variant of Nested Sampling, is a powerful tool for generating posterior samples and estimating marginal likelihoods. It is effective at solving complex problems including many where the posterior distribution is multimodal or has strong dependencies between variables. DNest4 is an open source (MIT licensed), multi-threaded implementation of this algorithm in C++11, along with associated utilities including: i) RJObject, a class template for finite mixture models, (ii) A Python package allowing basic use without C++ coding, and iii) Experimental support for models implemented in Julia. In this paper we demonstrate DNest4 usage through examples including simple Bayesian data analysis, finite mixture models, and Approximate Bayesian Computation.

The Kepler Mission has discovered thousands of exoplanets and revolutionized our understanding of their population. This large, homogeneous catalog of discoveries has enabled rigorous studies of the occurrence rate of exoplanets and planetary systems as a function of their physical properties. However, transit surveys like Kepler are most sensitive to planets with orbital periods much shorter than the orbital periods of Jupiter and Saturn, the most massive planets in our Solar System. To address this deficiency, we perform a fully automated search for long-period exoplanets with only one or two transits in the archival Kepler light curves. When applied to the $\sim 40,000$ brightest Sun-like target stars, this search produces 16 long-period exoplanet candidates. Of these candidates, 6 are novel discoveries and 5 are in systems with inner short-period transiting planets. Since our method involves no human intervention, we empirically characterize the detection efficiency of our search. Based on these results, we measure the average occurrence rate of exoplanets smaller than Jupiter with orbital periods in the range 2-25 years to be $2.0\pm0.7$ planets per Sun-like star.

Chemical tagging promises to use detailed abundance measurements to identify spatially separated stars that were in fact born together (in the same molecular cloud), long ago. This idea has not yielded much practical success, presumably because of the noise and incompleteness in chemical-abundance measurements. We have succeeded in substantially improving spectroscopic measurements with The Cannon, which has now delivered 15 individual abundances for ~100,000 stars observed as part of the APOGEE spectroscopic survey, with precisions around 0.04 dex. We test the chemical-tagging hypothesis by looking at clusters in abundance space and confirming that they are clustered in phase space. We identify (by the k-means algorithm) overdensities of stars in the 15-dimensional chemical-abundance space delivered by The Cannon, and plot the associated stars in phase space. We use only abundance-space information (no positional information) to identify stellar groups. We find that clusters in abundance space are indeed clusters in phase space. We recover some known phase-space clusters and find other interesting structures. This is the first-ever project to identify phase-space structures at survey-scale by blind search purely in abundance space; it verifies the precision of the abundance measurements delivered by The Cannon; the prospects for future data sets appear very good.

We present EVEREST, an open-source pipeline for removing instrumental noise from K2 light curves. EVEREST employs a variant of pixel level decorrelation (PLD) to remove systematics introduced by the spacecraft's pointing error and a Gaussian process (GP) to capture astrophysical variability. We apply EVEREST to all K2 targets in campaigns 0-7, yielding light curves with precision comparable to that of the original Kepler mission for stars brighter than $K_p \approx 13$, and within a factor of two of the Kepler precision for fainter targets. We perform cross-validation and transit injection and recovery tests to validate the pipeline, and compare our light curves to the other de-trended light curves available for download at the MAST High Level Science Products archive. We find that EVEREST achieves the highest average precision of any of these pipelines for unsaturated K2 stars. The improved precision of these light curves will aid in exoplanet detection and characterization, investigations of stellar variability, asteroseismology, and other photometric studies. The EVEREST pipeline can also easily be applied to future surveys, such as the TESS mission, to correct for instrumental systematics and enable the detection of low signal-to-noise transiting exoplanets. The EVEREST light curves and the source code used to generate them are freely available online.

Astronomical observations are affected by several kinds of noise, each with its own causal source; there is photon noise, stochastic source variability, and residuals coming from imperfect calibration of the detector or telescope. The precision of NASA Kepler photometry for exoplanet science---the most precise photometric measurements of stars ever made---appears to be limited by unknown or untracked variations in spacecraft pointing and temperature, and unmodeled stellar variability. Here we present the Causal Pixel Model (CPM) for Kepler data, a data-driven model intended to capture variability but preserve transit signals. The CPM works at the pixel level so that it can capture very fine-grained information about the variation of the spacecraft. The CPM predicts each target pixel value from a large number of pixels of other stars sharing the instrument variabilities while not containing any information on possible transits in the target star. In addition, we use the target star's future and past (auto-regression). By appropriately separating, for each data point, the data into training and test sets, we ensure that information about any transit will be perfectly isolated from the model. The method has four hyper-parameters (the number of predictor stars, the auto-regressive window size, and two L2-regularization amplitudes for model components), which we set by cross-validation. We determine a generic set of hyper-parameters that works well for most of the stars and apply the method to a corresponding set of target stars. We find that we can consistently outperform (for the purposes of exoplanet detection) the Kepler Pre-search Data Conditioning (PDC) method for exoplanet discovery.

$K2$'s Campaign 9 ($K2$C9) will conduct a $\sim$3.7 deg$^{2}$ survey toward the Galactic bulge from 7/April through 1/July of 2016 that will leverage the spatial separation between $K2$ and the Earth to facilitate measurement of the microlens parallax $\pi_{\rm E}$ for $\gtrsim$127 microlensing events. These will include several that are planetary in nature as well as many short-timescale microlensing events, which are potentially indicative of free-floating planets (FFPs). These satellite parallax measurements will in turn allow for the direct measurement of the masses of and distances to the lensing systems. In this white paper we provide an overview of the $K2$C9 space- and ground-based microlensing survey. Specifically, we detail the demographic questions that can be addressed by this program, including the frequency of FFPs and the Galactic distribution of exoplanets, the observational parameters of $K2$C9, and the array of resources dedicated to concurrent observations. Finally, we outline the avenues through which the larger community can become involved, and generally encourage participation in $K2$C9, which constitutes an important pathfinding mission and community exercise in anticipation of $WFIRST$.

The Second Workshop on Extreme Precision Radial Velocities defined circa 2015 the state of the art Doppler precision and identified the critical path challenges for reaching 10 cm/s measurement precision. The presentations and discussion of key issues for instrumentation and data analysis and the workshop recommendations for achieving this precision are summarized here. Beginning with the HARPS spectrograph, technological advances for precision radial velocity measurements have focused on building extremely stable instruments. To reach still higher precision, future spectrometers will need to produce even higher fidelity spectra. This should be possible with improved environmental control, greater stability in the illumination of the spectrometer optics, better detectors, more precise wavelength calibration, and broader bandwidth spectra. Key data analysis challenges for the precision radial velocity community include distinguishing center of mass Keplerian motion from photospheric velocities, and the proper treatment of telluric contamination. Success here is coupled to the instrument design, but also requires the implementation of robust statistical and modeling techniques. Center of mass velocities produce Doppler shifts that affect every line identically, while photospheric velocities produce line profile asymmetries with wavelength and temporal dependencies that are different from Keplerian signals. Exoplanets are an important subfield of astronomy and there has been an impressive rate of discovery over the past two decades. Higher precision radial velocity measurements are required to serve as a discovery technique for potentially habitable worlds and to characterize detections from transit missions. The future of exoplanet science has very different trajectories depending on the precision that can ultimately be achieved with Doppler measurements.

The extended Kepler mission, K2, is now providing photometry of new fields every three months in a search for transiting planets. In a recent study, Foreman-Mackey and collaborators presented a list of 36 planet candidates orbiting 31 stars in K2 Campaign 1. In this contribution, we present stellar and planetary properties for all systems. We combine ground-based seeing-limited survey data and adaptive optics imaging with an automated transit analysis scheme to validate 21 candidates as planets, 17 for the first time, and identify 6 candidates as likely false positives. Of particular interest is K2-18 (EPIC 201912552), a bright (K=8.9) M2.8 dwarf hosting a 2.23 \pm 0.25 R_Earth planet with T_eq = 272 \pm 15 K and an orbital period of 33 days. We also present two new open-source software packages which enable this analysis. The first, isochrones, is a flexible tool for fitting theoretical stellar models to observational data to determine stellar properties using a nested sampling scheme to capture the multimodal nature of the posterior distributions of the physical parameters of stars that may plausibly be evolved. The second is vespa, a new general-purpose procedure to calculate false positive probabilities and statistically validate transiting exoplanets.

Photometry of stars from the K2 extension of NASA's Kepler mission is afflicted by systematic effects caused by small (few-pixel) drifts in the telescope pointing and other spacecraft issues. We present a method for searching K2 light curves for evidence of exoplanets by simultaneously fitting for these systematics and the transit signals of interest. This method is more computationally expensive than standard search algorithms but we demonstrate that it can be efficiently implemented and used to discover transit signals. We apply this method to the full Campaign 1 dataset and report a list of 36 planet candidates transiting 31 stars, along with an analysis of the pipeline performance and detection efficiency based on artificial signal injections and recoveries. For all planet candidates, we present posterior distributions on the properties of each system based strictly on the transit observables.

From pulsating stars to transiting exoplanets, the search for periodic signals in K2 data, Kepler's 2-wheeled extension, is relevant to a long list of scientific goals. Systematics affecting K2 light curves due to the decreased spacecraft pointing precision inhibit the easy extraction of periodic signals from the data. We here develop a method for producing periodograms of K2 light curves that are insensitive to pointing-induced systematics; the Systematics-Insensitive Periodogram (SIP). Traditional sine-fitting periodograms use a generative model to find the frequency of a sinusoid that best describes the data. We extend this principle by including systematic trends, based on a set of 'Eigen light curves', following Foreman-Mackey et al. (2015), in our generative model as well as a sum of sine and cosine functions over a grid of frequencies. Using this method we are able to produce periodograms with vastly reduced systematic features. The quality of the resulting periodograms are such that we can recover acoustic oscillations in giant stars and measure stellar rotation periods without the need for any detrending. The algorithm is also applicable to the detection of other periodic phenomena such as variable stars, eclipsing binaries and short-period exoplanet candidates. The SIP code is available at https://github.com/RuthAngus/SIPK2.

We describe a method for removing the effect of confounders in order to reconstruct a latent quantity of interest. The method, referred to as half-sibling regression, is inspired by recent work in causal inference using additive noise models. We provide a theoretical justification and illustrate the potential of the method in a challenging astronomy application.

Kepler-296 is a binary star system with two M-dwarf components separated by 0.2 arcsec. Five transiting planets have been confirmed to be associated with the Kepler-296 system; given the evidence to date, however, the planets could in principle orbit either star. This ambiguity has made it difficult to constrain both the orbital and physical properties of the planets. Using both statistical and analytical arguments, this paper shows that all five planets are highly likely to orbit the primary star in this system. We performed a Markov-Chain Monte Carlo simulation using a five transiting planet model, leaving the stellar density and dilution with uniform priors. Using importance sampling, we compared the model probabilities under the priors of the planets orbiting either the brighter or the fainter component of the binary. A model where the planets orbit the brighter component, Kepler-296A, is strongly preferred by the data. Combined with our assertion that all five planets orbit the same star, the two outer planets in the system, Kepler-296 Ae and Kepler-296 Af, have radii of 1.53 +/- 0.26 and 1.80 +/- 0.31 R_earth, respectively, and receive incident stellar fluxes of 1.40 +/- 0.23 and 0.62 +/- 0.10 times the incident flux the Earth receives from the Sun. This level of irradiation places both planets within or close to the circumstellar habitable zone of their parent star.

A number of problems in probability and statistics can be addressed using the multivariate normal (Gaussian) distribution. In the one-dimensional case, computing the probability for a given mean and variance simply requires the evaluation of the corresponding Gaussian density. In the $n$-dimensional setting, however, it requires the inversion of an $n \times n$ covariance matrix, $C$, as well as the evaluation of its determinant, $\det(C)$. In many cases, such as regression using Gaussian processes, the covariance matrix is of the form $C = \sigma^2 I + K$, where $K$ is computed using a specified covariance kernel which depends on the data and additional parameters (hyperparameters). The matrix $C$ is typically dense, causing standard direct methods for inversion and determinant evaluation to require $\mathcal O(n^3)$ work. This cost is prohibitive for large-scale modeling. Here, we show that for the most commonly used covariance functions, the matrix $C$ can be hierarchically factored into a product of block low-rank updates of the identity matrix, yielding an $\mathcal O (n\log^2 n) $ algorithm for inversion. More importantly, we show that this factorization enables the evaluation of the determinant $\det(C)$, permitting the direct calculation of probabilities in high dimensions under fairly broad assumptions on the kernel defining $K$. Our fast algorithm brings many problems in marginalization and the adaptation of hyperparameters within practical reach using a single CPU core. The combination of nearly optimal scaling in terms of problem size with high-performance computing resources will permit the modeling of previously intractable problems. We illustrate the performance of the scheme on standard covariance kernels.

Among the available methods for dating stars, gyrochronology is a powerful one because it requires knowledge of only the star's mass and rotation period. Gyrochronology relations have previously been calibrated using young clusters, with the Sun providing the only age dependence, and are therefore poorly calibrated at late ages. We used rotation period measurements of 310 Kepler stars with asteroseismic ages, 50 stars from the Hyades and Coma Berenices clusters and 6 field stars (including the Sun) with precise age measurements to calibrate the gyrochronology relation, whilst fully accounting for measurement uncertainties in all observable quantities. We calibrated a relation of the form $P=A^n\times(B-V-c)^b$, where $P$ is rotation period in days, $A$ is age in Myr, $B$ and $V$ are magnitudes and $a$, $b$ and $n$ are the free parameters of our model. We found $a = 0.40^{+0.3}_{-0.05}$, $b = 0.31^{+0.05}_{-0.02}$ and $n = 0.55^{+0.02}_{-0.09}$. Markov Chain Monte Carlo methods were used to explore the posterior probability distribution functions of the gyrochronology parameters and we carefully checked the effects of leaving out parts of our sample, leading us to find that no single relation beween rotation period, colour and age can adequately describe all the subsets of our data. The Kepler asteroseismic stars, cluster stars and local field stars cannot all be described by the same gyrochronology relation. The Kepler asteroseismic stars may be subject to observational biases, however the clusters show unexpected deviations from the predicted behaviour, providing concerns for the overall reliability of gyrochronology as a dating method.

We have undertaken the largest systematic study of the high-mass stellar initial mass function (IMF) to date using the optical color-magnitude diagrams (CMDs) of 85 resolved, young (4 Myr < t < 25 Myr), intermediate mass star clusters (10^3-10^4 Msun), observed as part of the Panchromatic Hubble Andromeda Treasury (PHAT) program. We fit each cluster's CMD to measure its mass function (MF) slope for stars >2 Msun. For the ensemble of clusters, the distribution of stellar MF slopes is best described by $\Gamma=+1.45^{+0.03}_{-0.06}$ with a very small intrinsic scatter. The data also imply no significant dependencies of the MF slope on cluster age, mass, and size, providing direct observational evidence that the measured MF represents the IMF. This analysis implies that the high-mass IMF slope in M31 clusters is universal with a slope ($\Gamma=+1.45^{+0.03}_{-0.06}$) that is steeper than the canonical Kroupa (+1.30) and Salpeter (+1.35) values. Using our inference model on select Milky Way (MW) and LMC high-mass IMF studies from the literature, we find $\Gamma_{\rm MW} \sim+1.15\pm0.1$ and $\Gamma_{\rm LMC} \sim+1.3\pm0.1$, both with intrinsic scatter of ~0.3-0.4 dex. Thus, while the high-mass IMF in the Local Group may be universal, systematics in literature IMF studies preclude any definitive conclusions; homogenous investigations of the high-mass IMF in the local universe are needed to overcome this limitation. Consequently, the present study represents the most robust measurement of the high-mass IMF slope to date. We have grafted the M31 high-mass IMF slope onto widely used sub-solar mass Kroupa and Chabrier IMFs and show that commonly used UV- and Halpha-based star formation rates should be increased by a factor of ~1.3-1.5 and the number of stars with masses >8 Msun are ~25% fewer than expected for a Salpeter/Kroupa IMF. [abridged]

Kepler-91b is a rare example of a transiting hot Jupiter around a red giant star, providing the possibility to study the formation and composition of hot Jupiters under different conditions compared to main-sequence stars. However, the planetary nature of Kepler-91b, which was confirmed using phase-curve variations by Lillo-Box et al., was recently called into question based on a re-analysis of Kepler data. We have obtained ground-based radial velocity observations from the Hobby-Eberly Telescope and unambiguously confirm the planetary nature of Kepler-91b by simultaneously modeling the Kepler and radial velocity data. The star exhibits temporally correlated noise due to stellar granulation which we model as a Gaussian Process. We hypothesize that it is this noise component that led previous studies to suspect Kepler-91b to be a false positive. Our work confirms the conclusions presented by Lillo-Box et al. that Kepler-91b is a 0.73+/-0.13 Mjup planet orbiting a red giant star.

No true extrasolar Earth analog is known. Hundreds of planets have been found around Sun-like stars that are either Earth-sized but on shorter periods, or else on year-long orbits but somewhat larger. Under strong assumptions, exoplanet catalogs have been used to make an extrapolated estimate of the rate at which Sun-like stars host Earth analogs. These studies are complicated by the fact that every catalog is censored by non-trivial selection effects and detection efficiencies, and every property (period, radius, etc.) is measured noisily. Here we present a general hierarchical probabilistic framework for making justified inferences about the population of exoplanets, taking into account survey completeness and, for the first time, observational uncertainties. We are able to make fewer assumptions about the distribution than previous studies; we only require that the occurrence rate density be a smooth function of period and radius (employing a Gaussian process). By applying our method to synthetic catalogs, we demonstrate that it produces more accurate estimates of the whole population than standard procedures based on weighting by inverse detection efficiency. We apply the method to an existing catalog of small planet candidates around G dwarf stars (Petigura et al. 2013). We confirm a previous result that the radius distribution changes slope near Earth's radius. We find that the rate density of Earth analogs is about 0.02 (per star per natural logarithmic bin in period and radius) with large uncertainty. This number is much smaller than previous estimates made with the same data but stronger assumptions.

We establish the three-dimensional architecture of the Kepler-419 (previously KOI-1474) system to be eccentric yet with a low mutual inclination. Kepler-419b is a warm Jupiter at semi-major axis a = 0.370 +0.007/-0.006 AU with a large eccentricity e=0.85 +0.08/-0.07 measured via the "photoeccentric effect." It exhibits transit timing variations induced by the non-transiting Kepler-419c, which we uniquely constrain to be a moderately eccentric (e=0.184 +/- 0.002), hierarchically-separated (a=1.68 +/- 0.03 AU) giant planet (7.3 +/- 0.4 MJup). We combine sixteen quarters of Kepler photometry, radial-velocity (RV) measurements from the HIgh Resolution Echelle Spectrometer (HIRES) on Keck, and improved stellar parameters that we derive from spectroscopy and asteroseismology. From the RVs, we measure the mass of inner planet to be 2.5+/-0.3MJup and confirm its photometrically-measured eccentricity, refining the value to e=0.83+/-0.01. The RV acceleration is consistent with the properties of the outer planet derived from TTVs. We find that, despite their sizable eccentricities, the planets are coplanar to within 9 +8/-6 degrees, and therefore the inner planet's large eccentricity and close-in orbit are unlikely to be the result of Kozai migration. Moreover, even over many secular cycles, the inner planet's periapse is most likely never small enough for tidal circularization. Finally, we present and measure a transit time and impact parameter from four simultaneous ground-based light curves from 1m-class telescopes, demonstrating the feasibility of ground-based follow-up of Kepler giant planets exhibiting large TTVs.

We introduce a stable, well tested Python implementation of the affine-invariant ensemble sampler for Markov chain Monte Carlo (MCMC) proposed by Goodman & Weare (2010). The code is open source and has already been used in several published projects in the astrophysics literature. The algorithm behind emcee has several advantages over traditional MCMC sampling methods and it has excellent performance as measured by the autocorrelation time (or function calls per independent sample). One major advantage of the algorithm is that it requires hand-tuning of only 1 or 2 parameters compared to $\sim N^2$ for a traditional algorithm in an N-dimensional parameter space. In this document, we describe the algorithm and the details of our implementation and API. Exploiting the parallelism of the ensemble method, emcee permits any user to take advantage of multiple CPU cores without extra effort. The code is available online at http://dan.iel.fm/emcee under the MIT License.

We characterize the bulge, disk, and halo subcomponents in the Andromeda galaxy (M31) over the radial range 4 < R_proj < 225 kpc. The cospatial nature of these subcomponents renders them difficult to disentangle using surface brightness (SB) information alone, especially interior to ~20 kpc. Our new decomposition technique combines information from the luminosity function (LF) of over 1.5 million bright (20 < m_814W < 22) stars from the Panchromatic Hubble Andromeda Treasury (PHAT) survey, radial velocities of over 5000 red giant branch stars in the same magnitude range from the Spectroscopic and Photometric Landscape of Andromeda's Stellar Halo (SPLASH) survey, and integrated I-band SB profiles from various sources. We use an affine-invariant Markov chain Monte Carlo algorithm to fit an appropriate toy model to these three data sets. The bulge, disk, and halo SB profiles are modeled as a Sersic, exponential, and cored power-law, respectively, and the LFs are modeled as broken power-laws. We present probability distributions for each of 32 parameters describing the SB profiles and LFs of the three subcomponents. We find that the number of stars with a disk-like LF is ~5% larger than the the number with disk-like (dynamically cold) kinematics, suggesting that some stars born in the disk have been dynamically heated to the point that they are kinematically indistinguishable from halo members. This is the first kinematical evidence for a "kicked-up disk" halo population in M31. The fraction of kicked-up disk stars is consistent with that found in simulations. We also find evidence for a radially varying disk LF, consistent with a negative metallicity gradient in the stellar disk.

Kepler's immense photometric precision to date was maintained through satellite stability and precise pointing. In this white paper, we argue that image modeling--fitting the Kepler-downlinked raw pixel data--can vastly improve the precision of Kepler in pointing-degraded two-wheel mode. We argue that a non-trivial modeling effort may permit continuance of photometry at 10-ppm-level precision. We demonstrate some baby steps towards precise models in both data-driven (flexible) and physics-driven (interpretably parameterized) modes. We demonstrate that the expected drift or jitter in positions in the two-weel era will help with constraining calibration parameters. In particular, we show that we can infer the device flat-field at higher than pixel resolution; that is, we can infer pixel-to-pixel variations in intra-pixel sensitivity. These results are relevant to almost any scientific goal for the repurposed mission; image modeling ought to be a part of any two-wheel repurpose for the satellite. We make other recommendations for Kepler operations, but fundamentally advocate that the project stick with its core mission of finding and characterizing Earth analogs. [abridged]

In today's mailing, Hogg et al. propose image modeling techniques to maintain 10-ppm-level precision photometry in Kepler data with only two working reaction wheels. While these results are relevant to many scientific goals for the repurposed mission, all modeling efforts so far have used a toy model of the Kepler telescope. Because the two-wheel performance of Kepler remains to be determined, we advocate for the consideration of an alternate strategy for a >1 year program that maximizes the science return from the "low-torque" fields across the ecliptic plane. Assuming we can reach the precision of the original Kepler mission, we expect to detect 800 new planet candidates in the first year of such a mission. Our proposed strategy has benefits for transit timing variation and transit duration variation studies, especially when considered in concert with the future TESS mission. We also expect to help address the first key science goal of Kepler: the frequency of planets in the habitable zone as a function of spectral type.

We present and implement a probabilistic (Bayesian) method for producing catalogs from images of stellar fields. The method is capable of inferring the number of sources N in the image and can also handle the challenges introduced by noise, overlapping sources, and an unknown point spread function (PSF). The luminosity function of the stars can also be inferred even when the precise luminosity of each star is uncertain, via the use of a hierarchical Bayesian model. The computational feasibility of the method is demonstrated on two simulated images with different numbers of stars. We find that our method successfully recovers the input parameter values along with principled uncertainties even when the field is crowded. We also compare our results with those obtained from the SExtractor software. While the two approaches largely agree about the fluxes of the bright stars, the Bayesian approach provides more accurate inferences about the faint stars and the number of stars, particularly in the crowded case.

We present a probabilistic approach for inferring the parameters of the present day power-law stellar mass function (MF) of a resolved young star cluster. This technique (a) fully exploits the information content of a given dataset; (b) accounts for observational uncertainties in a straightforward way; (c) assigns meaningful uncertainties to the inferred parameters; (d) avoids the pitfalls associated with binning data; and (e) is applicable to virtually any resolved young cluster, laying the groundwork for a systematic study of the high mass stellar MF (M > 1 Msun). Using simulated clusters and Markov chain Monte Carlo sampling of the probability distribution functions, we show that estimates of the MF slope, {\alpha}, are unbiased and that the uncertainty, {\Delta}{\alpha}, depends primarily on the number of observed stars and stellar mass range they span, assuming that the uncertainties on individual masses and the completeness are well-characterized. Using idealized mock data, we compute the lower limit precision on {\alpha} and provide an analytic approximation for {\Delta}{\alpha} as a function of the observed number of stars and mass range. We find that ~ 3/4 of quoted literature uncertainties are smaller than the theoretical lower limit. By correcting these uncertainties to the theoretical lower limits, we find the literature studies yield <{\alpha}>=2.46 with a 1-{\sigma} dispersion of 0.35 dex. We verify that it is impossible for a power-law MF to obtain meaningful constraints on the upper mass limit of the IMF. We show that avoiding substantial biases in the MF slope requires: (1) including the MF as a prior when deriving individual stellar mass estimates; (2) modeling the uncertainties in the individual stellar masses; and (3) fully characterizing and then explicitly modeling the completeness for stars of a given mass. (abridged)