Constraining churning and blurring in the Milky Way using large spectroscopic surveys -- an exploratory study

We have investigated the possibilities to quantify how much stars move in the Milky Way stellar disk due to diffuse processes (i.e. so called blurring) and due to influences from spiral arms and the bar (i.e. so called churning). To this end we assume that it is possible to infer the formation radius of a star if we know their elemental abundances and age as well as the metallicity profile of the interstellar medium at the time of the formation of the star. Using this information, coupled with orbital information derived from Gaia DR2 data and radial velocities from large spectroscopic surveys, we show that it is possible to isolate stellar samples such that we can start to quantify how important the role of churning is. We use data from APOGEE DR14, parallaxes from Gaia and stellar ages based on C and N elemental abundances in the stars. In our sample, we find that about half of the stars have experienced some sort of radial migration (based solely on their orbital properties), 10 % have likely have suffered only from churning, whilst a modest 5-7 % of stars have never experienced either churning or blurring making them ideal tracers of the original properties of the cool stellar disk. Our investigation shows that it is possible to put up a framework where we can begin to quantify churning and blurring an important. Important aspects for future work would include to investigate how selection effects should be accounted for.


INTRODUCTION
It has long been understood that a star moving in the Galactic potential will be subject to transient gravitational interactions, such as passing close to another star or a giant molecular cloud or interacting with the density enhancement caused by a spiral arm (e.g., Wielen 1977). Many of these mechanisms have been studied in quite some detail and have allowed us to for example understand the structure in velocity space of the stars in the Hipparcos data as the result partially of interactions with spiral arms and the bar (Dehnen 2000). These interactions generally causes initially circular orbits to over time become less circular -an effect now known as "blurring".
The observation that the Sun is more metal-rich than E-mail: sofia@astro.lu.se (SF) the youngest stars in the solar neighbourhood has invoked the idea that the Sun potentially came from a region interior to the current position of the sun (Wielen et al. 1996). However, no reasonable set of physical interactions have been able to explain how the Sun could have undertaken this journey. Sellwood & Binney (2002) showed that another process, "churning", can move a star from a more or less circular orbit to another circular orbit thus erasing all memory of the past kinematic history of the star -this means that we can not integrate the orbits of stars backwards to figure out where they came from (an example can be found in Martínez-Barbosa et al. 2016).
Since the seminal paper by Sellwood & Binney (2002) several studies have shown that the effects of radial migration (churning and blurring combined) can be substantial (Roškar et al. 2008;Loebman et al. 2011). However, although there is no doubt that the mechanisms resulting in churn-ing and blurring are present a proper quantification of their respective importance in the Milky Way stellar disk is still lacking. This is partially due to a lack of a large enough and reliable enough data-set but also partially due to that we are still exploring efficient ways to study these effects. Frankel et al. (2018) provides a model that aims to quantify the global efficiency of radial migration amongst stars in the Milky Way. Using data from APOGEE DR12 (Alam et al. 2015) they find that with a radial orbit migration efficiency of 3.6±0.1 kpc. In this model the sun might have a formation radius of about 5.2 kpc. A quite substantial distance away from its present position.
In this study we explore a possible way to quantify the fraction of stars that have migrated radially in the Galaxy. We also attempt to derive the fractions of stars that have been blurred and churned, respectively. Our method borrows ideas from Minchev et al. (2018) and Grenon (1987) and is quite simple. What we do is to assume that there exists a model that describes how the radial metallicity gradient in the inter-stellar medium (ISM) in the Milky Way has evolved over time. Then if we know the age and the metallicity of a star it is straightforward to derive the Galacto-centric radius at which the star formed. This then allows us to calculate how far the star has moved radially in our Galaxy from when it formed till today. Combining such information with orbital data allows to further study how stars with different kinematic properties have moved -enabling us to put up a method to constrain both churning as well as blurring.
In this way we are able to identify the individual stars that have migrated and can, for example study the properties of blurred stars and contrast that with the properties of the churned stars. This allows, eventually, for a fine-grained understanding of the underlying stellar populations. It also allows a study of the properties of stars that have not moved from where they formed, providing further constraints on our understanding of how the stellar disk in the Milky Way formed.
This paper is organised as follows: Sect. 2 describes the data-set used in this exploratory study; Sect. 3 explains the method and describe the different sets of ISM profiles we use; in Sect. 4 we derive upper limits on how many stars have radially migrated and how many have been churned or blurred in our sample. We also discuss limitations of the sample (selection effects) and take a look at how cosmological simulations could potentially be used to provide the evolution of the ISM; Sect. 4.5 provides a summarising discussion of our results; Sect. 5 concludes the paper with a summary of our findings.

APOGEE and Gaia data
We use data from APOGEE (SDSS-IV data release 14, Majewski et al. 2017) and Gaia data release 2 (DR2, Gaia Collaboration et al. 2016Collaboration et al. , 2018a. From APOGEE DR14 we select all stars that fullfil the following four criteria • the surface gravity (log g) • the rotation • the overall iron abundance ([Fe/H] 1 ) (iii) Radial velocity could be determined for the star and the error in the radial velocity is < 0.5 km s −1 (iv) The uncertainty in the determination of [Fe/H] is < 0.05 dex Criterium (iii) ensures that we exclude potential binaries from the sample. The selected stars from APOGEE DR14 were cross-matched with Gaia DR2. We further require that the relative uncertainty in the parallax measured by Gaia t is less than 10 %. This cut allows us to do robust and straightforward calculations of the stellar orbits (compare, e.g., discussion in Gaia Collaboration et al. 2018b, about selecting the best stars for studying fine-structure in the observed Herzsprung-Russel diagram).
The combined sample drawn from APOGEE DR14 and Gaia DR2 fulfilling the criteria listed above amount to about 85 000 stars. Figure 1 shows [Mg/Fe] as a function of [Fe/H] for the sample.

Calculations of stellar orbits
Using the astrometric data from Gaia DR2 and radial velocities from APOGEE DR14 we used galpy 2 (Bovy 2015) to calculate velocities and orbital parameters for the stars. In parti such as the eccentricity of the orbit (ecc.).
Following Liu & van de Ven (2012) we also derive L z /L c for a more robust estimate of the circularity of the orbits. Here L z is the angular momentum in theẑ-direction (in cylindrical coordinates) while L c is the angular momentum in thê z-direction the star would have were it on a circular orbit characterised by the same energy as the current orbit. Thus where R c follows from solving where E is the orbital energy and v c is the circular velocity defined as v 2 c = R∂Φ/∂R for z = 0. The resulting distribution of L z /L c are shown in Fig. 2  d). Figure B1 shows a comparison of ecc. and L z /L c for all stars in our sample. Overall the two measures correlate well, but there are deviations. Following Liu & van de Ven (2012) we use the L z /L c for the majority of our investigations, however, we apply a stricter definition of orbital circularity than they do (see Sect. 3.4).

Stellar ages
We need stellar ages for our investigation. However, our stellar sample consists entirely of giant stars (most of them being on or near the red clump) and determining ages is difficult for such stars. We make use of the investigation by Martig et al. (2016) who used carbon and nitrogen to infer the mass of the stars and hence provide the possibility to place the star on the right track in the HR-diagram and derive its age. These ages are good to about 40 %, this means that we can not study ages of individual stars but the ensemble properties should be relatively robust.

Properties of the full sample
The final sample drawn from APOGEE DR14, Gaia DR2 (Gaia Collaboration et al. 2018a), andMartig et al. (2016) with stellar orbits calculated using galpy comprises some 18 000 stars. In this section we describe the overall characteristics of the sample as well as validate that the ages provide a reasonable description of the stars in the stellar disk. Figure 2 show the major characteristics of our sample. Figures 2 a) to d) show how the full sample covers Galactocentric distances from about 6 to 11 kpc and how their Galacto-centric distances relate to their other properties. The stars reach maximum distances above the plane (z max ) of several kpc, but the majority do not reach beyond 1.5 kpc. Metallicities range between -0.5 dex to super-solar and their distribution in Galacto-centric distances does not depend on metallicity. Stellar ages ranges from 0.5 to about 10 Gyr. The over-density of stars beyond the solar orbit have ages in the younger range whilst stars inside the solar orbit have somewhat larger ages. We will come back to these observations. Figure 3 provides a simple validation of the stellar ages. The stars have been divided in to high-and low-α stars, as indicated in Fig. 1. Figure 3 then shows the sample divided into the high-and low-α stars for four different ranges of z max . For each sample we also show the distributions of stellar ages for the high-and low-α stars. We find that the sample of high-[Mg/Fe] stars on average are older than the low-[Mg/Fe] stars for the sub-sample closest to the plane (plots in the top row). This is expected from solar neighbourhood studies of dwarf stars (e.g., Bensby et al. 2014;Fuhrmann et al. 2017) and thus confirms that our sample in all likelihood gives a reasonable description of the properties of Milky Way disk stars. The overall age of the low-α stars changes as a function of height above the Galactic plane such that the stars get older as we move to higher heights. The high-α stars, in contrast, has roughly the same median age and similar age spread at all z max . Thus resulting in an overall age-gradient for the whole sample as we move away from the plane. This is an expected behaviour.
These simple considerations serve as a validation that our ages give a reasonable representation of the stellar populations and we can use them with confidence in our investigation.
We note that that although the high-α stars are subdominant at low z max as expect, the relative number of lowα stars remains high also at large z max . At 1 kpc there are still more low-α stars than high-α stars. The difference is not very large and could be due to selection effects.
It is important to note that the method we use to infer formation radii for stars should only be applied to stars that are likely to have formed in a disc-like structure. It is very likely that so called high-α or thick disk stars have formed at a time when the interstellar medium was highly turbulent and/or in-homogenous (e.g., Bournaud et al. 2009). Only later did the interstellar medium settle enough to allow for disc formation (Kassin et al. 2012) where it is meaningful to study the effects of churning and blurring (compare Frankel et al. 2018). Figure 2 e) to h) show the properties of the sample when we restrict it to stars in the low-α trend as defined in Fig. 1 and with z max < 0.5 kpc.

ANALYSIS
Formation radii were derived in the same way as in Minchev et al. (2018) by assuming a model that describe how the radial metallicity gradient of the ISM changes over time and then simply referring each star to the relevant Galactocentric distance given its [Fe/H] and age. The difference between their work and ours is that we do not aim to derive the evolution of the ISM but assume that it is known. Figure 3 in Minchev et al. (2018) shows how Galacto-centric distances are assigned to stars based on their [Fe/H] and age.
Here we will only use the stars in the low-α sequence as defined in as defined in Fig. 1 (see discussion in Sect. 2.4).

Choice of ISM profiles
For our analysis we need a prescription of how the radial metallicity profile of the ISM varies over time. It is not our intention here to derive our own profiles, nor to test various profiles against observables (e.g., open cluster, O and B stars) but instead we are focusing on the method and exploring its strengths and weaknesses. Nevertheless, it is valuable to make use of a range of radial profiles to explore the method we develop. In the literature it is possible to find a number of radial metallicity profiles that describes how the ISM evolves over time.
We have selected four sets of radial profiles that describe how the ISM evolves. The profiles are shown in Fig. 3.1. Below we briefly summarise the major characteristics of each set, but we refer the reader to the original publications for further details (Minchev et al. 2018;Frankel et al. 2018 also obtained radial metallicity profiles from an on-going cosmological simulation (Agertz et al. in prep.). That analysis can be found in Sect. 4.4.

Minchev et al. (2018)
As described in their paper Minchev et al. (2018) derived the ISM profiles by forcing a (small) solar neighbourhood sample of stars to reproduce the distributions of formation radii for stars with different ages from the models by Minchev et al. (2013). Characteristic for these radial gradients is the steepening slope of the lines for older ages and the relatively narrow range of [Fe/H] at the smallest R gal . The radial gradients implemented in our study are shown in Fig.3.1 a).

Frankel et al. (2018)
Frankel et al. (2018) present a model that parametrises the star formation history, the chemical enrichment history and profiles that show where in the Galaxy a star forms given its age and metallicity. They assume that it is possible to describe the metallicity profile of the star forming gas as a product of a radial profile and a term describing how the chemical enrichment evolves over time. They then combine these with a model for stellar migration modelled as diffusive processes. This results in a model that can quantify the global efficiency of radial migration. For our purposes we essentially invert their prescription of where a star forms given its age and metallicity (see Appendix A1). The result is radial gradients that define the evolution of the ISM with time suitable for our purposes.
They are shown in Fig.3.1 b). The gradients are straight lines that spread out more and more for the oldest ages. This is similar to one of the rejected tests by Minchev et al. (2018).

Sanders & Binney (2015)
Sanders & Binney (2015) set up a model of the Galaxy based on analytic distribution functions. Of interest to us is their function that describes the relation between age and metallicity for each radius in the Galaxy. The relevant function is derived by fitting to the model results from Schönrich & Binney (2009), which includes full chemical evolution as well as gas accretion and outflow.
Here we invert their prescription to obtain a description of how the radial gradient of the ISM evolves with time (see Appendix A2). The result is shown in Fig.3.1 c). We note that for young ages there is hardly any differentiation at all -all gas share the same radial distribution at all times up and until 4 Gyr and for the oldest age (12 Gyr) the relation is flat at -1 dex since the model assumes that this is the metallicity of the ISM at the formation of the Galaxy (see Table 3  rates and yields. They use parametrised time-and radiusdependent diffusion coefficients to describe the radial migration. The parameterisation is based on N-body + SPH simulations. For further details see their paper. Although they do not provide a tabulated description of how the metallicity gradient in the ISM evolves with time, their Fig. 4 (bottom middle panel) show their results for three different ages. We reproduce these lines in Fig.3.1 d) and use them for our study. Because this model explicitly takes into account the results from the the N-body + SPH simulations the gradients have less idealised shapes. We use the three available lines and interpolate between them to find the formation radii of the stars. This is a simplification but one we deem reasonable as the lines appear relatively well-behaved.

Migratory distance -definition
In order to be able to quantitatively compare and analyse the results from the different model descriptions of how the radial metallicity gradient in the ISM has varied we define the concept of migratory distance. This measure is simply the difference between the star's current Galacto-centric distance and the Galacto-centric distance at its formation: (3)

Error on inferred formation radii
The principle to infer the formation radius of a star is simple, but we also need to consider the uncertainties in the properties used to derive the formation radius. In particular, we need to consider uncertainties in [Fe  radii we find that an error of 5% being the maximum of our derived distances taking parallax errors into account. For most stars the error is between 1 and 3%. We include these errors in our analysis in the following way. Iterating through each star in the full data set, we create a new star by assuming a normal error distribution and drawing new parameters for age, Galactic radius, and metallicity from a Gaussian where the mean is the original observed value and the standard deviation the associated uncertainty for each variable. We do this until we have generated 10 000 variants of each star and therefore 10 000 new and unique data sets.
We then proceed to analyse the data just as we would have done for the original sample but now with a sample where statistical uncertainties can be readily estimated. Figure 5 shows the resulting formation radii for stars of different ages for the four different sets of radial gradients. We note that Minchev et al. (2018), Frankel et al. (2018), and Kubryk et al. (2015a) show the same overall behaviour where older stars are on average formed further in towards the Galactic centre. This is largely in accordance with the criteria used in Minchev et al. (2018), i.e. an in-side out formation of the disk where older stars form in the inner parts of the disc. Using the radial gradients from Sanders & Binney (2015) on the other hand results in that the majority of the stars in our sample forming between roughly 5 and 8 kpc. Figure 6 shows the migratory distances as a function of age and formation radius of the stars. To create this figure we have used the approach to include error-estimates described in Sect. 3.3. Here we show the results for Frankel et al. (2018) but all four sets of radial gradients qualitatively show the same results.

Results
We also analyse a sub-sample of stars on highly circular  orbits. After some experimentation, we find that L z /L c > 0.99 is the best definition of a highly circular orbit for our sample. This is a somewhat arbitrary choice but is supported by Fig. B1. Figure 2 shows that such stars are present across the full range of Galacto-centric distances and metallicities. We believe that the choice of this cut does not significantly bias our investigation.
For all four age bins we find that stars on highly circular orbits (L z /L c > 0.99) have a smaller spread in migratory distance, on the order of 0.2 kpc. The median distance that a star has migrated depends on where it has formed -stars forming in the inner part of the disk have migrated significantly further distances than stars formed close to the sun or further out in the disk.

ANALYSIS, INTERPRETATION, AND DISCUSSION
In this section we show how combining the migratory distances with information about the stars' orbits can be used to constrain how many of the stars in a sample have, e.g., been churned. We discuss different ways of defining if a star has been churned or not and also look at stars that have not migrated at all. We also look at the formation radius of the sun (Sect. 4.3) and investigate the possibility to obtain the radial metallicity gradients for the ISM from a cosmological simulation (Sect. 4.4).

Stars that have migrated
In this section we will consider different ways to decide if a star or a stellar population has migrated in the Galaxy. We will provide some examples of how one can estimate the fractions of stars that have migrated and look at the fraction of stars that have only been churned. Figure 7 provides a summary of the results.

Estimating how many stars have radially migrated
The simplest definition of a star that has radially migrated in the Galaxy is a star for which its formation radius is not the radius the star formed at. As can be inferred from Sect. 3.4, the fraction of stars that have moved radially in the Galaxy is large. Adding orbital information allows for a more interesting analysis.
Here we consider a slightly more involved criterion for stars that have migrated in the Milky Way disk; namely stars that have a Galacto-centric formation radius that is outside the Galacto-centric radial range spanned by the apo-centre and the peri-centre of the current orbit of the star: Since there is no requirement set on the shape of the orbit this implies that the sample of such stars must include both stars that have been only churned as well as stars that have been both churned and blurred. Such stars are rather common, with about half of the stars in our sample fulfilling this criterium. The results are shown in Fig. 7, where we also show the results from the sections below. For the models by Minchev et al. (2018) and Kubryk et al. (2015a) we obtain a lower fraction of radially migrated stars (about 0.5) while for the other three models (Frankel et al. 2018;Sanders & Binney 2015, and, Agertz et al. in prep) we obtain a higher fraction.

Estimating how many stars have been churned
In this section we explore ways to estimate how many of our stars that have experienced churning (churning being the radial migration that causes a star to move radially without loosing the circularity of its orbit Sellwood & Binney 2002, see also Sect. 4.5). For this we need to define the sub-sample of stars that we think should have been subjected to just churning and not blurring. A star that has been blurred can also have been churned. This means that what we are trying to do here is to find a conservative lowest fraction of stars that have just been churned. We define a star to be a candidate as a churned star if it has a relatively circular orbit. We define such stars as those with L z /L c > 0.99.
We consider M D (migratory distance) as a function of stellar age and formation radius, compare Fig. 6. Table 1 summaries our results. The top part of the table lists first the fraction of stars with highly circular orbits defined as L z /L c > 0.99 (results for stars with L z /L c > 0.95 can be found in Appendix, B3). We find that younger stars have a larger fraction of stars on circular orbits. This fraction decreases from about 0.18 to 0.12 as the stellar age increases. An average over all age bins gives a fraction of about 0.15.
Although this is an interesting number we do not think that this gives a minimum fraction of stars on churned orbits. To be churned we further require that the the star's formation radius lies outside of the range of orbits their presentday orbit occupies. As expected this selection gives a smaller fraction of stars. Again, the fraction decreases with increasing age. The total fraction of stars that these criteria and hence are prime candidates for having been churned is about 0.1 (see lower part of Table 1).
The results are shown in Fig. 7 where we also include the same estimates but with a more relaxed criterion on circularity (0.95 instead of 0.99). With the more relaxed criterion on circularity for the orbits the minimum fraction of churned orbits increase to about 0.4. However, for reasons discussed below we do not believe that this gives an indication of the fraction of churned stars.

Estimating the size of churning and blurring
It is also interesting to attempt to estimate the size of churning and blurring. An analysis as the one provided in Fig. 6 allows us to do this. It is readily seen from this figure that the spread in M D for all stars is about 0.2 kpc larger than for stars on highly circular orbits. If we regard the stars on the highly circular orbits as essentially uninfluenced by blurring then this gap is an indication of the size of the increased orbital spread due to blurring. We note that the difference between the full sample and the sample on highly-circular orbits is the same for all formation radii for a given age. From our data we can not say if the difference between samples remains constant as age increases or not. As both churning and blurring can be modelled as diffusive processes it is possible that the difference remains the same over time. All five models investigated in this study show the same patterns.
If we instead consider the x-axis of the plot (M D) we find that stars in our sample that formed inside the solar circle have migrated on average more than those that formed outside the solar circle, with a monotonic change in M D as a function of formation radius. We also note a trend with age where older stars have a larger M D for the same formation radius as compared to younger stars. This is consistent with radial migration being a diffusive process. Although the youngest stars have smaller M D they still show substantial radial migration, indicating that churning must be a process that acts quickly on a stellar sample. We know it can not be blurring as blurring is a slow process and because even the circular sample shows significant migratory distances.
If we combine this finding with the finding that a minimum of about a tenth of the stars have been exclusively churned it indicates that about 40% or so of the stars in our sample have migrated either by being just blurred or blurred and churned.

Stars that have not moved
With our methodology it is possible to identify stars that have not moved radially in the Galaxy since they formed. We define such stars as those with |M D| < 200 pc. Table 2 summarizes the fraction of these stars for the different models. We find that between 5 and 7% of the stars in our sample have not moved radially in the Galaxy since they formed. The number decreases as the age of the stars increases for all models apart from those from Sanders & Binney (2015) where the fraction remains constant as a function of stellar age. This is consistent with what we saw in Sect. 3.4 and Fig. 5 where we show that this model produced roughly the same amount (and spread) in radial migration for all ages apart from the oldest stars. It appears natural that populations of older stars should have a smaller fraction of stars that have not moved radially in the Galaxy since they formed since over time a star's position in the Galaxy will be influenced by several phenomena, not the least blurring and churning. As time goes a star become more likely to have suffered from such phenomena and hence its orbit will start making excursions away from the original circular orbit. Figure 8 shows the properties of the sample of stars that have not moved. We note that at all radii there is a large spread in metallicity as well as age. Metallicity and age appear to correlate well at each radius, such that younger ages are associated with higher metallicities. In a slightly circular argument this could be taken as proof that indeed at a given radius there is a tight age-metallicity relation. When we look at all stars in our sample this is not the case, indeed, many studies of stars in the solar neighbourhood have shown that there is an acute lack of such a relation (e.g., Edvardsson et al. 1993;Feltzing et al. 2001;Casagrande et al. 2011). It should be noted for our method, that although per construction more metal-rich stars are younger at a given Galacto-centric radius they need not be on orbits that have not moved via churning and/or blurring. It would have been entirely possibly that there were only stars of one age or one metallicity at a given radius that were still on the same orbit they had when they formed. Thus, we believe that this does give observational support to the assumption that stars at a given radius in the Galaxy follow a tight age-metallicity relation.

The Galacto-centric formation radius of the Sun
Returning to the question wether or not the Sun has formed at the solar radius we find that using the radial metallicity gradients for the interstellar medium from Frankel et al. (2018), Sanders & Binney (2015), and Kubryk et al. (2015a) the sun formed at a distance from the Galactic centre of 5.7, 7.0, and 6.8 kpc, respectively. Clearly, there are uncertainties associated with these estimates, however, it is also clear that if we require an inside out-formation scenario in which the metallicity in the interstellar medium is enriched over time in such a fashion as to produce flatter and flatter radial gradients then the sun most likely formed between 5.5 and 7 kpc from the Galactic centre. This is largely consistent with other estimates of the sun's formation radius (Wielen et al. 1996;Minchev et al. 2013;Frankel et al. 2018). Haywood et al. (2019) discusses the possibility that the sun instead of coming from inside its current solar position in the Galaxy, it in fact has migrated in from the regions outside the sun's position. They find that this is supported by the numerical experiments carried out by Martínez-Barbosa et al. (2016), who uses backward integration in the Galactic potential to find out where a star comes from. As discussed earlier, with the inclusion of churning such integrations looses their ability to make such predictions, i.e. once a star has been churned you can not any more find out where it came from. As we have shown not all stars have been churned, some may just be blurred, others are untouched by dynamical encounters. So in theory the sun may have an un-churned orbit and you can retrace its orbit, however, we note that the fraction of stars in our sample fulfilling this criterium is small (compare Fig. 7). . Data based on the simulations of Agertz et al. (in prep.). a) Radial metallicity gradients in the ISM as a function of time. b) Distribution of formation radii derived using the gradients from the simulations by Agertz et al. (in prep.). See Sect. 4.4 for further details. Note that the simulation does not (yet) go to age zero.

Agertz et al. (in prep)
Agertz et al. (in prep) carried out a cosmological hydrodynamic+N-body zoom-in simulation of a Milky Waymass galaxy (M ∼ 6 × 10 10 M ) forming in a dark matter halo with virial mass M vir = 1.3 × 10 12 M at z = 0 . The simulation was carried out using the adaptive mesh refinement code RAMSES (Teyssier 2002), assuming a flat Λ-cold dark matter cosmology. We refer to Agertz et al. (in prep) for details, as well as Pehlivan Rhodin et al. (2019) for an extensive description of the included physics and simulation settings. The simulation reaches state-of-the-art resolution, with a mean resolution of ∼ 20 pc in the cold interstellar medium. For every simulation snapshot we identify the most massive progenitor to the z = 0 Milky Way-mass galaxy. We identify the disc plane from the stellar and gaseous angular momentum vector, and compute radially averaged [Fe/H]profiles from all neutral gas within a 2 kpc thick slab. The resulting radial ISM profiles are shown in Fig.9 a).
We note that insight drawn from abundance gradients in cosmological simulations of galaxy formation must be considered with care. While simulated galaxies can be selected to represent extended discs with global properties in line with the Milky Way's, their detailed assembly histories will not necessarily be compatible with that of the Milky Way galaxy. However, for this explorative study we find it informative to include these models as a part of our suite of diverse ISM models.
The radial gradients in Fig.9 a) share many of the overall characteristics of those derived in Minchev et al. (2018) and those used in Frankel et al. (2018) in as much as they steepen towards the inner galaxy and older profiles have lower iron content. Thus, we can expect that overall this description of the temporal evolution of the radial metallicity profile of the ISM should give rise to similar results as found for models explored previously. Indeed, that is also what we see, but with some modifications. Notably, Fig.9  b) shows that although the median radius moves to smaller and smaller radii as age increases (i.e. an inside out formation scenario as imposed in Minchev et al. 2018) there is significant spread at all ages. Also for this model we have derived the fractions of stars fulfilling the different criteria discussed in Sects. 4.1 and 4.2. We show these fractions in Fig. 7 together with the results from the other models. We note that the results from this model is largely the same as for the other models.

Discussion
One of the leading takeaways from our experimentation with this data-set is the apparently rapid onset of radial migration. We find that the median migration distances of the stars in each of our radial bins remain very constant in time. This extends out to our oldest age bin, whose stars generally show identical median migration distances to our youngest age bin. This implies that the bulk of radial migration must take place relatively early in a starâȂŹs lifetime. We suggest that this could be attributed primarily to the fact that the period of the spiral bar pattern, dominantly responsible for churning, is significantly shorter than the width of our youngest age bin. A considerable number of attempts have been made to constrain the pattern speed of the spiral arms using hydrodynamic simulations (see for example Chakrabarty 2007;Bissantz et al. 2003) and observations of nearby velocity fields (Fernández et al. 2001). Most have placed the period of the spiral arms between 250 and 350 Myr (Gerhard 2011). Thus, a spiral arm crossing is well-contained within the interval covered by our first age bin, and given that only a single transient interaction with a spiral arm is required to generate substantial changes to a stars angular momentum (Sellwood & Binney 2002), this should be sufficient to generate the displacements we observe.
This interpretation is also consistent with numerical Nbody simulations which find substantial migration distances within 1 Gyr due to churning. Indeed, the time-scale for churning to displace the stars can be as short as 0.5 Gyr (Grand & Kawata 2016). It has also been observed that the distribution of stellar radii occupied for stars formed in the same location spreads dramatically in the first few Gyr, and slows down considerably at later times (Kubryk et al. 2013). Our data provide empirical confirmation of these results.
Furthermore, the swift onset of radial migration, and subsequent invariability of median motions with time, indicates that churning is the dominant mechanism through which stars initially change their Galactic radii. It has long been understood that blurring alone is not sufficient to explain the distribution of stars observed in the local region and it is now understood that that churning plays an integral part in explaining the locally observed distributions (examples can be found in Schönrich & Binney 2009;Minchev et al. 2013). Though the relative strengths of churning and blurring, that is, how much of a stars displacement from the location it formed can be described by one process or the other, has not been much studied. This is due fundamentally to the fact that these processes act simultaneously, and are largely inseparable in observation.
It has been one of our aims in this work to investigate if it is possible to constrain the relative strengths of the two processes. Our ability to compare the migration distances of stars on circular orbits to the total sample enables us to observe how well the motions in our total sample are explained by a population that has not yet been blurred. In doing this, we see that churning immediately has a large effect on stellar radii, while blurring accounts for a comparatively small increase in spread at this time. This result is seen in the relatively small difference in spread between our total sample and the circular subset when compared to the overall spreads of these populations.
Blurring is traditionally modelled as a diffusive process, gradually increasing the spread of the radial distribution with time (Sellwood & Binney 2002;Schönrich & Binney 2009). As such, the deviations in radius are smaller and symmetric about the mean. In our data, this is represented by the equivalent median, though increased spread, in our total sample when compared to the population of stars on circular orbits. The difference in spread between these populations is presumably indicative of the added diffusion from blurring. Meanwhile, the expected changes to angular momentum caused by churning are quite large -several kpc from a single interaction with the spiral (Sellwood & Binney 2002) -when compared to the scales of the blurring found in our data: a spread only 0.2 kpc separating our total sample and circular subset. This leads us to believe that churning is a much stronger force than blurring in terms of the magnitude of displacement. This is coupled with the observation that our data increases its spread with time, though only slightly from age bin to subsequent age bin. We have noted that this consistent dispersion in time is present in equal magnitude for both our total sample and our highly-circular sub-sample. This implies that the aggregate spreading effect of churning and blurring in our total sample is on the same order as the effect purely from churning in our circular sample. Thus, beyond the initially large displacements caused by churning in the first hundreds of Myrs, churning begins to function indistinguishably from blurring, spreading only gradually with time. However, this is not to say that the distances covered by churning are less at later times. In fact, stars may be churned back and forth across the spiral patterns co-rotation radius many times in their lifetimes on so-called horseshoe orbits (Sellwood & Binney 2002). But this movement back and forth contributes little to any net migration when considering the average of the population (Halle et al. 2018), and serves only to broaden the distribution of present-day radii in a diffusive pattern reminiscent of blurring. Our results point to a net migration for a very large fraction of stars.
Previous attempts to constrain the strength of churning using observations from large spectroscopic surveys include Kordopatis et al. (2015); Hayden et al. (2018Hayden et al. ( , 2019. These studies have mainly identified potentially migrated stars via the eccentricities of their orbits. Kordopatis et al. (2015) found that for stars with super-solar metallicities about half of the stars in their RAVE sample had ecc. < 0.15 (which they took to imply a circular orbit). That such stars are present at the solar radius is interpreted as evidence for churning as the ISM is at solar metallicity today and hence, more metal-rich stars need to come from somewhere else than the solar neighbourhood. Hayden et al. (2018) used stars from the Gaia-ESO Survey. They selected stars with [Fe/H]>0.1 dex and found that about 20% of the stars in their sample do not reach the Galacto-centric radius at which they likely formed and can thus have been churned. Hayden et al. (2019) on the other hand, using a similar approach, find that for stars in the solar neighbourhood observed with GALAH and with ecc. < 0.2 as many as 70% has reached that location thanks to churning/radial migration. Compared with our more modest total conservative estimate of about 15% for our full sample this seems a very large number. On the other hand, Hayden et al. (2019) includes both high-and low-α stars, whilst we (in agreement with Frankel et al. 2018) argue that any attempts to constrain churning for α-enhanced stars is likely to fail due to the complex nature and formation channel of that stellar component which includes both accretion and turbulent formation scenarios in the early Universe (some examples include Bournaud et al. 2009;Agertz et al. 2009). We note that our definition of a circular orbit is much more stringent than used in these studies. Referring to Fig. B1 we can see that by instead using a cut in ecc. we would indeed have many more stars, but importantly, say if we cut at 0.15 then that would essentially include a very large portion of the stars with 0.95. We note that if we relax the criterion for circularity to 0.95 (from 0.99) then the difference between all stars and stars on circular orbits in Fig. 6 disappears. Meaning that then there is no difference between the two samples in terms of displacement, indeed it is not possible to distinguish between churning and blurring (see also discussion in Sect. 4). Hence, we conclude that although ecc. might appear as an easily understood measure of the orbital shape the more robust L z /L c is a better indicator of the orbital shape.
At this point in the discussion it is important to reiterate two things: 1) our study is of an exploratory nature, we wish to see if we can find means to constrain churning and blurring the stellar populations in the cold stellar disk in the Milky Way, 2) the sample we have used is far from perfect for the purposes. We believe that we have succeeded in showing that there are ways to effectively constrain the size and strength of churning using simple means to estimate the migratory distances for stars in the cold stellar disk. In this work we have made no attempt to account for the selection effects in our sample. Could some of our conclusions be influenced by this? Our main aim with this work is to establish ways to constrain the strength of churning and blurring. The results depends both on the models used (the radial metallicity gradients for the ISM) as well as the quality of stellar sample, including its physical distribution in the Galaxy. An inherent assumption is that there are no azimuthal changes in the stellar population in the Milky Way. This is likely incorrect, but current data does not allow to address this question. Future data releases from Gaia coupled with large spectroscopic surveys as well as dedicated follow-up of, e.g., Cepheids and A and B stars across the Milky Way disk will elevate this problem. Figure 2 f) and g) show that the stellar sample we use mainly is situated just outside the solar circle. Although there are stars between 5 -11 kpc there is a concentration around 9 kpc. There is also a slight trend between [Fe/H] and present day Galacto-centric radius such that more metalpoor stars are found further out in the Galaxy. We know that the Milky Way stellar disk is well populated also inside the solar circle but we know less about the metallicities of those stars. Taking the data at face-value, we can thus conclude that it is likely that we are missing stars on smaller radii meaning that we will not have such a good view of the radial migration experienced by the stars that are currently inside the solar circle. We think, however, that for the stars beyond the solar circle and inside about 10.5 kpc we have a pretty good sampling of the stellar population as of today. Thus our inferences about the strength of churning and blurring, as applied to these stars, should hold -churning is the stronger process and acts early on in the life of the stars. Later churning and blurring both acts as diffusive processes that grow slowly over time.
Future studies must still look in to how the selection function of the stellar sample influences the results. We have seen that for our sample of cold disk stars drawn from a combination of Gaia and APOGEE the results are largely model independent. That might not necessarily be the case with a sample defined in a different way and with a different selection function.

CONCLUSIONS
In this work we set out to explore the possibility to quantify how many stars have radially migrated in the stellar disk and, eventually, be able to put numbers on the strength and importance of the processes involved. We have taken some first steps on this path by utilising a sample of red giant stars from APOGEE DR14, Gaia parallaxes and proper motions, and stellar ages derived from C and N abundances in the stars. This sample has allowed us to quantify how large a fraction of the stars in the sample have had their orbits changed from initially circular to non-circular, how many re-main on circular orbits and how many are on circular orbits that might have been "churned".
We find that a conservative estimate is that about 10% of the stars in the sample have been churned. This is in contrast to recent studies that have much higher numbers, however, we note that those studies essentially look at stars with super-solar metallicities whilst we study stars of all metallicities. Furthermore, our definition of a highly circular orbit is deliberately conservative. If we instead select the stars have Galacto-centric radii at formation that lays outside of its apo-as well as peri-centre today, we find that about half of the stars have undergone some combination of churning and blurring. We estimate that a robust 5 − 7% of stars in our sample have not had their orbits blurred, nor have they been churned. These stars appear at all ages indicating that an individual star may escape these dynamical processes for quite a long time. Our study also provides tentative observational support to the assumption that stars at a given radius in the Galaxy follow a tight age-metallicity relation.
Looking towards the future, there are several on-going and upcoming large spectroscopic surveys that would be able to provide data to further explore the relative importance and strengths of churning and blurring in our Galaxy (e.g., WEAVE, 4MOST Balcells et al. 2010;de Jong et al. 2016, 2019, andJin et al. in prep.). From Gaia we will have the parallaxes and proper motions, which when combined with the radial velocities from the spectroscopic surveys will give the full 6D phase space information needed to calculate stellar orbits. The surveys will also provide the needed metallicity and elemental abundances. Stellar ages are difficult to derive. The best prospects are for turn-off stars (Sahlholdt et al. 2019) but to reach large volumes of the Milky Way require us to use red giants. In this study we have made use of stellar ages for red giants derived using elemental abundances of C and N and combinations thereof (Martig et al. 2016). It is essential that these and similar methods are further evaluated, developed and validated such that they may be used, at least in a statistical sense, in studies like the one presented here.   as compared to −0.07 dex kpc −1 in most other models (e.g., Minchev et al. 2018;Frankel et al. 2018).
Manipulating the equations in Sanders & Binney (2015) we arrive at the following equation that describes how the iron content in the ISM changes with radius (R) for a given age (τ):  Figure A2 b then shows the resulting radial gradients used in our work. For full references and discussion of the constants in Eq.(A2 and A3) we refer the reader to Sanders & Binney (2015). We note that as opposed to the other models used in our study this one hardly has any evolution of the ISM radial gradient with time for the time span that is of interest for the formation of the stars in the cold stellar disk, i.e. the last 6 to 8 Gyrs.

APPENDIX B: ORBITAL DATA
As described in Sect. 2.2 orbital data for the stars in our sample were calculated using the galpy-package (Bovy 2015). B1 Comparing L z /L c and ecc.
In this work, following Liu & van de Ven (2012), we have chosen to use L z /L c to characterise the circularity of the stellar orbits. Figure B1 shows a comparison L z /L c and ecc.. There is a clear correlation between the two properties. At a given L z /L c there is a substantial spread in ecc.. We note that our chosen cut for defining very circular orbits is 0.99 for L z /L c , this encompasses ecc. in the range 0 to 0.2.

B2 Including a bar in the Galactic potential
The galpy-package allows the user to include a bar when carrying out the orbital integrations (the reader is referred to Bovy 2015, for details about the potential) . Figure B2 shows the effect on ecc. when the bar is included in the potential. Although there are noticeable effects we note that for stars that we consider to be on circular orbits (which all have < 0.2, see Sect. B1) the effect is small and hence we have concluded that we did not need to include a bar in the Galactic potential for the purpose of this study.  B3 L z /L c < 0.95 For completeness we include here the resulting fraction of stars on circular orbits when the constraint has been reduced to 0.95 instead of the 0.99 we use in the final analysis (see Table 1 for the 0.99 results). The data are given in Table B1 and are also shown in Fig. 7.