Dynamical traceback age of the $\beta$ Pictoris moving group

Context: The $\beta$ Pictoris moving group is one of the most well-known young associations in the solar neighbourhood and several members are known to host circumstellar discs, planets, and comets. Measuring its age with precision is basic to study several astrophysical processes such as planet formation and disc evolution which are strongly age dependent. Aims: We aim to determine a precise and accurate dynamical traceback age for the $\beta$ Pictoris moving group. Methods: Our sample combines the extremely precise Gaia DR2 astrometry with ground-based radial velocities measured in an homogeneous manner. We use an updated version of our algorithm to determine dynamical ages. The new approach takes into account a robust estimate of the spatial and kinematic covariance matrices of the association to improve the sample selection process and to perform the traceback analysis. Results: We estimate a dynamical age of $18.5_{-2.4}^{+2.0}$ Myr for the $\beta$ Pictoris moving group. We investigated the spatial substructure of the association at birth time and we propose the existence of a core of stars more concentrated. We also provide precise radial velocity measurements for 81 members of $\beta$ Pic, including ten stars with the first determination of their radial velocities. Conclusions: Our dynamical traceback age is three times more precise than previous traceback age estimates and, more important, for the first time, reconciles the traceback age with the most recent estimates of other dynamical, lithium depletion boundary, and isochronal ages. This has been possible thanks to the excellent astrometric and spectroscopic precisions, the homogeneity of our sample, and the detailed analysis of binaries and membership.


Introduction
Young local associations and moving groups are fundamental structures to understand the stellar formation and evolution processes. They are small aggregates of stars (few dozens) sharing their dynamical properties. For this reason, it is implicitly assumed that they were born at the same time and at the same place (from the same molecular cloud), and therefore, they share the same chemical composition (de Zeeuw et al. 1999;Jayawardhana 2000). Most of the known associations are very nearby and allow a detailed study of their properties.
One of the most well-known associations is β Pictoris (β Pic). It was discovered a couple decades ago when Barrado y Navascués et al. (1999) identified the first two companions to the β Pic star and Zuckerman et al. (2001) identified an additional set of 17 co-moving stars. Since then, many studies have contributed to increase the number of members of this association (e. g. Torres et al. 2006, Malo et al. 2013, Binks et al. 2015, and Riedel et al. 2017b. Nowadays, there are a few hundreds of candidate members of the β Pic moving group, making it one of the richest associations. Its proximity (∼ 40 pc), and observational characteristics (it is visible both from the Southern and Northern hemispheres) facilitated the discovery of members with a large diversity of stellar masses and very interesting properties, such as discs, confirmed exoplanets, and exocomets (Kalas & Jewitt 1995, Kalas et al. 2004, Lagrange et al. 2010, Lagrange et al. 2019, Chauvin et al. 2012, and Kiefer et al. 2014. The age is one of the most fundamental parameters to study stellar formation and evolution. β Pic has an estimated age of ∼20 Myr (Barrado y Navascués et al. 1999, Barrado y Navascúes 2001, Mamajek & Bell 2014, Binks & Jeffries 2014 which is of particular interest to study several astrophysical processes such as disc evolution and planet formation. However, different methods lead to a relatively broad range of values and errors ranging from 10.8 Myr to 40 Myr (see Table 6 for a review of the literature age estimates of β Pic).
Among the few techniques available to determine stellar ages, dynamical ages 1 have the advantage that they are independent of stellar evolutionary models. The main assumption of this method is that the stars were formed together, in the past, at a time when the association was most concentrated. This assumption is supported by the lithium and isochronal ages where there is no evidence of a significant age spread (Mamajek & Bell 2014;Messina et al. 2016). Several authors in the literature used different techniques to traceback the positions and motions of the stars (linear trajectories, epicyclic approximation, orbital integration with a Galactic potential) and different definitions of the size of the association (e.g. standard deviation of the positions in a privileged direction, in 3D, maximum distance between members, pairwise encounters). Historically, the main limitations of the traceback analysis were the observational uncertainties in proper motions and the lack of trigonometric parallaxes and radial velocities to derive distances and spatial velocities (Ortega et al. 2002(Ortega et al. , 2004Song et al. 2003). After the Gaia Data Release 2 (DR2, Gaia Collaboration et al. 2016, 2018a we have a large, uniform sample of stars with extremely precise parallax and proper motions. Additionally, several authors measured radial velocities of β Pic members (e.g. Torres et al. 2006;Shkolnik et al. 2012; although these are highly inhomogeneous and systematic errors may exist between different studies. Currently, the main limitations of the traceback analysis are: 1) the availability of an homogeneous and precise dataset of radial velocities; 2) the design of a new strategy for the selection of kinematic members, adequate to the present high quality data, and 3) a statistically robust approach to analyse the orbits and to establish a dynamical traceback age. In this study, we made a special effort to prepare a clean sample with precise and uniform data.
This work is structured as follows. In Section 2 we present the spectroscopic observations we carried out and the process to measure precise radial velocities for new and archive data. We also describe the improvements of our method to select a bona fide sample of kinematic members from our initial list of candidates from the literature. In Section 3 we describe the algorithm used to derive the dynamical age and analyse in detail the orbits of the bona fide members. In Section 4 we discuss the results obtained and, finally we conclude in Sect. 5.

Data and sample selection
In this section we present a compilation of confirmed members and new candidates reported in the literature over the past decade. In order to have a sample with homogeneous stellar parameters, we use the 5D astrometric solution (positions, parallaxes, and proper motions) of the Gaia DR2 catalogue. We complement these data with a set of radial velocities (from our own observations plus archival data) analysed using the same methodology. In this study we use the radial velocities published in the literature and in Gaia DR2 only to compare with our own determinations. In Table 1 we review the selection process from the initial compilation to the final sample.
Our initial sample is based on Torres et al. (2008), Schlieder et al. (2012), Malo et al. (2013), Malo et al. (2014a), Gagné et al. (2015a,b), Alonso-Floriano et al. (2015), Messina et al. (2017),  and . This results in a sample of 236 stars after removing the sources in common among the several studies. Binaries and multiple systems are counted as one single object unless they have been resolved in previous studies. These authors used different algorithms based on the kinematics (and included the photometry in some cases) to identify new candidates or confirm members of β Pic. Most of these studies are pre-Gaia or were carried with partial information (missing parallaxes and/or radial velocities). For this reason, it has been mandatory to develop a tool to reject kinematic outliers with our homogeneous and precise astrometry and spectroscopy (see Sect. 2.3).

Proper motions and parallaxes
We use the proper motions and parallaxes of the Gaia DR2 catalogue which constitute the most recent and precise astrometric measurements available to date for our sample. To identify the Gaia DR2 counterparts of the stars in our sample we used the 2MASS source identifier (which are given in the original tables used to construct our initial sample) and the TMASS_BEST_NEIGHBOUR table available in the Gaia archive. For 42 sources we did not find a counter part with this procedure, so we manually refine the match considering position and magnitude. Finally, we find proper motions and parallaxes for 222 stars in our initial sample. There are 8 sources in Gaia DR2 with only the two-parameter solution and 6 not in Gaia DR2 (see App. A). The median of the uncertainties of this sample are ∼ 0.1 mas yr −1 in proper motions and 0.08 mas in parallax which lead to a median error in tangential velocity of 0.19 km s −1 , obtained taking into account the correlations among the astrometric parameters (see Table 1).

Radial velocities
The scarcity and quality of the radial velocities of β Pic stars are currently two of the main limitations to derive an accurate estimate of the dynamical age of the association. Even though many radial velocity measurements are available in the literature (e.g. Torres et al. 2006, Kharchenko et al. 2007, Shkolnik et al. 2012 we resort to reanalyse the spectra available in public archives to ensure that all the radial velocities are derived using the same methodology. The consistency and homogeneity of the individual measurements is indeed particularly important in dynamical traceback analysis (see e.g. Miret-Roig et al. 2018).

New spectroscopic observations
We performed spectroscopic observations of β Pic stars with three different instruments. The FEROS spectrograph (Kaufer et al. 1999) mounted on the ESO/MPG 2.2 m telescope operated at La Silla (Chile) was used to collect the spectra of 43 stars as part of programme 103.A-9009 (PI: W. Brandner). These observations were performed in OBJCAL mode that allows for simultaneous acquisition of the target spectrum and the calibration lamp during July and August 2019. We observed 8 stars with the CAFE spectrograph (Aceituno et al. 2013;Lillo-Box et al. 2020) mounted on the 2.2 m telescope of the Calar Alto Observatory (programme: H18-2.2-015, F19-2.2-002, PI: D. Barrado). The observations were carried out from July to October 2018, right after the upgrade of the instrument. The data were processed using the new instrument pipeline described in Lillo-Box et al. (2020), which performs the basic reduction and extracts the radial velocities. Finally, another 14 stars were observed with the SOPHIE spectrograph (Perruchot et al. 2008) mounted on the 1.93 m telescope of the Haute-Provence Observatory (programmes: 2018A-PNPS005, 2019A-PNPS008, PI: H. Bouy). These spectra were obtained in August 2018 and May 2019 and were processed with the instrument standard data reduction pipeline which measures radial velocities by numerical cross-correlation techniques. The median signal-to-noise ratio of our observations is 25.  Notes. Number of spectra analysed and number of new spectra obtained in this study with different instruments. The total number of spectra analysed is 723 and 141 of them are new. The (maximum) resolving power and spectral range of each spectrograph are indicated.

Spectroscopic archival data
In addition to the observations conducted by our team we did an exhaustive search for the spectra available in public archives. As shown in Table 2, a total of 582 spectra have been collected from the European Southern Observatory (ESO) and the ELODIE archives. We reanalysed all these data (see Sect. 2.2.3) and provide radial velocities for a larger number of stars. Table 2 shows the instruments that have been used in this study and the respective number of spectra analysed in each case. We specify the number of new spectroscopic observations presented in this work (see Sect. 2.2.1) which constitute a 20% of all the spectra. We note that some sources have been observed several times with the same or various instruments. In fact, the 723 spectra correspond to 81 different stars, 54% of which have been observed once, 18% twice, and the rest three or more times. We refer to Sect. 2.2.3 for a description of how we combined the different radial velocity measures for the same star.

Radial velocities determination
The observed and downloaded spectra were reduced using the official pipeline available for each instrument. We derived radial velocities by cross-correlating the reduced spectra of the stars with the closest mask to its spectral type. We used six different masks of spectral types A0, F0, G2, K0, K5, M5 and the iSpec routines for this purpose (Blanco-Cuaresma et al. 2014;Blanco-Cuaresma 2019). This procedure follows the cross-correlation technique (Baranne et al. 1996;Pepe et al. 2002) and fits a Gaussian profile to the cross-correlation function (CCF) to derive the radial velocity and associated uncertainty. We discard the radial velocity measurements resulting from a poor fit to the CCF due to, e.g., a low signal-to-noise ratio of the spectrum or a mismatch between the spectral type of the star and the adopted mask. We used the effective temperatures given in Gaia DR2 as a rough estimate of the spectral type of the star to choose the corresponding mask. For each star we compute the radial velocity scatter from the results obtained with 3 different masks: the closest mask (M) to the spectral type of the star, one before (M-1) and one after (M+1). We add this number in quadrature to the formal uncertainty returned from the iSpec routines. The later step accounts for the observed fluctuation on the radial velocity results derived from different masks 2 . We derive radial velocities for 81 stars of our initial sample of β Pic candidates by combining our own observations with archival spectra (see Table C.7). In the case of multiple radial velocity measurements for the same star we proceed as follows. For each radial velocity solution (for a given star) we generate a sample of 10 000 synthetic measurements from a Gaussian distribution where the mean and variance correspond to the radial velocity and its uncertainty. We repeat this process for all radial velocity measurements of the star. Then, we take the mean of the joint distribution of synthetic radial velocities as our final result for the radial velocity of the star. The uncertainties on the resulting radial velocity are computed from the 16% and 84% percentiles of the joint distribution of synthetic radial velocities. We note that for 10 stars (12% of the radial velocities we determine) our radial velocity is the first measurement ever made. This is an important product of our work since these data can be used to assess the membership and to study the dynamics of the association in 6D. Additionally, six of them are in our final bona fide sample of 26 stars (see Sect 2.3).
In Table 1 we compare the quality of our radial velocities with the Gaia DR2 catalogue and with previous ground-based spectroscopic surveys. We find a radial velocity in the literature for 137 sources in our initial sample. These measurements come from a variety of different surveys with different qualities and methods to determine the radial velocity. Our measurements are homogeneous and about 40% more precise than this compilation which is crucial for the success of our work. We see that the radial velocities of our sample are twice as precise as the Gaia DR2 radial velocities and we have a measurement for a larger number of sources. We identified and discarded 35 sources which have been classified as binaries in previous works. In order to include the binaries in our study we would need to determine the radial velocity motion of the centre of mass which is beyond the scope of this work. Figure 1 shows the comparison of the radial velocities derived in our study with the ones in Gaia DR2 and the ones in other spectroscopic surveys in the literature for the 42 single stars with 6D data in our sample. We found hints of binarity in two sources (2MASS J19312434-2134226 and 2MASS J22571130+3639451) and we discarded them from the analysis (see App. B). The median difference and root mean square error (RMSE) between the Gaia DR2 radial velocities and our measures are 0.7 km s −1 and 1.0 km s −1 , respectively. These values are computed disregarding the source with a radial velocity difference of about 5 km s −1 . The Gaia DR2 radial velocity of this star is based only in two transits which is probably the reason for its large uncertainty. If we compare the radial velocities from the literature and our sample we obtain a median difference and RMSE of 0.3 km s −1 and 0.9 km s −1 , respectively. Since we believe that the homogeneity and precision of our radial velocities is superior to any other sample we only use our measurements in the forthcoming analysis.

Kinematic sample selection
In this section we present the kinematic selection we designed to discard kinematic outliers in our sample. Kinematic outliers in the context of the present study refer to sources with a velocity significantly different than the group, either because they are non-members or because they are variable due to e.g. multiplicity. First, we introduce the notation we adopted to refer to position and velocity coordinate systems. We used the curvilinear heliocentric coordinates (ξ , η , ζ ) defined in Asiain et al. (1999) to place the stars in the configuration space. This coordinate system is centred at the current position of the Sun (R = 8.4 kpc) and rotates around the Galactic centre with a frequency of the circular velocity of ω = 28.81 km s −1 kpc −1 (Irrgang et al. 2013). It has the advantage that it minimises the variation in each component of the configuration space. The radial component ξ points towards the Galactic anti-centre, the azimuthal component η is measured along the circle of radius R and is positive in the sense of the galactic rotation, and the vertical component ζ is defined positive towards the north Galactic pole. We also refer to the corresponding velocities asξ ,η ,ζ . The second reference system considered in this work is the Cartesian heliocentric system. The spatial components X, Y, Z and velocity components  Table 1). Our final sample of 26 stars is represented by the filled markers. We note that our radial velocity uncertainties (horizontal error bars) are smaller than the markers in most of the cases. The largest uncertainty corresponds to the β Pic star (see footnote 2). data. The sample of 26 kinematically selected members is represented by the orange dots and the kinematically rejected sources are the black squares. The source 2MASS J11493184-7851011 (red dot) is retained by our kinematical criteria but is discarded due to its orbit (see text).
U, V, W, are defined with X, U pointing towards the Galactic centre, Y, V towards the direction of Galactic rotation and Z, W towards the north Galactic pole. We use a peculiar solar motion of (U , V , W ) = (11.1, 12.24, 7.25) km s −1 (Schönrich et al. 2010).
In Figure 2 we represent the velocity distribution of the 42 single stars with Gaia astrometry and radial velocities from this Table 3. Final bona fide sample of 26 members of β Pic selected to determine the dynamical age.   (2018); (6) Houk (1982) Notes. (a) Spectral types between parentheses were estimated from the absolute Gaia G−band magnitude.

2MASS ID
Article number, page 5 of 19 A&A proofs: manuscript no. aanda work. We see that a number of sources have a significant scatter. Most of them were classified as members of β Pic with pre-Gaia astrometry or with no radial velocity information and clearly appear to be kinematic outliers with our extremely precise data. We discard the kinematic outliers in the 3D velocity distribution (ξ ,η ,ζ ) in a similar way to what we did in Miret-Roig et al. (2018). The major improvement is that in this work we use a robust estimator of the covariance matrix (the Minimum Covariance Determinant from Sklearn, Pedregosa et al. 2011) to fit the central location ( #» µ ) and the covariance matrix (Σ) of the velocity ellipsoid of the association. Then, we compute the Mahalanobis distance of each object defined as: In Figure B.1 we show the distribution of Mahalanobis distances. We use the percentile p 65 to discard the kinematic outliers and retain 27 kinematic members (dots in Fig. 2). This threshold is empiric and represents the best compromise between rejecting kinematic outliers which hinder the traceback analysis and keeping kinematic members in the final sample. When we compute the orbits of our targets (see Sect. 3.1) we immediately see that one of them (2MASS J11493184-7851011, red circle in Fig. 2) has an orbit significantly different from the main group and thus, we discard this object. This star has a kinematics similar to β Pic but it is at > 3σ in positions with respect to β Pic. We also checked that this object has the largest Mahalanobis distance to the centre of the velocity distribution. We refer to App. B for a detailed discussion, source-by-source, of the kinematically rejected sources. The final sample contains 26 bona fide members of β Pic and their 3D positions and velocities are given in Table 3 (available at CDS).

Bona fide β Pic sample
In this paper we made a major effort to prepare a robust sample of β Pic members with the best precision possible in their determination of the positions in the 6D space phase. Then we used this valuable data to identify and remove kinematic outliers. In this section we review the main characteristics of our final sample.
The relative error in parallax of these members is less than 1% which allows us to compute the distance as the inverse of parallax. We note that four stars have a parallax error < 0.1% at distances up to 50 pc. The median relative errors in proper motions are of 0.3% in right ascension (µ α * ) and 0.09% in declination (µ δ ). The precision in µ α * and µ δ is similar but a few members have µ α * close to zero which increases the relative error. The β Pic star is the brightest source (G = 3.7 mag) and causes a fraction of the pixel used in the standard Gaia DR2 analysis to be saturated. Hence, measurements of its centroid position and the resulting astrometry are less precise than for fainter sources (G > 6 mag) (Lindegren et al. 2018).
In Table 4 we provide the median heliocentric position and velocity of β Pic. We see that the observational uncertainties in positions (σ err ) are of the order of tenths of parsecs and thus, the observed dispersion (σ obs ) can be interpreted as an intrinsic dispersion (σ int ). The dispersion in the Galactic plane (X, Y components) is twice the vertical dispersion (Z). When we look at the velocity dispersion we find that the median errors in velocity (σ err ) are significantly smaller than the velocity dispersion observed (σ obs ), indicating the presence of an intrinsic cosmic dispersion (σ int ). Therefore, the dispersion we observe in Fig. 2

Traceback analysis
In this section, we describe our methodology to perform the traceback analysis which is based on the work of Miret-Roig et al. (2018) with some improvements.

Towards a dynamical age estimate
We consider the same 3D Milky Way potential as in Miret-Roig et al. (2018), to integrate the equations of motion. This model is based on the Allen & Santillan (1991) potential which consists of a spherical central bulge, a disc, and a massive spherical halo, but with updated parameters taken from Irrgang et al. (2013, their Table 1). Hereafter we refer to this model as new A&S and we compare it with other axisymmetric models in Section 3.3. In Figure 3,   the time at which the members of the association were most concentrated in space. The algorithm to measure the degree of concentration, hereafter the size of the association, is of uttermost importance and different strategies to compute the size have been used in the literature. These different methodologies have significantly contributed to the large spread in the dynamical traceback ages determined. In general, the size of the association is estimated with the empiric standard deviation in the spatial coordinates. However, it is very sensitive to the presence of outliers, i.e. members which significantly deviate from the mean position of the association which are not necessarily contaminants. In this section we present three strategies to estimate the size of the association as function of time. Some of them are based on classical functions used in the literature (i.e. the variance) and others are novel, represent the overall variance of the association, and are independent of the coordinates chosen. In the following, we define the three functions we use to estimate the size of the association in a way that they all have units of length.
-The size in the radial, azimuthal, and vertical directions (S ξ , S η , S ζ ) are the squared root of the diagonal terms of the covariance matrix in each direction. -The Trace Covariance Matrix Size (S TCM ) is defined as: -The Determinant Covariance Matrix Size (S DCM ) is defined as: Each of these expressions are computed from the covariance matrix of the association in the configuration space. We used two  Fig. 5. Dynamical age distribution of the bona fide β Pic members, obtained with the robust estimate of the covariance matrix. The distribution obtained with the S TCM size estimator is colour coded with the 68%, 95%, and 99.7% highest density intervals. The distribution obtained with the S DCM size estimator is shown in dashed lines and the same highest density intervals are shown. The orbits were integrated using the new A&S potential and we computed 1 000 bootstrap repetitions. different algorithms to estimate the covariance matrix, namely the empirical covariance estimation, and the robust covariance estimation, both from the Sklearn packages (Pedregosa et al. 2011). Whereas the first corresponds to the classical maximum likelihood estimator, the second is less sensitive to outliers in the dataset.
The size estimators S ξ , S η , and S ζ , when computed with the empirical covariance estimation, correspond to the classical standard deviation in each direction. The other two size estimators (S TCM and S DCM ) can be interpreted from the eigenvalues of the covariance matrix. The trace of the association, often referred as the total variance of the covariance matrix, coincides with the sum of its eigenvalues. In Eq. 2, we introduced a factor 1/3 (in a 3D space) so that we can interpret the S TCM estimator as the arithmetic mean of the variances in the individual components. In any case, this multiplicative factor changes the absolute value of the size estimator but not the locus of the minimum, which is our main interest. The determinant of the covariance matrix, also known as the generalised variance, can be interpreted as the geometric mean of the eigenvalues of the covariance matrix. Then, the volume of the association is proportional to the squared root of the determinant of the covariance matrix. Finally, we define the diagonal of the Determinant Covariance Matrix Size (S dDCM ) analogously to the S DCM size but only considering the diagonal terms, i.e. neglecting the correlations among the three spatial components. This is not a good estimator of the size of the association since it neglects part of the information included in the covariance matrix. However, it can be understood as a geometric mean of the size estimators S ξ , S η , S ζ , so we include it only for comparison.
In Figure 4 we show the six parameters defining the size of β Pic (S ξ , S η , S ζ , S TCM , S DCM , S DV M ) computed with the empirical covariance estimate and the robust covariance estimate as function of time. It is remarkable that the minimum size obtained with the empirical covariance estimate (top panels) depends on the size estimator while we find a minimum at a similar times for all the size estimators considered with the robust covariance estimate (bottom panels). This is because the robust covariance estimate gives less weight to sources with a large dispersion, attenuating the impact of outliers.
From now on, we only consider the size estimates computed with the robust covariance estimates. In the left bottom panel of Fig. 4, we show the dispersion in the radial, azimuthal, and vertical direction independently. We see that the vertical component does not provide useful information to constrain the age of the association while the two components in the Galactic plane have a minimum at a similar time. In this panel we highlighted the azimuthal component (S η ) which is the size estimator we used in Miret-Roig et al. (2018).
In the middle bottom panel, we add the size from the determinant of the covariance matrix (S DCM ) and, for comparison, the inaccurate size using only the diagonal values of this matrix (S dDCM ), i.e. with and without correlations, respectively. Both curves have close minima with a time difference of ∼ 1 Myr, and are also similar to the age obtained with the S η size estimator. The correlations reduce the value of the determinant and in consequence, the absolute value of S DCM , estimating a birth size of the association of ∼ 5 pc. In the right bottom panel, we include the size estimator from the trace (S TCM ). As mentioned before, the S TCM and S DCM sizes correspond to the arithmetic and geometric mean of the eigenvalues of the covariance matrix, respectively. These two statistics are related by an inequality in which the arithmetic mean is always larger than the geometric mean, and they are only equal if all the individual values are the same. This corresponds to an isotropic covariance matrix, which is not the case in our study.
Currently, thanks to the excellent astrometric precision of Gaia and the homogeneous precise radial velocity sample derived in this work, the observational uncertainties are no longer what dominates the uncertainties in the dynamical age. We propagated the present uncertainties with an analytic approximation (Miret-Roig et al. 2018) and estimated that the dispersion due to observational uncertainties is 2 pc at the time of minimum size. At birth, the association had a S TCM size of ∼ 7 pc (see Fig. 4), similar to what has been observed in star forming regions such as Ophiuchus (Cánovas et al. 2019), Taurus (Galli et al. 2019), and Corona Australis (Galli et al. 2020).
As mentioned, the sample selection (i.e. the presence of contaminants or unidentified binaries) is extremely important. To estimate the impact of the sample selection on the age, we took 1 000 random samples of the 26 bona fide β Pic members and estimated the dynamical age with each. Then, we can determine a dynamical age and a robust uncertainty from the distribution of ages. In Figure 5 we report a kernel density estimate of the age distribution with a bandwidth of 1 Myr, this value is smaller and of the order of the age uncertainties. In Table 5 we report the mode and the 68%, 95%, and 99.7% highest density intervals 3 of the age distribution. Considering the S TCM size estimator and the 68% highest density interval we find a dynamical age of β Pic of 18.5 +2.0 −2.4 Myr (see Table 5). With the S DCM size estimator we obtain a similar age, 17.6 +3. 5 −2.9 Myr. We note that the two values agree within a 1 Myr difference which is significantly smaller than the age uncertainty.

Signs of substructure at birth time
When we look at Fig. 3 we see that at birth time some stars appeared to be more concentrated forming a core (filled dots) while a few members appeared to be more dispersed (empty dots). To identify these two populations we compute the Mahalanobis distance (Eq. 1) with the robust central location and covariance of the 3D spatial distribution (ξ , η , ζ ). In Figure C.1 we show the distribution of the Mahalanobis distances. We used the percentile p 68 to separate the core from the peripheral stars which results in 17 core stars and 9 peripheral stars (see Table 3). These stars were selected at birth time in the space of positions and in this space they appear most concentrated (see Fig. 3). Interestingly, in the present the stars forming the core appear more dispersed than those originally more dispersed. In the velocity space, both populations are mixed in the present and at birth time.
It is worth mentioning that if we use only the 17 core stars to study the dynamical age, we obtain an age very similar to the value we obtained in Sect. 3.1. With the S TCM size we find an age of 18.8 +1.7 −2.1 Myr and with the S DCM size an age of 17.6 +3.5 −1.2 Myr. As expected, in this case were all the stars are well concentrated at birth time, the age is independent of the covariance estimate used (empirical or robust). Additionally, the small bump we observe in Figure 6 at ∼ −15 Myr disappears with the age distribution obtained only with the 17 core stars. In short, if we use the core sample of 17 stars to trace back the age of β Pic we find variations of less than 1 Myr with respect to the value we obtained in Sect. 3.1, with all the covariance and size estimates considered in this study. This is the first time that the spatial distribution of β Pic is analysed in detail and these results should be revisited with a larger sample of members.

Effect of the Galactic potential
In this section we discuss the effect of considering different Galactic axisymmetric potentials and including nonaxisymmetric structures such as spiral arms, on the dynamical age. First, we considered two additional axisymmetric potentials, namely the one of McMillan (2017, hereafter McMillan17) and the one of Bovy (2015, hereafter MWPotential14). These two models together with the new A&S model have similar rotation curves in the range of radius relevant here with only slight differences in the mass distribution as can be seen in their respective rotation curves ( Figure  We also used a non-axisymmetric potential which accounts for the spiral arms in addition to the axisymmetric potential described in Sect. 3.1. The 3D spiral model is the PERLAS spiral arms from Pichardo et al. (2003, hereafter new A&S + spiral arms (P03)). The locus is the one following Drimmel & Spergel (2001) and has a pitch angle of 15.5 • . We take a pattern speed of Ω P = 21 km s −1 kpc −1 and a mass of 0.04% of the disc mass. These values are in agreement with the values proposed in Antoja et al. (2011). Recently, Eilers et al. (2020) estimated a density contrast at the solar radius of 20% which is similar to the amplitude of the arms used here which leads to a contrast of around 23% (Antoja et al. 2011).
In Figure 6 we present the dynamical age distribution obtained with different axisymmetric potentials and with the nonaxisymmetric potential with spiral arms. In Table 5 we report the percentiles of the dynamical age distribution for each of the potentials considered. The variations in the dynamical age due to the Galactic potential are minimal, and they are all compatible with the value we obtained in Sect. 3.1. Therefore, we conclude that the variations in the dynamical age produced by different Galactic potentials are much lower than our main source of uncertainty, i.e. the membership. This is valid for the potentials we have tested, the parameters of which are constrained by recent observations of the Milky Way and can be explained for the short integration times given the low age of the association. Given that the different Galactic potentials considered here lead to changes in the dynamical age smaller than the current uncertainties, we decided to keep the results obtained with the new A&S potential which has less parameters.

Discussion
In the previous section we discussed different strategies to determine the dynamical traceback age of β Pic. All of them are compatible, with differences of 1 Myr, significantly smaller than the age uncertainties. From now on, we adopt the an age of 18.5 +2.0 −2.4 Myr, obtained for the sample of 26 bona fide members, with the S TCM size, and with the axisymmetric potential. Our study, provides the first traceback age which conforms with other dynamical ages recently published in the literature such as the expansion or the forward modelling algorithms and with ages based on evolutionary models such as the lithium depletion or the isochronal ages. In Table 6 we present a compilation of previous age estimates published in the literature and we see that our determination is compatible with the majority of them. The first reliable age determination of the β Pic star and its moving group was an isochronal age presented in Barrado y Navascués et al. (1999), 20 ± 10 Myr, in full agreement with our current estimate.
The earliest traceback studies of β Pic obtained an age of 11-13 Myr (Ortega et al. 2002, Song et al. 2003, Ortega et al. 2004, younger that what we find here. These differences are most probably due to the large observational uncertainties of the pre-Gaia astrometry, and thus to the presence of a significant number of kinematic contaminants. Those authors used the maximum size between stars to determine the age of the association. We did not consider this size estimator in our study but it is clearly sensitive to the presence of outliers in the dataset. Miret-Roig et al. (2018) measured a dynamical age of β Pic of 13 +7 −0 Myr with a method very similar to the one presented in this work. We believe that the main differences between these two studies are 1) the precision of the 6D space phase positions, 2) the new sample selection based on a robust estimate of the 3D velocities covariance matrix, and 3) the new orbital analysis which uses an improved size estimator of the association. In our previous study, we used the Gaia DR1 astrometry and a compilation of radial velocities from the literature without any treatment. Now, we used the improved precision of Gaia DR2 and a uniform sample of radial velocities. The median uncertainty in the DR1 parallaxes and proper motions were 0.3 mas and 0.2 mas yr −1 , respectively, compared to the values 0.05 mas and < 0.1 mas yr −1 now available from the DR2. We discarded 10 objects of our previous sample for being classified as binaries and 5 others do not have a radial velocity measurement in our work. These leaves only 6 objects in common between the two works (23% of our new sample). The black solid line in Figure 4 (S η ) corresponds to the methodology used in Miret-Roig et al.
(2018) (their blue curve in Fig. 7). We see that with the empiric covariance (top panel) matrix we still find younger dynamical traceback ages. On the contrary, the size estimator S η , if we use a robust estimate of the covariance matrix, we recover a similar age to the one reported in Sect. 3.1. Another technique used in the literature to measure a kinematic age consists in studying if the association is under expansion. Torres et al. (2006) found a linear relation between the velocity and the position in the direction of the Galactic centre which results in an age of ∼ 18 Myr. In a similar approach, Mamajek & Bell (2014) found an age of 21 +10 −5 Myr taking into account the positions and velocities in the Galactic plane. Both results are compatible with what we obtain here. In addition, we used our new accurate sample to estimate an expansion age of β Pic. We fitted a line between the cartesian heliocentric positions XYZ and velocities UVW. We find evidence of expansion in the direction towards the Galactic centre and in the direction of Galactic rotation with slopes of κ X = 0.057 ± 0.006 km s −1 pc −1 and κ Y = 0.033 ± 0.008 km s −1 pc −1 , respectively. In the vertical direction we find a slope of κ Z = −0.02 ± 0.02 km s −1 pc −1 , slightly negative but compatible with zero. These coefficients result in an expansion age 4 of 17 ± 2 Myr and 29 ± 4 Myr in the radial and azimuthal directions. If we combine these measures with a weighted mean as done by Mamajek & Bell (2014), we obtain an expansion age of 20 ± 4 Myr, with excellent agreement with our traceback age.
Recently, Crundall et al. (2019) provided a new tool (Chronostar) to determine a dynamical age applying the forward-modelling technique, and obtained an age of 18.3 +1.3 −1.2 Myr. It is interesting to see how similar the results of their study are to ours despite the different sample of members (we have 15 members in common, 25% of their sample) and method used. These results prove that both methods are complementary. Their method allowed them to detect the β Pic members among a large catalogue of field stars while ours provides a deeper orbital analysis allowing us to discover, for example, the existence of a central core and a more dispersed structure at birth time.
Finally, it is important to mention that the age estimates in the literature based on the Li depletion or isochronal fitting obtained values very similar to the one obtained here and, in general, with a lower dispersion than the dynamical age estimates obtained up to now (see Table 6). If we exclude the work of Macdonald & Mullan (2010) which obtained an age of ∼ 40 Myr, twice the other works, we obtain a median value of 21 ± 4 Myr which is in good agreement with the age we measured. This is an important result since our method is independent of evolutionary models and these are two very different strategies.

Conclusions
In this work we measured a dynamical, traceback age of the β Pic moving group of 18.5 +2.0 −2.4 Myr which is compatible with ages based on evolutionary models. Our age estimate is the first traceback age that reconciles the ages determined by the traceback method with other dynamical ages (expansion, forward modelling), lithium depletion ages or isochronal ages.
The precision in the dynamical traceback age we achieved in this study is thanks to the combination of the Gaia DR2 astrometry and the uniform radial velocity sample of single stars that we produced in this work. We measured the radial velocity of 81 candidate members of β Pic in a uniform manner. For ten sources, our measure is the first radial velocity estimate. This is an important result of our work which allowed us to identify 15 kinematic outliers from our initial sample and two new potential spectroscopic binaries.
Our improved algorithm to determine the age (based on our previous work, Miret-Roig et al. 2018) provides a more rigorous kinematic sample selection and an improved orbital analysis. We showed the importance of using a robust estimate of the covariance matrix (instead of an empirical one) to minimise the impact of outliers (sources which deviate from the central locus of the association which are not necessarily contaminants). We explored different size estimators computed from the covariance matrix to determine the dynamical age (the standard deviation in different directions, the determinant, and the trace). All of them provide dynamical ages with differences of less than 1 Myr, i.e. compatible within the uncertainties, when computed from the robust covariance matrix. Our thorough orbital analysis allowed us to propose the existence of a central core of 17 stars which appeared more concentrate at birth time.
In this study we showed that different potentials (i.e. axysimmetrics and including the effect of spiral arms) lead to changes in the dynamical age that are within the current uncertainties. Nowadays, the major source of uncertainty in the dynamical, traceback age is the sample selection and the errors in the radial velocity estimates. For this reason, we stress the importance of choosing samples with accurate radial velocity data, with uncertainties comparable to the imminent eDR3 Gaia release. This is crucial to reject kinematic contaminants and binaries and for the success of a traceback analysis. ber of the Hercules Lyra association by López-Santiago et al. (2006) and  could not confirm its membership because they did not consider the Hercules Lyra association in their analysis. Our radial velocity estimate (12.01 ± 0.12 km s −1 ) is similar to a recent value from the literature (11.658 ± 0.006 km s −1 , Soubiran et al. 2018).

2MASS J02232663+2244069
Our radial velocity measurement (12.60 ± 0.15 km s −1 ) is consistent with the one from the Gaia DR2 catalogue (12.1 ± 0.6 km s −1 ). This source was listed as a high probability (99%) member of β Pic by Malo et al. (2013) based on a radial velocity and a proper motion which differ by 2 km s −1 and 6 mas yr −1 , respectively, from Gaia DR2. The different data could explain why this source was discarded by our kinematic selection and the membership of this source has been revised with our data.

2MASS J03573393+2445106
We have three spectra for this source with radial velocity measures of 13.46 ± 0.18 km s −1 (2018-08-12), 13.44 ± 0.18 km s −1 (2018-08-14), and 15.30±0.14 km s −1 (2019-11-30). This source is rotationally variable (0.86 days, Hartman et al. 2011) which could explain the variations in the radial velocity that we measure. This source is a candidate of spectroscopic binary which requires more follow-up observations to confirm it. We also note that Crundall et al. (2019) classified it as a field contaminant.

2MASS J05004928+1527006
This source was classified as a member by Schlieder et al. (2012) based on a predicted radial velocity of 13.70 ± 2.03 km s −1 . We measure a radial velocity of 18.4 ± 0.3 km s −1 , similar to what Article number, page 12 of 19 N. Miret-Roig et al.: Dynamical traceback age of the β Pictoris moving group has been reported in the literature (White et al. 2007), and significantly different to the predicted value used in the previous membership analysis. Additionally, this source has been classified as a member of the Taurus-Auriga complex , and therefore is a likely contaminant in β Pic.

2MASS J08475676-7854532
This source was classified as a candidate of β Pic based on a predicted radial velocity of 13.4 ± 1.5 km s −1 (Malo et al. 2013). This value is significantly different from our measurement of 23.1 ± 0.3 km s −1 and with the literature (23.4 ± 0.3 km s −1 from Malo et al. 2014a). Using Gaia, it has been proposed as a member of η Chamaeleontis (Cantat-Gaudin et al. 2018).

2MASS J11493184-7851011
This source was classified as a β Pic candidate based on a predicted radial velocity of 10.8 ± 1.6 km s −1 and a predicted distance of 68 pc ( = 14.7 mas) by Malo et al. (2014a). Our two radial velocity measurements differ by about 1.3 km s −1 between them but have a mean value of 14.5 ± 0.8 km s −1 , which is not compatible with the predicted radial velocity in that study. The Gaia DR2 parallax of this source is 9.92 ± 0.03 mas, indicating this source is probably a contaminant. A recent study classified this source as a Chamaeleontis (Schneider et al. 2019).

2MASS J13545390-7121476
This source was classified as a candidate member of β Pic by Malo et al. (2014a) based on proper motions values which differ of about 20 mas yr −1 from the values of Gaia DR2. This source is probably a contaminant.

2MASS J19312434-2134226
Our radial velocity measurement (−36.6 ± 1.8 km s −1 ) is not consistent with the literature (e.g. Shkolnik et al. 2012 measured a radial velocity of −26.0 ± 1.8 km s −1 and Malo et al. 2014a −25.6 ± 1.5 km s −1 ) with a difference of about 10 km s −1 . We checked the CCF and there are hints it might be a spectroscopic binary. In addition, a recent study classified this as a member of the Argus association (Janson et al. 2017).

2MASS J21212873-6655063
This source was classified by Malo et al. (2014a) as a high probability (99.9%) member of β Pic. However, their analysis was based on pre-Gaia astrometry and the proper motions they used differ about 20 mas yr −1 from the one of Gaia DR2, indicating the membership should be revised.

2MASS J23314492-0244395
This source was classified as a β Pic candidate member by Malo et al. (2013). However, their analysis was based on pre-Gaia astrometry and the proper motions they used differ about 10 mas yr −1 to the ones of Gaia DR2.

2MASS J23512227+2344207
Our radial velocity measurement (−1.0 ± 0.3 km s −1 ) differs by about 1 km s −1 from the measurement of Shkolnik et al. (2012) (−2.1 ± 0.5 km s −1 ). Binks & Jeffries (2016) provided another radial velocity measure for this source (38.6 ± 1.6 km s −1 ), with a discrepancy of several tens of km s −1 . Based on their measurement, they rejected this source as a β Pic member and also suggested the possibility of a binary system to explain the differences observed. Messina et al. (2016) classified this source as a single star based on a study of photometric variability. Further work is required to confirm the binarity of this source. Additionally, other authors have classified this source as member of other moving groups (e.g. Shkolnik et al. 2012, Klutsch et al. 2014.

2MASS J21183375+3014346
This source was classified as a candidate member of β Pic by Schlieder et al. (2012) with a predicted radial velocity of −15.1 ± 0.9 km s −1 . This value is significantly different from our radial velocity measurement (−22.0 ± 0.3 km s −1 ). Additionally, Shkolnik et al. (2017) recently measured a radial velocity similar to ours (−22.5 ± 0.8 km s −1 ) and rejected the β Pic membership of this source.

2MASS J22571130+3639451
This source was classified as a candidate member of β Pic by Schlieder et al. (2012) with a predicted radial velocity of −10.0 ± 0.9 km s −1 although their measured radial velocity was −20 ± 1.2 km s −1 . We have analysed 8 spectra of this source and obtained a variable radial velocity between −10 km s −1 and −20 km s −1 , indicating this is probably an unresolved spectroscopic binary.
2MASS J16120516-4556242, 2MASS J17092947-5235197, 2MASS J18430597-4058047, and 2MASS J20105054-3844326 These sources were classified as new members of the β Pic moving group by  with no radial velocity measurements. The first estimation of their radial velocity provided for the first time in the present work, shows that their velocity is not compatible with the velocity distribution of β Pic, suggesting they might be contaminants.   (1) 2MASS source ID; (2) number of spectra used to measure the radial velocity; (3) spectrograph; (4-6) radial velocity, radial velocity error and julian day of individual measures; (7) radial velocity; (8) binary flag (from literature), the asterisk indicates the binaries identified in this work; (9) contamination flag (see App. B). This file will be available at the CDS.