Dynamical Mass of the Exoplanet Host Star HR 8799

HR 8799 is a young A5/F0 star hosting four directly imaged giant planets at wide separations ($\sim$16-78 au) which are undergoing orbital motion and have been continuously monitored with adaptive optics imaging since their discovery over a decade ago. We present a dynamical mass of HR 8799 using 130 epochs of relative astrometry of its planets, which include both published measurements and new medium-band 3.1 $\mu$m observations that we acquired with NIRC2 at Keck Observatory. For the purpose of measuring the host star mass, each orbiting planet is treated as a massless particle and is fit with a Keplerian orbit using Markov chain Monte Carlo. We then use a Bayesian framework to combine each independent total mass measurement into a cumulative dynamical mass using all four planets. The dynamical mass of HR 8799 is 1.47$^{+0.12}_{-0.17}$ \Msun assuming a uniform stellar mass prior, or 1.46$^{+0.11}_{-0.15}$ \Msun with a weakly informative prior based on spectroscopy. There is a strong covariance between the planets' eccentricities and the total system mass; when the constraint is limited to low eccentricity solutions of $e<0.1$, which is motivated by dynamical stability, our mass measurement improves to 1.43$^{+0.06}_{-0.07}$ \Msun. Our dynamical mass and other fundamental measured parameters of HR 8799 together with MESA Isochrones&Stellar Tracks grids yields a bulk metallicity most consistent with [Fe/H]$\sim$ -0.25-0.00 dex and an age of 10-23 Myr for the system. This implies hot start masses of 2.7-4.9 \Mjup for HR 8799 b and 4.1-7.0 \Mjup for HR 8799 c, d, and e, assuming they formed at the same time as the host star.


INTRODUCTION
HR 8799 is a A5/F0 V star (Gray & Kaye 1999) that is most well known for hosting four exoplanets that have been directly imaged within a cold extended debris disk and beyond a warmer interior dust belt (Marois et al. 2008Su et al. 2009;Booth et al. 2016;Wilner et al. 2018;Faramaz et al. 2021; see also review by Konopacky & Barman 2018). Among directly imaged systems, HR 8799 stands out as an unusual case of protoplanetary disk evolution due to its multiple planets orbiting at wide separations. In addition to HR 8799, the only other currently known systems with multiple imaged exoplanets are the young stars PDS 70, which hosts two accreting protoplanets within its transition strained the orbits of the companions (e.g., Soummer et al. 2011;Currie et al. 2012;Esposito et al. 2013;Goździewski & Migaszewski 2014;Pueyo et al. 2015;Maire et al. 2015;Zurlo et al. 2016;Konopacky et al. 2016;Wertz et al. 2017;Wang et al. 2018b;Gravity Collaboration et al. 2019), but the host star mass has been assumed to be fixed at ∼1.5 M in the majority of these analyses based on estimates from surface gravity measurements and high-mass evolutionary models (e.g., Gray & Kaye 1999;Moya et al. 2010a;Baines et al. 2012). For example, Maire et al. (2015), Zurlo et al. (2016) and Wertz et al. (2017) adopt a mass of 1.51 M based on the results of Baines et al. (2012), who inferred the mass by comparing their stellar radius measurement and effective temperature constraint to high-mass stellar evolution grids. Similarly, Konopacky et al. (2016) and Wang et al. (2018b) use a mass prior of 1.52±0.15 M for their orbit fitting analysis, which is also based on the results of Baines et al. (2012) coupled with an additional uncertainty to account for model systematics. Dynamical mass measurements, on the other hand, are "gold standards" because they rely only on Kepler's laws (e.g., Hillenbrand & White 2004, see also review by Serenelli et al. 2021). While host star dynamical masses have been measured in other systems with a directly imaged substellar companion (e.g., β Pic, Nielsen et al. 2014;Wang et al. 2016;Dupuy et al. 2019;Gravity Collaboration et al. 2020;Nielsen et al. 2020;Lagrange et al. 2020;Vandal et al. 2020;Brandt et al. 2021a), it is more often the case that the modest orbital coverage from long-period planets cannot admit a well constrained dynamical mass (e.g., Bowler et al. 2020).
The precise age of HR 8799 remains unsettled but is important because it is necessary to infer mass estimates of the orbiting companions using substellar evolutionary models (e.g., Marley et al. 2007;Fortney et al. 2008;Marleau & Cumming 2014). If HR 8799 is a member of the Columba Association (Doyon et al. 2010;Zuckerman et al. 2011;Malo et al. 2013) then its age is 42 +6 −4 Myr (Bell et al. 2015). However, Hinz et al. (2010) calculated the UVW space motion of HR 8799 using revised Hipparcos astrometry and called into question whether it is an unambiguous kinematic member of this group. They used their revised space motion together with modeling of the epicyclic orbit of HR 8799 and found that its distance from the center of Columba is implausibly large. Moya et al. (2010a) conducted an asteroseismic analysis of HR 8799 and estimated plausible ages as young as 26 Myr and as old as ∼1.6 Gyr; however, their interpretation depends on the unknown inclination of the host star 1 . Faramaz et al. (2021) reevaluated the moving group membership status of HR 8799 by calculating the UVW space motion with the latest Gaia EDR3 astrometry (Gaia Collaboration et al. 2021) and using the BANYAN Σ tool (Gagné et al. 2018). When using the stellar radial velocity constraints derived from Wang et al. (2018a) and Ruffio et al. (2019), they found that HR 8799 is more likely to be a field star rather than a member of Columba. Brandt et al. (2021b) recently presented a dynamical mass for the planet HR 8799 e and used it to derive a hot-start cooling age of 42 +24
In this study, we independently constrain the dynamical mass of HR 8799 by fitting orbits to relative astrometry of each of the four companions utilizing all previously published astrometry together with a new epoch we obtained of planets b, c, and d using Keck/NIRC2 in 2015. We then jointly derive the dynamical mass distribution of HR 8799 in a Bayesian framework by assuming the total mass interior to each orbit is simply the stellar mass (M total ≈ M ). HR 8799's dynamical mass, measured radius, bolometric luminosity, and effective temperature together with high-mass stellar evolutionary models offer a way to determine the system age independent of kinematic constraints, which in turn can be used to refine the inferred masses of the four planets using substellar evolution models.
In §2 we describe our new observations of the outer three planets with Keck/NIRC2 and summarize the additional astrometry that we use in our analysis. We detail the orbit fitting procedure and our Bayesian framework to jointly constain the host star mass in §3. Results are summarized in §4 and we explore the implications of our dynamical mass in §5.  Wright et al. (2011) conducted a spectroscopic asteroseismic analysis and found the stellar inclination to be i 40 • . However, asteroseismic studies of HR 8799 are generally difficult to carry out due in part to the complex nature of its pulsation frequencies (e.g., Sódor et al. 2014). during this time the sequence was interrupted for about 20 min during which time we obtained nearby 3 µm sky frames to serve as flats and 15 unsaturated images of the host star (each with an integration time of 0.1 s and 10 coadds) with the mask removed to photometrically calibrate the deeper coronagraphic dataset. Natural seeing in the visible (0.5 µm) ranged from 0. 5-0. 6 as recorded by the CFHT's Differential Image Motion Monitor. Bad pixels and cosmic rays were corrected using a nearest neighbor averaging algorithm, and each image was flat-fielded using a median-combined and normalized master sky frame. The bright sky background at 3 µm is removed by subtracting the master sky frame for each image. Images are registered following the description in Bowler et al. (2015); a 2D elliptical Gaussian is fit to the host star, which is visible behind the coronagraph, and an ADI data cube centered on HR 8799 is assembled. Each image is corrected for optical distortions with the distortion solution from Service et al. (2016). PSF subtraction is carried out using LOCI (Lafrenière et al. 2007) with the following parameter values that control the angular tolerance as well as optimization and subtraction zone sizes and shapes: W =7, N A =300, g=1.0, N δ =0.5, and dr=2.0. Circular masks are placed at the locations of HR 8799 b, c, d, and e to avoid biasing the optimization regions during PSF subtraction.
The final processed image and signal-to-noise ratio (S/N) map are shown in Figure 1. Planets b, c, and d are clearly visible, but planet e is not evident at a significant level in the processed frame. S/Ns for each planet are computed using non-overlapping circular apertures of 6 pix at 36 azimuthal angles centered on HR 8799; we measure S/Ns of 9.8, 7.3, and 7.2 for planets b, c, and d, respectively. Astrometry and relative photometry for each source is calculated using the negative PSF injection approach described in Bowler et al. (2018). The median-combined unsaturated image of HR 8799 is used as the PSF model for each planet. The PSF model is inserted in the raw images with a starting guess position and (negative) amplitude. PSF subtraction is carried out in the same fashion as for our final processed image, and the RMS of the residuals at the location of the planet is used to assess the quality of subtraction. The amoeba downhill simplex algorithm (Nelder & Mead 1965) is used to find the best-fitting separation, P.A., and flux ratio of each planet ( Figure 2). The uncertainties in these quantities are computed using the standard deviation of the last ten iterations of the amoeba algorithm as it settles in its final value that minimizes the

Previous Astrometry
We make use of all astrometry of the HR 8799 planets published prior to December 2020 2 in our orbit fitting analysis in addition to our new epoch from 2015. Observations of this system have been taken with a wide range of telescopes and instruments, including HST, Subaru, Keck, Gemini, MMT, VLT, Palomar, and LBT. Astrometry published prior to 2016 are presented in Bowler (2016, Table 4 in their Appendix) 3 . Additionally, we incorporate published astrometry obtained with Gem-ini/GPI (Wang et al. 2018b) and with VLTI/GRAVITY (Gravity Collaboration et al. 2019). For astrometry in ∆RA and ∆Dec, we convert to separation and position angle and propagate the uncertainties in a Monte Carlo fashion. For multiple reductions of the same dataset, or where more than one dataset exists for the same epoch, we use the measurements that appear to have the most realistic astrometric errors. Specifically, for epoch 1998.  .58, 2007.81, 2008.72, 2009.58, 2009.83, 2010.53, and 2010Wertz et al. (2017)

Accounting for Potential Systematic Errors
Using all available astrometry in the orbit fitting could introduce systematic errors from a variety of sources, including the use of differently calibrated instruments, PSF-subtraction algorithms, and strategies to measure relative astrometry. One approach to avoid this problem is to adopt a uniform dataset from a small number of well-calibrated instruments. For example, Konopacky et al. (2016) used data obtained solely from the NIRC2 camera at Keck Observatory to create self-consistent measurements for their orbit fitting. Several other recent studies adopted similar strategies (e.g., Wang et al. 2018b;Gravity Collaboration et al. 2019). However, while this approach may mitigate systematic errors, it also limits the total time baseline and number of astrometric epochs for the orbit fits. Here we use all available astrometry of the HR 8799 planets as of December 2020 with the goal of maximizing the information content and time baseline of the data. Given that potential systematic uncertainties may affect our dynamical mass measurement, we consider an additional uncertainty term that we refer to as astrometric "jitter" (σ jit ) to be added in quadrature to each set of raw astrometric errors following the approach in Bowler et al. (2020). For each planet, we assume a lin-ear model for both separation (ρ) and position angle (θ) as a function of time because the astrometric time baseline relative to the orbital period is small ( 0.16 for all planets). We adopt a reduced χ 2 value (χ 2 ν ≡ χ 2 /ν, where ν is the number of degrees of freedom) as a metric to assess reliability of the quoted uncertainties. The χ 2 ν values for all four planets' separations as a function of time are < 1.0 (0.37, 0.59, 0.61, and 0.69 for HR 8799 b, c, d, and e, respectively), indicating that the . Orbit fits for HR 8799 b, c, d, and e using orbitize!. Each panel shows 150 orbits drawn from the converged MCMC chains as described in §3.2 together with the measured astrometry. For each planet, the left panel displays the orbits as they appear on-sky and the right panels compare the same orbits to the astrometry. raw uncertainties appear to be reasonable. These are left unchanged in this study. However, all four planets' position angles as a function of time have χ 2 ν > 1.0 (5.4, 2.4, 2.3, and 2.3 for HR 8799 b, c, d, and e, respectively), suggesting that the quoted uncertainties are underestimated or that systematic errors are present. The PA jitter term (σ jit,PA ) for each planet is found by iteratively increasing a starting value of 0 in steps of 10 −5• until χ 2 ν reaches unity. This results in σ jit,PA values of 0.37 • , 0.26 • , 0.62 • , and 0.69 • for HR 8799 b, c, d and e, respectively. These are added in quadrature to the original uncertainties: σ 2 total,PA = σ 2 jit,PA + σ 2 PA . The original and adjusted astrometry are displayed in Figure 3.

Orbit Fitting
We fit Keplerian orbits to the astrometric data of each planet using the orbitize! . In the Keplerian model used in orbitize! 4 , the following parameters are varied directly when fitting with relative astrometry: a (semi-major axis), e (eccentricity), i (inclination), ω (argument of periapsis), Ω (position angle of the ascending node), τ (time of last periapsis), π (parallax), and M tot (total mass of the system). We treat each planet independently, corresponding to four separate orbital constraints.
A Gaussian prior of 24.46 mas with a standard deviation of 0.05 mas from Gaia EDR3 (Gaia Collaboration et al. 2016 is used for π. A log-uniform prior distribution is adopted for a (i.e., P (a) ∝ a −1 ) with a range of 1 to 150 au. An isotropic inclination distribution is adopted for i such that P (i) ∝ sin i (which spans 0 to π rad). Linearly uniform priors are chosen for the remaining free parameters. Their ranges are: 0 to 1 for e 5 ; 0 to 2π rad for both ω and Ω; 0 to 1 for τ ; and 0.5 to 4.0 M for M tot . To avoid artificially reducing the precision of the final dynamical mass measurement, we make no further assumptions about the nature of the planets' orbits such as coplanarity, resonances, stability, or orbit crossing events.
Each ptemcee ensemble is initialized with 16 temperatures and 1500 walkers per temperature. Each walker samples from the parameter space 1.2×10 5 times and a burn-in length of the first 6.0×10 4 steps is discarded. To assist their convergence, the walkers are initialized at positions near the most probable values found from conducting preliminary orbit fits. Saving only the samples from the lowest-temperature set of chains and applying a thinning factor of 20 to curtail the effect of correlation results in 4.5×10 6 posterior samples for each orbit fit.

Bayesian Framework to Constrain the Mass of HR 8799
We employ a Bayesian inference framework to derive the joint probability distribution for the host star mass based on our individual orbit fits for each of the four planets. Bayes' theorem naturally updates prior information with new data to improve posterior probabilities: (1) Here P (θ) represents prior knowledge of the parameters θ, P (D | θ) represents the likelihood function of the data D given the model, and P (θ | D) represents the posterior probability of the model parameters given the data. For this study of HR 8799, our orbit fits are carried out for a uniform prior distribution in total (stellar) mass, so each of the four total mass posterior distributions from the independent MCMC fits also represent the likelihood function for the dynamical mass of HR 8799 (to within a constant of proportionality) given the set of astrometry for that planet. This arises by treating the planets as massless "test" particles such that the total mass of the system is assumed to be confined to the host star (M p M , or M ≈ M tot ). If the HR 8799 planets have masses 10 M Jup , as expected from evolutionary models (e.g., §5.2), then this approximation is accurate at the 2% level. We use the following notation to apply Bayes' Theorem to each of the four independent total mass distributions for planets b, c, d, and e that resulted from orbit fitting: (2) P i is therefore the likelihood function of the astrometry for planet i conditioned on the stellar (total) mass M . The prior probability for M , P (M ), is informed by prior knowledge of the stellar mass (e.g., based on spectral type or evolutionary models).  planets b and c this is The jointly constrained mass from all four planets is then: This setup allows us to only run the orbit fits once for each planet with a uniform mass prior, but then easily adjust the prior P (M ) as desired for the joint mass constraint.

RESULTS
For each planet, the resulting M tot posterior distributions are linearly interpolated onto a fine grid spanning 0.5-4.0 M ( Figure 5). These are then multiplied together following Equation 3 to arrive at the jointly constrained dynamical mass distribution of HR 8799. Given that the individual total mass constraints are non-Gaussian and right skewed, we summarize the results using the posterior distribution median and highestdensity credible interval (HCI) metrics.
For the initializing stellar mass prior P (M ), we adopt a linearly uniform distribution to let the astrometry drive the dynamical mass. The final result for this baseline case is shown as Figure 6. The median and 68% HCI of the final joint mass posterior is 1.47 +0.12 −0.17 M 6 . The 95% and 99.7% HCI span 1.23-1.86 M and 1.16-2.32 M , respectively.
Our choice of a linearly uniform prior may not be especially realistic given that we know HR 8799 is a late-A/early-F star. Our joint stellar mass measurement may be further constrained by using a more informative prior that better captures knowledge of the host star before our experiment takes place while at the same time not being overly restrictive to ensure the astrometry still dominate the posterior. Gray & Kaye (1999) conducted a spectroscopic analysis of HR 8799 and estimated its mass to be 1.47 ± 0.30 M based on their surface gravity measurement combined with an estimate of the radius. Following the same analysis detailed above but now using this Gaussian prior results in a dynamical mass of 1.46 +0.11 −0.15 M for HR 8799. Here the 95% and 99.7% HCI span 1.25-1.75 M , and 1.18-1.97 M , respectively. The effect of the informed prior is to further reduce the uncertainties on the final joint stellar mass by ≈10%, which is expected given that this prior is in good agreement with the original dynamical mass (which used a uniform prior).
While not the focus of this study, we summarize results for the orbital parameters from our fits with cornerplots and a table of best-fit values in Appendix A. At the 95% HCI levels, the eccentricities of the planets are constrained to e 0.5, with the exception of HR 8799 d at e 0.6. Similarly, at the 68% HCI levels, the eccentricities of planets b, c, and e are constrained to e 0.3 while HR 8799 d has e 0.4. This is consistent with previous studies that also found the allowed eccentricities of HR 8799 d to be marginally greater than those of the other three planets (e.g., Pueyo et al. 2015;Wertz et al. 2017). The most probable inclinations for the four planets span i ≈25-35 • , which are similar to each other as well as with that of the extended debris disk (i ≈ 31-33 • , Wilner et al. 2018;Faramaz et al. 2021). However, the longitude of the ascending node (Ω) 7 must be accounted for in addition to i to assess coplanarity. For any two orbital planes, these parameters are related to the true mutual inclination ψ by the geometric relation ψ = arccos(cos i 1 cos i 2 + sin i 1 sin i 2 cos(Ω 1 − Ω 2 )). We place a 95% confidence upper limit of ψ < 55 − 75 • for all permutations of the HR 8799 planet orbital planes. Orbital parameters in this study are generally consistent with results from prior published orbit fits, although our posteriors tend to be larger. This is an expected consequence of our methodology: we solved for the total system mass using agnostic priors, made no constraining assumptions on allowed orbits based on stability, and inflated the astrometric uncertainties to account for potential systematic errors.

Low-Eccentricity Scenario
The posteriors in Figure 8 show that there is strong covariance between planet eccentricities and total system mass; higher planet eccentricities imply a substantially larger total mass. To assess how the cumulative dynamical mass might change if near-circular orbits are considered, we reconduct the analysis in an identical manner but we now limit the possible orbit solutions for each planet to e ≤ 0.1. Although low eccentricities on their own are not necessarily indicative of orbital stability (e.g., Fabrycky & Murray-Clay 2010), this exercise gives some insight into how the stellar mass constraint depends on eccentricity and what we might expect if long-term orbital stability is taken into account 8 . The joint dynamical mass in this scenario of near-circular orbits is 1.43 +0.06 −0.07 M , with 95% and 99.7% HCI intervals spanning 1.31-1.57 M and 1.28-1.65 M , respectively. The smaller uncertainties are expected given the strong correlation between eccentricity and total mass in our baseline orbit fit. Using the informed Gaussian mass prior of 1.47 ± 0.30 M does not significantly influence the joint stellar mass constraint (M = 1.43 +0.06 −0.07 M in this scenario). The orbit fitting results for this low eccentricity scenario are presented in Appendix B. We conclude that a restricted eccentricity significantly impacts the dynamical mass by reducing the overall uncertainty by a factor of ≈2. For this study, however, we adopt the unrestricted solution in which eccentricity is allowed to vary from 0 to 1 to avoid assumptions about mutual planet interactions, which are not yet robustly validated by the astrometry.

DISCUSSION
We fit orbits to all the available astrometry of the HR 8799 planets with minimal dynamical assumptions, then combined the independent total mass constraints into a joint dynamical mass measurement. Table 1 from Baines et al. (2012) compiles many of the existing mass estimates of HR 8799, which have been determined with a variety of methods and generally span the range of ∼1.2-1.8 M . Our dynamical mass measurement is consistent with all these prior estimates at its 95% HCI level, including the estimates of 1.513 +0.023  (Yi et al. 2001) evolutionary models. Our dynamical mass measurement thus serves as a reassuring validation of the ≈1.5 M stellar mass commonly adopted in previous studies of the HR 8799 system. Wang et al. (2018b) used a self-consistent astrometric dataset from NIRC2 and GPI (55 epochs in total) to simultaneously measure the companion orbital parameters and total mass of HR 8799 in an MCMC framework. In their "unconstrained" fits, they found a total mass of 1.48 +0.05 −0.04 M which, under the approximation M ≈ M tot , can be considered another dynamical mass measurement of HR 8799 to compare with ours. The larger uncertainties in this study are likely influenced by our adoption of a linearly uniform total mass prior compared to their adoption of a well-informed Gaussian prior of 1.52±0.15 M . Furthermore, by using a dataset only from NIRC2 and GPI without the inflated astrometric uncertainty that we accounted for, their astrometric measurement uncertainties were overall smaller compared to the astrometry used in this study. However, despite some differences in procedure and data used, both results agree well. The dynamical mass presented in this work should be considered a conservative, model-independent measurement that prioritized astrometric baseline over per-epoch precision and instrument uniformity.

Age of HR 8799
Several stellar parameters of HR 8799 (radius, luminosity, and effective temperature) were directly constrained by Baines et al. (2012) using a Center for High Angular Resolution Astronomy (CHARA) Array stellar angular diameter measurement together with photometry and the parallax from Hipparcos. They found a radius of R = 1.44±0.06 R , a bolometric luminosity of L = 5.05±0.29 L , and an effective temperature of T eff = 7193±87 K; these are based on a bolometric flux of F bol = (1.043±0.012)×10 −7 erg s −1 cm −2 , an angular diameter (accounting for limb darkening) of θ LD = 0.342±0.008 mas, and a distance of D = 39. They also updated the radius to R = 1.502±0.035 R when using their Gaia EDR3 distance together with the stellar angular diameter from Baines et al. (2012).
The measured mass, luminosity, radius, and effective temperature of HR 8799 (Table 2) can be used to constrain its bulk metallicity and age using stellar evolutionary models. We downloaded a custom grid of MESA Isochrones & Stellar Tracks (MIST) evolutionary models (Paxton et al. 2011(Paxton et al. , 2013(Paxton et al. , 2015(Paxton et al. , 2018 References-(1) This work (assumes a uniform stellar mass prior); (2)  like HR 8799 follow a meandering trajectory through the HR diagram as they quickly reach the main sequence and then evolve off it. In doing so their luminosity and effective temperature will repeatedly rise and fall. As a result, a given luminosity and temperature can correspond to two or more degenerate ages during the preand post-main sequence evolutionary phases. The additional mass and radius measurements help by limiting the ages consistent with all of these constrained parameters.
We begin our assessment of the bulk metallicity and age by randomly drawing a stellar mass rounded to the nearest 0.01 M with a probability proportional to the dynamical mass posterior distribution of HR 8799 found in §4. We also draw a luminosity from a normal distribution with a mean of 5.441 L and a standard deviation of 0.066 L . The iso-mass track will either be inconsistent with the given luminosity throughout the entire evolution of that star, or it will intersect the given luminosity at least once. If the former occurs, that trial is ignored and another mass is drawn; for the latter, ages are interpolated across the grid sampling at the luminosity crossing point. If the corresponding model radius and temperature at that point are within 3σ of the measured values of 1.502 ± 0.035 R and 7193 ± 87 K, respectively, then the ages are retained. This process of drawing a mass and luminosity and saving ages consistent with the models is repeated 10 4 times to build up a distribution of ages consistent with the measured values. Results for the linearly uniform mass prior are shown in Figure 7. It is shown in black for the pre-main sequence phase of evolution and then becomes gray for the post-main sequence phase. Similarly, the dashed black and gray lines correspond to the 68% HCI. Each column of panels corresponds to a different metallicity from [Fe/H] = -0.5 to +0.25 dex. The orange ellipse represents the measurement and 1-, 2-, and 3-σ uncertainties of L and T eff for HR 8799. Middle Row : The inferred age distributions of HR 8799 following the method described in §5.1. Bottom Row : Corresponding hot-start masses for the planets derived directly from the age distribution and planet luminosities.
Our use of four different [Fe/H] values enables a qualitative constraint of the stellar bulk metallicity: the stellar parameters are most consistent with [Fe/H]∼ -0.25-0.00 dex but only marginally consistent with [Fe/H]= -0.5 and +0.25 dex. The λ Boötis nature of HR 8799 has made previous attempts to characterize the bulk metallcity difficult because the surface abundances may not be an accurate tracer of the internal composition (e.g., Moya et al. 2010b;Murphy et al. 2021). In spite of the coarse sampling in the metallicities as well as potential systematic uncertainties in the MESA grids themselves, these qualitative constraints support a picture in which HR 8799 may have a solar-like or slightly sub-solar bulk composition.
There are two general solutions for the age of HR 8799 that correspond to an ambiguity related to whether it is moving toward or evolving from the zero age main sequence. For the [Fe/H]= 0.00 dex grids, these correspond to ages of ∼10 7−7.4 yr and ∼10 8−9.1 yr. Fabrycky & Murray-Clay (2010) found that the outer three HR 8799 planets must be 20 M Jup to remain stable on timescales of tens of Myr. Motivated by this dynamical stability requirement, we ignore the older age scenario because this results in implausibly high planet masses up to the hydrogen-burning limit. For each metallicity grid, we discard ages corresponding to a planet mass cutoff of >20 M Jup using hot-start evolutionary models (see §5.2).
Excluding the [Fe/H]= +0.25 dex grids, which are least consistent with the measured properties of HR 8799, an age of 10-23 Myr best agrees with our dynamical mass. This is younger than the age of 38-48 Myr (Bell et al. 2015) that is associated with HR 8799 if it is a member of the Columba moving group. This is also somewhat younger than the hot-start cooling age of 42 +24 −16 Myr for HR 8799 e recently derived by Brandt et al. (2021b), although at twice their quoted uncertainties, both ages are consistent. Baines et al. (2012) use Yonsei-Yale grids fixed at solar metallicity to derive an age consistent with either ∼33 Myr or ∼90 Myr. Regardless of kinematic association with a young moving group, these model-based constraints all support the picture that the HR 8799 system is <100 Myr and probably younger than 50 Myr.
To assess the impact of mass priors and restricted eccentricities on the results, we repeat this exercise for the other three scenarios considered in §4: the dynamical mass with the informed Gaussian prior, and the two low-eccentricity cases. Summary statistics for the resulting trimmed distributions are tabulated in Table 3. Figures showing the ages and corresponding hot-start masses using our other dynamical masses under different assumptions for the mass prior and range of eccentricities can be found in Appendix C. The values of the ages and hot-start masses do not significantly change depending on which dynamical mass is used, and the constraints of [Fe/H] also remain consistent.

Hot-Start Mass Constraints of the Planets
In light of these revised age estimates, we use hotstart evolutionary models from Burrows et al. (1997) to infer the masses of the imaged companions using their bolometric luminosities (Figure 7). Here we assume the age of the star to be the age of the planets. For planet b, we adopt a bolometric luminosity of log(L bol /L ) = -5.1 ± 0.1 dex, and for planets c and d we adopt log(L bol /L ) = -4.7 ± 0.1 dex (Marois et al. 2008;Rajan et al. 2015). The spectrum and near-infrared contrast of planet e is similar to that of planets c and d (e.g., Marois et al. 2010;Greenbaum et al. 2018), thus we adopt the same value of log(L bol /L ) = -4.7 ± 0.1 dex for this inner planet.
Excluding the [Fe/H]= +0.25 dex grids and the discarded solutions corresponding to old (∼Gyr) ages, the resulting mass estimates are 2.7-4.9 M Jup for HR 8799 b and 4.1-7.0 M Jup for HR 8799 c, d, and e (Table 3). These masses are within the giant planet regime assuming a boundary of ≈13 M Jup and are lower than, but still broadly consistent with, many of the previous estimates in the literature. For example, Marois et al.  M Jup for HR 8799 e by using a Gaia EDR3 host star proper motion anomaly together with results from Wang et al. (2018b) for orbital solutions and planet mass ratios. This is somewhat higher than the values we found and could be caused by an underestimated age, compositional difference between the star and planet, higher planet luminosity, or systematic errors in the models themselves. Altogether, the mass estimates presented in this work imply a minimum-mass extrasolar nebula of M ≥15-26 M Jup for HR 8799.

CONCLUSION
We have measured the model-independent dynamical mass of HR 8799 in a Bayesian framework using all available astrometry of its four planets and treating them as massless independent particles. The joint dynamical mass of HR 8799 is 1.47 +0.12 −0.17 M assuming a uniform prior, validating regularly used values of ≈1.5 M in previous studies. The modest but realistic uncertainties improve to 1.46 +0.11 −0.15 M when an informative prior based on the stellar spectroscopy is used. When only near-circular orbits of e <0.1 are allowed, the joint dynamical mass measurement becomes 1.43 +0.06 −0.07 M . We used our stellar mass constraints together with other previously measured parameters of the host star to investigate the age and bulk metallicity of HR 8799. The age constraint from MESA models is 10-23 Myr after excluding ages corresponding to hot-start masses of ≥20 M Jup . This supports an intermediate age for this system independent of kinematic membership status to young moving groups. We also find that the favored bulk metallicity is [Fe/H]∼ -0.25-0.00 dex.
Finally, using these inferred ages, we derived hot start masses of 2.7-4.9 M Jup for HR 8799 b and 4.1-7.0 M Jup for HR 8799 c, d, and e. These are somewhat lower than typical estimates, as expected from the younger inferred age. Continued astrometric monitoring of this system will enable future studies to progressively reduce the uncertainty in the orbital elements of the planets and the mass of HR 8799. However, eventually the assumption that M ≈ M tot will break down as the total mass precision becomes comparable to the planet masses. At that point, a more complete dynamical modeling that includes mutual planet interactions will be needed to properly account for the total system mass interior to each orbit (e.g., Lacour et al. 2021). for constructive suggestions as well as Henry Ngo and Sarah Blunt for their consulting with the orbitize! package.
A.G.S acknowledges support from the Texas Astronomy Undergraduate Research experience for Underrepresented Students (TAURUS) Scholars program. TAURUS graciously acknowledges funding from the National Science Foundation (NSF), National Aeronautics and Space Administration (NASA), the Cox Board of Visitors for the University of Texas' Department of Astronomy and McDonald Observatory, and many generous donors who gave to the TAURUS HornRaiser campaign. TAURUS is also grateful for the countless time commitments of students, staff, and faculty at UT to build a healthier, happier, and more representative astronomy community. This material is This work received computational support from UTSA's HPC cluster SHAMU, operated by the Office of Information Technology. A.G.S. thanks Eric Schlegel for the sponsorship that enabled the use of this resource.
Some of the data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain. 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 institu-tions, in particular the institutions participating in the Gaia Multilateral Agreement. This research has made use of NASA's Astrophysics Data System Bibliographic Services.
Facility: Keck:II (NIRC2) Software: orbitize! , ptemcee (Foreman-Mackey et al. 2013;Vousden et al. 2016), matplotlib (Hunter 2007), numpy (van der Walt et al. 2011;Harris et al. 2020 The results of our MCMC orbit fits for the loweccentricity scenario are summarized in Figure 9 and Table 5. We also display a sample of the sky-projected orbits drawn from the MCMC chains in Figure 10.

C. MESA AGES AND HOT-START PLANET MASSES
This Appendix includes a summary of results from stellar and substellar evolutionary models following our procedure described in §5.1 but using the other three dynamical masses from §4. Figure 11 shows results using the dynamical mass with the informed Gaussian prior. Results with the low-eccentricity assumption using the dynamical mass with a uniform prior and with the informed Gaussian prior are displayed in Figures 12 and  13, respectively.  a Maximum a posteriori probability.  a Maximum a posteriori probability.