We investigate whether the galaxy and star formation model used for the Auriga simulations can produce a realistic globular cluster (GC) population. We compare statistics of GC candidate star particles in the Auriga haloes with catalogues of the Milky Way (MW) and Andromeda (M31) GC populations. We find that the Auriga simulations do produce sufficient stellar mass for GC candidates at radii and metallicities that are typical for the MW GC system (GCS). We also find varying mass-ratios of the simulated GC candidates relative to the observed mass in the MW and M31 GC systems for different bins of galactocentric radius-metallicity (r$_{\text{gal}}$ -[Fe/H]). Overall, the Auriga simulations produce GC candidates with higher metallicities than the MW and M31 GCS and they are found at larger radii than observed. The Auriga simulations would require bound cluster formation efficiencies higher than ten percent for the metal-poor GC candidates, and those within the Solar radius should experience negligible destruction rates to be consistent with observations. GC candidates in the outer halo, on the other hand, should either have low formation efficiencies, or experience high mass loss for the Auriga simulations to produce a GCS that is consistent with that of the MW or M31. Finally, the scatter in the metallicity as well as in the radial distribution between different Auriga runs is considerably smaller than the differences between that of the MW and M31 GCSs. The Auriga model is unlikely to give rise to a GCS that can be consistent with both galaxies.

Stars born at the same time in the same place should have formed from gas of the same element composition. But most stars subsequently disperse from their birth siblings, in orbit and orbital phase, becoming 'field stars'. Here we explore and provide direct observational evidence for this process in the Milky Way disc, by quantifying the probability that orbit-similarity among stars implies indistinguishable metallicity. We define the orbit similarity among stars through their distance in action-angle space, $\Delta (J,\theta)$, and their abundance similarity simply by $\Delta$[Fe/H]. Analyzing a sample of main sequence stars from Gaia DR2 and LAMOST, we find an excess of pairs with the same metallicity ($\Delta\mathrm{[Fe/H]}<0.1$) that extends to remarkably large separations in $\Delta (J,\theta)$ that correspond to nearly 1 kpc distances. We assess the significance of this effect through a mock sample, drawn from a smooth and phase-mixed orbit distribution. Through grouping such star pairs into associations with a friend-of-friends algorithm linked by $\Delta (J,\theta)$, we find 100s of mono-abundance groups with $\ge 3$ (to $\gtrsim 20$) members; these groups -- some clusters, some spread across the sky -- are over an order-of-magnitude more abundant than expected for a smooth phase-space distribution, suggesting that we are witnessing the 'dissolution' of stellar birth associations into the field.

Many astrophysical and galaxy-scale cosmological problems require a well determined gravitational potential which is often modeled by observers under strong assumptions. Globular clusters (GCs) surrounding galaxies can be used as dynamical tracers of the luminous and dark matter distribution at large (kpc) scales. A natural assumption for modeling the gravitational potential is that GCs accreted in the same dwarf galaxy merger event move at the present time on similar orbits in the host galaxy and should therefore have similar actions. We investigate this idea in one realistic Milky Way like galaxy of the cosmological N-body simulation suite Auriga. We show how the actions of accreted stellar particles in the simulation evolve and that minimizing the standard deviation of GCs in action space, however, cannot constrain the true potential. This approach known as adaptive dynamics does therefore not work for accreted GCs.

We explore the origin of phase-space substructures revealed by the second Gaia data release in the disc of the Milky Way, such as the ridges in the $V_{\phi}$-$r$ plane, the undulations in the $V_{\phi}$-$r$-$V_r$ space and the streams in the $V_{\phi}$-$V_r$ plane. We use a collisionless N-body simulation with co-spatial thin and thick discs, along with orbit integration, to study the orbital structure close to the Outer Lindblad Resonance (OLR) of the bar. We find that a prominent, long-lived ridge is formed in the $V_{\phi}$-$r$ plane due to the OLR which translates to streams in the $V_{\phi}$-$V_r$ plane and examine which closed periodic and trapped librating orbits are responsible for these features. We find that orbits which carry out small librations around the $x_1(1)$ family are preferentially found at negative $V_r$, giving rise to a `horn'-like feature, while orbits with larger libration amplitudes, trapped around the $x_1(2)$ and $x_1(1)$ families, constitute the positive $V_r$ substructure, i.e. the Hercules-like feature. This changing libration amplitude of orbits will translate to a changing ratio of thin/thick disc stars, which could have implications on the metallicity distribution in this plane. We find that a scenario in which the Sun is placed close to the OLR gives rise to a strong asymmetry in $V_r$ in the $V_{\phi}$-$V_r$ plane (i.e. Hercules vs. `the horn') and subsequently to undulations in the $V_{\phi}$-$r$-$V_r$ space. We also explore a scenario in which the Sun is placed closer to the bar corotation and find that the bar perturbation $alone$ cannot give rise to the these features.

Action space synthesizes the orbital information of stars and is well-suited to analyse the rich amount of kinematic disk substructure in the Gaia DR2 radial velocity sample (RVS). In this work, we revisit one of the strongest perturbers in the Milky Way (MW) disk: the m=2 bar. We investigate how its resonances affect the actions of individual test particle stars, i.e., (JR,Lz,Jz) estimated in an axisymmetric MW potential. We confirm that the stars' behaviour is well approximated by scattering and oscillation along a slope Delta JR / Delta Lz = l/m centered on the l:m resonance lines. The Outer Lindblad Resonance (OLR, l=+1,m=2) creates signatures in the stellar action space that can be used to identify the Galactic bar's OLR in the Gaia DR2 RVS data: (a) The JR dependence of the oscillation causes an overdensity ridge (underdensity region) at Lz larger (smaller) than the resonance line in the (Lz,JR) plane. (b) For the first time, we demonstrate that the OLR is expected to cause a gradient in average Jz with Lz across the resonance. (c) We show that the change of predominantly outward to inward motions at the OLR occurs along the resonance line in action space. The latter signature allows us to identify three candidates for the bar's OLR - and therefore its pattern speed Omega_bar - in the Gaia data within 3 kpc from the Sun: 1.85 Omega0, 1.2 Omega0, and 1.6 Omega0 (with ~0.1 Omega0 uncertainty). This demonstrates that (i) the local Gaia action data is consistent with both the short-fast and long-slow bar models in the literature, and that (ii) axisymmetrically estimated actions are a powerful diagnostic even in non-axisymmetric systems.

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 DR2 has provided an unprecedented wealth of information about the kinematics of stars in the Solar neighbourhood, and has highlighted the degree of features in the Galactic disc. We confront the data with a range of bar and spiral models in both action-angle space, and the $R_{\mathrm{G}}-v_{\phi}$ plane. We find that the phase mixing induced by transient spiral structure creates ridges and arches in the local kinematics which are consistent with the Gaia data. We are able to produce a qualitatively good match to the data when combined with a bar with a variety of pattern speeds, and show that it is non trivial to decouple the effects of the bar and the spiral structure.

The quality and quantity of 6D stellar position-velocity measurements in the second Gaia data release (DR2) allows to study small-scale structure in the orbit distribution of the Galactic disc beyond the immediate Solar neighborhood. We investigate the distribution of orbital actions $(J_R,J_\phi=L_z,J_z)$ of $\sim 3.5$ million stars within $1.5~\text{kpc}$ of the Sun, for which precise actions can be calculated from Gaia DR2 alone. This distribution $n(J_R,L_z)$ reveals a remarkable amount of sub-structure. The known moving groups in the $(U,V)$-plane of the Solar neighborhood correspond to overdensities in $(J_R,L_z)$, as expected. But $n(J_R,L_z)$ also exhibits a wealth of density clumps and ridges that extend towards higher $J_R$. These $n(J_R,L_z)$ features are most prominent among orbits that stay close to the Galactic plane and remain consistently visible out to $\sim1.5~\text{kpc}$, as opposed to the sub-structure in velocity space. Some of these $n(J_R,L_z)$ ridges resemble features expected from rapid orbit diffusion along particular ($J_R,L_z$)-directions in the presence of various resonances. Several of these $n(J_R,L_z)$ structures show a dramatic imbalance of stars moving in or out, suggesting that stars are not phase-mixed along orbits or on resonant orbits. Orbital action and angle space of stars in Gaia DR2 is therefore highly structured over kpc-scales, and appears to be very informative for modeling studies of non-axisymmetric structure and resonances in the Galactic disc.

We compare the distribution in position and velocity of nearby stars from the Gaia DR2 radial velocity sample with predictions of current theories for spirals in disc galaxies. Although the rich substructure in velocity space contains the same information, we find it more revealing to reproject the data into action-angle variables, and we describe why resonant scattering would be more readily identifiable in these variables. We compute the predicted changes to the phase space density, in multiple different projections, that would be caused by a simplified isolated spiral pattern, finding widely differing predictions from each theory. We conclude that the phase space structure present in the Gaia data shares many of the qualitative features expected in the transient spiral mode model. We argue that the popular picture of apparently swing-amplified spirals results from the superposition of a few underlying spiral modes.

Spectrophotometric distances to stars observed by large spectroscopic surveys offer a crucial complement to parallax distances that remain very important also after the future Gaia data releases. Here we present a probabilistic approach to modeling spectroscopic information for a subset of 4,000 main sequence stars with good parallaxes ($\sigma_\varpi/\varpi<0.1$) from the LAMOST $\times$ TGAS $\times$ 2MASS cross-match, yielding a precise spectroscopic distance estimator with uncertainties of $\sim$6% for single stars. Unlike previous approaches to this problem, we explicitly account for the individual parallax uncertainties in the model building and fully incorporate the fraction of near-equal binaries of main sequence stars, which would lead to biased distance estimates if neglected. Using this model, we estimate the distance for all (150,000) main sequence stars from LAMOST Data Release 5, without parallax information. As an application, we compute their orbital actions, where our more precise distances result in 5 times smaller action uncertainties. This illustrates how future studies of the Milky Way's orbital structure can benefit from using our model. For the fainter and more distant stars of most current spectroscopic surveys, an approach such as the one presented in this work will deliver better distances than Gaia Data Release 2.

Investigations of the origin and evolution of the Milky Way disk have long relied on chemical and kinematic identification of its components to reconstruct our Galactic past. Difficulties in determining precise stellar ages have restricted most studies to small samples, normally confined to the solar neighbourhood. Here we break this impasse with the help of asteroseismic inference and perform a chronology of the evolution of the disk throughout the age of the Galaxy. We chemically dissect the Milky Way disk population using a sample of red giant stars spanning out to 2~kpc in the solar annulus observed by the {\it Kepler} satellite, with the added dimension of asteroseismic ages. Our results reveal a clear difference in age between the low- and high-$\alpha$ populations, which also show distinct velocity dispersions in the $V$ and $W$ components. We find no tight correlation between age and metallicity nor [$\alpha$/Fe] for the high-$\alpha$ disk stars. Our results indicate that this component formed over a period of more than 2~Gyr with a wide range of [M/H] and [$\alpha$/Fe] independent of time. Our findings show that the kinematic properties of young $\alpha$-rich stars are consistent with the rest of the high-$\alpha$ population and different from the low-$\alpha$ stars of similar age, rendering support to their origin being old stars that went through a mass transfer or stellar merger event, making them appear younger, instead of migration of truly young stars formed close to the Galactic bar.

RoadMapping is a dynamical modeling machinery developed to constrain the Milky Way's (MW) gravitational potential by simultaneously fitting an axisymmetric parametrized potential and an action-based orbit distribution function (DF) to discrete 6D phase-space measurements of stars in the Galactic disk. In this work we demonstrate RoadMapping's robustness in the presence of spiral arms by modeling data drawn from an N-body simulation snapshot of a disk-dominated galaxy of MW mass with strong spiral arms (but no bar), exploring survey volumes with radii 500pc<=r_max<=5kpc. The potential constraints are very robust, even though we use a simple action-based DF, the quasi-isothermal DF (qDF). The best-fit RoadMapping model always recovers the correct gravitational forces where most of the stars that entered the analysis are located, even for small volumes. For data from large survey volumes, RoadMapping finds axisymmetric models that average well over the spiral arms. Unsurprisingly, the models are slightly biased by the excess of stars in the spiral arms. Gravitational potential models derived from survey volumes with at least r_max=3kpc can be reliably extrapolated to larger volumes. However, a large radial survey extent, r_max~5kpc, is needed to correctly recover the halo scale length. In general, the recovery and extrapolability of potentials inferred from data sets which were drawn from inter-arm regions appear to be better than those of data sets drawn from spiral arms. Our analysis implies that building axisymmetric models for the Galaxy with upcoming Gaia data will lead to sensible and robust approximations of the MW's potential.

We investigate the matter distribution of a spiral galaxy with a counter-rotating stellar core, SDSS J1331+3628 (J1331), independently with gravitational lensing and stellar dynamical modelling. By fitting a gravitational potential model to a quadruplet of lensing images around J1331's bulge, we tightly constrain the mass inside the Einstein radius R_ein = (0.91 +/- 0.02)'' (~= 1.83 +/- 0.04 kpc) to within 4%: M_ein = (7.8 +/- 0.3) x 10^10 M_Sun. We model observed long-slit major axis stellar kinematics in J1331's central regions by finding Multi-Gaussian Expansion (MGE) models for the stellar and dark matter distribution that solve the axisymmetric Jeans equations. The lens and dynamical model are independently derived, but in very good agreement with each other around ~R_ein. We find that J1331's center requires a steep total mass-to-light ratio gradient. A dynamical model including an NFW halo (with virial velocity v_200 ~= 240 +/- 40 km/s and concentration c_200 ~= 8 +/- 2) and moderate tangential velocity anisotropy (beta_z ~= -0.4 +/- 0.1) can reproduce the signatures of J1331's counter-rotating core and predict the stellar and gas rotation curve at larger radii. However, our models do not agree with the observed velocity dispersion at large radii. We speculate that the reason could be a non-trivial change in structure and kinematics due to a possible merger event in J1331's recent past.

We present RoadMapping, a full-likelihood dynamical modelling machinery that aims to recover the Milky Way's (MW) gravitational potential from large samples of stars in the Galactic disk. RoadMapping models the observed positions and velocities of stars with a parametrized, three-integral distribution function (DF) in a parametrized axisymmetric potential. We investigate through differential test cases with idealized mock data how the breakdown of model assumptions and data properties affect constraints on the potential and DF. Our key results are: (i) If the MW's true potential is not included in the assumed model potential family, we can - in the axisymmetric case - still find a robust estimate for the potential, with only <~ 10% difference in surface density within |z| <= 1.1 kpc inside the observed volume. (ii) Modest systematic differences between the true and model DF are inconsequential. E.g, when binning stars to define sub-populations with simple DFs, binning errors do not affect the modelling as long as the DF parameters of neighbouring bins differ by < 20%. In addition, RoadMapping ensures unbiased potential estimates for either (iii) small misjudgements of the spatial selection function (i.e., <~ 15% at the survey volume's edge), (iv) if distances are known to within 10%, or (v) if proper motion uncertainties are known within 10% or are smaller than delta mu <~ 1 mas/yr. Challenges are the rapidly increasing computational costs for large sample sizes. Overall, RoadMapping is well suited to making precise new measurements of the MW's potential with data from the upcoming Gaia releases.