---
abstract: 'A star in the Milky Way’s disk can now be at a Galactocentric radius quite distant from its birth radius for two reasons: either its orbit has become eccentric through radial heating, which increases its radial action $J_R$ (‘blurring’); or merely its angular momentum $L_z$ has changed and thereby its guiding radius (‘churning’). We know that radial orbit migration is strong in the Galactic low-$\alpha$ disk and set out to quantify the relative importance of these two effects, by devising and applying a parameterized model (${\ensuremath{ \mathbf{p_m}  }}$) for the distribution $p(L_z, J_R, \tau, {\ensuremath{\mathrm{[Fe/H]}}}| {\ensuremath{ \mathbf{p_m}  }})$ in the stellar disk. This model describes the orbit evolution for stars of age $\tau$ and metallicity ${\ensuremath{\mathrm{[Fe/H]}}}$, presuming coeval stars were initially born on (near-)circular orbits, and with a unique ${\ensuremath{\mathrm{[Fe/H]}}}$ at a given birth angular momentum and age. We fit this model to APOGEE red clump stars, accounting for the complex selection function of the survey. The best fit model implies changes of angular momentum of $\sqrt{\langle \Delta L_z \rangle^2} \approx 619\, \mathrm{kpc~km/s~}(\tau/\mathrm{6~Gyr})^{0.5}$, and changes of radial action as $\sqrt{\langle \Delta J_R \rangle^2} \approx 63\, \mathrm{kpc~km/s~} (\tau/\mathrm{6~Gyr})^{0.6}$ at 8 kpc. This suggests that the secular orbit evolution of the disk is dominated by diffusion in angular momentum, with radial heating being an order of magnitude lower.'
author:
- Neige Frankel
- Jason Sanders
- 'Yuan-Sen Ting(丁源森)'
- 'Hans-Walter Rix'
bibliography:
- 'lit.bib'
title: 'Keeping it Cool: Much Orbit Migration, yet Little Heating, in the Galactic Disk'
---

[UTF8]{}[gbsn]{}

Introduction {#sec:introduction}
============

What dominates the secular orbit evolution of nearly isolated disk galaxies? The Milky Way’s last major merger is thought to have occurred before 7-8 Gyr ago [e.g. @rix_bovy_2013; @bland-hawthorn_gerhard_2016; @belokurov_etal_2018; @Helmi_etal2018], leaving a long subsequent period for internal processes to dominate the dynamical evolution of the disk: non-axisymmetries (bar, spiral arms) and giant molecular clouds can rearrange stellar orbits in different ways [e.g. @sellwood2014; @sellwood_binney_2002], as they resonantly interact or scatter.

The change in a star’s orbit can be decomposed into (1) ‘cold processes’: the orbit’s size or angular momentum changes, but remains circular, and (2) ‘heating processes’: the orbit’s eccentricity and vertical extent change. [@sellwood_binney_2002] dubbed the first process “radial migration" and postulated that it could be important for restructuring stellar discs. Pure radial migration occurs when stars are near corotation of a non axisymmetry and can change angular momentum ($L_z$) by some amount $\Delta L_z$, without significant change in their radial action ($J_R$): $\Delta J_R = 0 \Delta L_z$. Such secular change of angular momentum can be very substantial – at least in dynamical simulations – indeed of order unity [@roskar_etal_2008a; @minchev_2011; @kubryk_etal_2013; @halle_2015; @loebman_etal_2016]. Nonetheless, [@minchev_2011] and [@daniel_2019] used simulations to argue that such interactions should come with non-zero changes in $J_R$ due to possible resonance overlaps between, for example, co-rotation and Lindblad resonances of different non-axisymmetries, or higher order resonances. But this effect could also increase the angular momentum changes. Heating processes could arise from interactions with non-axisymmetric perturbations through other resonances, or other heating agents such as satellites [e.g. @VelazquezWhite_1999]. The different components of the velocity dispersion have been measured to increase with stellar ages both in the Solar neighbourhood and over the disk [@wielen_1977; @soubiran_2008; @sanders_das_2018; @ting_rix_2019; @mackereth_2019].

There is now persuasive evidence that [*radius*]{} migration in the Galactic disk is strong. In [@frankel_etal_2018] we measured this radius migration, which must be the combined effect of radial heating and angular momentum diffusion. We modeled radial migration as a global process, fitting the distributions $p(R|R_0, \tau)$ of present-day Galactocentric radii ($R$) as a function of their birth radii ($R_0$) and age ($\tau$) over a wide range of radii. We assumed that stars were born on initially tight metallicity-birth radius relations, and it is the the radial diffusion that introduces the observed present-day scatter in this relation [@edvardsson_etal_1993; @casagrande_etal_2011]. Inferring radial migration this way was originally proposed in the seminal paper by [@schonrich_binney_2009a]. Our global model fits turned out to constrain the overall migration scale well, to $\sigma_{RM} = 3.6 \mathrm{kpc}\sqrt{\tau/\mathrm{8 ~Gyr}}$: the typical star migrates about by a scale-length over the age of the disk.

However, this work did not disentangle diffusion in angular momentum from the increase in radial action, but only measured the combined effect of ‘heating’ and ‘cold $L_z$ diffusion’ as ‘orbit migration’. Therefore, the relative contributions of $\Delta J_R$ [or ‘blurring’ in the terminology of @sellwood_binney_2002; @schonrich_binney_2009a] and $\Delta L_z$ (or ‘churning’) to the evolution of the stars’ orbits in the Milky Way has yet to be quantified. This is what we set out to do here.

We now set out to build on that model and disentangle the strength of diffusion in angular momentum from that of increase in radial action. We will do this by generalizing the model of F18 from constraining the radius migration, i.e. the ‘diffusion rate’ in $R$, to constraining the secular orbit evolution in the disk plane, by quantifying the diffusion rates of both actions $L_z$ and $J_R$. In Section \[section\_data\], we describe our data set. We then construct the model and data likelihood in sections \[sec:model\] and \[sec:likelihood\] and present the result best fit model in section \[sec:results\]. We discuss the implications in Section \[sec:discussion\].

![image](2019-08-27_data_distributions.pdf){width="\textwidth"}

Data: APOGEE-DR14 Red Clump Giants {#section_data}
==================================

Since the disk build up and its secular dynamical evolution involve processes occurring on large spatial and time scales, our science case requires a sample of stars with extensive coverage of the Galactic disk, a wide range of ages, and accurate and reliable 6D phase-space coordinates to calculate orbits. In practice, this requires spectroscopic and astrometric information, precise distances (as these are the dominant uncertainty in the action determination), and a way to limit the impact of dust extinction at low Galactic latitudes. A crossmatch of the APOGEE-RC catalogue [@bovy_etal_2014; @majewski_etal_2017; @Abolfathi_2018_apogeedr14] with Gaia DR2 [@gaia_dr2_2018; @Lindegren_2018_gaiadr2] satisfies these desiderata: red clump stars are bright standard candles and APOGEE collected their spectra in near-infrared (limiting the effect of extinction).

Data Selection and Catalogs
---------------------------

We start with the stars from the 14th data release of the APOGEE near-infrared spectroscopic survey and restrict to the low-$\alpha$ sample as in [@frankel_etal_2018], with stars mostly born between $\sim$ 8 Gyr ago and now. We further select high-fidelity red clump stars derived by [@ting_hawkins_rix_2018] cross-matched with Gaia DR2.

Red clump stars are low-mass stars in the core helium burning stage. They are good standard candles [to $\sim 0.1$ mag, @girardi_2016; @hawkins_etal_2017; @hall_etal_2019; @sanders_smith_evans_2019; @chan_bovy_2019], allowing the determination of a precise photometric distance. The identification of the red clump stage, done in [@ting_hawkins_rix_2018], relies on spectroscopic estimates of the asteroseismic parameters $\Delta \nu$ and $\Delta P$ which contain information on the evolutionary stage and stellar mass.

We further restrict our sample to the low latitude ($|b|< 25$ deg) “short cohort” [as defined in @zasowski_2013] fields of APOGEE, which consists of the brightest stars with $H$ band apparent magnitude $7 \leq H \leq 12$ to ease our modeling and reduce the fraction of stars with large uncertainties. Stars in longer cohorts are fainter, so at larger distances, have greater distance uncertainties and are more extinguished. Including them would increase computational expenses and complicate the selection function without commensurately increasing the information content. So we restrict the analysis to the short cohort stars for which we can work out the probability that they were selected for targeting in APOGEE ([*i.e.*]{} we can determine the selection function). We reject all APOGEE “special targets” in our data set as well as all the stars that are in APOGEE fields for which we could not work out a probability of selection [see @frankel_etal_2019]. Since the number of crossmatch failures between the red clump sample and Gaia DR2 is negligible ($ < 1\%$), we assume in the following that the selection of our data is purely determined by APOGEE selection function, and that Gaia is complete within the APOGEE’s short cohort selection cuts.

Basic Data and their Uncertainties
----------------------------------

We use eight basic pieces of information: Galactic longitude ($l$), Galactic latitude ($b$), distance ($D$), line of sight velocity ($v_\mathrm{los}$), metallicity ([$\mathrm{[Fe/H]}$]{}) from APOGEE, proper motion in right ascension ($\mu_\alpha$) and declination ($\mu_\delta$) from Gaia, and age ($\tau$) derived from the full spectrum in [@ting_rix_2019], which are calibrated to asteroseismology.

The photometric distances ($D$) were determined to about 7% in [@ting_rix_2019], using near-infrared and Gaia G photometry for red clump stars as standard candles, exploiting the fact that interstellar extinction is weaker at longer wavelengths [e.g., @indebetouw_2005; @wang_2019_extinction]. The line of sight velocity ($v_\mathrm{los}$) is taken from the APOGEE-DR14 catalog, and the proper motions ($\mu_\alpha, \mu_\delta$) from Gaia DR2. Spectroscopic age estimates $\tau$ for this sample were obtained by [@ting_rix_2019] from a data-driven method built to determine ages from the APOGEE spectra, trained on the APOKASC2 red clump sample [@pinsonneault_2018]; this approach has a precision of about 0.15 dex, with possible systematics at large ages, because the C and N spectral features tracing the age, [e.g. @Martig2016; @ness_etal_2016] vary more weakly at large ages and hence contains less information. A more extensive discussion on the possible implications of such systematics on the modeling, and comparisons of different age estimates can be found in [@frankel_etal_2019]. The metallicity estimates ([$\mathrm{[Fe/H]}$]{}) are taken from the ASPCAP pipeline with typical uncertainties below 0.05 dex [@holtzman_2018A_dr14].

Galactocentric Rest Frame and Orbital Actions
---------------------------------------------

From these basic data, we extract and pre-compute the quantities that are more directly used in our model, and propagate uncertainties via Monte Carlo sampling of 80 points. In particular, we compute the Galactocentric radius ($R$), height above the plane ($z$), azimuthal velocity ($v_\phi$) and radial velocity ($v_R$) in Galactocentric coordinates, assuming the distance between the Sun and the Galactic center $R_\odot = 8.2 $ kpc [@gravity_2019], and the Sun’s height above the Galacctic plane $z_\odot \approx 20.8$ pc [@bennett_bovy_2019]. We also assume the Solar velocity with respect to the Local Standard of Rest $v_\odot = [-11.1, 12.24,~7.25]~\mathrm{km.s^{-1}}$[@schoenrich_2010] and the tangential velocity of the Sun 247.4 km/s [@gravity_2019; @reid_2004].

We further compute the orbital parameters relevant to our modeling of the Galactic disk; going from phase-space coordinates $(\vec{x},\vec{v})$ to orbits, as quantified [*e.g.* ]{} by their actions, requires the adoption of a gravitational potential. We compute the stars’ angular momenta ($L_z$, which are the azimuthal actions $J_\phi$ in an axissymmetric potential) and radial actions ($J_R$) using the Python package [*Galpy*]{} [@bovy_2015_galpy] based on the algorithm of [@binney_2012]. The orbits are integrated over the potential [*MWPotential2014*]{} of this package, which is also the potential used in our model, scaled to the Galactocentric rest-frame as described above. The 1D marginal distributions of angular momentum, radial action, age, [$\mathrm{[Fe/H]}$]{}  and Galactocentric radius are illustrated in Figure \[fig:distributions\_data\]. The distributions in angular momentum and Galactocentric radius are directly affected by the APOGEE on-sky footprint and implicit distance limit; else we would expect these distributions to be approximately exponential. The age distribution shows a prominent peak at at $\sim$2 Gyr, reflecting the mass dependence of the RC lifetime [@girardi_2016; @bovy_etal_2014].

Chemo-dynamical Model for the Evolution of the Galaxy’s Low-$\alpha$ Disk \[sec:model\]
=======================================================================================

We now present a global model for the formation and evolution of the Galaxy’s low-$\alpha$ disk, building directly on [@sanders_binney_2015], [@frankel_etal_2018] and [@frankel_etal_2019]. The most important astrophysical assumption of the model is that stars form on near-circular orbits from chemically well mixed cold gas with an inside-out star formation history. Over time, the orbits of stars evolve, and the gas is enriched in metals. The combination of age and [$\mathrm{[Fe/H]}$]{} in this model implies a birth angular momentum [$ L_{z0}$]{}  or a birth radius. [*Radial heating*]{} of stellar orbits is modeled as an increase of their mean radial action $J_R$ and [*radial migration*]{} is modeled as a global diffusion process in angular momentum $L_z$. In our modeling, vertical heating is only implicit: we ignore the (weak) coupling between in-plane and vertical motions, and model an age-dependent vertical profile for the disk to incorporate the 3D spatial selection function (see Section \[sec:likelihood\]). The overall model for the dataset $p_\mathrm{dataset}(l, b, D, {\ensuremath{\mathrm{[Fe/H]}}}, \tau, v_R, v_\phi|{\ensuremath{ \mathbf{p_m}  }})$, and how it is combined with the Galactic disk model $p_{MW}(L_z, J_R{\ensuremath{\mathrm{[Fe/H]}}}, \tau |{\ensuremath{ \mathbf{p_m}  }})$ is summarized in Figure \[fig\_graphical\_model\] and in Table \[table:model\_aspects\]. The model aspects are then combined together in the following section in a likelihood function, used to constrain its parameters. In this way we can disentangle the strength of radial migration and radial heating and compare them, but only under a set of (physically sensible) assumptions that we will now lay out in some detail.

Model Assumptions \[section\_assumptions\]
------------------------------------------

Our modelling assumptions are as follows:

1.  Stars are born with a tight relation between [$\mathrm{[Fe/H]}$]{} and [$ L_{z0}$]{} at any given $\tau$. This is a ‘weak chemical tagging’ assumption, supported by [@krumholz_ting_2018] and [@ness_etal_2019], who showed that in the low-$\alpha$ disk, [$\mathrm{[Fe/H]}$]{}and age can predict present-day $L_z$ precisely (and we here attribute the scatter in that relation to radial migration). . We parametrize this tight relation (Section \[section\_enrichment\]) and fit for the parameters.

2.  Secular processes (processes happening on time scales longer than a typical orbital time scale) have dominated the evolution of the Galaxy’s low-$\alpha$ disk. As the last major merger presumably occurred before the formation of the low-$\alpha$ disk [e.g., @rix_bovy_2013; @bland-hawthorn_gerhard_2016; @belokurov_etal_2018; @Helmi_etal2018], the Milky Way is thought to have evolved in near isolation (i.e. interacting mainly with itself and occasionally with lower mass satellites, e.g., the Sagittarius dwarf galaxy) for the past 7-8 Gyr, which leaves a lot of time for slower, more gradual processes to occur and affect the Galaxy’s evolution. We assume that secular processes cause orbits to diffuse, and set out to measure the strength of this diffusion, but we do not make assumptions on the nature of this process nor try to identify the agents driving them.

3.  We assume that the present-day potential of the Milky Way is sufficiently well approximated by the axisymmetric potential [*MWPotential2014*]{} in the [*Galpy*]{} package [@bovy_2015_galpy], and that the present equilibrium state of the Milky Way disk can be described by separable distribution functions $p({\bf J}) = p(L_z)p(J_R|L_z)p(J_z|L_z)$ where the vertical motion is independent of the radial action [the so-called adiabatic approximation, @binney_2010]. This is manifestly an approximation since we know the disk is not axisymmetric, nor in equilibrium: there are spiral perturbations, a bar and a warp [@bland-hawthorn_gerhard_2016; @beane_etal_2019].

4.  When computing the actions, we use the Stäckel approximation [@binney_2012] as implemented in the [*Galpy*]{} package [@bovy_2015_galpy]. Furthermore, we implicitly assume in the models that the height above the plane $z$ and vertical velocity $v_z$ are independent of the radial action [the adiabatic approximation, @binney_2010]. We use this approximation because we wish to focus on the in-plane distributions $p(J_R|L_z)p(L_z)$ and for our considered orbits, which predominately lie close to the Galactic plane, the approximation is valid.

Modeling the Gradual Build-up of the Stellar Disk \[subsection:pLzo\]
---------------------------------------------------------------------

We now describe the parameterized version of our modeling for the successive build-up of the Galactic stellar disk, encompassing the star-formation history with inside-out growth.

We model the time-integrated distribution of angular momenta at birth as $$\label{eq:Lz0_distribution}
    p({\ensuremath{ L_{z0}}}|{\ensuremath{ \mathbf{p_m}  }}) = \frac{{\ensuremath{ L_{z0}}}}{\langle {\ensuremath{ L_{z0}}}\rangle ^2} \exp\left(-\frac{{\ensuremath{ L_{z0}}}}{\langle {\ensuremath{ L_{z0}}}\rangle}\right),$$ where we define the mean angular momentum at birth as $\langle {\ensuremath{ L_{z0}}}\rangle = {\ensuremath{R_{d0}}}\times 235 \mathrm{km.s^{-1}}$, to fit for and interpret the parameter ${\ensuremath{R_{d0}}}$ as a global scale-length. This scale-length is time-integrated and reflects the global profile of the disk after all stars are born. For a cold disk at birth, this model is approximately equivalent to an exponential surface density profile ($\Sigma(R_0) \propto \exp(-{\ensuremath{R_0}}/ {\ensuremath{R_{d0}}})$) with a scale-length ${\ensuremath{R_{d0}}}$. But the actual spatial scale-length of the disk at birth may be different since the spatial distribution of stars will depend on gradual changes of the potential due to the on-going build up of the disk. The possible inside-out growth (illustrated in Fig. \[fig:model\_aspects\]) is modeled through an ${\ensuremath{ L_{z0}}}$-dependent star formation, where the star formation time-scale depends linearly on birth angular momentum, causing the inner disk to form stars on shorter time-scales than the outer disk. The star formation history is adapted from [@frankel_etal_2019], but now taken to be a function of birth angular momentum (rather than birth radius): $$\label{eq:SFH_def}
    \begin{split}
            \mathrm{SFH}(&\tau ~|~{\ensuremath{ L_{z0}}}, {\ensuremath{ \mathbf{p_m}  }}) = c({\ensuremath{ L_{z0}}}, {\ensuremath{ \mathbf{p_m}  }}) \\ 
            &\times \exp\left[\frac{1}{{\ensuremath{\tau_{\tiny{SFR}}}}}\left((1-{\ensuremath{ x_{io}}}\frac{{\ensuremath{ L_{z0}}}/235\mathrm{~km\,s}^{-1}}{\mathrm{8.2~kpc}})\tau - \tau_m \right) \right].
    \end{split}$$ Here, $c({\ensuremath{ L_{z0}}}, {\ensuremath{ \mathbf{p_m}  }})$ is a normalization constant such that $\int \mathrm{SFH}d\tau = 1$ at any given [$ L_{z0}$]{}. Since the star formation history is not the primary focus of this work, we treat it as a nuisance aspect of the model. We fit for and marginalize over the parameters ${\ensuremath{ x_{io}}}$ and [$\tau_{\tiny{SFR}}$]{}. The parameter $\tau_m$ corresponds to the maximum stellar age we consider in the low-$\alpha$ disk, fixed to 6 Gyr.

Present-day Gravitational Potential \[subsection:gravitational\_potential\]
---------------------------------------------------------------------------

We assume that the present-day gravitational potential of the Milky Way disk is well described by the [*MWPotential2014*]{} in the [*Galpy*]{} Package. In the present work, we only use the present-day potential of the Milky Way disk and its derived quantities (circular velocity $v_\mathrm{circ}$, epicyclic frequencies $\kappa$ and $\nu$ etc.) and make no assumptions about its past evolution.

As the Milky Way’s stellar disk has gradually grown from inside-out over the past 7-8 Gyr (Subsection \[subsection:pLzo\]), the mass distribution of the stellar and gas disks have changed, and the potential $\Phi_\mathrm{pot}(\tau)$ has evolved accordingly. Linking stellar birth radii ${\ensuremath{R_0}}$ to their birth angular momenta ${\ensuremath{ L_{z0}}}$ would require to know exactly how the potential has evolved. This could be done by modeling the mass distribution in different Milky Way components, including a growth for the stellar disk. But the present-day contributions of each component of the Milky Way are already fairly unconstrained today [e.g., @de_salas_2019; @eilers_etal_2019], so we do not attempt to infer them in the past. Instead, our modeling is fully based on linking birth actions to present-day actions rather than birth positions. In an axisymmetric potential, if the gradual build up of the disk is slow and adiabatic, then the actions of stars should be conserved: stars on the same orbit but at different phases will, if the change in potential is sufficiently slow, experience the same changes of potential (averaged over a period) and thus conserve their actions. However, stellar velocities and positions should change: as the disk mass increases, stars will on average sink to orbits closer to the Galactic center. Therefore, measuring a change of orbital action $\Delta \mathrm{J}$ bypasses other orbital changes of the stars and gives direct insights into the secular processes in the disk.

Modeling the Angular Momentum Evolution \[subsection:p\_Lz\]
------------------------------------------------------------

We follow the argument of [@sellwood_binney_2002] that the radial orbit redistribution of stars is caused by a sequence of stochastic processes, of some nature that we do not determine (which could be, for example, short-lived spiral perturbations). In this limit, stars follow a random walk in angular momentum so that radial migration can be modeled as a diffusion process in angular momentum, dubbed “churning” by [@sellwood_binney_2002]. Following [@sanders_binney_2015] we adopt the parameterized angular momentum diffusion equation $$\label{eq:diffusion}
    \frac{\partial f}{\partial t} = \frac{\partial}{\partial L_z}\left(-D^{(1)}
    f+ \frac{\sigma ^2}{2}\frac{\partial  f}{\partial L_z } \right).$$ The nature of the radial migration process determines the diffusion coefficients of this diffusion equation. For simplicity, we assume that the diffusion coefficients in the equation above do not depend on $L_z$ and that $\sigma$ is independent of time. These assumptions imply that we are only constraining an effective global mean (in radius and time) of the overall $L_z$ diffusion.

[@herpich_etal_2017] showed that if radial migration were asymptotically efficient, the angular momentum distribution in a disk should go as $f(L_z) = \exp(-L_z/\langle L_z \rangle)/\langle L_z \rangle$, and their model predictions match the observed angular momentum profile of external galaxies well. On the other hand, surface brightness profiles of disk galaxies have been observed to have exponential or de Vaucouleur profiles [@deVaucouleur_1948], which would correspond to (in the limit of circular orbits) an angular momentum distribution $f(L_z) \propto L_z  \exp(-L_z/\langle L_z \rangle)$ to which a bulge could be added at the center. Since we are not modeling a bulge in the present work, our model will in any case be inadequate in the inner few kpc of the Milky Way.

Drawing on these considerations, we impose that the steady state solution for the $L_z$ diffusion equation is either the exponential distribution in $L_z$, or the exponential surface density profile. At $L_z >> \langle L_z \rangle$[^1], these distribution are similar and constrain the diffusion coefficient $D^{(1)}$ to $D^{(1)} =  \frac{-{\ensuremath{\sigma_\mathrm{Lz}}}^2}{2 \langle {\ensuremath{ L_{z0}}}\rangle \tau_m}$, where we pose the mean specific angular momentum $\langle {\ensuremath{ L_{z0}}}\rangle = 235 \mathrm{km.s^{-1}} {\ensuremath{R_{d0}}}$. This ensures approximate conservation of angular momentum. In the microscopic limit, [@schonrich_binney_2009a] model this process with each star having a probability to move to a radius $r_i$ proportional to the stellar mass $m_i$ at that radius. This ensures that the number of stars migrating from radius $i$ to radius $j$ is proportional to $m_i m_j$, which is equal to the number of stars migrating from $j$ to $i$, thereby conserving the total disk profile. However, in reality the total disk angular momentum need not be conserved, as external torques, e.g. exchanges with the halo, could change it; here we do neglect this effect, since it is not well quantified or understood even in simulations: [@buck_2019] show that some Milky Way-like simulated galaxies see their scale-length increase over time (net outward migration), while others don’t. Here, with our fixed choice of $D^{(1)}$ disk profile remains constant with a mean angular momentum $\langle {\ensuremath{ L_{z0}}}\rangle$. But individual populations of age $\tau$ born on profiles with $\langle {\ensuremath{ L_{z0}}}\rangle (\tau) < \langle {\ensuremath{ L_{z0}}}\rangle$ will, on average, broaden, and those born with $\langle {\ensuremath{ L_{z0}}}\rangle (\tau) > \langle {\ensuremath{ L_{z0}}}\rangle$ will, on average, shrink.

The present day angular momentum of a star of age $\tau$ would then be related to its birth angular momentum by $$\label{eq:p_radial_migration}
\begin{split}
p({\ensuremath{ L_{z}}}~|~ & {\ensuremath{ L_{z0}}}, \tau,~ {\ensuremath{ \mathbf{p_m}  }}) = N({\ensuremath{ \mathbf{p_m}  }})\\
& \times \exp{\left(-\frac{({\ensuremath{ L_{z}}}- {\ensuremath{ L_{z0}}}- D^{(1)}\tau)^2}{2\sigma^2(\tau)}\right)}
\end{split}$$ where $\sigma(\tau) = {\ensuremath{\sigma_\mathrm{Lz}}}\sqrt{\frac{\tau}{\tau_m}}$ is the radial migration strength in angular momentum units (kpc km/s$^{-1}$). The parameter ${\ensuremath{\sigma_\mathrm{Lz}}}$ is to be fit and $N({\ensuremath{ \mathbf{p_m}  }})$ is a normalizing constant such that $\int_0 ^\infty p({\ensuremath{ L_{z}}}~|~ {\ensuremath{ L_{z0}}}, \tau, {\ensuremath{ \mathbf{p_m}  }}) d L_z= 1$ (there are no counter rotating stars in the disk). $\tau_m$ is the maximum age of the disk we consider, which is fixed to 6 Gyr here since there are so few constraining stars in the data at larger ages (see 3rd panel in Fig. \[fig:distributions\_data\]).

Radial Heating \[subsection:p\_heating\]
----------------------------------------

We assume that stars are born on near-circular orbits, and that their mean radial action increases with time as their orbits are kinematically heated (dubbed “blurring” by [@sellwood_binney_2002]). We adopt the isothermal disk model from [@sanders_binney_2015; @binney_2010], and we fit for the increase of mean radial action as a function of age and final position. The radial action distribution is then written as $$\label{eq:p_radial_action}
p(J_R~|~\tau, {\ensuremath{ L_{z}}},~ {\ensuremath{ \mathbf{p_m}  }}) = \frac{1}{2\pi}\frac{\kappa}{\sigma_R^2}
\exp\bigl (-\frac{\kappa J_R}{\sigma_R^2(\tau,{\ensuremath{ L_{z}}})}\bigr ),$$ with $\kappa = \kappa(R_\mathrm{circ}({\ensuremath{ L_{z}}}))$ the frequency of radial motion (the epicycle frequency), which depends on the guiding radius $R_\mathrm{circ}$, the radius of circular orbit of angular momentum [$ L_{z}$]{}. The velocity dispersion $\sigma_R = \sigma_R(\tau, R_\mathrm{circ}({\ensuremath{ L_{z}}}))$ traces the heating history of stars in the disk: $$\label{eq:heating}
\begin{split}
    \sigma_R(\tau, R_\mathrm{circ}({\ensuremath{ L_{z}}})) &= \sigma_\mathrm{vR0} \left(\frac{\tau + \tau_1}{\tau_m + \tau_1} \right)^\beta\\ 
    & \times \exp\left(\frac{8\,\mathrm{kpc} - R_\mathrm{circ}}{R_{\sigma_R}}\right),
    \end{split}$$ where $\sigma_\mathrm{vR0}$ is the velocity dispersion of stars of age $\tau_m = 6$ Gyr in the Solar neighbourhood (to be fitted), $\tau_1$ is set to reflect a dynamical time, $\sim$ 110 Myr, as in [@sanders_binney_2015] and allows stars to be born with small but non-zero eccentricity. The scale $R_{\sigma_R}$ represents a possible decay of the velocity dispersion with Galactocentric radius, as in [@sanders_binney_2015] (to be fitted). Note again that the radial actions are adiabatic invariants, and so the gradual build up of the disk should not influence the radial action of stars after their birth.

![image](2020-01-29_feht.pdf)

Weak Chemical Tagging: \[Fe/H\]-$\tau$-${\ensuremath{ L_{z0}}}$ Relation \[section\_enrichment\]
------------------------------------------------------------------------------------------------

For the current modeling, we use [$\mathrm{[Fe/H]}$]{}, and assume that it is approximately only a function of time and ${\ensuremath{ L_{z0}}}$; in some sense, this is the weakest form of “chemical tagging”.

We model the ${\ensuremath{\mathrm{[Fe/H]}}}$-$\tau$-${\ensuremath{ L_{z0}}}$ relation as a a power law in time and a broken line in the variable ${\ensuremath{ L_{z0}}}$ with an inner disk radial gradient of $-0.03~\mathrm{dex/kpc}$ for ${\ensuremath{ L_{z0}}}<3\times235\,\mathrm{kpc\,km/s}$. Since we have no constraining data on the inner metallicity profile in the disk, we use two external arguments to motivate this imposed flattening in the inner disk: (1) the observed metallicity profile in the stars is flatter in the inner disk/bulge [@pietrukowicz_2015], and (2) chemical evolution models predict a flatter metallicity profile in the gas in the inner disk [e.g., @schonrich_binney_2009a]. Namely, $$\label{eq:eq_Fe_Ro_t}
\begin{split}
{\ensuremath{\mathrm{[Fe/H]}}}= & {\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{max}f(\tau ) + b_{\ensuremath{\mathrm{[Fe/H]}}}\\ 
&+ {\ensuremath{\nabla\mathrm{[Fe/H]}}}({\ensuremath{ L_{z0}}})  \frac{{\ensuremath{ L_{z0}}}}{235~\mathrm{kms^{-1}}}
\end{split}$$ where the gradient $${\ensuremath{\nabla\mathrm{[Fe/H]}}}({\ensuremath{ L_{z0}}}) =\begin{cases}
    \nabla\mathrm{inner} &\text{$\frac{{\ensuremath{ L_{z0}}}}{235 \mathrm{km.s^{-1}}} < 3$ kpc}\\
    {\ensuremath{\nabla\mathrm{[Fe/H]}}}& \text{otherwise}.
  \end{cases}$$ with the inner metallicity gradient $\nabla\mathrm{inner}$ fixed at $-0.03~\mathrm{dex.kpc^{-1}}$ and the outer metallicity gradient ${\ensuremath{\nabla\mathrm{[Fe/H]}}}$ is to be fitted. $b_{\ensuremath{\mathrm{[Fe/H]}}}$ is introduced so the overall profile is continuous at ${\ensuremath{ L_{z0}}}= 3\times 235~\mathrm{kpc~km/s}$: $$b_{\ensuremath{\mathrm{[Fe/H]}}}=\begin{cases}
     0 &\text{$\frac{{\ensuremath{ L_{z0}}}}{235 \mathrm{km.s^{-1}}} < 3$ kpc}\\
    \frac{(\nabla_\mathrm{inner}-{\ensuremath{\nabla\mathrm{[Fe/H]}}}) {\ensuremath{ L_{z0}}}}{235~\mathrm{km~s^{-1}}} & \text{otherwise}.
  \end{cases}$$ The central metallicity is governed by the parameter ${\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{max}$ and we use a time dependence of $$\label{eq:enrichment_time}
f(\tau ) = \left( 1 - \frac{\tau}{12~\mathrm{Gyr}} \right)^{\ensuremath{{\gamma_\mathrm{[Fe/H]}}}}.$$ The set of model parameters we fit for are [$\nabla\mathrm{[Fe/H]}$]{}, [${\gamma_\mathrm{[Fe/H]}}$]{}, and ${\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{max}$. These functions are plotted in Fig. \[fig:tagging\_birth\_feh\].

Vertical Distribution of Stars \[subsection:p\_z\]
--------------------------------------------------

The vertical distribution of stars in the Milky Way disk is not the primary focus of this work (since we are essentially interested in the in-plane motions), but we must model it because the spatial selection of our data is in three spatial dimensions and so we require a 3D model for the disk. For simplicity, we model the vertical distribution of stars in the disk with the best fit model of [@ting_rix_2019] in the regime of the isothermal disk: $$\begin{split}
p(z~|~L_z, {\ensuremath{ L_{z0}}}, \tau, {\ensuremath{ \mathbf{p_m}  }}) &=  \frac{1}{2h_z(L_z, {\ensuremath{ L_{z0}}}, \tau)}\\ & \times \mathrm{sech}^2\left(\frac{z}{h_z(L_z, {\ensuremath{ L_{z0}}}, \tau)}\right),
\end{split}$$ where the scale-height $$h_z = a_z \sqrt{\frac{2 \overline{J_z}(L_z, {\ensuremath{ L_{z0}}}, \tau)}{\nu (R)}},$$ with the vertical frequency $\nu$ defined such that $\nu ^2 = \frac{\partial^2\Phi}{\partial z^2}$. $\overline{J_z}(L_z, {\ensuremath{ L_{z0}}}, \tau)$ is the mean vertical action of stars of angular present momentum $L_z$, birth angular momentum ${\ensuremath{ L_{z0}}}$ and of age $\tau$. [@ting_rix_2019] studied the vertical heating history of the Galactic disk using an APOGEE red clump data set, and published an analytic fit for $ \overline{J_z}(L_z, {\ensuremath{ L_{z0}}}, \tau)$. We model the age- and radial-dependent vertical distribution of red clump stars according to this form. Since the model of [@ting_rix_2019] used different age dependencies and birth radii scales than those in our work, we allow for an overall scaling of their heating law parametrized by $a_z$, which we anticipate will be near unity (see resulting fits in Table \[table:model\_aspects\]).

  Model aspect                                                                                                  Functional Family                                                                                                                                Model Parameters                                 Max. Likelihood Value
  -------------- ---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- ------------------------------------------------------------------------------ -----------------------
                                                                                       [ [$p({\ensuremath{ L_{z0}}}| {\ensuremath{ \mathbf{p_m}  }})$]{}]{}                                                                                                   [Birth scale-length]{}                             
                                                              $\propto {\ensuremath{ L_{z0}}}\exp{(-{\ensuremath{ L_{z0}}}/ {\ensuremath{v_\mathrm{circ}}}{\ensuremath{R_{d0}}})} $                                                                      [[${\ensuremath{R_{d0}}}$]{}]{}                               [2.8 kpc]{}
                                                                                     [ $ p(\tau | {\ensuremath{ L_{z0}}}, {\ensuremath{ \mathbf{p_m}  }})$]{}                                                                                  [Inside-out coefficient ${\ensuremath{ x_{io}}}$]{}                        0.65
                  [$ \propto \exp\left[\frac{1}{{\ensuremath{\tau_{\tiny{SFR}}}}}\left((1-{\ensuremath{ x_{io}}}\frac{{\ensuremath{ L_{z0}}}/{\ensuremath{v_\mathrm{circ}}}}{\mathrm{8.2~kpc}})\tau - \tau_m \right) \right]$]{}           [[SFR time-scale $ {\ensuremath{\tau_{\tiny{SFR}}}}$]{} ]{}                    1 Gyr
                                                                                                                                                                                                                                                               Maximum age $\tau_m$                                  Fixed to 6 Gyr
                                                                                   [$ p(L_z | {\ensuremath{ L_{z0}}}, \tau, {\ensuremath{ \mathbf{p_m}  }})$]{}                                                                                               [Diffusion strength]{}                             
                                               [$\propto \exp{\left(-\frac{({\ensuremath{ L_{z}}}- {\ensuremath{ L_{z0}}}- D\tau)^2)}{2{\ensuremath{\sigma_\mathrm{Lz}}}^2\tau/\tau_m} \right)}$]{}                                                        [[$\sigma_\mathrm{Lz}$]{}]{}                                582 kpckm/s
                                                                                  [ $ p(J_R | {\ensuremath{ L_{z}}}, \tau,  {\ensuremath{ \mathbf{p_m}  }})$]{}                                                                                            [Time dependence $\beta$]{}                                     0.3
                                                                                                   [$\propto \exp(-\kappa J_R / \sigma_R^2)$]{}                                                                                                     Velocity dispersion $\sigma_\mathrm{vR0}$                            49 km/s
                                            $    \sigma_R = \sigma_\mathrm{vR0} \left(\frac{\tau + \tau_1}{\tau_m + \tau_1} \right)^\beta \exp\left(\frac{8\,\mathrm{kpc} - R_c}{R_{\sigma_R}}\right)$                                                   Radial dependence $R_{\sigma_R}$                                19 kpc
                                                                                   $ p(z | L_z, {\ensuremath{ L_{z0}}}, \tau, {\ensuremath{ \mathbf{p_m}  }})$                                                                                                   $a_z$ (scaling)                                          1.16
                                                                        $ p({\ensuremath{\mathrm{[Fe/H]}}}| {\ensuremath{ L_{z0}}}, \tau, {\ensuremath{ \mathbf{p_m}  }})$                                                                        Time dependence [${\gamma_\mathrm{[Fe/H]}}$]{}                          0.456
                                    [$ {\ensuremath{\mathrm{[Fe/H]}}}= {\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{max} \left( 1 - \frac{\tau}{12~\mathrm{Gyr}} \right)^{\ensuremath{{\gamma_\mathrm{[Fe/H]}}}}$]{}                                   [Radial gradient [$\nabla\mathrm{[Fe/H]}$]{}]{}                     -0.0936 dex/kpc
                                 [$ + {\ensuremath{\nabla\mathrm{[Fe/H]}}}({\ensuremath{ L_{z0}}})  \frac{{\ensuremath{ L_{z0}}}}{v_\mathrm{circ}({\ensuremath{ L_{z0}}})} + b_{\ensuremath{\mathrm{[Fe/H]}}}$]{}                  [Max [$\mathrm{[Fe/H]}$]{}  ${\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{max}$]{}          0.7 dex

 The Likelihood of the APOGEE-RC $\times$ Gaia DR2 Data \[sec:likelihood\]
==========================================================================

To determine how this model is constrained by our data, we must construct the data’s likelihood for any given set of model parameters, and on this basis sample the parameter’s posterior probability distribution, [*pdf*]{}. We now lay out how to implement this.

For each star, our ‘data’ are $$\mathcal{D} = \{l, b, D, \tau, {\ensuremath{\mathrm{[Fe/H]}}}, v_R, v_\phi\}$$ with their associated uncertainties (see Section \[section\_data\]). As noted in Section \[section\_data\], these data are not sampled directly from the underlying Milky Way disk distribution: the sample very much reflects both the distribution of stars in the Milky Way $p_{MW}(\mathcal{D})$, and the selection process $S(l, b, D, \tau)\equiv p(\mathrm{select}| \mathcal{D})$. Therefore, the probability distribution of the data given our model parameters is $$\begin{split}
p_\mathrm{dataset}(\mathcal{D}~&|~ {\ensuremath{ \mathbf{p_m}  }}) = C\cdot p_\mathrm{MW}(\mathcal{D} ~|~{\ensuremath{ \mathbf{p_m}  }})\\
&\times S(l,b,D)f_{RC}(\tau)
\end{split}$$ where $p_\mathrm{MW}(\mathcal{D} ~|~{\ensuremath{ \mathbf{p_m}  }})$ is the model of the Galactic disk, which is a combination of the model aspects described above as detailed in Appendix \[sec:appendix\_model\]. $C$ is a normalizing term specified in Eq. \[eq:survey\_volume\] and Section \[sec:survey\_volume\]. $S(l,b,D) = p(\mathrm{select}|~l, b, D)$ is the selection function, or the probability that a star ends up in the catalog given its observable properties. The observables for each star are a combination of its intrinsic properties and its position and velocity with respect to us, and the evolutionary stage is determined using the spectrum [@ting_hawkins_rix_2018]. Typically the selection function is a strong function of apparent magnitude. In the case of standard candles where the apparent magnitude is a function of distance (and extinction) only, such as RC stars, the selection function essentially reduces to a function of $(l,b,D)$. Finally, the term $f_{RC}(\tau)$ is the probability of a star to be on the red clump evolutionary stage given its age $\tau$ [@bovy_etal_2014].

We determine the spatial selection function of APOGEE following the methods laid out in [@rix_bovy_2013; @bovy_etal_2014] and [@frankel_etal_2019], and extend the methodology to the DR14 data which includes APOGEE-2. The main difference in the APOGEE-2 target selection for the main disk fields with respect to APOGEE-1 is the inclusion of two color bins: blue ($(J-K_s) < 0.8$ mag) and red ($(J-K_s) > 0.8$ mag) [@zasowski_2017]. This makes the selection function more complex because the fraction of selected stars of different stellar types differs between APOGEE fields. However, we are not affected by this complexity because we are working with RC stars, assumed perfect standard candles [with $(J-K)_0 = 0.68$ mag, @hawkins_etal_2017] lying fully in the blue color bin, such that we need not account for the selection fractions of the red color bin. In practice, the spatial selection function, $S(l,b,D)$, is a piece-wise function in each APOGEE field, such that for a field $i$ centred at $(l,b)_i$ it can be expressed solely as a function of distance $S_i(D)$. The dependence on distance is more complex than a constant between $D_{min}$ and $D_{max}$ set by the stars absolute magnitude and the survey’s magnitude cuts, due to the 3D spatial distribution of dust, which limits the fraction of stars seen at a given absolute magnitude and distance. We model this using the 3D dust map of [@green_etal_2019]. A detailed description of how this is incorporated into the model can be found in [@frankel_etal_2019][^2].

Normalization of the PDF: Survey Volume \[sec:survey\_volume\] 
---------------------------------------------------------------

To make $p_\mathrm{dataset}$ a probability density function, it must be normalized by $$\label{eq:survey_volume}
C^{-1} 
= V_S({\ensuremath{ \mathbf{p_m}  }}) 
= \int _\mathcal{D} p_\mathrm{MW}(\mathcal{D} ~|~{\ensuremath{ \mathbf{p_m}  }}) S(l,b,D) \mathrm{d}\mathcal{D}$$ which is a 7-dimensional integral over all the physical properties of the data We refer to $V_S$ as the ‘survey volume’. Two of the integrals $(l, b)$ can be transformed into a sum over APOGEE fields, if we assume that the properties of stars in the sky within a single APOGEE pointing are uniform. This is a valid assumption as the APOGEE fields are typically 3 deg across. We compute $V_S$ in the next subsection, and assemble the data likelihood in Subsection \[subsection:likelihood\].

Expanding Eq. \[eq:survey\_volume\], the survey volume can be written $$\begin{split}
    V_S({\ensuremath{ \mathbf{p_m}  }}) &= \sum_{\mathrm{field}~i} \int_\mathrm{Data~space} p_\mathrm{MW}(\tau, {\ensuremath{ L_{z0}}}, J_R, {\ensuremath{ L_{z}}}, z, ~|~{\ensuremath{ \mathbf{p_m}  }})\\
    & \times S_i(D) f_{RC}(\tau) \Omega_i D^2  \mathrm{d}{\ensuremath{ L_{z0}}}\mathrm{d}D \mathrm{d}v_R
    \mathrm{d}v_\phi \mathrm{d}\tau,
\end{split}$$ This integral is not straightforward to compute. The integrand is proportional to a probability density that we cannot fully compute, but that we can sample. Therefore, we compute this integral through iterative importance sampling (details in Appendix \[sec:appendix\_Vs\]). We start with a ‘blind’ optimization of the data likelihood by normalizing the likelihood with samples generated from a model with initial plausible guesses for the model parameters $\mathbf{p_{m,\mathrm{prop}}}$. We then use these best fit parameters to generate new samples, which better reflect the shape of the likelihood and can be used normalize the likelihood for further calculations.

![image](2020-01-17_chains_10000_betasrm_1p49.pdf)

Data Likelihood Function and Parameters Posterior\[subsection:likelihood\]
--------------------------------------------------------------------------

The overall Milky Way disk model combined with the selection function predicts the likelihood of the data for any star in the sample. Assuming all measurements are independent, we write the total likelihood of the entire data set given our model with parameters [$ \mathbf{p_m}  $]{} as $$\label{eq:likelihood}
\begin{split}
    &p_\mathcal{L}( \{{\ensuremath{\mathrm{[Fe/H]}}}, \tau, l,b,D,v_R, v_\phi\}|{\ensuremath{ \mathbf{p_m}  }}) \\
    &= \prod _{i=1}^{N_\mathrm{stars}} p_\mathrm{dataset}({\ensuremath{\mathrm{[Fe/H]}}}_i, \tau_i, l_i,b_i,D_i,v_{R,i}, v_{\phi,i}|{\ensuremath{ \mathbf{p_m}  }}).
\end{split}$$ We use uniform priors with wide ranges in the parameter space, and enforce distances, spatial scales and time-scales to be positive. We first maximize the likelihood [@NeldMead65], which gives the results in Table \[table:model\_aspects\]. We start the optimizer from different initial conditions to lower the chances of the optimizer becoming stuck in local maxima. Using the MCMC sampler [*emcee*]{} [@foreman-mackey_2013], we then sample the posterior, $$\begin{split}
    &p_\mathrm{pos}({\ensuremath{ \mathbf{p_m}  }}~|~ \{{\ensuremath{\mathrm{[Fe/H]}}}, \tau, l,b,D,v_R, v_\phi\})\\
    & \propto p_\mathrm{prior}({\ensuremath{ \mathbf{p_m}  }}) p_\mathcal{L}(\{{\ensuremath{\mathrm{[Fe/H]}}}, \tau, l,b,D,v_R, v_\phi\}~|~{\ensuremath{ \mathbf{p_m}  }})
\end{split}$$ using 12,000 iterations ($>$ 50 times the auto-correlation time) and 52 walkers. We initialize the MCMC sampling uniformly in a hypercube (of size greater than 8 times the error bars quoted in Fig \[fig:results\_corner\]) centered on the maximum likelihood estimates.

The posteriors of the model parameters are illustrated in Fig. \[fig:results\_corner\]. The parameters of interest are the dynamical parameters: $\{{\ensuremath{\sigma_\mathrm{Lz}}}, \sigma_{vR0}, \beta, R_{\sigma_R} \}$, in red and orange in Fig. \[fig:results\_corner\]. We treat the other model parameters as nuisance parameters and marginalize over them, although we comment briefly on them in Section \[sec:results\]. There are no degeneracies and only weak covariances. In particular, the estimate of radial migration (or spread in angular momentum) [$\sigma_\mathrm{Lz}$]{} is slightly correlated with the metallicity gradient. This is expected because the information on radial migration comes from the scatter in metallicity. Stars born in a galaxy with a shallower metallicity gradient will need to migrate larger distances to produce the same metallicity spread. The other notable covariance is between $\sigma_{vR0}$ and $\beta$, which together quantify heating in the radial direction: $\sigma_{vR0}$ is the strength, and $\beta$ the time dependency. If most of the sample is dominated by stars of a given age (here 2 Gyr), then the two ways to reproduce the local velocity dispersion of stars of that age is to either have a greater strength, and a faster time evolution (small $\beta$), or a smaller strength and a more linear time evolution (larger $\beta$).

The other model aspects also seem well fitted since the overall distributions of the data are well recovered by the model (Fig. \[fig:distributions\_data\]) and their estimate is robust to the tests we have performed (Subsection \[subsection:tests\]). This shows that the best fit model is also qualitatively a good fit. The parameter analogous to the scale-length of the disk at birth is about 2.9 kpc. The covariances are more important for the inside-out growth model aspects, which are treated as nuisance parameters in the present work. The measurements of inside-out growth are consistent with our previous estimate in [@frankel_etal_2019] with ${\ensuremath{ x_{io}}}\approx 0.7$, where the implications and shortcomings of fitting data with large age uncertainties are discussed extensively.

![image](2020-01-05_model_aspects_polar.pdf)

Best fit Milky Way Disk Model \[sec:results\]
=============================================

Fig. \[fig:results\_corner\] illustrates that all model parameters are well constrained, that the [*pdfs*]{} are approximately Gaussian, and that for most parameter combinations the covariances are small. The model seems to be well-posed for and well-constrained by the data set. This holds true for its dynamical aspects (brown in Fig. \[fig:results\_corner\]), inside-out growth (green), and ISM enrichment (blue). The best fit model is illustrated in a schema in Fig. \[fig:model\_aspects\].

In the following subsections, we look at these different aspects more closely.

Migration strength and Age-Radial Velocity Dispersion
-----------------------------------------------------

The strength of radial migration is encompassed in the model aspect ${\ensuremath{\sigma_\mathrm{Lz}}}(\tau)$, the width of the distribution of stars of age $\tau$ about their mean angular momentum. We find ${\ensuremath{\sigma_\mathrm{Lz}}}(\tau) = 567~\mathrm{kpc~km/s} \sqrt{\tau / 6~\mathrm{Gyr}}$. With a circular velocity of about $235$ km/s, this corresponds to a migration scale of $2.4$ kpc for the 6 Gyr stars.

Radial heating leads stars to increase their random motion in the radial direction. Using Eq. \[eq:heating\], we find that the velocity dispersion at the Sun is about 43 km/s for the $\sim$ 6 Gyr stars, in line with [@nordstrom_etal_2004] and [@mackereth_2019], and that the age dependency is $\tau^{0.3}$. In the epicycle approximation, radial heating leads to epicycle amplitudes typically of $A_\mathrm{epi} = \sqrt{ \sigma _r^2(\tau)/\kappa ^2 }\approx 1.5 \mathrm{kpc}$ for the older stars, as illustrated in Fig. \[fig:model\_aspects\]. The parameter $R_{\sigma_R}$ best fit value is $\approx 27$ kpc, whereas it has been commonly assumed to be of order 2 $\times$ disk scale length ($\sim6-7$ kpc). The original motivation for the exponential decay of the radial velocity dispersion was to keep a disk with a constant scale height $h_z$ and a constant ratio $\sigma_r/\sigma_z$ [@van_der_kruit_1982], with $\sigma_z \propto h_z\sqrt{\rho} \propto \exp(-R/2R_d)$. But (1) the agents driving the evolution of the vertical and radial motion are likely different [e.g @sellwood2014] so $\sigma_r$ and $\sigma_z$ need not be related, and (2) the Galactic disk is now known and expected to flare [@ting_rix_2019; @minchev_2015; @bland-hawthorn_gerhard_2016; @bovy_etal_2016; @kawata_2017; @sanders_das_2018; @mackereth_etal_2017; @mackereth_2019]: the scale-height of coeval stellar populations increases with radius, hence there is no need for $R_{\sigma_R}$ to be small.

Inside-out Growth and the Metallicity Profile
---------------------------------------------

The model fit favors a global inside-out growth of the disk, where stars formed first from low angular momentum gas, and star formation moved gradually to higher angular momentum ($x = 0.68$). This is in accord with the results of [@frankel_etal_2019], who used a similar model in radius on the APOGEE-RC data set of the 12th data release. Our model is illustrated for different values of $x$, as well as for the best fit, in Fig. \[fig:model\_aspects\] in green (top right quadrants). As the disk forms from inside-out, the enrichment in metals (here iron, \[Fe/H\]) proceeds with a radial gradient and still on-going enrichment (Figure \[fig:tagging\_birth\_feh\]).

The Orbit-Age-Abundance Distributions
-------------------------------------

The distributions most directly affected by radial mixing are the metallicity distribution functions of stars at given radius $p({\ensuremath{\mathrm{[Fe/H]}}}|R)$ [@hayden_2015; @loebman_etal_2016], whose shape is influenced by the amount of metal-rich stars incoming from the inner disk the metal-poorer stars coming from the outer disk, and the initial metallicity profile of the gas from which stars formed, set by the inside-out star formation history [@schoenrich_mcmillan_2017; @schonrich_binney_2009a].

The available data, and the framework we developed in this work, allow us to make comparisons in more dimensions, reproducing the entire data set’s 5-dimensional distribution $p({\ensuremath{\mathrm{[Fe/H]}}}, \tau, J_R, L_z, R)$ resulting from the global evolution of the Milky Way disk (and selection procedures). Figure \[fig:data\_space\_comparisons\] illustrates the distributions $p({\ensuremath{\mathrm{[Fe/H]}}}, \tau, J_R, L_z, R)$ observed in the data set (brown) and predicted by the best fit model for the data set (green). Overall, the observed distributions are globally well recovered by the model fit. The metallicity radial profile (bottom left panel of the corner plot in Fig \[fig:data\_space\_comparisons\]), which arises from (1) chemical evolution and (2) subsequent orbit evolution is well described, and illustrates that the metallicity distribution functions $p({\ensuremath{\mathrm{[Fe/H]}}}, R)$ and the age-metallicity distributions (4th panel, bottom), usually the main diagnostic used in the literature, are all well reproduced. The distributions in $L_z$ and $J_R$ reflect both the secular evolution of the disk and APOGEE’s spatial selection function.

![image](2019-12-25_dataspace.pdf){width="\textwidth"}

Verifications: Model Variants and Parameter Recovery \[subsection:tests\]
-------------------------------------------------------------------------

During the construction of this model, we tested a set of model variants. We first tested the parameter recovery on noised mock data, and rejected all models for which we could not recover the true parameters. For example, we could not add a parameter to quantify how much angular momentum is actually conserved while stars diffuse in angular momentum (by e.g., fitting for a simple form of diffusion coefficient $D^{(1)}$ from Eq. ) because there is a degeneracy with the enrichment model.

Once the best model candidate (which is presented in Section \[sec:model\]) was identified, the best fit parameters found with MLE, and the posterior sampled with MCMC, we performed a series of additional tests to verify different aspects of our results. We investigated whether uncertainties in our integral calculations introduced systematics or biases in the estimation of the parameters, by computing it using different model realizations (by changing the model parameters) and Monte Carlo samples of different sizes, and found no change in the results within 1$\sigma$ (for both the data used in the present paper and mock samples generated from our model). We generated and fitted mock data with different noise levels (increasing or decreasing the formal uncertainties by a few percent), and the best fit parameters are well recovered too.

We studied the effect of varying the model of the potential on the estimates of the dynamical parameters. Using the default [*MWPotential2014*]{} in [*Galpy*]{} with a circular velocity of 220 km/s at 8 kpc, we found small changes in the estimates of [$\sigma_\mathrm{Lz}$]{}. These changes are however expected and quantifiable. The estimate for radial migration dropped by about $\sigma_{Lz220} \approx {\ensuremath{\sigma_\mathrm{Lz}}}\times \frac{220}{235}$, which is what one expect since the metallicity scatter, and the radial metallicity gradient in the stars are the same, and ${\ensuremath{\sigma_\mathrm{Lz}}}\approx \sigma_{\ensuremath{\mathrm{[Fe/H]}}}v_\mathrm{circ}/{{\ensuremath{\nabla\mathrm{[Fe/H]}}}}$. But a full exploration of alternate gravitational potentials is beyond the scope of the present work. The potential model we are currently using is well constrained by external data [@bovy_rix_2013; @bovy_2015_galpy].

Model Limitations {#subsection:model_limitations}
-----------------

Forward modeling the orbit-age-\[Fe/H\] structure of the APOGEE$\times$Gaia dataset, with an interpretable model that accounts simultaneously for diverse aspects of Galaxy evolution, data uncertainties and the survey selection function, provides a framework with a great potential for Galactic archaeology. The best fit model is well tested and reproduces the observed trends and distributions of the data set well. However, at present our model lacks some features that one might desire from a full physical model of the Galaxy:

1.  Radial migration of stars in the Milky Way is measured indirectly through its impact on the age-metallicity distributions. The strongest assumption we have made in that direction is that the birth age-metallicity relation was tight and monotonic, and that our choice of functional forms to model the evolution of \[Fe/H\] were flexible enough and adequate. However, any inadequacy or inappropriate rigidity in that model will be measured as radial migration in this context, so ${\ensuremath{\sigma_\mathrm{Lz}}}$ could be a lower limit on the strength of radial migration.

2.  Our description of $L_z$ change is only an approximation of the solution to the diffusion equation (Eq. \[eq:diffusion\]) that is valid far from the Galactic center, and is not self-consistent. It is likely that the strength of migration is a function of radius and time [e.g. @kubryk_etal_2013; @toyouchi_chiba_2018], which in our simple diffusion picture is ignored. Additionally, Eq. \[eq:diffusion\] should contain a source term for star formation, which we have modeled separately. This should not impact the results drastically, as the two extreme regimes are recovered: in the limit where radial migration is asymptotically inefficient, stars’ current $L_z$ distribution is a Dirac function of their birth distribution, and the overall $L_z$ distribution is the exponential profile used for initialization in Eq. \[eq:Lz0\_distribution\]. In the limit where radial migration is asymptotically efficient, our imposed steady state solution for the $L_z$ distribution recovers the results of [@herpich_etal_2017] for a flat circular velocity curve, with the same scale length [$R_{d0}$]{}.

3.  The secular evolution processes in the disk (diffusion in $L_z$ and increase in $J_R$) were treated independently. However, they should be covariant depending on the heating agents, (e.g. as in Eq. \[eq:sellwood\_binney\_delta\_Jr\]) . Ideally, we should treat the entire Fokker-Planck equation in action space. However, we argue that this should have only a small impact on the present results, because the disk remains relatively cold ($J_R$ does not increase much) and $L_z$ diffusion is stronger by an order of magnitude . This may blur out such covariances on large spatial and time-scales. However, in the Solar neighbourhood, there are indeed clear over-densities in the $L_z-J_R$ plane, some of them arising along Lindblad resonances with the bar or other non-axisymmetries [@sellwood_binney_2002; @trick_coronado_rix_2019; @trick_2019]). These portray the impact of the most recent set of non-axisymmetries.

4.  The treatment of radial heating, with the time dependence $\beta$ could depend on angular momentum. Indeed, as shown in simulations, the heating time dependence $\beta$ could depend on spiral arms and Galactocentric radius [@aumer_2016; @BT08]. The radial dependence of $\beta$ in the age-velocity dispersion relation in the Milky Way disk was confirmed in [@mackereth_2019]. The present model fit leads to $\beta \approx 0.3$, which is typical and agrees with the largest part of the disk values in [@mackereth_2019], but not in the outer disk where their $\beta $ decreases, possibly due to the weakening of spiral strength.

![Secular evolution of the Milky Way disk implied by the best fit model, as diffusion in angular momentum (top) and radial action (bottom) drawn from the MCMC samples illustrated in Fig \[fig:results\_corner\] and using the model equations, Galactocentric radii 6, 8 and 12 kpc. The right-hand side y axis shows the equivalent in distance units. Top: we assumed a constant circular velocity curve $R_\mathrm{circ} = L_z/v_\mathrm{circ}$ with $v_\mathrm{circ} = 235$ km/s. Bottom, shown for the Solar Galactocentric radius: we assume the radial spread about the guiding radius due to the epicycle motion corresponding to a radial action $J_R$ at $R$ is $\sqrt{\langle (R-R_g)^2 \rangle} = \sqrt{J_R/\kappa} = A/\sqrt{2}$ with $\kappa$ the epicycle frequency and $A$ the epicycle amplitude.[]{data-label="fig:secular_evolution"}](2020-01-20_Lz_Jr_r.pdf){width="\columnwidth"}

Astrophysical Implications and Discussion \[sec:discussion\]
============================================================

Secular processes drive stellar orbits to evolve, which can have a large impact on a disk galaxy’s evolution. If strong, such diffusion processes lead to a near complete dynamical memory loss, challenging efforts in Galactic archaeology to try to infer a galaxy’s history from its present-day properties. Radial migration is such a process and was shown to be strong in simulations [e.g. @minchev_etal_2013], but had not been well quantified across the Milky Way disk.

We have presented the first global model relating stellar ages, chemistry, and dynamics, which we have fitted to high-quality data over a large extent of the Galactic’s low-$\alpha$ disk, accounting rigorously for data uncertainties and selection effects. The model builds on and extends [@sanders_binney_2015]’s pioneering framework. The previous modeling to [@frankel_etal_2018] was based on more local data and did not incorporate stellar ages. In this work, we have forward-modelled the distributions of stellar ages, \[Fe/H\], $L_z$ and $J_R$ with a radially-dependent star formation history accounting for the inside-out growth of the disk, a parametric chemical evolution model, and birth and action-based radial heating and radial migration orbit-evolution model. We have fitted all aspects of this model simultaneously, but focus on the orbit evolution, and treated the other model aspects as a nuisance. We have used the APOGEE$\times$Gaia red clump sample, a dataset that covers a large part of the Galactic disk, from $R \sim 4$ kpc to 13 kpc, a scale larger than the typical radial migration scale of a few kpc.

Large efforts have been made to measure the strength of radial migration in the Milky Way, but previous works lacked data on large spatial scales, making quantifying a large scale diffusive process a hard task. Furthermore, previous approaches did not use stellar ages [e.g. @sanders_binney_2015] instead relying more tangentially on the known local relationships between stellar age and kinematics. Other methods measured a scatter in the \[Fe/H\]-\[$\alpha$/Fe\]-$R$ directly, but without accounting simultaneously for the galactic evolution processes that could contribute to it, thus lacking a framework that could use and describe the entire dataset [e.g. @hayden_2015]. In [@frankel_etal_2018; @frankel_etal_2019], we developed a framework accounting for the main evolution aspects of the Milky Way disk, with a radially-dependent star formation history, chemical evolution of the disk, and evolution of the stars’ Galactocentric radius. This model constrained a global orbit migration scale of about 3 kpc $\sqrt{\tau/6~\mathrm{Gyr}}$, implying that radial mixing happens on scales comparable to the scale-length of the Milky Way disk. However, this description of the disk only measured diffusion in Galactocentric radius, and not angular momentum, so failed to disentangle the two major processes causing stars to change Galactocentric radius (‘churning’ and ‘blurring’).

Secular Dynamical Evolution
---------------------------

Our model describes the in-plane secular evolution of the Milky Way’s low-$\alpha$ disk and disentangles the contributing processes: diffusion in angular momentum (‘churning’) and increase in radial action (‘blurring’, or ‘radial heating’). As both processes are diffusion in action space, we can quantitatively compare their strengths in a meaningful way, and here we choose to inspect the root-mean-squared (rms) deviation in the actions.

We first work out the rms deviations of the actions expected in the Solar neighbourhood from external data, and will then show that our more global model recovers this particular case. From a simplistic perspective, the rms deviation in the radial action is related to the disc properties as $$\sqrt{\langle(\Delta J_R)^2\rangle}\approx\frac{\sqrt{2}\sigma_R^2}{\kappa}.$$ Taking the ‘textbook’ quantities for the radial epicyclic frequency $\kappa\approx37\,\mathrm{km/s/kpc}$ and $\sigma_R\approx38\,\mathrm{km/s}$ for old stars from [@BT08], we find $\sqrt{\langle(\Delta J_R)^2\rangle}\approx55\,\mathrm{kpc\,km/s}$. Likewise, the rms deviation in the angular momentum can be simply expressed as $$\sqrt{\langle(\Delta L_z)^2\rangle}\approx\frac{v_\mathrm{circ}\sigma_\mathrm{[Fe/H]}}{|\mathrm{d[Fe/H]}/\mathrm{d}R|}.
\label{eqn::lzsimple}$$ Using approximate values for the solar neighbourhood of $\sigma_\mathrm{[Fe/H]}=0.2\,\mathrm{dex}$ for the metallicity dispersion [@nordstrom_etal_2004], $|\mathrm{d[Fe/H]}/\mathrm{d}R|=0.062\,\mathrm{dex/kpc}$ for the radial metallicity of young stars [@lucklambert2011] and $v_\mathrm{circ}=235\,\mathrm{km/s}$ we find $\sqrt{\langle(\Delta L_z)^2\rangle}\approx750\,\mathrm{kpc\,km/s}$. These simple calculations confirm that plausibly $\sqrt{\langle(\Delta L_z)^2\rangle}$ is an order of magnitude larger than $\sqrt{\langle(\Delta J_R)^2\rangle}$.

Expanding to the greater extent of the disk with our full model, from Eq. \[eq:p\_radial\_migration\], the variance of the angular momentum distribution of a stellar population of age $\tau$ is $$\langle (L_z - L_{z0})^2 \rangle = \sigma(\tau)^2 + (D\tau)^2.$$ where the drift term $(D \tau)^2$ is subdominant, contributing only $\sim20$% to $\langle (L_z - L_{z0})^2 \rangle $.

Similarly, the variance of the radial action of a population of age $\tau$ is $$\langle (J_R - J_{R0})^2 \rangle = 2 \langle J_R \rangle^2 + 2J_{R0}^2 -2 J_{R0} \langle{J_R}\rangle,$$ with $J_{R0}$ the radial action at birth, which we assume here is zero since in the model stars are born on near-circular orbits. $\langle J_R \rangle = \sigma_R^2/\kappa$ is the mean radial action as defined in Eq. \[eq:p\_radial\_action\]. Both of these quantities are plotted as functions of $\tau$ in Fig. \[fig:secular\_evolution\]. Using a reference age of $\tau_m=6\,\mathrm{Gyr}$, we find at 8 kpc (see Fig. \[fig:secular\_evolution\] for the spatial variations) $$\begin{split}
\sqrt{\langle (L_z - L_{z0})^2 \rangle} &\approx (619~\mathrm{kpc~km/s})
\left( \frac{\tau}{\mathrm{6~Gyr}}\right) ^{0.5}
,\\
\sqrt{\langle (J_R - J_{R0})^2 \rangle} &\approx (63~\mathrm{kpc~km/s})
\left( \frac{\tau}{\mathrm{6~Gyr}}\right) ^{0.6}.
\end{split}$$ We note that up to a factor between 1.2 and $\sqrt{2}$, these quantities are very close to [$\sigma_\mathrm{Lz}$]{}=572 kpc km/s  and the mean radial action $\langle J_R \rangle=\sigma_R^2/\kappa = 45$kpc km/s respectively, so our general conclusions do not depend much on the details of our choice of reference quantities (i.e. [$\sigma_\mathrm{Lz}$]{}versus $\sqrt{(\Delta L_z)^2}$). A spatial representation of the diffusion in angular momentum and increase of radial action is illustrated in Fig. \[fig:model\_aspects\]. We show two examples for the secular evolution of the disk for the first two panels (weak mixing and heating, and strong mixing and heating), and the best fit in the third panel. The second $y$ axis of Fig. \[fig:secular\_evolution\] also illustrates this more quantitatively.

Using the same APOGEE RC dataset, [@frankel_etal_2018; @frankel_etal_2019] measured a migration strength in Galactocentric radius of $3.1~ \mathrm{kpc} \sqrt{\tau / \mathrm{6~Gyr}}$, which is slightly larger than we have found here. These models were purely spatial and ignored the dynamics. The spreads in the metallicity distributions in these models are wholly accounted for by the radius migration and its strength is more simply linked to the radial metallicity gradient (as in Eq. ). In the new dynamical model presented here, spreads in the metallicity distribution are due to a combination of both migration and heating, the latter of which introduces more extreme metallicity stars from the inner and outer disc at their apo- and pericentres respectively, and so further broadens the metallicity distributions. This reduces our measured radial migration strength to $\sqrt{\langle(\Delta L_z)^2\rangle}/v_\mathrm{circ}=2.6~ \mathrm{kpc} \sqrt{\tau / \mathrm{6~Gyr}}$ with the difference coming from the radial heating.

Implications of a Strong $L_z$ Diffusion Process
------------------------------------------------

If the strength of angular momentum diffusion, of the order of the mean angular momentum of the Galaxy, is typical to all disk galaxies, this redistribution has important implications for galactic archaeology for external galaxies. For the Milky Way, the strength of radial migration can be measured through a physical \[Fe/H\]-age-$L_z$ scatter, obtained by data for individual stars. Such a framework may not be applicable in external galaxies beyond the Local Group, where all properties are integrated. Studying stellar populations in external galaxies may lead to good present-day age histograms and present-day mass-weighted age gradients, but reflect only mildly the formation of galaxies due to important dynamical memory loss: i.e. age and metallicity radial gradients weaken [@frankel_etal_2019].

As argued in [@herpich_etal_2017], an asymptotically strong redistribution of stellar angular momenta in cold disks could naturally lead disk profiles to follow exponential distributions, as is observed in disk galaxies [e.g., @deVaucouleur_1948; @freeman_1970ApJ]. Since disks are not always expected to form with exponential profiles [e.g. @roskar_etal_2008a], a strong diffusive process that leads to an exponential profile irrespective of the initial conditions could reconcile the observations with simulations of galaxy formation.

Even though mixing processes are strong, they are not strong enough to erase all gradients in which case even using chemical-age information would not rewind stars back to their birth conditions because the final state of the system would be independent from its initial state. In the Milky Way, metallicity and age radial gradients are weakened, but not erased.

However, any modelling of the Milky Way’s chemical evolution requires us to account for the strong radial orbit redistribution. For instance, a local age distribution might reflect better the global star formation history of the disk rather than the local star formation history, as a local sample of stars, even those on circular orbits, may contain stars born kiloparsecs away and lack stars born locally.

Disentangling $L_z$-Diffusion from Heating
------------------------------------------

Dynamical processes produce correlated changes in the actions of stars. Assuming $\Delta J_R = f(\Delta L_z)$, the general function $f$ will depend on the specifics of the dynamical interaction. In the present work, we do not explicitly model the possible interactions leading to radial heating and angular momentum diffusion; we only measure their effect over 6 Gyr of evolution with an effective model, and find that across the disc $$\sqrt{\langle (J_R - J_{R0})^2 \rangle} \approx 0.1 \sqrt{\langle (L_z - L_{z0})^2 \rangle}.$$ This result can already provide some global constraints on the nature of dynamical processes across the Milky Way, but not on the details of the secular interactions. The net changes are not directly comparable to the changes expected over single migration events $\Delta J_R = f(\Delta L_z)$. To zeroth order, near the main resonances of a non-axisymmetry rotating at a constant pattern speed $\Omega_P$ (e.g. the bar or a spiral wave), the change of radial action is related to the change in angular momentum through [e.g. @sellwood_binney_2002] $$\label{eq:sellwood_binney_delta_Jr}
 \Delta J_R = \frac{\Omega - \Omega_P}{\kappa} \Delta L_z.$$ Here $\kappa$ and $\Omega$ are the radial and azimuthal frequencies. [@sellwood_binney_2002] point out that near corotation ($\Omega = \Omega_P$), $\Delta J_R$ should be very small even though $\Delta L_z$ can be large – a star can move from circular orbit to circular orbit. Therefore, there is no *dynamical* evidence that a star found on a near circular orbit at radius $R$ today was not born on a different circular orbit. Around the Lindblad resonances, where $\kappa = \pm m(\Omega - \Omega_P)$, interactions with non-axisymmetries tend to heat the disk $\Delta J_R=\pm\Delta L_z/m$. More recent works argue that $\Delta J_R = f(\Delta L_z)$ is not necessarily linear, and that angular momentum redistribution at corotation might not always occur without changes in $J_R$: resonances can overlap, leading to non-linear effects and stochastic motions of the stars [e.g @minchev_2011; @minchev_2012; @daniel_2019]. In this model context, our findings suggest that migration near corotation was important.

In addition to these non-linear effects, $\Delta L_z$ and $\Delta J_R$ may deviate from Eq. \[eq:sellwood\_binney\_delta\_Jr\] if spiral perturbations do not rotate as solid bodies as seen in simulations [i.e. with a pattern speed that changes with Galactocentric radius @quillen_2011; @grand_etal_2012] or indirectly with extragalactic observations [@merrifield_2006; @masters_2019], or if their pattern speed is a function of time.

The simulations described above, as well as those of [@brunetti_etal_2011] and [@loebman_etal_2016], have brought understanding of the processes involved in the secular evolution of disk galaxies and the processes at play in radial orbit migration, and have pioneered qualitative comparisons with Milky Way data. However, they are not directly comparable to observed data in the Milky Way because observed data are noisy and do not represent the full Milky Way disk. More importantly, the present view of the Milky Way only represents the equivalent of the final snapshot of a simulation, as argued in [@aumer_2016]. This means age-kinematic relations differ from heating histories (as much as time differs from age) and our effective models cannot be used to recover robustly the full evolutionary history of the Galaxy. Simulations are necessary to guide the construction of realistic and physically motivated forward models and to make the link between the global measure of $\sqrt{\langle (J_R - J_{R0})^2 \rangle}$, $\sqrt{\langle (L_z - L_{z0})^2 \rangle}$ and the instantaneous changes $\Delta J_R$ and $\Delta L_z$.

Implications for the Sun and the Solar System
---------------------------------------------

In a model with significant radial migration, the Sun potentially formed quite far from its present Galactocentric radius. Here we analyse the most likely history of the Sun using our model.

Using Eq. \[eq:eq\_Fe\_Ro\_t\] with the Sun’s age of $\tau_\odot = 4.6$ Gyr [@bonanno_2002_agesun] and ${\ensuremath{\mathrm{[Fe/H]}}}_0 = 0\pm 0.05$ dex [@asplund_2009], we find the birth angular momentum of the Sun was $L_{z0\odot} \approx 1824 \pm 127$ kpc km/s . If, at the time of its formation, the circular velocity corresponding to this angular momentum were 235 kpc km/s, this would correspond to a birth Galactocentric radius of $7.8 \pm 0.6$ kpc, which is 5% closer to the Galactic center than today. This is quite different from previous estimates of the Solar birth location in [@frankel_etal_2018] (5.3 kpc), but in better agreement with [@minchev_2018] (7.3 kpc) and [@haywood_etal_2019]. We interpret this significant change of Solar birth location from our previous estimates as a consequence of two model modifications. Firstly, the introduction of the drift term $D$ towards the inner disk in the diffusion equation \[eq:diffusion\]. In [@frankel_etal_2018], $D$ was set to zero, which resulted in global outwards migration due to the negative density gradient, with a disk profile that broadens with time. Here, we approximately conserve angular momentum, with a disk profile remaining approximately constant over time. As a result, stars have a higher probability to migrate inwards than outwards. Secondly, our chemical enrichment description (Section \[section\_enrichment\]) is different: it is a function of birth angular momentum and not birth Galactocentric radius, and we have imposed a flattening of the \[Fe/H\] profile in the inner disk, which is more physically and observationally motivated.

The Solar birth Galactocentric radius is still widely debated and not well constrained. Most chemical evolution arguments lead to birth radii estimates between its present-day radius and 3 kpc closer to the Galactic center [e.g. @wielen_1977; @nieva_przybilla_2012; @minchev_2018; @sanders_binney_2015; @feltzing_2019; @kubryk_etal_2015a; @frankel_etal_2018], except for [@haywood_etal_2019] who argue that the Sun is a typical outer disk star. Even though different models infer different birth radii for individual stars, the overall radial migration rate estimate remains similar in all models, as shown in [@feltzing_2019]. Finally, [@martinez-barbosa_etal_2015] use backward integration over the Sun’s lifetime, concluding that the Sun was born in the outer disk. However, without knowledge of the past evolution of the Milky Way’s potential, such an exercise is not trivial. Better estimates of the Solar birth place may additionally allow to put tighter constraints on the environment in which the solar system has evolved [e.g. encounters with Giant Molecular clouds @kokaia_davies_2019] .

Application to the Solar Siblings’ Orbit Distributions
------------------------------------------------------

We estimate the possible present-day $L_z$ and $J_R$ ranges occupied by stars that were born with the same ${\ensuremath{ L_{z0}}}$, the same [$\mathrm{[Fe/H]}$]{}, and at the same time as the Sun with Eq. \[eq:p\_radial\_migration\]. Assuming solar siblings undergo phase mixing rapidly (the Sun has undergone $\sim$ 20 Galactic orbits), there is then no dynamically noticeable difference between ‘sharing the same birth cluster as the Sun’ and ‘being born with the same [$ L_{z0}$]{}, and time’. From our model fit, 95% of these stars should currently have $550 \leq L_z \leq 2770$ kpc km/s and $J_R \leq 130$ kpc km/s. This is roughly consistent with the results of [@webb_etal_2019], who used simulations to investigate the present-day positions of solar siblings in the $(L_z, J_R, J_z)$ space in different possible potentials and constrained present-day solar siblings angular momenta to $353 \leq L_z \leq 2110$ kpc km/s and $J_R \leq 116$ kpc km/s. The exact values of these bounds should depend on the detailed history of the Milky Way disk, but their model gives an angular momentum range of about 2000 kpc km/s, which is close to our 2$\sigma(\tau)$ value.

However, if the abundance profile of the gas in the Galactic disk is really axisymmetric, then there is no chemical information on the phase to disentangle whether a star is born in the same birth cluster as the Sun, or just at the same Galactocentric radius (or [$ L_{z0}$]{}). . The recent analysis of [@ness_etal_2019] shows that stellar orbits and abundances can be well predicted with only \[Fe/H\] and age, implying that our present analysis contains the most essential elements for chemical tagging.

Limitations and caveats
-----------------------

The physical limitations of our modeling were discussed extensively while presenting the best fit Milky Way disk model in Subsection \[subsection:model\_limitations\]: the model could improve by allowing a time- and radius-dependent strength of radial migration and radial heating.

We now discuss another approximation we have made while constructing the model: the model for the population selection of the red clump stars. In practice, the red clump selection is based on a neural network trained to classify stellar evolutionary stage from their spectra, trained on asteroseismic data [@ting_hawkins_rix_2018]. In the model, we have approximated the selection of red clump stars as a cut in logg-Teff-color space as in [@bovy_etal_2014]. This approximation is well motivated because (1) the classifier is currently one of the best methods to obtain a pure and complete sample of red clump stars given their spectra [@ting_hawkins_rix_2018], and (2) the modeled cut in logg-Teff-color space is, in theory, a good approximation to selecting core helium burning stars.

However, this approximation is conceptually not satisfying (as for other methods based on data driven selection functions), and with no known applicable and rigorous solution: our treatment of the population selection function in the model is conceptually inconsistent with the actual selection of stars. The actual data-driven selection of stars is not trivial to forward model: this method takes a star’s spectrum as input and returns asteroseismic parameters $\Delta \nu$ and $\Delta P$, thereby disentangling RC and secondary RC. Therefore, to assess selection effects in our forward model properly, we should generate a set of theoretical spectra of various ages, $T_\mathrm{eff}$, $\log g$ and various abundances, add noise and instrumental effects, and then pass these spectra to the neural networks that selected the red clump stars to evaluate the fraction of generated stars that becomes classified as red clump, as a function of age and metallicity (or any stellar quantity that we wish to model). Hence, if for example the neural networks were to fail in some areas of the parameter space (where e.g. the training data are sparse), the effects would be fully modelled.

As machine learning-based selection functions are (justifiably) becoming important in the field of astrophysics, they will also need come with a consistent modeling framework in the future.

Instead, we have approximated this population selection effect as a function of age from [@bovy_etal_2014] as in [@frankel_etal_2018; @frankel_etal_2019] and assumed it was sufficient for the purpose of our work. We argue qualitatively that this should have only a small impact on the present work: such a selection should mainly affect the distribution of ages (at large ages, were the training sample was sparse) of our sample, that is most closely linked to the star formation history in our modeling context. However, we treat this star formation history as a nuisance aspect and marginalize over it, and assume the red clump selection does not bias significantly the dynamical parameters of the stars (beyond the covariances between age and kinematics, which we do model), and we focus on radial migration and diffusion in action space. As can be seen in Fig. \[fig:results\_corner\], the dynamical parameters are not correlated with the star formation history.

Summary
=======

We have presented and applied a global model for the secular evolution of the Milky Way’s low-$\alpha$ disk. The model accounts for an inside-out star formation history, \[Fe/H\] gradual enrichment, and the subsequent evolution of a stars’ orbits as diffusion in action space. Applying this model to the APOGEE red clump stars, we have fully accounted for the selection function of the survey and for data uncertainties. The data are constraining, and the model fit with MCMC implies, in this context, that

1.  $L_z$ redistribution evolves as $\sqrt{\langle (L_z - L_{z0})^2 \rangle} \approx (619~\mathrm{kpc~km/s})\left(\frac{\tau}{\mathrm{6~Gyr}}\right)^{0.5}$, which corresponds to a migration distance of about 2.6 kpc for the 6 Gyr old stars (see Fig. \[fig:secular\_evolution\]). In other words, for a coeval population of stars, 68% of them will be within $ 2.6~\mathrm{kpc}\sqrt{\tau{/\mathrm{6~Gyr}}}$ of their birth radius and the remaining 32% will have migrated further;

2.  $J_R$ evolves as $\sqrt{\langle (J_R - J_{R0})^2 \rangle} \approx$\
    $ (63~\mathrm{kpc~km/s})
    \left( \frac{\tau}{\mathrm{6~Gyr}}\right) ^{0.6}$. This shows that redistribution in angular momentum is stronger than increase of radial action, by a fator $\sim 10$, leading us to conclude that radial migration dominates the evolution of the Galaxy’s low-$\alpha$ disk. This leads the disk to remain kinematically cold, but with a strong dynamical memory loss, making it necessary to use chemical and age information to recover the birth conditions;

3.  the Sun’s birth angular momentum inferred from the best fit is relatively close to its present-day angular momentum with $L_{z0,\odot} \approx 1824 \pm 127$ kpc km$s^{-1}$, but its siblings may have a large distribution in action space with a width of 2000 kpc km/s in $L_z$ and 130 kpc km/s in $J_R$.

We have demonstrated that our approach can disentangle the diversity of dynamical phenomena that have shaped the Milky Way’s disk. However, our parametric model is purely effective so does not capture the real complexities of individual dynamical processes, but only their average effect on the Milky Way. We hope that this may be solved in the future, by applying and coupling this model to more detailed simulations of galaxy evolution.

Acknowledgements {#acknowledgements .unnumbered}
================

It is a pleasure to thank Michele and Georges Laillet for hosting the writing of early drafts of this work. We thank Paul McMillan, James Binney, Gregory Green, Paola Di Matteo and Misha Haywood for interesting discussions.

N.F. acknowledges support from the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). J.L.S acknowledges the support of the Royal Society, and the Leverhulme and Newton Trusts. H.-W.R. received support from the European Research Council under the European Union’s Seventh Framework Programme (FP 7) ERC Grant Agreement n. \[321035\]. YST is grateful to be supported by the NASA Hubble Fellowship grant HST-HF2-51425 awarded by the Space Telescope Science Institute.

The following softwares were used during this research: Astropy [@astropy], Matplotlib [@matplotlib], [*Galpy*]{} [@bovy_2015_galpy], Emcee [@foreman-mackey_2013]. Figure \[fig:results\_corner\] was produced using the package Corner [@corner].

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is <https://www.cosmos.esa.int/gaia>. The Gaia archive website is <https://archives.esac.esa.int/gaia>.

Combining the Model Aspects into a Global PDF \[sec:appendix\_model\]
=====================================================================

Global Milky Way Disk model
---------------------------

We show here how the different model aspects presented in subsections (\[subsection:pLzo\], \[subsection:p\_Lz\], \[section\_enrichment\], \[subsection:p\_heating\], \[subsection:p\_z\]) are combined together to build the overall model for the Milky Way disk. Applying the probabilistic chain rule, and marginalizing over the dummy variable [$ L_{z0}$]{}, the different aspects of the model appear: $$\label{eq_p_model}
\begin{split}
p_
\mathrm{MW}(\tau, {\ensuremath{\mathrm{[Fe/H]}}}, J_R, {\ensuremath{ L_{z}}}, z |{\ensuremath{ \mathbf{p_m}  }}) &= 
 \int p({\ensuremath{ L_{z0}}}| {\ensuremath{ \mathbf{p_m}  }}) p(\tau|{\ensuremath{ L_{z0}}}, {\ensuremath{ \mathbf{p_m}  }})
 p({\ensuremath{\mathrm{[Fe/H]}}}, J_R, {\ensuremath{ L_{z}}}, z |{\ensuremath{R_0}}, \tau, {\ensuremath{ \mathbf{p_m}  }}) \mathrm{d}{\ensuremath{ L_{z0}}}\\
& = \int p({\ensuremath{ L_{z0}}}|{\ensuremath{ \mathbf{p_m}  }})p(\tau|{\ensuremath{ L_{z0}}},{\ensuremath{ \mathbf{p_m}  }})
 p({\ensuremath{\mathrm{[Fe/H]}}}|{\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }})p(z |{\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }})
 p(J_R, {\ensuremath{ L_{z}}}|{\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }}) \mathrm{d}{\ensuremath{ L_{z0}}}.
\end{split}$$ The first term on the right hand side is the stars’ birth angular momentum distribution (the first part of Subsection \[subsection:pLzo\]). The second is the star formation history conditioned on birth angular momentum, resulting from an inside-out star formation history (second part of Subsection \[subsection:pLzo\]). The third term is the distribution of metallicity in the star forming disk in function of time (modeled here as a Dirac function, since we are assuming a tight ${\ensuremath{ L_{z0}}}-\tau-{\ensuremath{\mathrm{[Fe/H]}}}$ relation, Subsection \[section\_enrichment\]), the fourth term is the vertical distribution of stars in the disk, and the last one is the joint distribution of in-plane orbital properties (which we take as the azimuthal action, or angular momentum, and the radial action). It can be split $$\begin{split}
p(J_R,& {\ensuremath{ L_{z}}}~|~{\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }}) =  p(J_R ~|~ {\ensuremath{ L_{z}}}, {\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }})p({\ensuremath{ L_{z}}}~|~{\ensuremath{ L_{z0}}}, \tau,{\ensuremath{ \mathbf{p_m}  }}),
\end{split}$$ where the first part corresponds to radial heating, and is conditioned on both birth angular momentum and present-day angular momentum. Heating through scattering should happen over the entire trajectory of the star, so in some sense at an average of the birth and final angular momenta, but we will drop the dependence on birth angular momentum as an approximation (see Subsection \[subsection:p\_heating\]). The second term corresponds to radial migration, modelled as diffusion in angular momentum. Here it shows the probability of a star to be at angular momentum $L_z$ given it was born at ${\ensuremath{ L_{z0}}}$ a time $\tau$ ago (see Subsection \[subsection:p\_Lz\]). Each of these model aspects are presented in Section \[sec:model\] and assembled together to form the Milky Way model in Section \[sec:likelihood\].

Modeling the Dataset: Noise Model, Selection Function, and the Observables
--------------------------------------------------------------------------

The Milky Way model described above cannot be directly applied to our dataset, since the stars were selected in a given survey and the data are noisy. We therefore write the model for the data set, in the space of the noisy observables (with subscript ’obs’). We (1) marginalize over uncertainties, and (2) apply the selection function to the model: $$\begin{split}
    p&_\mathrm{dataset} = p_\mathrm{dataset}(l, b, D_\mathrm{obs}, v_{X,\mathrm{obs}}, v_{Y,\mathrm{obs}}, {\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{obs},
    \tau_\mathrm{obs} | {\ensuremath{ \mathbf{p_m}  }}, \mathbf{\sigma}) \\
    & = \int p_\mathrm{dataset}(l, b, D_\mathrm{true}, v_{X,\mathrm{true}}, v_{Y,\mathrm{true}}, {\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{true},
    \tau_\mathrm{true}|{\ensuremath{ \mathbf{p_m}  }}) p_\mathrm{noise}(\mathrm{obs}~|~\mathrm{true}, \mathbf{\sigma}) \mathrm{d}^n\mathrm{true}\\
    & = \frac{1}{V_s({\ensuremath{ \mathbf{p_m}  }})} \int p_\mathrm{MW}(l, b, D_\mathrm{true}, v_{X,\mathrm{true}}, v_{Y,\mathrm{true}}, {\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{true},
    \tau_\mathrm{true}|{\ensuremath{ \mathbf{p_m}  }}) S(l, b, D_\mathrm{true})f_{RC}(\tau_\mathrm{true}) p_\mathrm{noise}(\mathrm{obs}~|~\mathrm{true}, \mathbf{\sigma}) \mathrm{d}^n\mathrm{true}\\
    & = \frac{1}{V_s({\ensuremath{ \mathbf{p_m}  }})}\frac{1}{(2\pi)^2}\int p_\mathrm{MW}(L_z, J_R, z,{\ensuremath{\mathrm{[Fe/H]}}}_\mathrm{true},
    \tau_\mathrm{true}|{\ensuremath{ \mathbf{p_m}  }}) D_\mathrm{true}^2 \cos(b)S(l, b, D_\mathrm{true})f_{RC}( \tau_\mathrm{true})p_\mathrm{noise}(\mathrm{obs}~|~\mathrm{true}, \mathbf{\sigma}) \mathrm{d}^n\mathrm{true}.
\end{split}$$ From the first to the second line, we marginalize over data uncertainties with a noise model $p_\mathrm{noise}(\mathrm{obs}~|~\mathrm{true}, \mathbf{\sigma})$ where ‘noise’ denotes all the noisy variables used here (those which have subscript ‘obs’), and the uncertainty parameter array $\sigma$ reflects the uncertainties described in the Section \[section\_data\]. From the second to the third line, we split the dataset model as the product of the Milky Way model extensively described in Section \[sec:model\] and the selection function $S(l, b, D, \tau)$. From the third to the fourth line, we describe the disk in the 2D space of actions ($J_R, L_z$) instead of the 4D phase space ($x, y, v_x, v_y$). We note the slight inconsistency in the actual action calculation based on the Staeckel approximation [@binney_2012] which uses the full 6D phase space information, and our simplified modelling assumption based on the adiabatic distribution functions [@binney_2010]. This assumption should not matter much since we restrict to the thin ($|b|<25$ deg), young (red clump), low-$\alpha$ disk where vertical excursions are very limited. The $D_\mathrm{true}^2 \cos(b)$ term is the Jacobian to change from Galactic to Cartesian coordinates. In practice, we perform this integral by sampling the true values (in phase space and propagating directly to action space) from an approximate noise model $p(\mathrm{true}|\mathrm{obs})$, and then we use importance sampling (weighting the integrand with $p(\mathrm{obs}|\mathrm{true})/p(\mathrm{true}|\mathrm{obs})$) and Monte Carlo integrate by summing the remaining terms over these samples. The $p(\mathrm{obs}|\mathrm{true})$ term is the uncertainty model as described in Section \[section\_data\].

Accounting for the Survey Volume $V_S({\ensuremath{ \mathbf{p_m}  }})$ \[sec:appendix\_Vs\]
===========================================================================================

The survey volume in Eq. \[eq:survey\_volume\] is a 7 dimensional integral (or 5, after having integrated over $(l, b)$ by assuming the distribution function does not vary over an individual APOGEE field). We choose to compute this integral by importance sampling. This method works best if the proposal distribution $p_\mathrm{proposal}$ is similar to the target distribution. Ideally, we would like to generate from our full models using a set of realistic parameters $\mathbf{p_\mathrm{prop}}$. We could then systematically use these samples in a Monte Carlo integration of the survey volume corresponding to a new set of model parameters. The advantage of this method is that our samples approximately trace the best fit model so provide an accurate computation of the normalization with a minimal number of samples. It also provides a tractable way to handle the selection function in the survey volume which is automatically incorporated in our sampling distribution.

However, one cannot sample points directly from our full models as because they are constructed in a complex and un-normalized way. As a workaround, we choose to sample from the full model in two steps, using a simpler Galaxy model as an intermediate distribution. First, we sample stars from the simple proposal distribution that is easy to normalize. Then, we down-sample these data through importance sampling using our proposal model. The simple Galaxy model $p_s$ is an exponential disk of constant scale length $R_{ds}$ and scale-height $h_{zs}$ chosen close to the analogous parameters in the overall model (e.g. $R_{d,\mathrm{prop}}$). Given their positions in the disk, stars velocities are sampled from a Gaussian centered on $(v_R,v_\phi)=(0,v_\mathrm{circ}(R))$ with large standard deviations that envelope the known velocity dispersion.

To generate $N_\mathrm{prop}$ samples from our proposal distribution, we use the following procedure:

1.  For every APOGEE field $\mathrm{i}$, we sample on-sky positions $(l, b)$ using boundaries defined by the selection function. We then sample distances $D$ using the cumulative distribution function of $$\label{eq:appendix_simple_proposal}
        p_s(D~|~l, b) \sim D^2 \cos(b) \exp(-R(l, b, D)/R_{ds}) \mathrm{sech}^2(z(l,b,D)/h_{zs}).$$

2.  We down-sample from these positions using the relative normalization of the simple model in each field. The field that contains the greatest number of stars is not down-sampled, and the other fields are down-sampled by accepting the points with probability $$P_\mathrm{accept}(\mathrm{point}~|~ \mathrm{field~i}) =F_\mathrm{proposal}(l_\mathrm{i}, b_\mathrm{i})/F_\mathrm{max},$$ where $$F_\mathrm{proposal}(l_\mathrm{i}, b_\mathrm{i}) = \int_{D_{\mathrm{min,i}}}^{D_{\mathrm{max,i}}} p_s(D~|~l_\mathrm{i}, b_\mathrm{i})S(l_\mathrm{i}, b_\mathrm{i}, D) \mathrm{d}D$$ and $F_\mathrm{max}=\mathrm{max}(F_\mathrm{proposal}(l_\mathrm{i}, b_\mathrm{i}))$. After this step, the samples have a distribution that follows our simple Galaxy model but now accounts fully for the spatial selection function of APOGEE and the 3D extinction.

3.  We now sample the remaining variables from our simple model before a further downsampling. Ages are sampled from a uniform distribution. Birth angular momenta [$ L_{z0}$]{} are sampled from a normal distribution centered on $R\times 235$ kpc km/s with a standard deviation that increases as $\sigma \sim \sqrt{\tau}$ (to mimic radial migration). The strength $\sigma$ is chosen larger than the analogous parameter in the global model ([$\sigma_\mathrm{Lz}$]{}) such that this distribution envelopes the radial migration model. Radial and azimuthal velocities ($v_R$, $v_\phi$) are sampled from normal distributions centered on $0$ and $v_\mathrm{circ}(R)$ respectively.

4.  We further down-sample these points by accepting them with a probability $$P_\mathrm{accept} =\frac{p_\mathrm{proposal}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z~|~\mathbf{p_\mathrm{prop}})f_{RC}(\tau)}{p_s({\ensuremath{ L_{z0}}}, R, z, v_r, v_\phi)},$$ where $p_\mathrm{proposal}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z~|~\mathbf{p_\mathrm{prop}})$ is our full proposal model with fixed parameters $\mathbf{p_\mathrm{prop}}$. The term $f_{RC}(\tau)$ contains the age-dependent fraction of stars on the red clump evolutionary stage. This down-sampling leads to a set of $N_\mathrm{prop}$ points.

Our $N_\mathrm{prop}$ samples are now drawn from the pdf $S(l,b,D)f_{RC}(\tau) p_\mathrm{proposal}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z~|~\mathbf{p_\mathrm{prop}})$, so we can calculate the survey volume as $$\label{eq:survey_volume_appendix}
    \begin{split}
    V_S({\ensuremath{ \mathbf{p_m}  }}) &= \int _\mathcal{D} p_\mathrm{model}(\mathcal{D} ~|~{\ensuremath{ \mathbf{p_m}  }}) S(l,b,D)f_{RC}(\tau) \mathrm{d}\mathcal{D}    \\    
    & = \int_\mathcal{D} \frac{p_\mathrm{model}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z ~|~{\ensuremath{ \mathbf{p_m}  }}) }{p_\mathrm{proposal}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z~|~\mathbf{p_\mathrm{prop}})} \times S(l,b,D) f_{RC}(\tau) p_\mathrm{proposal}(L_z, J_R, \tau, {\ensuremath{ L_{z0}}}, z~|~\mathbf{p_\mathrm{prop}})  \mathrm{d}\mathcal{D}\\
    & \approx \frac{1}{N_\mathrm{prop}}\sum_i^{N_\mathrm{prop}} \frac{p_\mathrm{model}(L_{z\mathrm{prop}i}, J_{R\mathrm{prop}i}, \tau_{\mathrm{prop}i}, L_{z0\mathrm{prop}i}, z_{\mathrm{prop}i} ~|~{\ensuremath{ \mathbf{p_m}  }}) }{p_\mathrm{proposal}(L_{z\mathrm{prop}i}, J_{R\mathrm{prop}i}, \tau_{i\mathrm{prop}}, L_{z0\mathrm{prop}i}, z_{\mathrm{prop}i}~|~\mathbf{p_\mathrm{prop}})}. 
    \end{split}$$ This integration has several advantages over using a regular grid (which is inefficient as the number of dimensions in the data increases) or re-sampling the normalization sample each time with new parameters. First, we only need to produce Monte Carlo samples once, not each time we need to evaluate the model, which is computationally more efficient. Secondly, since the initially generated samples are fixed, we need not recompute actions from these samples each time the model is evaluated, which saves additional computation time. Thirdly, as highlighted by [@mcmillanbinney2013], for fixed samples the stochastic noise is limited and the overall normalized model is a smooth function of the model parameters which is a desirable property when we want to optimize the likelihood to fit the model parameters).

Step 2 of our procedure (down-sampling from the points generated at step 1) is in principle not mandatory. But in practice, the more alike the two distributions in the integral are, the greater the effective sample size.

We have tested this integration method on mock data, and the results remain robust as long as the proposal distribution is broad enough to envelope the distribution that we want to normalize. Additionally, we have tested the overall optimization scheme (MCMC ran on a model using this integration method) on mock data, and recovered the true parameters largely within the uncertainties.

[^1]: The region $L_z \leq \langle L_z \rangle $ corresponds to the inner 3 kpc of the disk, where we currently do not have data, see the right most panel of Fig. \[fig:distributions\_data\]

[^2]: The selection function is published with a tutorial of its use at <https://github.com/NeigeF/apogee_selection_function>
