Validation of the Gaia Early Data Release 3 Parallax Zero-point Model with Asteroseismology

The Gaia mission has provided astrometric information for over 1.5 billion stars as part of Gaia Early Data Release 3 (EDR3), which is largely complete down to G ∼ 21 in uncrowded regions (Gaia Collaboration et al. 2012). This release successfully builds upon the previous Data Release 2 (DR2; Gaia Collaboration et al. 2016, 2018), with improvements to the mission's angular resolution, completeness, and astrometric precision. In particular, the reported parallax precision in EDR3 has improved by a factor of 2 from ≈40 μas for sources with G < 15 to ≈20 μas.

Because of the nominal increase in parallax precision, it is all the more important to understand systematic uncertainties in the Gaia parallaxes. Following indications of systematic errors in the Gaia Data Release 1 parallaxes from the Gaia team and from subsequent independent investigations (Michalik et al. 2015; De Ridder et al. 2016; Gaia Collaboration et al. 2016; Jao et al. 2016; Lindegren et al. 2016; Stassun & Torres 2016; Davies et al. 2017; Huber et al. 2017; Zinn et al. 2017), several studies investigated the Gaia zero point in Gaia DR2. Many of the studies investigated the level of the global offset, which is due to degeneracies between variations in the angular separation of the spacecraft's fields of view and the parallax zero point of the astrometric solution (Butkevich et al. 2017; Lindegren et al. 2018). These studies required independent parallax estimates to compare the Gaia parallaxes against and ranged from independent trigonometric parallaxes (Leggett et al. 2018), classical Cepheid photometric parallaxes (Groenewegen 2018; Riess et al. 2018), RR Lyrae photometric parallaxes (Muraveva et al. 2018; Layden et al. 2019; Marconi et al. 2021), open cluster/OB association isochronal parallaxes (Yalyalieva et al. 2018; Melnik & Dambis 2020; Sun et al. 2020), very long baseline interferometric parallaxes (Kounkel et al. 2018; Bobylev 2019), statistical parallaxes (Schönrich et al. 2019; Muhie et al. 2021), and red giant asteroseismic parallaxes (Hall et al. 2019; Khan et al. 2019).

Apart from the 10–100 μas global offset inferred by the aforementioned studies, the Gaia team also identified position-, color- and magnitude-dependent trends in the Gaia zero point, which were thought to result from Gaia's scanning pattern and CCD response (Arenou et al. 2018; Lindegren et al. 2018). Independent studies subsequently confirmed similar trends (Leung & Bovy 2019; Zinn et al. 2019a, 2019b; Chan & Bovy 2020; Fardal et al. 2021).

The global offset and position-, color-, and magnitude-dependent parallax systematics quantified in DR2 are also present, though to a lesser extent, in EDR3 (Lindegren et al. 2020a, 2020b). The reduction in systematics is a result of EDR3 benefitting from a new astrometric solution using 34 months of data as opposed to DR2's 22 months; improvements in EDR3 to the Velocity error and effective Basic Angle Calibration (VBAC) model, which models the basic angle variations that contribute to the global offset component of the parallax zero point; as well as improvements in the photometric image parameter determination (Rowell et al. 2020) and its iterative inclusion in the astrometric solution (Lindegren et al. 2020b).

Whereas the Gaia team did not recommend a specific parallax zero-point correction in DR2, a model for the parallax zero point has been provided for EDR3. This model is fit according to an iterative solution based ultimately on a sample of quasars distributed across the sky (Lindegren et al. 2020a), for which the EDR3 parallaxes should be effectively zero. The resulting model, denoted Z₅ or Z₆, depending on whether the parallax is from a five-parameter or six-parameter astrometric solution, ² once subtracted from the raw EDR3 parallaxes, is meant to remove most magnitude-, color-, and position-dependent errors.

The Gaia team has checked the Z₅ model using stars in the Large Magellanic Cloud (LMC), open cluster members, and wide binaries (Fabricius et al. 2020; Lindegren et al. 2020a), finding good performance across a range of magnitude and color parameter space.

Since the EDR3 release, independent analyses have begun to quantify what adjustments may be required of the Z₅ model. Riess et al. (2021), using a sample of classical Cepheids (G ≲ 11, 1.35 μm⁻¹ ≲ ν_eff ≲ 1.6 μm⁻¹), infer an offset of −14 ± 6 μas such that corrected Gaia parallaxes are larger than Cepheid photometric parallaxes. Bhardwaj et al. (2021), appealing to blue RR Lyrae (ν_eff > 1.5 μm⁻¹), find an offset with corrected Gaia parallaxes of −25 ± 5 μas in the same direction. Stassun & Torres (2021), using 76 bright, blue eclipsing binaries (5 ≲ G ≲ 12, ν_eff > 1.5 μm⁻¹), show no statistically significant parallax residuals after correction according to the Gaia parallax model (+15 ± 18 μas). Moreover, an analysis using photometric parallaxes of red clump stars has recently indicated parallax residuals for G < 10.8 and as a function of ecliptic longitude (Huang et al. 2021).

There is, however, nearly a complete absence of stars with 1.4 μm⁻¹ < ν_eff < 1.5 μm⁻¹, G < 13 in either the LMC sample used for validation of the Z₅ solution (Lindegren et al. 2020a) or any of the independent validation test samples thus far. This range in magnitude space is especially interesting to explore given that it links the G > 13 checks at similar colors with the LMC by the Gaia team to the G < 11 study from Riess et al. (2021). I therefore consider in this paper the evidence for adjustments to the Z₅ model by appealing to asteroseismic parallaxes of more than 2000 first-ascent red giant branch stars, which occupy this region of magnitude and color parameter space. Specifically, I look for errors in Gaia parallaxes corrected according to the Gaia Z₅ parallax zero-point model as a function of color and magnitude. In so doing, I also consider the evidence for corrections to the formal statistical uncertainties in Gaia parallaxes, as well as evidence for spatially correlated uncertainties in Gaia parallaxes.

Asteroseismic data are adopted from APOKASC-2 (Pinsonneault et al. 2018), which consists of asteroseismic parameters for red giant branch stars, ${\nu }_{\max }$ and Δν, derived from light curves acquired by the Kepler mission (Borucki et al. 2010), as well as spectroscopic temperatures and metallicities from the Apache Point Observatory Galactic Evolution Experiment (APOGEE; Majewski et al. 2010), a survey within Sloan Digital Sky Survey IV (Blanton et al. 2017). I also appeal to an independent analysis of Kepler data from Yu et al. (2018) using the SYD asteroseismic pipeline (Huber et al. 2009). Whereas the APOKASC-2 asteroseismic measurements are average measurements based on results from five asteroseismic pipelines, using data from SYD alone provides a check on the sensitivity of the result to the asteroseismic data.

I make use of APOGEE DR16 (Ahumada et al. 2020) parameters in this work, updated from the APOGEE DR14 (Holtzman et al. 2018) parameters provided in Pinsonneault et al. (2018). Both APOGEE DR14 and DR16 data are taken using the R = 22,500 APOGEE spectrograph on the 2.5 m telescope of the Sloan Digital Sky Survey (Gunn et al. 2006; Wilson et al. 2019). The data reduction pipeline for APOGEE is described in Nidever et al. (2015), and the spectroscopic analysis is performed using the APOGEE Stellar Parameter and Chemical Abundance Pipeline (ASPCAP; García Pérez et al. 2016). The APOGEE temperatures are calibrated to be on the infrared flux method scale of González Hernández & Bonifacio (2009; Holtzman et al. 2015), which, for the red giant sample I work with here, implies temperatures are on an absolute scale to within ≈20 K (Zinn et al. 2019b); I adopt statistical uncertainties of 30 K.

I cross-matched the APOKASC-2 sample with EDR3 by first matching the APOKASC-2 stars to Gaia DR2 sources using the DR2-2MASS cross-match (Marrese et al. 2019), and then using the Gaia DR2 designations to match to EDR3 sources using the EDR3-DR2 cross-match (Torra et al. 2020). In so doing, I only kept stars that have angular separations between DR2 and EDR3 sources of less than 100 mas.

Most of the APOKASC-2 stars have five-parameter solutions instead of six-parameter solutions (44 versus 2130 first-ascent red giant branch stars after quality cuts), so the following analysis only uses stars with five-parameter astrometry solutions and thus concentrates on the validation of the Z₅ model.

To ensure that the astrometric solutions are not affected by binarity, I followed the quality cuts of Fabricius et al. (2020), rejecting stars with ruwe ≥ 1.4, ipd_frac_multi_peak > 2, and ipd_gof_harmonic_amplitude ≥ 0.1. ipd_frac_multi_peak corresponds to the fraction of observations for which the source was identified as having two resolved peaks in the image, and therefore is indicative of resolved binaries; ipd_gof_harmonic_amplitude indicates the level of variation in the image goodness-of-fit as a function of scan direction, which, if large, would imply the image is asymmetric and therefore suggestive of an unresolved binary; ruwe indicates that the astrometric solution does not completely describe the motion of the source, and so can therefore identify unresolved binaries with variable photocenters (Lindegren et al. 2020b).

I further restricted the sample to those with G < 13, which are observed within window classes WC0a and WC0b (Lindegren et al. 2020b). These magnitude-defined windows define what pixel mask the source is read out with, and WC0 was divided into two for EDR3: sources with 11 ≲ G ≲ 13 are observed with the WC0b window and those with G ≲ 11 with the WC0a window. Magnitude-dependent systematics in the astrometry may therefore be related to the different behavior of photometry in the different window classes (see, e.g., Figure A.4 in Lindegren et al. 2020b). The G < 13 regime is further complicated, however, by the use of differing integration times through the use of "gates," and which themselves depend on magnitude. I will investigate the possibility of magnitude-dependent errors in corrected Gaia parallaxes in what follows.

I adopt nu_eff_used_in_astrometry (denoted ν_eff here) as a proxy for color, following Lindegren et al. (2020a), and which is defined based on DR2 photometry for the five-parameter astrometry I use in this analysis (Lindegren et al. 2020b); bluer stars have larger ν_eff, and redder stars have smaller ν_eff. I did not make cuts in color, because the sample lies well within the regime 1.24 μm⁻¹ < ν_eff < 1.72 μm⁻¹; outside of this range, the image chromaticity point-spread function and line-spread function corrections are not calibrated, and so there is reason to believe they would suffer from a stronger color-dependent zero-point error (Lindegren et al. 2020a). I will also fit for color-dependent terms to describe residual differences between asteroseismic and Gaia parallaxes in what follows.

The parallax of a star can be derived from the asteroseismic radius via the Stefan–Boltzmann law, written in the following form:

$$\begin{eqnarray}&&\begin{array}{l}{\varpi }_{\mathrm{seis}}({T}_{\mathrm{eff}},F,R)={F}^{1/2}{\sigma }_{\mathrm{SB}}^{-1/2}{T}_{\mathrm{eff}}^{-2}{R}^{-1}\\ \quad =\ {f}_{0}^{1/2}{10}^{-1/5(m+{{BC}}_{b}({T}_{\mathrm{eff}})-{A}_{b})}{\sigma }_{\mathrm{SB}}^{-1/2}{T}_{\mathrm{eff}}^{-2}{R}^{-1},\end{array}\end{eqnarray} \tag{ 1 }$$

where R is the asteroseismic radius (see below), T_eff is the effective temperature, σ_SB is the Stefan–Boltzmann constant, and F is the stellar bolometric flux. Here, I have rewritten the flux in terms of a bolometric correction in the photometric band b, BC_b ; the extinction in that passband, A_b ; and a magnitude-flux conversion factor that assumes the solar irradiance of f₀ = 1.361 × 10⁶ erg s⁻¹ cm⁻² (Mamajek et al. 2015), and an apparent solar bolometric magnitude of m_bol = −26.82 (Torres 2010). I work here with 2MASS K_s photometry (Skrutskie et al. 2006) to reduce extinction effects. I use a bolometric correction based on MIST models (Paxton et al. 2011, 2013, 2015, 2018, 2019) as implemented in isoclassify (Huber et al. 2017; Berger et al. 2020). Visual extinctions were adopted from Rodrigues et al. (2014) and converted into ${A}_{{K}_{{\rm{s}}}}$ assuming ${A}_{{K}_{{\rm{s}}}}=0.113\,{A}_{V}$ (Schlafly & Finkbeiner 2011).

The asteroseismic radii required to yield parallaxes according to Equation (1) can be computed according to a scaling relation,

$$\begin{eqnarray}&&\displaystyle \frac{R}{{R}_{\odot }}\approx \left(\displaystyle \frac{{\nu }_{\max }}{{\nu }_{\max ,\odot }}\right){\left(\displaystyle \frac{{\rm{\Delta }}\nu }{{f}_{{\rm{\Delta }}\nu }{\rm{\Delta }}{\nu }_{\odot }}\right)}^{-2}\ {\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{T}_{\mathrm{eff},\odot }}\right)}^{1/2},\end{eqnarray} \tag{ 2 }$$

where ${\nu }_{\max }$ is an asteroseismic observable that corresponds to the frequency at which stochastically excited oscillations in the stellar envelopes have their highest amplitude and Δν is an asteroseismic observable that describes the approximately constant separation in frequency between two oscillation modes of adjacent radial order but the same spherical degree. The former is theoretically and empirically related to stellar surface gravity (Brown et al. 1991; Kjeldsen & Bedding 1995; Chaplin et al. 2008; Belkacem et al. 2011). The latter is related to the mean stellar density (Ulrich 1986; Kjeldsen & Bedding 1995), though there is evidence that the observed Δν requires a correction (e.g., White et al. 2011; Guggenberger et al. 2016; Sharma et al. 2016), which is represented in the above by f_Δν . Δν_⊙ and ${\nu }_{\max ,\odot }$ represent measurements of these quantities for the Sun, to which the scaling relations are tied. The solar reference values of Δν_⊙ = 135.146 μHz and ${\nu }_{\max ,\odot }=3076\,\mu \mathrm{Hz}$ that I adopt from Pinsonneault et al. (2018) are calibrated such that asteroseismic masses agree with dynamical masses from eclipsing binaries in NGC 6791 (Grundahl et al. 2008; Brogaard et al. 2011, 2012) and NGC 6819 (Jeffries et al. 2013; Sandquist et al. 2013; Brewer et al. 2016). f_Δν , also adopted from Pinsonneault et al. (2018), amounts to ≈3% corrections of Δν, and which vary on a star-by-star basis according to mass, surface gravity, temperature, and metallicity.

In this work, I also consider effectively correcting both Δν and ${\nu }_{\max }$ using nonlinear scaling relations from Kallinger et al. (2018), which are derived from fits to dynamical surface gravities and mean stellar densities in addition to the open clusters NGC 6791 and NGC 6819:

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{R}{{R}_{\odot }}\approx {\left(\displaystyle \frac{{\nu }_{\max }}{{\nu }_{\max ,\mathrm{ref}}}\right)}^{\kappa }{\left(\displaystyle \frac{{\rm{\Delta }}\nu }{{\rm{\Delta }}{\nu }_{\mathrm{ref}}}\right)}^{-2}\\ \quad \times \,{\left[1.0-\gamma {\left({\mathrm{log}}_{10}\displaystyle \frac{{\rm{\Delta }}\nu }{{\rm{\Delta }}{\nu }_{\mathrm{ref}}}\right)}^{2}\right]}^{2}\ {\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{T}_{\mathrm{eff},\odot }}\right)}^{1/2},\end{array}\end{eqnarray} \tag{ 3 }$$

where κ = 1.0075 ± 0.0021, γ = 0.0043 ± ± 0.0025, and I use the Δν_ref value from Kallinger et al. (2018) appropriate for the method of deriving APOKASC-2 Δν values, Δν_ref = 133.1 ± 1.3 μHz. I adopt ${\nu }_{\max ,\mathrm{ref}}=3076\,\mu \mathrm{Hz}$ (Pinsonneault et al. 2018). In what follows, I marginalize over errors in the asteroseismic radii due to, e.g., the choice of ${\nu }_{\max ,\mathrm{ref}}$ or Δν corrections, by fitting for a multiplicative factor that brings the radii into agreement with the Gaia parallaxes.

I conservatively restricted the analysis to stars with R < 30R_⊙, based on findings from Zinn et al. (2019b) that asteroseismic radii are likely inflated in the evolved red giant regime. I also restricted the analysis to the first-ascent red giant branch stars with $15\,\mu \mathrm{Hz}\lt {\nu }_{\max }\lt 200\,\mu \mathrm{Hz}$ and 2 μHz < Δν < 10 μHz, which is the range occupied by the majority of the open cluster calibrators used by Pinsonneault et al. (2018) and Kallinger et al. (2018) to set the ${\nu }_{\max }$/Δν solar reference values/corrections. For parts of the following analysis using the Yu et al. (2018) data set, I use the nonlinear scaling relations of Equation (3), because the APOKASC-2 Δν corrections (f_Δν in Equation (2)) assume the ${\nu }_{\max }$ and Δν of the APOKASC-2 catalog.

The approach I take to investigate the need, if any, for adjustments to the Z₅ five-parameter solution parallax zero-point model is to compare Gaia parallaxes corrected using Z₅ to asteroseismic parallaxes. By taking the difference between the corrected Gaia parallaxes and the asteroseismic parallaxes, which I refer to as the parallax residuals in what follows, one can simultaneously constrain asteroseismic parallax problems due to asteroseismic radius and any additive problems in the Gaia parallax that may remain after correction according to Z₅. The simultaneous calibration is possible because errors in the asteroseismic parallax due to the asteroseismic radius will be fractional and therefore dependent on parallax, given the multiplicative dependence of asteroseismic parallax on the asteroseismic radius (Equation (1)). By contrast, Gaia parallax zero-point problems are found to be additive and can be described by terms that depend on magnitude, position, and color (Lindegren et al. 2020a).

Formally speaking, the model for describing any adjustments to the Z₅ model can therefore be described as a likelihood of the form

$$\begin{eqnarray}&&\begin{array}{l}{ \mathcal L }(c,d,e,a,k,g,s| {\hat{\varpi }}_{\mathrm{Gaia}},{T}_{\mathrm{eff}},{\rm{\Delta }}\nu ,{\nu }_{\max },{A}_{{K}_{{\rm{s}}}},\\ \quad {K}_{{\rm{s}}},{BC},{\nu }_{\mathrm{eff}},G,\sin \beta )\\ \quad \propto \,\displaystyle \frac{1}{\sqrt{{\left(2\pi \right)}^{N}| {C}^{{\prime} }| }}\exp \left[-\displaystyle \frac{1}{2}{\left({\boldsymbol{y}}-{\boldsymbol{x}}\right)}^{{\rm{T}}}{C}^{{\prime} -1}({\boldsymbol{y}}-{\boldsymbol{x}})\right].\end{array}\end{eqnarray} \tag{ 4 }$$

In the above,

$$\begin{eqnarray*}&&{\boldsymbol{y}}\equiv Y(a,{d}_{a},{\nu }_{\mathrm{eff}}){\varpi }_{\mathrm{seis}}({T}_{\mathrm{eff}},{\rm{\Delta }}\nu ,{\nu }_{\max },{A}_{{K}_{{\rm{s}}}},{K}_{{\rm{s}}},{BC})\end{eqnarray*}$$

and

$$\begin{eqnarray*}&&{\boldsymbol{x}}\equiv {\hat{\varpi }}_{\mathrm{Gaia}}-{Z}_{5}({\nu }_{\mathrm{eff}},G,\sin \beta )+{\rm{\Delta }}Z(c,d,e,{\nu }_{\mathrm{eff}},G),\end{eqnarray*}$$

where

$$\begin{eqnarray}&&Y(a,{d}_{a},{\nu }_{\mathrm{eff}})\equiv 1-a+{d}_{a}({\nu }_{\mathrm{eff}}-1.48)\end{eqnarray} \tag{ 5 }$$

describes fractional errors in asteroseismic parallaxes: a describes errors in the asteroseismic radius scale due to, e.g., solar reference value choice, and would be unity in the absence of any. As I discuss in Section 4.1, the term d_a describes color-dependent corrections to the asteroseismic parallaxes due to, e.g., bolometric correction systematics that are a function of temperature/color.

ΔZ is the model for the residual parallax difference left unexplained by the Z₅ model. ${\hat{\varpi }}_{\mathrm{Gaia}}$ is the raw EDR3 parallax, $\sin \beta$ is the sine of ecliptic latitude, and Z₅ is evaluated using the Python implementation of the correction, zero_point. ³ I will only be using Gaia parallaxes corrected in this way in the analysis, which I denote ϖ_Gaia, as opposed to the raw Gaia parallaxes, denoted ${\hat{\varpi }}_{\mathrm{Gaia}}$.

I consider several different forms for ΔZ, allowing for a local offset to the corrected Gaia parallaxes as well as color- and magnitude-dependent terms. Because I only analyze stars in the Kepler field, it is not possible to constrain ecliptic latitude-dependent terms as the Gaia team has done for the Z₅ model (though it is possible to statistically infer the level of small-scale variations as a function of angular scale, as I discuss below and in Section 4.7). At its most complicated, the model for the parallax residuals takes the form

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}Z(c,d,e,{\nu }_{\mathrm{eff}},G)=\left\{\begin{array}{l}c+{c}_{2}+{d}_{1}({\nu }_{\mathrm{eff}}-1.48)\\ +{d}_{3}{\left(1.48-{\nu }_{\mathrm{eff}}\right)}^{3}+{e}_{1}(G-12.2)\\ -{\rm{\Delta }}{Z}_{G=13},\\ \mathrm{for}\ G\lt 10.8\\ c+{d}_{2}({\nu }_{\mathrm{eff}}-1.48)\\ +{d}_{3}{\left(1.48-{\nu }_{\mathrm{eff}}\right)}^{3}+{e}_{2}(G-12.2)\\ -{\rm{\Delta }}{Z}_{G=13},\\ \mathrm{for}\ G\geqslant 10.8.\end{array}\right.\end{array}\end{eqnarray} \tag{ 6 }$$

I consider nested models under various permutations of the parameters, c, d, e, etc., as listed in Table 1; a blank entry indicates an unused parameter, such that it can be considered to be zero in Equation (6). In some models, there is a single color term across all magnitudes, i.e., d₁ ≡ d₂ ≡ d. Models otherwise differ by removing any number of terms by setting them to zero (e.g., the preferred model, Model 0, has a single color term and no magnitude term, i.e., d₁ ≡ d₂ ≡ d and e₁ ≡ e₂ ≡ e = 0). The term ΔZ_G=13 is not a fitted parameter, but rather a constant defined such that c describes the mean offset at G = 13 and ν_eff = 1.48 μm⁻¹ remaining after correction with d and/or d₃, viz., ΔZ_G=13 ≡ 0.8e₂ μas. c therefore can be interpreted as the average local offset of EDR3 parallaxes with respect to the rest of the sky/quasar frame of reference that the Z₅ model was calibrated to. The pivot point of ν_eff = 1.48 μm⁻¹ is chosen to be consistent with the q_1k and q_2k color terms in the Z₅ model, which take the form (ν_eff − 1.48) and ${(1.48-{\nu }_{\mathrm{eff}})}^{3}$ (Lindegren et al. 2020a), respectively. The pivot point of G = 10.8 is motivated by the piecewise functions in magnitude used to describe the magnitude dependence of Z₅ and which have a breakpoint at 10.8. This is also in the transition region of G ≈ 11 between the window classes WC0a and WC0b. The faintest breakpoint in the Z₅ model that overlaps with the sample's magnitude range occurs at G = 12.2, which sets the pivot point for the e term.

Note that one is able to simultaneously constrain (1) fractional errors in the asteroseismic radii, which would naturally arise from problems in the radius-scaling relation (Equation (3)) entering into the asteroseismic parallaxes (Equation (1)) as parameterized by the a and d_a terms in Equation (5), as well as (2) additive errors in the corrected Gaia parallaxes, as parameterized by the c, d, and e terms in Equation (6).

The elements of the covariance matrix in Equation (4) are given by

$$\begin{eqnarray}&&\begin{array}{l}{C}_{{ij}}^{{\prime} }={C}_{{ij}}+{k}^{2}{\delta }_{{ij}}{\sigma }_{{\varpi }_{\mathrm{Gaia}}}^{2}+\left[{g}^{2}+{\left(\displaystyle \frac{\partial Y}{\partial {\varpi }_{\mathrm{seis}}}\right)}^{2}\right]{\delta }_{{ij}}{\sigma }_{{\varpi }_{\mathrm{seis}}}^{2}\\ \quad +\,\left[{\left(\displaystyle \frac{\partial Y}{\partial {\nu }_{\mathrm{eff}}}\right)}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Delta }}Z}{\partial {\nu }_{\mathrm{eff}}}\right)}^{2}\right]{\delta }_{{ij}}{\sigma }_{{\nu }_{\mathrm{eff}}}^{2}+{\left(\displaystyle \frac{\partial {\rm{\Delta }}Z}{\partial G}\right)}^{2}{\delta }_{{ij}}{\sigma }_{G}^{2}\\ \quad +\,{\delta }_{{ij}}{s}^{2},\end{array}\end{eqnarray} \tag{ 7 }$$

and describe the covariance between the asteroseismic–Gaia parallax difference for two stars, i and j, where δ_ij is the Kronecker delta function. In the above, k and g describe corrections to the formal statistical uncertainties of Gaia and asteroseismic parallaxes, which can be constrained because the two uncertainties are not strongly correlated (see Section 4.6). s can be thought of as an intrinsic scatter in the parallax residuals that capture any variation not described by the model, or, alternatively, as an additive correction to the Gaia and/or asteroseismic parallax uncertainties. σ_G is the uncertainty on G, and ${\sigma }_{{\nu }_{\mathrm{eff}}}$, the uncertainty in nu_eff_used_in_astrometry, is adopted to be 1% of ν_eff, because the latter is a fixed and not fitted quantity for five-parameter astrometric solutions—this is ≈2 times larger than the uncertainty in astrometric_pseudo_colour, which is a fitted color as part of the six-parameter solutions. The statistical uncertainties in the asteroseismic parallaxes are denoted ${\sigma }_{{\varpi }_{\mathrm{seis}}}$ and include contributions from Δν, ${\nu }_{\max }$, T_eff, ${A}_{{K}_{{\rm{s}}}}$, K_s, and BC via linear propagation of uncertainty. The Gaia spatially correlated errors are encapsulated in C_ij , which describes the spatial covariance in Gaia parallax between stars i and j separated by an angular distance, θ_ij .

By appealing to quasars in EDR3, Lindegren et al. (2020b) fit an exponential function to the parallax spatial covariance for θ_ij ≳ 05 of ${C}_{{ij},\mathrm{QSO}}={\rho }_{\mathrm{QSO}}{e}^{\left(-\mathrm{ln}2{\theta }_{{ij}}/{\theta }_{1/2,\mathrm{QSO}}\right)}$, with ρ_QSO = 142 μas² and θ_1/2,QSO ≈ 11°.

In this work, I derive an estimate of the spatial covariance matrix (Section 4.7), modeling it with an equation of the same form as C_ij,QSO:

$$\begin{eqnarray}{C}_{{ij}}=\left\{\begin{array}{ll}\rho {e}^{\left(-\mathrm{ln}2{\theta }_{{ij}}/{\theta }_{1/2}\right)}, & \mathrm{for}\ 0^\circ \lt {\theta }_{{ij}}\leqslant 10^\circ \\ 0, & \mathrm{otherwise},\end{array}\right.\end{eqnarray} \tag{ 8 }$$

with best-fitting parameters according to Table 2. Note that θ_1/2 is defined analogously to the half-angle in the Gaia team description of the covariance, θ_1/2,QSO. I tested the impact of spatial covariance in the likelihood analysis using an approximation described in Zinn et al. (2017). Briefly, I divided the Kepler field into ∼25 × 25 squares corresponding to the Kepler modules, treating each one as independent from the other, given that spatially correlated Gaia parallax errors are ≈1% the level of statistical uncertainties on scales larger than 3° (see Section 4.7). The results thus accounting for spatial correlations are not significantly different from those without spatial correlations, and so I ignore the spatial covariance term in Equation (7), C_ij , in what follows. I describe how I infer the spatial covariance in Gaia parallaxes and present the fit to Equation (8) in Section 4.7.

Table 2. Gaia EDR3 Parallax Spatially Correlated Systematic Uncertainties

ρ	θ1/2	χ2/dof	N
${15}_{-9}^{+9}\,\mu {\mathrm{as}}^{2}$	${0.49}_{-0.17}^{+0.17}\deg$	2.9	37

$\sqrt{{C}_{{ij}}({\theta }_{{ij}}={\theta }_{1/2})}$ a	$\sqrt{{C}_{{ij}}({\theta }_{{ij}}=0^\circ )}$ a	$\sqrt{{C}_{{ij}}({\theta }_{{ij}}=0\buildrel{\circ}\over{.} 5)}$ a	$\sqrt{{C}_{{ij}}({\theta }_{{ij}}=1^\circ )}$ a	$\sqrt{{C}_{{ij}}({\theta }_{{ij}}=5^\circ )}$ a
${2.5}_{-0.8}^{+0.6}\,\mu \mathrm{as}$	${3.9}_{-1.4}^{+1.1}\,\mu \mathrm{as}$	${2.5}_{-0.8}^{+0.6}\,\mu \mathrm{as}$	${1.6}_{-0.6}^{+0.5}\,\mu \mathrm{as}$	${0.1}_{-0.1}^{+0.1}\,\mu \mathrm{as}$

Notes. Best-fitting parameters for spatially correlated uncertainties of Gaia parallaxes, as inferred in the Kepler field and parameterized by Equation (8). ρ describes the variance at small scales. θ_1/2 is the characteristic angular scale for the spatially correlated uncertainties.

^a The square root of the covariance, C, is a measure of the systematic uncertainty in the Gaia parallaxes due to spatial correlations.

Download table as:  ASCIITypeset image

I fit the model of Equation (4) using Markov Chain Monte Carlo (MCMC), rejecting from analysis stars whose Gaia and asteroseismic parallaxes disagree by more than 2.5σ. Figure 1 shows the resulting posterior distributions of the parameters for the preferred model (Model 0 in Table 1). The best-fitting parameters are provided in Table 1, which are taken to be the means of the posteriors; the uncertainties are taken to be the standard deviations of the posteriors. I also provide the Bayesian Information Criterion difference (ΔBIC; Schwarz 1978) for all models that I considered compared to the preferred model, where a smaller value indicates stronger evidence for the model, and a difference in 6 between models is taken to be a strong preference (Kass & Raftery 1995). In what follows, I take a conservative approach to potential adjustments to Z₅, by default assuming terms in ΔZ are null unless there is strong evidence for them.

Zoom In Zoom Out Reset image size

Before correction by Z₅, the raw Gaia and asteroseismic parallaxes have a mean difference of +22 μas (scatter of 23 μas), in the sense that Gaia parallaxes are smaller than asteroseismic parallaxes. Even before correction according to Z₅, this is a significant improvement over the ≈+50 μas underestimation of DR2 parallaxes compared to asteroseismic parallaxes for this sample of red giant branch stars (Zinn et al. 2019a). For stars with G > 10.8, correcting the Gaia parallaxes according to Z₅, adjusting the asteroseismic radii with a, and removing the local offset unique to the Kepler field c, yields a residual of +9 μas, which is due to a color trend. I note that a nonzero color term, d, is strongly preferred to fit the parallax residuals. However, I take caution in interpreting this term as solely due to adjustments needed of Z₅. As I explain in Section 4.1, the color term as fitted in Model 0, d, may have contributions due to color-dependent asteroseismic parallax errors (d_a in Equation (5)), for which there is not strong enough evidence to confirm. Under a conservative assumption, removing finally the color term implies the Z₅ model leaves an insignificant residual of +0.3 μas, with uncertainty in the mean of this residual of 0.4 μas. For stars with G ≲ 11, the Z₅ model appears to overcorrect the parallaxes by −15 ± 3 μas, as I discuss in Section 4.4.

The parallax residuals (asteroseismic – Z₅-corrected Gaia) modeled here are shown as a function of asteroseismic parallax and various other parameters in Figures 2 and 3. The observed difference between corrected Gaia parallaxes and asteroseismic parallaxes is shown as error bars, and the best-fitting, preferred model to describe these residuals, ΔZ (Equation (6) and parameters from Model 0 in Table 1), is shown as the yellow band. The purple band shows the model for what the parallax residuals would look like if the asteroseismic radii were exactly on the Gaia parallactic scale (i.e., Y(a) = 0), assuming there is no color term in ΔZ (e = 0) or in asteroseismic parallax (d_a = 0), and in the absence of a local offset for the Kepler field's part of the sky (c = 0). As such, it represents the c₂ term of Equation (6), and which is the adjustment to the Z₅ model for which there is the strongest evidence (Section 4.4). The purple band thus demonstrates that the Z₅ model performs very well for G ≳ 11, apart from a constant, local offset for the Kepler field not shown by the purple band (c = −15 ± 2 μas from Model 0 in Table 1).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Below, I discuss aspects of the parallax residuals as a function of the different terms in Equation (6), considering evidence for not only refinements to the Z₅ model with a Kepler field-specific local offset, c; a color term, d; and a magnitude term, e; but also evidence for corrections to the asteroseismic parallaxes with the color term, d_a , and the asteroseismic radius-rescaling term, a. I also discuss the uncertainty budget in the parallax difference, including evidence for corrections to the fractional parallax uncertainties and quantifying systematic spatial variations in Gaia parallaxes.

4.1. Corrections to Gaia and Asteroseismic Parallaxes: The Color Term

It should first be noted that the color term in this analysis, d, describes a smaller range in ν_eff (1.4 μm⁻¹ ≲ ν_eff ≲ 1.5 μm⁻¹) than the linear color terms q_1k of the Z₅ model attempt to explain (1.24 μm⁻¹ ≲ ν_eff ≲ 1.72 μm⁻¹). This means that stronger local gradients in the parallax systematics may still remain after correction by Z₅, which the sample could be sensitive to. The sample will be most sensitive to color trends between 12 ≲ G ≲ 13, where most of the data are, meaning that any color-dependent refinement to Z₅ I find are not necessarily valid in other magnitude regimes and may not describe color trends within the brighter range of the sample, 9 ≲ G ≲ 12.

With these caveats in mind, a significant color term of d = −300 ± 25 μas μm is strongly favored to describe the parallax residuals (Model 11 versus Model 0). This trend is also clearly visible in Figure 4(e).

Zoom In Zoom Out Reset image size

Here, I take care to consider the possibility that d may not only describe an additive correction to EDR3 parallaxes but also may partially describe small problems in the asteroseismic parallax via the bolometric correction and/or Δν correction.

Regarding possible contributions to d from the asteroseismic parallaxes, consider, for example, that bolometric corrections used to calculate the asteroseismic parallax (Equation (1)) are found to vary as a function of temperature/color depending on the prescription by ≈2%–4% (Zinn et al. 2019b; Tayar et al. 2020). This is enough to explain some of the color dependence of the parallax difference: for a median parallax of the sample of ≈0.7 mas, the resulting effect on the asteroseismic parallaxes would be 7–14 μas, corresponding to differences in ≈ 70–100 μas μm in d, which is comparable to the best-fitting color term of d₂ = −300 ± 25 μas μm.

Additional color residuals may come about through the choice of the Δν correction, which depends on temperature and therefore color. I test the sensitivity of the color trend to Δν corrections by considering the nonlinear scaling relations of Kallinger et al. (2018) instead of the APOKASC-2 Δν corrections (Model 16). The resulting color term of −210 ± 25 μas μm reveals a ≈100 μas μm level variation due to the Δν correction choice. Similarly, the Yu et al. (2018) data using the same nonlinear scaling relations yield d = −180 ± 30 μas μm (Model 14). Taken together, these differences of ∼100 μas μm are suggestive of the level of color-dependent uncertainties in Δν corrections used to calibrate asteroseismic radii.

Due to these indications of contributions to d from asteroseismic parallax systematics instead of EDR3 parallax systematics, I considered a color-dependent correction explicitly for the asteroseismic parallaxes of the form d_a ϖ_seis(ν_eff − 1.48) (Equation (5)). This term is designed to capture color systematics in the asteroseismic parallax, which will tend to be fractional (dependent on ϖ_seis), given that bolometric corrections and Δν enter multiplicatively in Equation (1). (The same rationale motivates the fractional correction to asteroseismic parallaxes, a, in Equation (5).) Simultaneously allowing for an additive color term, d, and this fractional term, d_a , yields d_a = −0.14 ± 0.06 μm and a less substantial d = −220 ± 40 μas μm (Model 2). For the typical star in the sample with ϖ_seis ≈ 700 μas, this value of d_a corresponds to a color term due to systematics in the asteroseismic parallax of ≈ −100 μas μm, which is consistent with expectations as outlined above. Nevertheless, there is no strong evidence (i.e., ΔBIC < −6) for adding d_a either as a replacement for the additive color correction (Model 9 versus Model 0) or in addition to d (Model 2 versus Model 0). I therefore conservatively correct for a color term in what follows, without definitely attributing it fully to Gaia parallax systematics. To be sure, a color correction of some sort is required to explain the parallax residuals, though it is likely partially due to an additive adjustment needed of Z₅ as well as a properly parallax-dependent color term arising from bolometric correction and Δν systematics.

Non-Gaia contributions notwithstanding, there would appear to remain a ≈ −200 μas μm residual that would be attributable to color residuals in the Z₅ model. In this regard, one can helpfully refer to the consistency check performed by the Gaia team in the LMC. On average, there is not a large gradient in the LMC parallax across the whole range 1.1 μm⁻¹ ≲ ν_eff ≲ 1.9 μm⁻¹ after correction by the Gaia zero-point model (Figure 23 of Fabricius et al. 2020). However, it is clear that gradients exist on smaller scales, which, for 1.4 μm⁻¹ ≲ ν_eff ≲ 1.5 μm⁻¹, is approximately −140 μas μm, in the sense that bluer stars have too-small parallaxes. Although most of the LMC members used by the Gaia team have G > 13, this nonetheless indicates that linear residuals in parallax as a function of color exist after correction by Z₅. Moreover, this gradient is similar in magnitude but opposite in sign to the gradient I identify. Recalling that the sines of the ecliptic latitudes of the LMC and Kepler field are approximately the same magnitude (≈1) but of opposite signs, this is suggestive of a refinement to the q₀₁ term for G = 12.2 in the Z₅ model: q₀₁ sets the magnitude of a correction of the form $\sin \beta ({\nu }_{\mathrm{eff}}-1.48)$, and an increase from ≈40 μas μm adopted for the Z₅ model to at least 100 μas μm could seemingly be accommodated, given the uncertainties in q₀₁ (Figure 11 of Lindegren et al. 2020a). Keeping in mind both the potential errors in the Z₅ coefficients due to the bootstrap fitting approach in the bright regime (G < 10.8; Lindegren et al. 2020a) and the two times larger range in color explained by the q₀₁ term as opposed to d, I believe this is a plausible scenario.

I attempted to model the color-dependent residuals with the addition of a cubic term in the color dependence (d₃ in Equation (6)), which was motivated by the Gaia team's cubic term in Z₅, q₂₀. However, d₃ is strongly disfavored compared to having no cubic term, according to the ΔBIC in Table 1 (Model 5 versus Model 0). The data also prefer having a single color term instead of one to describe G < 10.8 and one to describe G ≥ 10.8, as is seen from the difference in ΔBIC between Model 4 and Model 0.

4.2. Corrections to Asteroseismic Parallaxes: The Radius-rescaling Term

4.3. Corrections to Gaia Parallaxes: The Constant Offset Term

4.4. Corrections to Gaia Parallaxes: The Magnitude Term

4.5. Comparison to Gaia DR2

4.6. Statistical Uncertainty in Gaia Parallaxes

4.7. Spatial Correlations in Gaia Parallaxes

In addition to pointing out the underestimation of statistical uncertainties in EDR3 parallaxes, the Gaia team quantifies systematic variations in parallax as a function of position on the sky induced by the scanning pattern of the satellite (Lindegren et al. 2020b). The present analysis extends the Gaia team's analysis to smaller scales, redder colors, and brighter magnitudes compared to their estimate with faint, blue quasars. The operating principle behind the estimate provided here is that asteroseismic parallaxes are astrophysically expected to have no intrinsic correlation with position in the 10 × 10 deg _² Kepler field, and, thus any spatial systematics in EDR3 parallaxes would manifest as spatially correlated differences between asteroseismic and EDR3 parallaxes.

In detail, the estimate of the spatial correlations in Gaia parallaxes follows Zinn et al. (2017). ⁴ I correct the Gaia parallaxes according to the Z₅ prescription and also according to Model 0, further subtracting off any residuals by forcing the Gaia parallaxes to have the same mean parallax as the asteroseismic parallaxes. I then compute the covariance of the parallax difference ϖ_seis − ϖ_Gaia, ${C}_{{ij}}=\langle {({\varpi }_{\mathrm{seis}}-{\varpi }_{\mathrm{Gaia}})}_{i}{({\varpi }_{\mathrm{seis}}-{\varpi }_{\mathrm{Gaia}})}_{j}\rangle$ for pairs of stars, i and j, in bins of angular separation, θ_ij , estimating the uncertainties in the covariance using from bootstrap sampling (Zinn et al. 2017). ⁵

The resulting covariance at given angular scales are shown as error bars in Figure 5. A positive value indicates a correlation between parallaxes at a given angular separation, and a negative value indicates an anticorrelation. If there were no spatial correlations in the Gaia parallaxes, the covariance would be zero.

Zoom In Zoom Out Reset image size

I then fit the covariance with the model of Equation (8) via MCMC, adopting uncertainties on each binned covariance point as shown in Figure 5. The resulting model is shown as a solid black curve in Figure 5. The gray curves are random draws of the best-fitting model, according to the covariance among the best-fitting parameters, as estimated from the MCMC chains.

I show for reference the Gaia team covariance model from Lindegren et al. (2020b) as the solid brown curve. I note that the Gaia model is fit to large scales θ_ij ≳ 05, and so I only show the Lindegren et al. (2020b) model in that regime. I have also subtracted the variance inferred on scales larger than the Kepler field of 76 μas² (Lindegren et al. 2020b) to make it more comparable to the results of this analysis, though I expect residual differences to remain (see below).

As can be seen in Figure 5 and from the best-fitting parameters in Table 2, the systematic covariance floor at the smallest scales is 4.0 ± 1 μas, which corresponds to 35% ± 10% of the statistical uncertainty of ϖ_Gaia in the sample. Already at a low level on small scales, the covariance becomes negligible on scales larger than a few degrees.

In the overlap region of 05 ≲ θ_ij ≲ 10° between the covariance estimated here and that estimated by Lindegren et al. (2020b) using quasars, the quasar covariance model is markedly larger than the covariances of the Kepler field (brown versus black curves in Figure 5). This may be indicative of smaller spatial systematics in the G < 13 regime compared to the G > 17 regime probed by quasars or may signal variation in the small-scale systematics across the sky. Note that there will be an offset between the covariance estimated here using the Kepler field and the covariance inferred from quasars because the latter has contributions from scales larger than 10°. As noted above, I attempt to correct for this difference by removing 8.7 μas in quadrature from the quasar covariance, which is the estimated rms scatter in quasar parallaxes on scales larger than 7° (Lindegren et al. 2020b); a large-scale rms scatter of ≈11 μas would be required to reconcile the quasar covariance and the Kepler field covariance at 05, which is also consistent with the Kepler field local offset I find of c = −15 ± 2 μas (Section 4.3). I also note that the quasar parallaxes were not corrected for the Z₅ model, which will leave spatial systematics due to the ecliptic latitude-dependent zero-point terms modeled with Z₅, and which could help explain the ≈2 μas larger rms scatter required to match the quasar and Kepler field covariances. I note that Fardal et al. (2021) find that the level of spatial correlation in Gaia DR2 parallaxes increases by a factor of 6 between G = 13 and G = 20. Under the assumption that this magnitude dependence is preserved in EDR3 spatial correlations, this would be another reason to expect a larger spatial covariance among quasars than I find here.

The spatial covariance found in this work would be expected to be in better agreement with the Gaia team's estimates of the covariance using the LMC than those using quasars because there are no contributions from large-scale variations. Additionally, the LMC members are more comparable in magnitude to the sample here. Indeed, the estimate from Lindegren et al. (2020b) of a systematic uncertainty of 6.9 μas for 01 ≲ θ_ij ≲ 45 accords better with the present estimate of 4 ± 1 μas at the smallest scales, keeping in mind the larger systematic expected of the LMC estimate due to not correcting the parallaxes using Z₅.

I stress that the spatial covariance estimation procedure used here will not be sensitive to any average parallax systematic that persists over the entire Kepler field of view, because I correct the parallaxes according to Equation (6). However, the local offset term that I find, c = −15 ± 2 μas, is a conservative estimate of what sort of average offset the asteroseismic and Gaia parallaxes have, averaged over the whole Kepler field (see Section 4.3). This is consistent with the expected rms of 8.1 μas on scales larger than 10° using quasars (Lindegren et al. 2020b). Studies making use of the small-scale covariance as described by Equation (8) and Table 2 would therefore be advised that the large-scale contributions to the covariance may well be larger than the small-scale ones that are available using the Kepler field.

I have compared Gaia EDR3 parallaxes with five-parameter solutions to parallaxes derived from knowledge of the temperature, bolometric flux, and asteroseismic radius for more than 2000 first-ascent red giant branch stars in the Kepler field, and conclude the following:

• 1.  
  The Z₅ zero-point model (Lindegren et al. 2020a) brings Gaia EDR3 parallaxes into agreement with asteroseismic parallaxes to an average of ≈2 μas before any adjustments to Z₅ or asteroseismic parallaxes for the regime 10.8 ≤ G < 13, 1.4 μm⁻¹ ≲ ν_eff ≲ 1.5 μm⁻¹ (purple band in Figure 2).
• 2.  
  I find strong evidence that the Z₅ zero-point model overcorrects parallaxes for 9 ≲ G ≲ 11, such that corrected parallaxes are too large by 15 μas ± 3 μas (purple band in Figure 4(f)), consistent with the findings for even brighter Cepheids distributed across the sky (6 ≲ G ≲ 10; Riess et al. 2021).
• 3.  
  I identify a significant color dependence to the asteroseismic−Gaia parallax residuals of −300 ± 25 μas μm for G ≳ 11 and 1.4 μm⁻¹ ≲ ν_eff ≲ 1.5 μm⁻¹, of which ∼100 μas μm may plausibly be attributed to asteroseismic parallax systematic uncertainties. If this color term were due to unmodeled behavior in Z₅, there is reason to believe it may have an ecliptic latitude dependence.
• 4.  
  After removing the aforementioned color term from the asteroseismic−Gaia parallax residuals and after adjusting Z₅ according to an additive −15 ± 3 μas offset for G < 10.8, the Gaia EDR3 and asteroseismic parallaxes agree to within a constant offset of c = −15 ± 2 μas. The latter may be interpreted as a conservative estimate of the large-scale variation in the parallax zero point, given its dependence on color trends, and is broadly consistent with the Gaia team's estimates of spatially correlated parallax variations on scales larger than the Kepler field.
• 5.  
  Gaia EDR3 parallax statistical uncertainties in the 9 ≲ G ≲ 13, 1.4 μm⁻¹ ≲ ν_eff ≲ 1.5 μm⁻¹ regime probed by the sample are underestimated by 22% ± 6%, consistent with estimates from the Gaia team (Fabricius et al. 2020) and in good agreement with independent analysis using wide binaries (El-Badry et al. 2021).
• 6.  
  The spatial correlations in Gaia EDR3 parallaxes on the scales probed by the Kepler field (≲10°)—described by Equation (8) with parameters in Table 2—agree broadly with those inferred by Lindegren et al. (2020b) to within an additive constant due to correlations on scales larger than the Kepler field. These spatial correlations incur a ≈4 μas systematic parallax uncertainty on angular scales less than ≈01, falling to below a ≈0.1 μas level on scales larger than ≈5° (Figure 5).

The zero point in EDR3 is not expected to be significantly different in the next installment of Data Release 3 scheduled in 2022, because the next complete astrometric solution will only be done for Gaia DR4. At that time, further reduction in the parallax zero point is expected in part due to the reversal of the Gaia satellite's precession (Lindegren et al. 2020b). Until such time, further characterization of the EDR3 Z₅ model would be fruitful, particularly to definitively constrain any color adjustments needed of it and to validate it as a function of ecliptic latitude.

I am grateful to Adam Riess for encouraging, productive conversations and instructive comments. I would also like to thank Ruth Angus and the Center for Computational Astrophysics for providing office space during the COVID-19 pandemic. Thanks go especially to the Gaia team for their revolutionary work. I am supported by an NSF Astronomy and Astrophysics Postdoctoral Fellowship under award AST-2001869.

This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions.

SDSS-IV acknowledges support and resources from the Center for High Performance Computing at the University of Utah. The SDSS website is www.sdss.org.

SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, Center for Astrophysics—Harvard & Smithsonian, the Chilean Participation Group, the French Participation Group, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, the Korean Participation Group, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astronomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatário Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronomical Observatory, United Kingdom Participation Group, Universidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, University of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation.

Software: asfgrid (Sharma & Stello 2016), NumPy (Walt et al. 2011), pandas (McKinney 2010), Matplotlib (Hunter 2007), IPython (Pérez & Granger 2007), SciPy (Virtanen et al. 2020), isoclassify (Huber et al. 2017; Berger et al. 2020), corner (Foreman-Mackey 2016).