When and under what conditions did the first planet form? The oldest planetary system that we know of is ∼11.2 Gyr old (Kepler-444, Campante et al. 2015;Mack et al. 2018), with a host star metallicity of [Fe/H] = -0.52. The first generation of stars and planets is expected to have metallicities similar to or below that of Kepler-444. Such a metal-poor environment poses challenges for planet formation. For gas giant planets, core mass growth via accretion halts in the extremely metal-poor regime, resulting in a strong dependence of planet occurrence rate on stellar metallicity-known as the planet-metallicity correlation (e.g., Fischer & Valenti 2005). For terrestrial planets, the planet-metallicity correlation is weaker (Wang & Fischer 2015) and perioddependent (Wilson et al. 2018;Petigura et al. 2018). Therefore, metallicity is a key to regulate planet for- * NSF Astronomy and Astrophysics Postdoctoral Fellow.
mation. More interestingly, some elements (e.g., αelements) may play a bigger role than others (Adibekyan et al. 2012;Bashi & Zucker 2019).
The abundance of the elements increases as the Milky Way evolves. Among the stellar populations in the Milky Way, halo stars are the most metal-poor on average and therefore serve as an ideal test ground for planet formation in metal-poor environments. However, identifying halo stars is difficult because halo stars are relatively rare in the solar neighborhood, representing of order 1% of the local population; discovering them has historically been laborious; and detailed characterization via high-resolution spectroscopy has been timeconsuming (see Carney & Latham (1987) and associated works).
The Gaia mission (Gaia Collaboration et al. 2016) significantly changed the landscape of the study of halo stars, because we can use it to identify metal-poor stars based on their kinematics, thanks to the well-known correlation between the two (Eggen et al. 1962). Gaia's astrometric and radial velocity measurements of one bil-lion stars in the Milky Way provide an efficient way of identifying halo stars. As such, we set out to identify stars with halo kinematics using the Gaia DR2 data (Gaia Collaboration et al. 2018). The kinematically selected halo stars are then validated using data from the APOGEE spectroscopic survey (Majewski et al. 2017) and other photometric surveys to check the purity of the sample.
Many of the halo stars are in the field of view of the Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015). Searching for transiting planets around halo stars are therefore made possible by the Gaia and TESS missions. Our halo sample is much larger than previous samples used for planetary searches via radial velocity (Sozzetti et al. 2006;Faria et al. 2016). Since > 70% stars in our sample have [Fe/H]< -1, our search will likely result in planetary systems that are more metal-poor than the current record holder at [Fe/H]=-0.84 (Anglada-Escudé et al. 2014). The methods and results of the planet search are described in a companion paper (Boley et al. 2021).
It is of interest to determine the absolute amount of metal deficiency of the protoplanetary disk as input into simulations of planet formation in disks (e.g., Johnson & Li 2012). However, halo stars in the literature can have values that differ by ∼ a factor of two, even when all studies use high-resolution/high signal-to-noise data (e.g., Jofré et al. 2014). Among the sources of systematic differences in the studies, the validity of using ionization equilibrium in LTE to determine gravity (Thévenin & Idiart 1999) is another area where Gaia information is crucial. Gaia parallaxes resolve the question of whether a star is a subgiant or main-sequence star. We therefore consider a selection of well-studied halo stars from the literature with new spectra and analyses in this work.
The paper is organized as follows: Section 2 outlines the process and criteria used to select a sufficiently pure sample of halo dwarf stars, while Section 3 offers verification of the success of the sample selection. Section 4 details the observation and data reduction process used to generate the spectra used in Section 5 to analyze the iron abundances of selected halo stars from the literature. The results of this analysis are discussed and compared with other papers in Section 6.
Stars were selected preferentially from the CTLv08 01 (Stassun et al. 2018(Stassun et al. , 2019) ) to have dwarf-like radii (raddflag = 1); have Gaia DR2 parallaxes and proper motions (PMFlag == 'gaia2'; PARFlag == 'gaia2'; gaiaqflag == 1); be within 1kpc according to the TICv8, which uses the Bailer-Jones et al. ( 2018) method to calculate distance (Stassun et al. 2019) (this is to ensure that the kinematic cut described below is valid, which was designed based on a simulation of the nearest 1kpc stars); and be outside of the ecliptic plane (|lat| > 6) to ensure candidates are observable by TESS. Note that, although the majority of targets from the CTL have T < 13, the Cool Dwarf list extends to T = 16 (for a description of this special list, see Muirhead et al. 2018).
Additional stars were selected by an equivalent search in the TICv8 (Stassun et al. 2018(Stassun et al. , 2019) ) for stars with 13 ≤ T < 15 (otherwise the same cuts as for the CTL are applied), effectively extending the magnitude limit of the halo dwarf search to T < 15.
To select stars with halo-like kinematics, we ran a Milky Way Galaxia (Sharma et al. 2011) simulation of stars within 1kpc. We designed a kinematic cut to choose halo stars based on only proper motion using this simulation, so that a Gaia radial velocity was not required. Our Galaxia simulations were run using the default parameters, which are described in Sharma et al. (2011). In particular, the Galaxia stellar halo geometry follows the Besançon model of an oblate spheroid, with density proportional to Max(ac,a) R n , where a 2 = R 2 + z 2 2 , a c = 500 pc, = 0.76, and n = -2.44 (Robin et al. 2003), with R and z denoting Galactocentric cylindrical radius and height above the plane in pc. Following Bond et al. (2010), the Galaxia halo velocity distribution is assumed to be an ellipsoid with σ R = 141 km s -1 and σ θ = σ φ = 75 km s -1 , which describes well the observed SDSS halo velocity distribution within 10kpc (Smith et al. 2009) as well that of Gaia halo stars (Posti et al. 2018).
The resulting kinematic selection was applied to the CTL/TIC sample, such that v δ < 0 km/s and either
Finally, stars were selected to fall below the binary sequence of a MIST 13 Gyr [Fe/H] = -1.0 (Dotter 2016;Choi et al. 2016) isochrone in J -K s -L space, where the luminosity, L, is adopted from the TIC. The final number of halo dwarf candidates selected in this way is 16940. The stars passing the kinematic cuts versus those passing the additional binary main sequence cut are shown in Figure 1. An excerpt from the final sample is displayed in Table 1. The entire catalog is available electronically in csv format.
Our sample selection can be tested by looking at the subset of stars with full 3-D kinematics or with metallicities. While the goal of the sample is to locate metal-poor stars, confirming that the proper-motion selected stars fall into kinematically hot populations when radial velocities are added is also valuable, given the correlation between the two properties.
We calculated the Galactic space velocity (Johnson & Soderblom 1987) for the 780 stars in the halo sample that had available Gaia DR2 radial velocity data and the 61 stars in the sample with APOGEE radial velocity data. Because radial velocity data are not used our selection procedure, we can further verify the validity of our selection with the galactic space velocity distribution of the subsample with radial velocity data. The Toomre diagram (Sandage & Fouts 1987) in Figure 2 shows that the galactic space velocities of these stars are all above 200 km/s. We note that the distribution is consistent with that from the literature that studied metal-poor halo stars, which have V tot > 180 km s -1 (e.g., Nissen & Schuster 2010;Schuster et al. 2012).
In addition to the 3-D kinematic verification of the sample, we can also directly measure the degree to which it is indeed metal-poor by appealing to independent metallicities from the literature, which have been compiled by the TIC.
We analyzed the distribution of the TIC metallicities, which are available for roughly 11% (n = 1,897) of the sample, yielded an overall purity estimate of 70%, the proportion of TIC metallicities which are below the low metallicity threshold of [Fe/H] = 1.0. We note that while the fraction of stars in the halo sample which have TIC metallicities is relatively small, the roughly 2,000 data points available should provide a metallicity distribution which is sufficiently representative of that of the entire halo sample. The distribution of these TIC metallicities is shown in Figure 3.
It should be noted the metallicities in the TIC/CTL are compiled from literature spectroscopic values, and are therefore heterogeneous. In order of preference, metallicities for the catalogue are adopted from: SPOCS (Brewer et al. 2016); PASTEL (Soubiran et al. 2016); Gaia-ESO DR3 (Gilmore et al. 2012); TESS-HERMES DR1 (Sharma et al. 2018); GALAH DR2 (Buder et al. 2018); APOGEE-2 DR14 (Abolfathi et al. 2018); LAM-OST DR4 (Luo et al. 2015); RAVE DR5 (Kunder et al. 2017); and Geneva-Copenhagen DR3 (Holmberg et al. 2009), as described in Stassun et al. (2019). While the TIC metallicities may not be self-consistent and may have a few tenths dex uncertainties, it provides a representative distribution of the halo dwarf candidate sample metallicities.
Some of our kinematically-selected halo stars have been previously known and studied (Reddy et al. 2006;Sozzetti et al. 2009;Boesgaard et al. 2011). However, there are discrepancies in derived stellar properties such as effective temperature and surface gravity, both among the literature, and when the literature is compared to the results of this paper. Table 2 shows effective tem-Figure 2. A projection of the stars' 3-dimensional space velocity onto a 2-dimensional plane, where the x-axis represents velocity in the direction of galactic rotation (V), and the y-axis represents the magnitude of the vector sum of radial velocity (U) and velocity perpendicular to the galactic plane (W). A star's total velocity is represented by its distance from the origin point (0,0). For this reason, iso-velocity lines are shown in gray for readability. All velocities are reported relative to the local standard of rest. Error bars are approximately 3 km/s along either axis and thus are not visible on the scale of the plot. • BD+51 1696, the star measured by the greatest number of studies we compared with, shows significant spread in both T eff and log(g) between studies.
• In many cases, Boesgaard et al. (2011) derives surface gravities characteristic of post-main-sequence stars, despite the selection process of this paper being designed to limit the sample to main sequence dwarfs.
This motivates us to observe 13 previously-known halo stars using the PEPSI spectrograph (Strassmeier et al. 2015) at the Large Binocular Telescope (LBT, 2 × 8.4 m on Mt. Graham, Arizona, USA). The higher spectral resolution (R=120,000) and a homogeneous abundance analysis, in addition to the external constraints from the Gaia DR2 data, help to reconcile the discrepancies and lead to a more accurate determination of stellar properties. We detail the PEPSI observations and our abundance measurements below.
LBT PEPSI observations were made on March 25, May 15, June 23 and 25, 2019 UT with 200 µm fiber (R = 130 000 and 1.75 on sky) in two spectral regions 4800 -5441 Å and 6278 -7419 Å (cross-dispersers (CD) 3 and 5) with 10-60 min integration time depending on the star brightness. A typical signal-to-noise ratio achieved is 350 in CD5 and 260 in CD3 for a star V=9.7 in one hour integration time with two LBT mirrors. The whole sample has signal-to-noise ratio ranging from 100 to 600.
The data reduction is done using the Spectroscopic Data Systems (SDS) with its pipeline adapted to the PEPSI data calibration flow and image specific content. Its description is given in Strassmeier et al. (2018).
The specific steps of image processing include bias subtraction and variance estimation of the source images, super-master flat field correction for the CCD spatial noise, échelle orders definition from the tracing flats, scattered light subtraction, wavelength solution using ThAr images, optimal extraction of image slices and cosmic ray removal, wavelength calibration and merging slices in each order, normalization to the master flat field spectrum to remove CCD fringes and blaze function, a global 2D fit to the continuum of the normalized image, and rectification of all spectral orders in the image to a 1D spectrum for a given cross-disperser.
The spectra from two sides of the telescope are averaged with weights into one spectrum and corrected for the barycentric velocity of the Solar system. The wavelength scale is preserved for each pixel as given by the wavelength solution without rebinning. The wavelength solution uses about 3000 ThAr lines and has an uncertainty of the fit at the image center of 4 m/s.
To determine the metallicities and stellar parameters of the stars observed with PEPSI, we first calculated the equivalent widths (EW) of absorption lines in each star's spectrum. This was done using an automated program which displays graphs and measurements from fits of a Voigt and Gaussian distribution, along with its direct measurement of the observed data. From this we chose the method which most closely fit the data trend. This manual screening allows for the best measurements to be kept for all lines, mitigating the effects of noise, contamination, and improper line fitting. If these effects are too great, the program allows us to discard the measurement from the final data.
To streamline this process, if all three EW measurements (Gaussian, Voigt, and raw) are within 20% of each other, usually the result of a clean, isolated line feature, the program automatically takes the average of the three and uses the result as that line's EW. On the other hand, if both curve fitting functions fail, the line is automatically discarded without prompting the user. This latter scenario usually occurs where a line is weak enough that it is completely drowned out by noise and is indistinguishable from the continuum.
Line information was compiled from the NIST Atomic Spectra Database (Kramida et al. 2020). A small section of the PEPSI spectrum for the star HD 160693 is shown in Figure 4 with line features annotated with their corresponding species and equivalent width.
For each star, we chose to keep only lines for which EW > 5m Å, to minimize the effects of noise on weak lines in our measurements. Furthermore, manual sigmaclipping based on each individual line's derived abundance was performed after the initial analysis to remove lingering outlying measurements, after which the process was repeated. This was done twice. A complete list of our final EW measurements can be found in Table 3.
To derive abundances and stellar atmospheric parameters, we used a program written by the authors which combines PyMOOGifoot_0 (Adamow 2017), which is a Python implementation of the Fortran code MOOGfoot_1 (Sneden 1973), with several Python functions designed to iteratively derive stellar parameters, abundances, and the uncertainties thereof.
We used the ATLAS9 model atmospheres computed for the APOGEE survey (Mészáros et al. 2012), and the PyKMOD atmosphere interpolatorfoot_2 to create the model stellar atmospheres for our analysis. To convert the resulting abundance values into quantities relative to solar, we took the solar iron abundance [Fe/H] = 7.48 as stated in Palme et al. (2014).
We used Gaia DR2 (Gaia Collaboration et al. 2018), 2MASS (Skrutskie et al. 2006), and WISE (Wright et al. 2010) photometry compared against Dartmouth theoretical isochronesfoot_3 (Dotter et al. 2007) to derive T eff and log(g). The photometry was queried from its respective databases using Astroquery, a module of the Astropy software library (Price-Whelan et al. 2018). We used distances computed by Bailer-Jones et al. ( 2018) to convert these values from apparent to absolute magnitude.
The process is as follows: the star's photometric profile is compared against the synthetic photometry for each age and equivalent evolutionary phase (EEP) at a fixed [Fe/H] from the Dartmouth isochrones. From this comparison, a 2-dimensional grid of residuals is created in the age-EEP plane. Thus, each point on the grid represents the sum of the differences of magnitudes at each pass band between those observed of the star and those from the given point on the isochrone.
From the absolute minimum point on this grid, we define the area around this point for which the residuals are within 10% of the minimum as valid points in the grid from which to extract T eff and log(g) values.
The averages of the resulting temperature and gravity distributions are taken as the parameters derived from the given metallicity, with the standard deviation of the mean being taken as the uncertainty.
To determine a value for the microturbulence parameter (ξ), we sought to remove the correlation between reduced equivalent width (REW) and abundance for Fe I lines. This was done by calculating the slope of this correlation for a fixed set of microturbulence values (ranging from 0 to 4 km/s), and linearly interpolating the microturbulence value at which the resulting slope is 0.
In the case where this process did not succeed (for example, if there is no point where the slope value crosses 0 along the range of microturbulence values), we chose to follow the convention used by Boesgaard et al. (2011), which is to simply set microturbulence equal to 1.5 km/s.
The process first assumes a metallicity of [Fe/H] = -1.0. From here, it gets T eff and log(g) values according to Section 5.2.1, and interpolates a model atmosphere from the grid with these photometric parameters and a microturbulence set to 1.5 km/s. Then, it runs MOOG with this model atmosphere and the line list for a given star.
It then reads the MOOG output, gathering the metallicity data and a value for [Fe/H], and then repeats, now using isochrones of the output [Fe/H] to derive new parameters for the model atmosphere. Convergence is reached and the process ends when the metallicity output by MOOG is the same (to within 0.1 dex, the precision of the isochrone grid) as that used to derive that run's stellar parameters.
Note that the initial metallicity choice is arbitrary. It is simply used as an initial guess for the iteration process, which lasts a variable amount of time based on how accurate this initial guess is.
After the process converges on a solution for [Fe/H], we calculate the microturbulence as outlined in Section 5.2.2, re-deriving stellar parameters as necessary as the microturbulence affects the metallicity. At this point, the program proceeds with the uncertainty analysis.
Uncertainty in the microturbulence parameter was determined by perturbing the microturbulence until the slope of the REW-abundance correlation fell outside of the uncertainty range.
All other uncertainties, including those for T eff , log(g), and abundances were calculated using iteration of the process outlined in Section 3.2 of Epstein et al. (2010). The method used by those authors creates a matrix of partial derivatives of parameters with respect to one another and uses various equations to account for the effects of uncertainties in each parameter on the uncertainties of all the others.
Since, by the methodology outlined in Section 5.2.1, the uncertainty in metallicity can affect the uncertainty in T eff and log(g), and by the methodology outlined here, the reverse is also true, we are required to iterate these calculations repeatedly, using the output uncertainties from one run as the input uncertainties of the next.
Initial uncertainties were simply taken as the standard deviation of the mean value derived for each parameter, and the calculations taken from Epstein et al. were repeated until the [Fe/H] uncertainty was changed by less than 0.01 dex from one iteration to the next.
We found the stars to fall within a temperature range of 5600 ≤ T eff ≤ 6800 and a surface gravity range of 4.1 ≤ log(g) ≤ 4.7. Detailed stellar parameter and abundance results for each star can be found in Table 4. The stars also fall well within the expected metallicity values given their halo classification, as we found them to fall within the range -2.0 < [Fe/H] < -0.6. One method of deriving stellar parameters involves using the spectroscopic abundance analysis, adjusting them until the following conditions are met:
1. Effective Temperature: The star should be in excitation equilibrium, i.e. the correlation between abundance derived from each line and excitation potential (EP) should be removed.
2. Surface Gravity: The star should be in ionization balance, i.e. the abundances of Fe I and Fe II should be within 1σ of each other.
3. Microturbulence: There should be no correlation between the abundance derived from each line and reduced equivalent width (log( EW λ )). We followed this method for microturbulence, but our T eff and log(g) values were derived photometrically from isochrones. However, for every star except BD+75 839, the resulting log(g) led to ionization balance.
Also, the resulting T eff led to excitation equilibrium for many of the stars, such that the slope of the EPabundance regression line was equal to 0 at the 1-sigma level for four stars, and at the 3-sigma for seven, where sigma is defined as the standard deviation of the linear regression fit assuming normal distribution of residuals.
Two stars (HD 64090 and HD 160693) featured slopes which were 4-sigma away from zero. Notably however, these two stars were also unable to achieve microturbulence convergence.
For the six stars in common with Reddy et al. (2006), the mutual discrepancy between metallicities is on average 0.06 dex, which we consider to be reasonably identical within a margin of error. Their process of determining stellar parameters was similar to ours in that it also does not rely on spectroscopy for T eff and log(g). They use (b -y) and (V -K s ) to derive T eff , whereas we use magnitudes, rather than colors, and use G, Bp, Rp, J, H, K s , W 1, W 2, W 3, and W 4.
For surface gravities, they make use of Hipparcos astrometry whereas we use the same photometric method for log(g) as was used for T eff . And lastly, for microturbulence, they use a previously-calculated relation between other parameters and ξ for metal-poor dwarfs, where we in this case use the spectroscopic approach for microturbulence.
We also compared our results with those for the seven stars our sample shares with Boesgaard et al. (2011). In the paper, they determined the stellar parameters spectroscopically according to the conditions in Section 6.1. Their metallicities are systematically lower than ours by an average of 0.19 dex, except in the case of BD+51 1696, where theirs is higher by 0.23 dex.
Although similar stellar parameter convergence conditions (ionization balance, excitation equilibrium, microturbulence convergence) were met to a degree of uncertainty in both papers, it can be seen in Tables 2 and4 that the resulting parameters are in many cases quite different, especially surface gravity.
From these discrepancies we can assume the importance of the additional photometric constraints on the parameters.
Based on the available information, it seems that the spectroscopic method can lead to cases where there are multiple points in the T eff -log(g)-ξ parameter space where the convergence conditions are met. This leaves the chance that, if additional constraints are not used, any given analysis may not necessarily converge on the correct point, leading to parameters which conflict with other stellar data.
We investigated the magnitude of the effects of the LTE assumption on our final abundance measurements. We accessed data from Bergemann et al. (2012) via a web tool by M. Kovalev et al. (2018), and discovered that all corrections for lines we tested were within our metallicity uncertainty thresholds, thus making them negligible corrections to our measurements.
In this paper we discuss our selection process for metal-poor halo stars, detailing kinematic and photometric criteria resulting in a final sample of ∼ 16, 940 stars, where [Fe/H] < -1.0 for roughly 70% of the stars included based on comparison to literature metallicities.
We also present a re-analysis of 13 halo stars from the literature using new observations taken with the highresolution PEPSI spectrograph, attempting to rectify discrepancies of stellar properties in the literature. We measure the metallicity of these stars to an accuracy of roughly σ([Fe/H]) = 0.1 dex. In this process, we also used Gaia, 2MASS, and WISE photometry to derive accurate effective temperature and surface gravity values using Dartmouth theoretical isochrones.
In summary, given the overall fidelity of our sample, the halo dwarf candidates presented here will prove to be useful targets for planet studies in metal-poor host systems with TESS. An analysis of the sample in the context of planet occurrence rates in this metal-poor regime is described in Boley et al. 2021, the companion to this paper.
https://github.com/madamow/pymoogi
http://www.as.utexas.edu/ ∼ chris/moog.html
https://github.com/kolecki4/PyKMOD
http://stellar.dartmouth.edu/models/