This paper documents the sixteenth data release (DR16) from the Sloan Digital Sky Surveys; the fourth and penultimate from the fourth phase (SDSS-IV). This is the first release of data from the southern hemisphere survey of the Apache Point Observatory Galactic Evolution Experiment 2 (APOGEE-2); new data from APOGEE-2 North are also included. DR16 is also notable as the final data release for the main cosmological program of the Extended Baryon Oscillation Spectroscopic Survey (eBOSS), and all raw and reduced spectra from that project are released here. DR16 also includes all the data from the Time Domain Spectroscopic Survey (TDSS) and new data from the SPectroscopic IDentification of ERosita Survey (SPIDERS) programs, both of which were co-observed on eBOSS plates. DR16 has no new data from the Mapping Nearby Galaxies at Apache Point Observatory (MaNGA) survey (or the MaNGA Stellar Library "MaStar"). We also preview future SDSS-V operations (due to start in 2020), and summarize plans for the final SDSS-IV data release (DR17).

We make use of APOGEE and $Gaia$ data to identify stars that are consistent with being born in the same association or star cluster as the Sun. We limit our analysis to stars that match solar abundances within their uncertainties, as they could have formed from the same Giant Molecular Cloud (GMC) as the Sun. We constrain the range of orbital actions that solar siblings can have with a suite of simulations of solar birth clusters evolved in static and time-dependent tidal fields. The static components of each galaxy model are the bulge, disk, and halo, while the various time-dependent components include a bar, spiral arms, and GMCs. In galaxy models without GMCs, simulated solar siblings all have $J_R < 122$ km $\rm s^{-1}$ kpc, $990 < L_z < 1986$ km $\rm s^{-1}$ kpc, and $0.15 < J_z < 0.58$ km $\rm s^{-1}$ kpc. Given the actions of stars in APOGEE and $Gaia$, we find 104 stars that fall within this range. One candidate in particular, Solar Sibling 1, has both chemistry and actions similar enough to the solar values that strong interactions with the bar or spiral arms are not required for it to be dynamically associated with the Sun. Adding GMCs to the potential can eject solar siblings out of the plane of the disk and increase their $J_z$, resulting in a final candidate list of 296 stars. The entire suite of simulations indicate that solar siblings should have $J_R < 122$ km $\rm s^{-1}$ kpc, $353 < L_z < 2110$ km $\rm s^{-1}$ kpc, and $J_z < 0.8$ km $\rm s^{-1}$ kpc. Given these criteria, it is most likely that the association or cluster that the Sun was born in has reached dissolution and is not the commonly cited open cluster M67.

Studies of the ages, abundances, and motions of individual stars in the Milky Way provide one of the best ways to study the evolution of disk galaxies over cosmic time. The formation of the Milky Way's barred inner region in particular is a crucial piece of the puzzle of disk galaxy evolution. Using data from APOGEE and Gaia, we present maps of the kinematics, elemental abundances, and age of the Milky Way bulge and disk that show the barred structure of the inner Milky Way in unprecedented detail. The kinematic maps allow a direct, purely kinematic determination of the bar's pattern speed of 41+/-3 km/s/kpc and of its shape and radial profile. We find the bar's age, metallicity, and abundance ratios to be the same as those of the oldest stars in the disk that are formed in its turbulent beginnings, while stars in the bulge outside of the bar are younger and more metal-rich. This implies that the bar likely formed ~8 Gyr ago, when the decrease in turbulence in the gas disk allowed a thin disk to form that quickly became bar-unstable. The bar's formation therefore stands as a crucial epoch in the evolution of the Milky Way, a picture that is in line with the evolutionary path that emerges from observations of the gas kinematics in external disk galaxies over the last ~10 Gyr.

The kinematics of the Milky Way disc as a function of age are well measured at the solar radius, but have not been studied over a wider range of Galactocentric radii. Here, we measure the kinematics of mono-age, mono-$\mathrm{[Fe/H]}$ populations in the low and high $\mathrm{[\alpha/Fe]}$ discs between $4 \lesssim R \lesssim 13$ kpc and $|z| \lesssim 2$ kpc using 65,719 stars in common between APOGEE DR14 and $\it{Gaia}$ DR2 for which we estimate ages using a Bayesian neural network model trained on asteroseismic ages. We determine the vertical and radial velocity dispersions, finding that the low and high $\mathrm{[\alpha/Fe]}$ discs display markedly different age--velocity-dispersion relations (AVRs) and shapes $\sigma_z/\sigma_R$. The high $\mathrm{[\alpha/Fe]}$ disc has roughly flat AVRs and constant $\sigma_z/\sigma_R = 0.64\pm 0.04$, whereas the low $\mathrm{[\alpha/Fe]}$ disc has large variations in this ratio which positively correlate with the mean orbital radius of the population at fixed age. The high $\mathrm{[\alpha/Fe]}$ disc component's flat AVRs and constant $\sigma_z/\sigma_R$ clearly indicates an entirely different heating history. Outer disc populations also have flatter radial AVRs than those in the inner disc, likely due to the waning effect of spiral arms. Our detailed measurements of AVRs and $\sigma_z/\sigma_R$ across the disc indicate that low $\mathrm{[\alpha/Fe]}$, inner disc ($R \lesssim 10\,\mathrm{kpc}$) stellar populations are likely dynamically heated by both giant molecular clouds and spiral arms, while the observed trends for outer disc populations require a significant contribution from another heating mechanism such as satellite perturbations. We also find that outer disc populations have slightly positive mean vertical and radial velocities, likely because they are part of the warped disc.

Gaia measures the five astrometric parameters for stars in the Milky Way, but only four of them (positions and proper motion, but not parallax) are well measured beyond a few kpc from the Sun. Modern spectroscopic surveys such as APOGEE cover a large area of the Milky Way disk and we can use the relation between spectra and luminosity to determine distances to stars beyond Gaia's parallax reach. Here, we design a deep neural network trained on stars in common between Gaia and APOGEE that determines spectro-photometric distances to APOGEE stars, while including a flexible model to calibrate parallax zero-point biases in Gaia DR2. We determine the zero-point offset to be $-52.3 \pm 2.0uas$ when modeling it as a global constant, but also train a multivariate zero-point offset model that depends on $G$, $G_{BP} - G_{RP}$ color, and $T_\mathrm{eff}$ and that can be applied to all 139 million stars in Gaia DR2 within APOGEE's color--magnitude range. Our spectro-photometric distances are more precise than Gaia at distances $\approx 2kpc$ from the Sun. We release a catalog of spectro-photometric distances for the entire APOGEE DR14 data set which covers Galactocentric radii $2kpc\lesssim R \lesssim19kpc$; $\approx 150,000$ stars have <10% uncertainty, making this a powerful sample to study the chemo-dynamical structure of the disk. We use this sample to map the mean [Fe/H] and 15 abundance ratios [X/Fe] from the Galactic center to the edge of the disk. Among many interesting trends, we find that the bulge and bar region at $R \lesssim 5kpc$ clearly stands out in [Fe/H] and most abundance ratios.

Deep learning with artificial neural networks is increasingly gaining attention, because of its potential for data-driven astronomy. However, this methodology usually does not provide uncertainties and does not deal with incompleteness and noise in the training data. In this work, we design a neural network for high-resolution spectroscopic analysis using APOGEE data that mimics the methodology of standard spectroscopic analyses: stellar parameters are determined using the full wavelength range, but individual element abundances use censored portions of the spectrum. We train this network with a customized objective function that deals with incomplete and noisy training data and apply dropout variational inference to derive uncertainties on our predictions. We determine parameters and abundances for 18 individual elements at the $\approx 0.03$ dex level, even at low signal-to-noise ratio. We demonstrate that the uncertainties returned by our method are a realistic estimate of the precision and they automatically blow up when inputs or outputs outside of the training set are encountered, thus shielding users from unwanted extrapolation. By using standard deep-learning tools for GPU acceleration, our method is extremely fast, allowing analysis of the entire APOGEE data set of $\approx250,000$ spectra in ten minutes on a single, low-cost GPU. We release the stellar parameters and 18 individual-element abundances with associated uncertainty for the entire APOGEE DR14 dataset. Simultaneously, we release astroNN, a well-tested, open-source python package developed for this work, but that is also designed to be a general package for deep learning in astronomy. astroNN is available at https://github.com/henrysky/astroNN with extensive documentation at http://astroNN.readthedocs.io.

We have used AAVSO visual and photoelectric V data, and the AAVSO time-series package VSTAR and the Lomb-Scargle time-series algorithm to determine improved pulsation periods, "long secondary periods" (LSPs). and their amplitudes in 51 shorter-period pulsating red giants in the AAVSO photoelectric photometry program, and the AAVSO LPV (long period variable) binocular program. As is well known, radial pulsation becomes detectable in red giants at about spectral type M0, with periods of about 20 days. We find that the LSP phenomenon is also first detectable at about M0. Pulsation and LSP amplitudes increase from near zero to about 0.1 at pulsation periods of about 100 days. At longer periods, the pulsation amplitudes continue to increase, but the LSP amplitudes are generally between 0.1 and 0.2 on average. The ratios of LSP to pulsation period cluster around 5 and 10, presumably depending on whether the pulsation period is the fundamental or first overtone. The pulsation and LSP phase curves are generally close to sinusoidal, except when the amplitude is small, in which case they may be distorted by observational scatter or, in the case of the LSP amplitude, by the pulsational variability. As with longer-period stars, the LSP amplitude increases and decreases by a factor of two or more, for unknown reasons, on a median time scale of about 20 LSPs. The LSP phenomenon is thus present and simiklar in radially pulsating red giants of all periods. Its cause remains unknown.