Wolf 503 b: Characterization of a Sub-Neptune Orbiting a Metal-Poor K Dwarf

Using radial velocity measurements from four instruments, we report the mass and density of a $2.043\pm0.069 ~\rm{R}_{\oplus}$ sub-Neptune orbiting the quiet K-dwarf Wolf 503 (HIP 67285). In addition, we present improved orbital and transit parameters by analyzing previously unused short-cadence $K2$ campaign 17 photometry and conduct a joint radial velocity-transit fit to constrain the eccentricity at $0.41\pm0.05$. The addition of a transit observation by $Spitzer$ also allows us to refine the orbital ephemeris in anticipation of further follow-up. Our mass determination, $6.26^{+0.69}_{-0.70}~\rm{M}_{\odot}$, in combination with the updated radius measurements, gives Wolf 503 b a bulk density of $\rho = 2.92\pm ^{+0.50}_{-0.44}$ $\rm{g}~\rm{cm}^{-3}$. Using interior composition models, we find this density is consistent with an Earth-like core with either a substantial $\rm{H}_2\rm{O}$ mass fraction ($45^{+19.12}_{-16.15}\%$) or a modest H/He envelope ($0.5\pm0.28\%$). The low H/He mass fraction, along with the old age of Wolf 503 ($11\pm2$ Gyrs), makes this sub-Neptune an opportune subject for testing theories of XUV-driven mass loss while the brightness of its host ($J=8.3$ mag) makes it an attractive target for transmission spectroscopy.


INTRODUCTION
One of the most notable discoveries in the exoplanet field is the ubiquity of not one, but two new classes of planet frequently found orbiting late-type stars with periods less than 100 days; the super Earths and sub-Neptunes. This planet sub-population, first discovered over a decade ago through Doppler surveys of the southern sky (Mayor & Udry 2008;Lovis et al. 2009), has expanded dramatically with the Kepler/K2 and TESS missions (Fulton et al. 2017). With over 2,000 confirmed planets, these missions have presented us with a diversity of worlds that we previously had not anticipated. As with any discovery, these planets have forced us to rethink and reformulate not only our theories of planet formation, but also the evolution of exoplanets and their atmospheres as well as how they are affected by their host stars.
By combining both radius and mass measurements with models of planetary interiors it appears that these short period sub-Neptunes (1.6-3.2 R ⊕ ) potentially range from volatile-rich worlds with hydrogen/helium envelopes constituting nearly a third of their mass to rocky cores stripped of their atmosphere by their host star (Rogers & Seager 2010;Owen & Wu 2013;Lopez & Fortney 2014). Interior models even suggest some sub-Neptunes could host hydrospheres of super-critical water blanketed by steam-dominated envelopes (Zeng & Sasselov 2014;Thomas & Madhusudhan 2016;Mousis et al. 2020). With core-accretion as the prevailing theory of planet formation, the frequency of these planets was pre- viously thought to be, at best, rare. The low density, high temperature environment of the proto-planetary disk within the snowline makes building planets larger than Earth through core-accretion inefficient, let alone planets with substantial gaseous envelopes. However, in-situ formation of sub-Neptunes can still be possible as large dust grains drift from the outer disk inward and accumulate at the inner-most regions of the disk in a pebble-accretion scenario. In contrast, formation beyond the snowline and subsequent migration inward can also lead to sub-Neptunes (Bodenheimer & Lissauer 2014;Venturini & Helled 2017). The added complexity required in our formation theories gives rise to other intriguing questions: How has planetary migration affected the exoplanet populations we observe today? Are the properties of host stars reflected in the planets that orbit them?
In the post-Kepler era, efforts have now shifted from discovery to characterization and the answers to those questions seem to be on the horizon. Radial velocity (RV) surveys to measure the masses of previously discovered planets have been essential in placing these planets in context. Precise masses are especially important in modeling the potential atmospheric conditions on these planets. The degeneracy between the mean molecular weight of the atmosphere (a potential indicator of metal content) and the surface gravity can have a serious impact on a planet's potential for follow-up observations since these parameters can have similar effects on transmission spectra. The best way to break the degeneracy is to have mass uncertainties ≤ 20% (Batalha et al. 2019). With the James Webb Space Telescope (JWST) set to be the premier facility for studying exoplanet atmospheres, it is essential to not only measure precise masses but also constrain the transit-timing uncertainties beyond the K2 values; Spitzer has been instrumental in assuring valuable telescope time will not be wasted on missed transit events.
In this paper we characterize the sub-Neptune Wolf 503 b (Peterson et al. 2018). This planet orbits a bright (J = 8.3 mag) K-dwarf star making it an intriguing candidate for future atmospheric follow-up. This work is outlined as follows: We begin by describing the properties of the host star as well as re-deriving key stellar parameters with new Gaia Early Data Release 3 (EDR3) parallax values in Section 2. In Section 3 we present the analysis of K2 photometry using a photo-eccentric model which implies and eccentric orbit for Wolf 503 b. Section 4 describes the observations made with the Spitzer space telescope and uses the results to further constrain the orbital ephemeris. In Section 5 we use radial velocity data to further confirm the planet's eccentricity and conduct a joint RV-transit fit to better constrain the orbital parameters. Section 6 discusses possible interior composition Wolf 503 b along with the potential for atmosphere characterization with JWST.

TARGET SYSTEM PARAMETERS
Wolf 503 (EPIC 212779563, HIP 67285) is a bright (J = 8.3 mag) K3.5V main sequence dwarf. At 44.630 ±0.033 pc, this nearby star is currently known to host one planet, Wolf 503 b, that was discovered in K2 campaign 17 photometry in 2018 (Peterson et al. 2018). Wolf 503 b was found to be a 2.043 R ⊕ planet that completes one orbit roughly every 6 days. At 0.06 AU from its star, Wolf 503 b has an equilibrium temperature ∼ 800 K; an intermediate temperature compared to other sub-Neptunes discovered. In order to increase the accuracy of derived parameters such as mass and radius, we re-derived key stellar parameters for Wolf 503 using the Gaia mission's new parallax measurements the uncertainty on which has been reduced by a factor of three (EDR3 Gaia Collaboration et al. 2020). With values for spectroscopic parameters T eff , [Fe/H], and log(g) from Peterson et al. (2018) and the photometric magnitude in the K-band, we use Isoclassify (Huber et al. 2017;Berger et al. 2020) to obtain the distance, R * , and L * of Wolf 503. Isoclassify determines stellar parameters using a sample of 2200 Kepler stars in combination with Gaia data with uncertainties on those parameters based on MIST data. We use the direct method described in Huber et al. (2017) to determine these parameters which are listed in Table 1. The values found for distance, R * , and L * agree with values previously found by Peterson et al. (2018) but the uncertainties see slight reductions (≤ 1% in the case of radius and luminosity).

A Metal-Poor Host
The age and metallicity of its host star sets Wolf 503 b apart from the majority of sub-Neptunes. Wolf 503's age is estimated to be between 9-13 Gyrs and the star has a metallicity of [Fe/H] = −0.47 ± 0.08 (Peterson et al. 2018) making it one of the more metal poorer stars to host a sub-Neptune ( Figure 1). It has been well established that Jupiter-class planets are frequently found orbiting stars of increasing metallicity with a correlation between close-in giant planet occurrence rate and host star metal enrichment (Gonzalez 1997;Santos et al. 2004;Thorngren et al. 2016). This correlation can be understood in the context of core accretion; massive planets need more solid material in order to trigger a runaway accretion of gas. Although this trend weakens with decreasing planetary size , warm sub-Neptune occurrence is still correlated with host star metallicity . However, the formation of sub-Neptune planets require specific disk conditions that balance the build up of a massive core while also preventing the runaway accretion that results in a gas giant. Venturini & Helled (2017) found that formation scenarios with low solid accretion rates ( 10 −6 M ⊕ yr −1 ) resulted in the highest sub-Neptune occurrence rate. This accretion rate is compatible with disks of low metallicites but is also possible in low mass disks as well.

K2 SHORT CADENCE PHOTOMETRY
While the detection of eccentric orbits is usually done with radial velocity observations, through the photoeccentric effect (Dawson & Johnson 2012), one can obtain broad constraints on a planet's eccentricity from its lightcurve if an independent measurement of the stellar density can be made. In this section we extract previously unused K2 short cadence photometry and, with a stellar density obtained in Section 2, use a photoeccentric transit model to determine if Wolf 503 b is on a circular or eccentric orbit.

Lightcurve Extraction
Wolf 503 was observed by Kepler from 2018 March to 2018 May. We extract photometry from K2 's target pixel file (TPF) using the Lightkurve package (Lightkurve Collaboration et al. 2018). TPF's are the main data product of the Kepler/K2 and TESS missions consisting of stacks of "postage stamp" frames centered on the target star. Each frame represents one timestamp (or cadence) in which data was taken. For Kepler/K2 short cadence, the sampling rate is about a minute between exposures whereas long cadence only samples every 30 minutes.
After the failure of two of Kepler's reaction wheels, the solution that allowed K2 to be possible resulted in the target stars drifting across the detector over the length of the campaign. This drift causes changes in flux levels and needs to be corrected for. Lightkurve implements the Self Flat-Fielding (SFF) technique introduced by Vanderburg & Johnson (2014) to account for the motion of the Kepler spacecraft. Aperture photometry was performed on the TPF using the a circular pixel mask of radius 5 pixels centered on the star. We experimented with various aperture sizes ranging from 4 to 6 pixels. The 5 pixel radius produced the lowest outof-transit spread in the data after SFF was applied. The result of the self flat-fielding technique is show in Figure 2 with red tick marks indicating clear transit events with the exception of the first and tenth transits which suffered from thruster burns. There remains occasional decreases in flux between transits that are not periodic and are likely due to extreme differences in pixel sensitivity across the detector. Since these points are not explicitly used in fitting process they have no impact on the parameters derived in the following section.

Photo-eccentric Model
We fit the short-cadence data using the exoplanet package (Foreman-Mackey et al. 2021) which uses a Hamiltonian Monte Carlo (HMC) routine to explore the posterior probability distribution. We minimize a negative log-likelihood function using the period of the orbit (P ), time of inferior conjunction (T conj ), impact parameter (b), scaled planet radius (R p /R * ), and stellar density assuming a circular orbit (ρ * ,circ ) as free parameters and use a quadratic limb darkening law with the parameters held at u 0 = 0.5916 and u 1 = 0.1322 obtained from Claret & Bloemen (2011). Using the values obtained from the minimization algorithm, we initialized the HMC sampler with four parallel chains running 8,000 tuning steps and 6,000 sampling steps. Loose Gaussian priors were placed on P and T conj and instead of sampling directly in ρ * ,circ , we reparameterize according to Sandford & Kipping (2017) using log 10 (ρ * ,circ ) with a uniform prior.
With a median R p /R * value of 2.534 ± 0.020%, our work agrees reasonably well with the previous analysis of this system. However, we obtain ρ * = 16±1 g cm −3 ; higher than what is expected using values found in Section 2. This large stellar density can be explained in terms of the photo-eccentric effect and indicates that Wolf 503b is not on a circular orbit but an eccentric one.
While a planet on a circular orbit has a constant velocity, a planet in an eccentric orbit has a maximum velocity at periapsis and a minimum at apoapsis. This creates a dependence of the transit length on the argument of periapse (ω) of the orbit. An observer viewing the transit at periapsis (ω = 90 • ) would record a shorter transit than an observer at apoapsis. From Winn (2010) the length of the transit (T ) can be given as, where g is expressed as, By employing Kepler's 3rd Law, we can substitute a/R * in favor of the stellar density and obtain a key equation from Kipping (2010), A planet on a circular orbit will give ρ circ = ρ * . However, a planet on an eccentric orbit, transiting near periapsis, would give a larger ρ circ compared to an independent measurement of the star's density.
Using the methodology of Dawson & Johnson (2012), we sample values of e and ω from uniform distributions on the interval [0, 1] and [−π, π], respectively and compute the g parameter using Equation 2. These values of g are then used together with the ρ * ,circ to calculate what true density is implied from Equation 3 and is then compared to the value found with isoclassify using a likelihood function. The likelihoods are then used to reweight the samples of e and ω in order to obtain a an estimate of the eccentricity yielding a 1σ range of 0.59 to 0.82.

SPITZER PHOTOMETRY
After the discovery of Wolf 503 b, we were awarded Director's Discretionary Time (Crossfield et al. 2019) to observe the planet's transit with Spitzer.
On 2019/11/10 we observed one transit using the 4.5 µm channel (IRAC2, Fazio et al. 2004) with 2.0 s integrations taken in subarray mode; the transit observation encompassed 208 frames and spanned 7 hr 27 min. In addition, we acquired a short observation before and after the transit to check for bad pixels. Our observations were scheduled following standard best practices for precise Spitzer photometry, including using Peak-Up mode to place the star as closely as possible to the well-characterized "sweet spot" of the IRAC2 detector.

POET Reduction Pipeline
To extract photometry from the Spitzer observations, we use the Photometry for Orbits Eclipses and Transits (POET 1 ) package (Cubillos et al. 2013;May & Stevenson 2020). In summary, POET creates a bad pixel mask and discards bad pixels based on the Spitzer Basic Calibrated Data (BCD). Outlier pixels are also discarded using sigma-rejection. Then, the center of the point spread function (PSF) is determined. POET provides multiple routines to determine the PSF center and since we see no evidence of any source near Wolf 503 (Peterson et al. (2018)) we opt for a simple 2-D Gaussian fitting technique. After the center of the PSF is found, interpolated aperture photometry is used to extract the lightcurve. The resulting data is then fit with a model that accounts for both the lightcurve itself in addition to a ramp-like trend attributed to "charge trapping" (discussed in 4.3) and the sub-pixel sensitivity of the detector. The posterior distribution is sampled using an MCMC algorithm with chains initialized at the best fit values.

Interpolated Aperture Photometry
The quality of the fit is dependent not only on the model and the aperture size, but also the method of interpolation and the bin-size used. To find the best result, we tested various aperture sizes (ranging from 2 -6 pixels in increments of 1 pixel) with both nearest neighbor (NNI) and bilinear (BLI) interpolation using different bin sizes (0.1, 0.03, 0.01, and 0.003 square). For each case, the standard deviation of the normalized residuals (SDNR) was calculated and compared. The method resulting in the lowest SDNR was an aperture size of 5 pixels with interpolated photometry performed with bilinear interpolation and a bin size of 0.03 x 0.03 pixels.

Spitzer Systematics
At 4.5 µm, the primary systematic effect is the subpixel sensitivity variation causing the measured flux to be dependent on the target's position on the array (Stevenson et al. 2012;Charbonneau et al. 2005;Cubillos et al. 2013). To mitigate this variability, we utilize the BiLinearly-Interpolated Subpixel Sensitivity (BLISS) mapping described in Stevenson et al. (2012) which has been shown to be a more effective method Figure 2. The top panel shows the raw lightcurve extracted from the short cadence target pixel file using aperture photometry while the bottom panel shows the detrended result after employing the Self Flat-Fielding (SFF) technique to account for Kepler's motion during the K2 campaigns. Red tick marks indicate the transits of Wolf 503 b and the blue overlay shows which points were used in the fitting process. The 1st and 10th transits were omitted as they coincided with thruster burns.
of mapping the sub-pixel sensitivity of the detector as compared to polynomial fits or the weighted sensitivity function of Ballard et al. (2010).
There is also a temporal systematic that induces a ramp-like trend in the extracted lightcurve (see Figure  4a). The cause of this is thought to be charge trapping (Agol et al. 2010) and is an issue especially for brighter targets. During read-out of the detector, not all electrons are drained from the pixel leaving to be "trapped" in the pixel. As the observation progresses, the electrons build up increasing the effective gain of the detector. This manifests itself as a ramp-like trend in the lightcurve. To account for this, our fitted model uses a linear trend of the form: Where r 1 is the slope of the linear model and a free parameter which is fit for. t 0 is a constant term approximated as the mid-point phase of the transit and, in our case, set to be 1.0.

Lightcurve Fitting
The model lightcurve is generated using the batman package. Since we observed only one transit event with Spitzer (as opposed to the eight analyzed from K2 ) the uncertainties on transit parameters will necessarily be larger. The primary advantage of observing Wolf 503 b with Spitzer is making accurate predictions of future transits (Section 4.5).
After performing a least squares fit, an MCMC routine is initialized on the best-fit values and allowed to go though 100,000 iterations taking the first 3,000 as burn in. The parameters involved the analysis are R P /R * , T conj , a/R * , cos(i), and the parameters in Eq. 5. A Gaussian prior was placed on a/R * informed by the value found from our radial velocity analysis in Section 5.6. We use a quadratic limb darkening law with the parameters held at u 0 = 0.0973 and u 1 = 0.1276 obtained from Claret & Bloemen (2011). To increase computational efficiency, the eccentricity and argument of periastron are held at values found in 5.6. The median posterior value for the transit depth is given in table 2 and the lightcurve fit to the photometry is shown in Figure 4b.

Ephemeris Improvements
With K2 photometry alone, the uncertainties on a planet's period (P ) and the time of inferior conjunction (T conj ) degrade our ability to predict transits in the future. The uncertainty in mid-transit time (T n ) scales linearly with the number of orbits (n) since the initial observation (Beichman et al. 2016):  By the time JWST is operational (2022), the 3-sigma uncertainty in the transit time, calculated from the longcadence period alone, would be nearly 2 hours. For a planet whose total transit time lasts little over an hour, there is a likely chance we would only observe a partial transit, or in the worst case, miss the transit entirely. Both of these situations are unacceptable uses of valuable telescope time making precise knowledge of when a transit will occur crucial for future follow up studies.
In order to tighten our constraint on the mid-transit time we use the conjunction time obtained with both K2 and Spitzer and use a weighted least squares routine to obtain a more precise value of the period. A weighted least squares method is used instead of a joint fit to both the K2 and Spitzer photometry as the systematics in the Spitzer data tend to be so strong that it is difficult to model the systematics independently of the transit itself. We obtain a new period of 6.001274 ± 2.1e − 05 days. With the period obtained from the short cadence K2 photometry analyzed in this work, the precision improves to 34 minutes and with the addition of the Spitzer transit we ultimately come to a mid-transit time precision of just 21 minutes; a 5-fold improvement from the K2 long cadence prediction.

RADIAL VELOCITY ANALYSIS
We obtained radial velocity measurements of Wolf 503 from four instruments: the Keck Observatory's High Resolution Echelle Spectrometer (HIRES, Vogt et al. 1994), the Calar Alto high-Resolution search for M dwarfs with Exoearths with Near-infrared and optical Echelle Spectrographs (CARMENES, Quirrenbach et al. 2014Quirrenbach et al. , 2018, the High Accuracy Radial velocity Planet Searcher North (HARPS-N, Cosentino et al. 2012), and the Planet Finder Spectrograph (PFS, Crane et al. 2006Crane et al. , 2008Crane et al. , 2010. Observations were taken from May 2018 to March 2020 totalling 110 data points. In the following sections we describe the observations and reductions performed for each instrument and the subsequent analysis. The radial velocity points are available in Table 3 and displayed in Figure 7.

PFS Spectroscopy
We observed Wolf 503 with PFS from UT 24 May 2018 to UT 03 August 2018 with each exposure totaling 20 minutes producing 42 velocity measurements. The mean internal uncertainty is 1.53 m s −1 .
PFS is an iodine cell-based precision RV spectrograph installed on the 6.5m Magellan Clay telescope with an average resolution of R 130,000. RV values are measured by placing a cell of gaseous I 2 , which has been scanned with the NIST FTS spectrometer (Nave 2017) at a resolution of 1 million, in the converging beam of the telescope. This cell imprints the 5000-6200Å region of the incoming stellar spectra with a dense forest of I 2 lines that act as a wavelength calibrator and provide a proxy for the point spread function (PSF) of the spectrometer (Marcy & Butler 1992). The resulting spectra are split into 2Å chunks, each of which is analyzed using the spectral synthesis technique described in Butler et al. (1996), which deconvolves the stellar spectrum from the I 2 absorption lines and produces an independent measure of the wavelength, instrument PSF, and Doppler shift. The final Doppler velocity from a given observation is the weighted mean of the velocities of all the individual chunks (∼800 for PFS). The final internal uncertainty of each velocity is the standard deviation of all 800 chunk velocities about that mean.

HIRES Spectroscopy
A total of 27 radial velocity observations of Wolf 503 were obtained from the HIRES spectrograph during the period of 2018 May to 2019 April. HIRES is an iodine (I 2 ) cell-based spectrograph installed on the 10-m Keck I telescope capable of resolutions of R 50,000 operating between 360-800 nm. Observations were made in collaboration with the California Planet Search (CPS). Spectra were taken with the 14" by 0.861" "C2" decker with exposures averaging 17 minutes in order to reach the requisite signal to noise ratio of 200 per pixel. We obtained an average signal-to-noise ratio of 223 at 550 nm and an average internal velocity error of 1.08 m s −1 . Spectra were reduced and radial velocities calculated as described in Howard et al. (2010).

CARMENES Spectroscopy
We obtained 21 high-resolution spectra of Wolf 503 between June 2018 and July 2018 with the CARMENES instrument mounted on the 3.5 m telescope at the Calar Alto Observatory, Almería, Spain, under the observing program S18-3.5-021 (PI Pallé). The CARMENES spectrograph has two arms, the visible (VIS) arm covering the spectral range 0.52-0.96µm and a near-infrared (NIR) arm covering the spectral range 0.96-1.71µm.
Here we use only the VIS channel observations to derive radial velocity measurements. All observations were taken with exposure times of 1200 s resulting in SNR per pixel at 745 nm of CARMENES VIS spectra in the range 41-131. CARMENES performance, data reduction and wavelength calibration are described in Trifonov et al. (2018) and Kaminski et al. (2018).
Relative radial velocity values, chromatic index (CRX), differential line width (dLW), and Hα index values were obtained using serval 2 . For each spectrum, we also computed the crosscorrelation function (CCF) and its full width half maximum (FWHM), contrast (CTR) and bisector velocity span (BVS) values, following Lafarga et al. (2020). The RV measurements were corrected for barycentric motion, secular acceleration and nightly zero-points. Due to the low declination of the star (δ = −6.14 deg), Wolf 503 was observed from Calar Alto at relatively high airmasses (ranging from 1.5 to 2.1), which has a high impact on the telluric contamination of the spectra. Therefore to achieve the highest RV precision, we correct the spectra from telluric absorption using Molecfit Kausch et al. 2015) following the method presented in Nortmann et al. (2018) and Salz et al. (2018).

HARPS-N Spectroscopy
We collected a total of 20 radial velocity observations of Wolf 503 between June 2018 and March 2020 with the HARPS-N spectrograph installed on the 3.6-m Telescopio Nazionale Galileo (TNG) at the Observatorio del Roque de los Muchachos in La Palma, Spain. These observations were part of the HARPS-N Collaboration's Guaranteed Time Observations (GTO) program. Wolf 503 has an apparent magnitude V = 10.28, so we obtained spectra with signal-to-noise ratios in the range SNR = 41 -128 (average SNR = 83), at 550 nm in 30 minute exposures, depending on the seeing and sky transparency. A summary of the observations is provided in Table 3. The average internal RV error of the observations is 1.19. ± 0.46 m s −1 .
The spectra were reduced with version 3.7 of the HARPS-N Data Reduction Software (DRS), which includes corrections for color systematics introduced by variations in seeing (Cosentino et al. 2014). The radial velocities were computed using a numerical weighted mask following the methodology outlined by Baranne et al. (1996). Masks are chosen based on the closest spectral type of the star and in this case the K5 mask was chosen.

Stellar Activity & Rotation
The Mount Wilson S HK index is a commonly used metric of chromospheric activity defined as the ratio of flux in the Ca II H & K line cores (3968.5Å and 3933.7Å, respectively) to the flux in the nearby continuum (Wilson 1963;Duncan et al. 1991). As part of the California Planet Search (CPS), Isaacson & Fischer (2010) compiled a catalogue of S HK values for over 2,000 stars. A key finding was that K-dwarfs with a color index 1.0 < B-V < 1.3 produce the lowest levels of velocity noise that tend to mimic the radial velocity signals of a planet (known as "jitter"). Wolf 503 finds itself in this color range with a B-V color of 1.02 suggesting that it is a particularly good radial velocity target. This is confirmed by comparing Wolf 503's S HK index to those of Isaacson and Fischer's sample within that color range. High Resolution Echelle Spectrometer (HIRES) measurements give a median S HK for Wolf 503 of 0.246, much lower than the sample average of 0.536 indicating that Wolf 503 is chromospherically quiet.
Using Generalized Lomb-Scargle (GLS) periodograms (Zechmeister & Kürster 2009), we verify that the 6 day signal is present in the radial velocity data (Appendix, Figure 9). The planetary signal is prominent in both NUV-V colors and effective temperatures for nearby (within 50pc) K-dwarfs from the TESS Input Catalog. Most stars are in a quiescent group signifying low NUV emission while some stars are outliers. A number of these outliers (e.g HIP 106335, HD 8049) are known to have higher stellar activity. Wolf 503 (red star) is a member of the quiescent group which agrees with measurements of stellar activity such as SHK .
HIRES and HARPS-N data, but is not clearly seen in the data from PFS or CARMENES.
Due to Wolf 503's low S HK , we do not expect stellar activity to impact radial velocity measurements. To verify this, Figure 5 compares Generalized Lomb-Scargle (GLS) periodograms of S HK measurements against the full radial velocity data set (Table 3) to search for stellar activity on the timescale of the orbital period. Data acquisition is described in Section 5. The S HK periodogram contains a number of low frequency peaks (below 0.1 d −1 ). The presence of features such as starspots on the stellar surface have the ability to mimic radial velocity signals with periods that can reflect that of the star's rotational period. Peterson et al. (2018) reported a projected rotational velocity of v sin i * = 0.8 ± 0.5 km s −1 . Combining this with the stellar radius of 0.690 ± 0.02 R we obtain a maximum rotation period of 43 ± 27 days. This wide window of possible rotational periods coincides with the low frequency peak structures seen in both the S HK periodogram and RV periodogram ( Figure 5) but these peaks do not coincide.
Another method of assessing the activity of K-dwarf stars is by comparing their high energy flux, in particular near ultra-violet (NUV). At birth, stars have strong magnetic fields and large high energy emissions. As the star ages, a decay in the rotation rate causes a subsequent decrease in this high energy emission. Since this decrease is thought to begin rather quickly, approximately 100 Myrs after formation (Richey-Yowell et al. 2019), many K-dwarfs should fall into a quiescent group with low NUV flux. We test this by selecting TESS Input Catalog (TIC) K-dwarfs within 50 pc in the effective temperature range 3850 K < T ef f < 5340 K corresponding to spectral types K9V-K0V (Pecaut & Mamajek 2013). This list is then cross-referenced with the stars in the Galaxy Evolution Explorer (GALEX ) catalog (Bianchi et al. 2017) to obtain their NUV magnitudes.
In Figure 6 we take the NUV-V colors versus T ef f which shows the majority of K-dwarfs have large NUV-V colors forming a quiescent group. NUV-V increases with decreasing stellar temperature, a result that is consistent with a study of M-dwarfs made by Lépine et al. (2013). We note that Wolf 503, represented as a red star in Fig. 6, has NUV-V = 7.94 placing it within the quiescent group. In contrast, we find a number of K-dwarfs with lower NUV-V colors that could be considered active. Among these few stars, some are known to be active such as HIP 106335, identified to be an "active/fast rotator" by Santos et al. (2011). Additionally, HD 8049 has a high (relative to K-dwarfs, e.g Isaacson & Fischer 2010) S HK value of 0.678 (Arriagada 2011) and is also found in the active group with an NUV-V = 4.76. Interestingly, some members of the active group, TYC 422-1303-2 among them with the lowest NUV-V, have gone largely unstudied.

RV-Only Analysis
The RV measurements are analyzed using the open source, orbit-fitting tool-kit RadVel (Fulton et al. 2018).
With RadVel, a model orbit is fit to the data with the orbital parameters being period (P , with a Gaussian prior informed by the value found in 4.5), time of inferior conjunction (T conj ), RV amplitude (K), eccentricity (e), and argument of periastron (ω). Other parameters that are fit include an RV offset (γ) and jitter (σ) terms for all instruments. During the fitting process, √ e cos ω and √ e sin ω are used in lieu of e and ω alone in order to avoid biasing the eccentricity.
Our analysis consists of comparing a simple model of a circular orbit to models with additional parameters such as eccentricity and a linear trend. An MCMC routine is initialized on best-fit values and used to determine the median value of the posterior distribution as well as obtaining an uncertainty for each parameters. As discussed in Section 5.5, neither short-term stellar activity nor rotation are expected to affect our results and so methods of mitigating those effects (e.g. Gaussian processes) are not implemented. We also consider the potential for the Rossiter-Mclaughlin (RM) effect to bias any RV measurements taken during transit. Using Equation 4 from Winn (2010), with the best case scenario of an impact parameter of zero, the maximum amplitude of the RM effect would be ∼ 0.6 m s −1 ; smaller than the average uncertainty for each instrument. In reality, Wolf 503 b likely has a high impact parameter (see Sec. 5.7) which renders any bias due to the RM effect even more negligible.
In order to measure the justification of any added parameters we utilized the Akaike information criterion (AIC). An AIC score allows us to compare the goodness of fit of different models while also taking into account over-fitting. The model that minimizes the AIC is considered optimum. The difference between the lowest AIC and the AIC of a model in question (∆AIC = AIC model − AIC min ) allows us to reject models that either poorly describe the data or contain too many parameters: ∆AIC < 2 shows little difference between the two models, 2 < ∆AIC < 10 indicates less support for the model, and a ∆AIC > 10 means the model is strongly disfavored.
When comparing models of circular and eccentric orbits both with and without acceleration terms we find that an eccentric orbit with a linear trend is by far the preferred model with a circular orbit being disfavored (∆AIC = 9.57) and any model without a trend included  Table 2 are the median values of the posterior distributions. We add in quadrature the RV jitter terms listed in Table 2 with the measurement uncertainties for all RVs. b) Residuals to the best fit 1-planet model. c) RVs phase-folded to the ephemeris of planet b. The small point colors and symbols are the same as in panel a. Red circles are the same velocities binned in 0.08 units of orbital phase.
being entirely ruled out (∆AIC = 17). Our analysis of the radial velocity data alone reveals an orbital eccentricity of 0.35±0.09. The discrepancy between this value and the one found from the K2 photometry is addressed with a joint RV-transit fit in Section 5.7.
A trend in RV data can indicate the presence of a long-period, massive companion however they can also be caused by long term stellar activity. We observe positive trends in both the S-index and Full Width at Half Maximum values from HIRES and HARPS-N suggesting that this trend is stellar in origin rather than evidence of planet 'c'. However, we also note that there exists only a slight correlation between the S-index values and the radial velocity measurements with a Pearson correlation coefficient of 0.26. Further monitoring of this system is likely needed to determine the nature of this trend.

Joint RV-Transit Analysis
The discrepancy between the eccentricity values predicted from the photo-eccentric modeling of the K2 photometry and from the RV data alone suggests that a joint RV-transit analysis may be necessary for Wolf 503 b. Often, the degeneracy between the impact parameter and eccentricity can result in small estimates of b  (Dawson & Johnson 2012). We attempt to resolve this discrepancy by modeling the photometry and radial velocity measurements simultaneously. Our joint model is constructed using exoplanet using the same parameters from the photo-eccentric and radial velocity models. Priors were placed on ρ * using the values in Table 1. Without the orbital information we gain from the RV analysis, the impact parameter derived from photometry is both small and unconstrained at b = 0.18 ± 0.11 but our joint model revises this value to b = 0.65±0.06 and produces a new, slightly higher eccentricity estimate of e = 0.41 ± 0.05. The scaled planet radius is also affected, due to the dependence on both b and e, increasing to 2.79 ± 0.05%. All other parameters remained consistent with the values found with either the photo-eccentric model or RV only model. A sum-mary of model and derived parameters is provided in Table 2.

DISCUSSION
From our radial velocity analysis, we find Wolf 503 b has a mass of 6.26 +0.69 −0.70 M ⊕ and, combining this with a radius of 2.043±0.069 R ⊕ , has a bulk density of 2.92 +0.50 −0.44 g cm −3 . These measurements allow us to place this planet in context and investigate its viability as a target for atmospheric characterization.

Interior Models and Formation Theories
Sub-Neptunes are typically described as low-density planets with modest H-He envelopes making up anywhere between 0.1-10% of the planet's mass. Super-Earths, on the other hand, are thought to be smaller planets with higher densities stripped bare of any enve- lope. With the newly acquired mass of Wolf 503 b we employ the Structure Model Interpolator tool (smint, Piaulet et al. 2020) to determine the envelope mass fractions for both a H 2 O dominated planet and one with an H-He envelope. smint uses the model grids of Fortney 2014 andZeng et al. 2016 to determine the mass fraction for H-He and water, respectively. The former model grids assume a core composed of a 2:1 mix of rock and Iron while the latter employs a 2-layer model reflective of Earth's core and mantle. We find that Wolf 503 b is entirely consistent with an Earth-like core of rock and iron with either an H 2 O mass fraction of 45 +19 −16 % or an H/He mass fraction (f H,He ) of 0.49±0.28%. These values are consistent with common definitions of sub-Neptunes.
The lack of sub-Neptunes orbiting the Sun means that we still have much to discover about their origins. Early investigations into planet formation focused on replicating the system architecture of our own Solar System and, even though the planets in our system can be formed at their current positions, it is generally accepted that this is not feasible for hot/warm sub-Neptunes through classic core-accretion (Bodenheimer & Lissauer 2014;Inamdar & Schlichting 2015;Venturini et al. 2020). Schlichting (2014) calculated the required enhancement of the minimum mass solar nebula (MMSN) needed for in-situ formation of sub-Neptunes for various masses and at varying distances from its host star. For a roughly 5 M ⊕ planet forming at 0.05 AU, the MMSN would need to contain 90 times more solid material in the inner-disk. If the metallicity of the host star is reflective of the refractory content of the proto-planetary disk then this enhancement requirement is especially unreasonable for Wolf 503 whose metal content is only 30% that of the Sun. Pebble accretion could offer an in-situ formation pathway for Wolf 503 b, however, pebble accretion may tend to form systematically drier planets since the pebbles should lose most of their volatiles during their journey to the inner disk (Oka et al. 2011;Ida 2019). Planetesimals that form beyond the snow line are more likely to retain their volatiles and can contain 10-50% water by mass (Izidoro et al. 2019) resulting in vastly different compositions for migrating and in-situ planets. Although planetary compositions derived from bulk density alone are degenerate, the high bulk water composition of Wolf 503 b could imply a formation beyond the snowline and subsequent migration inwards.
When considering planets of high T eq , it is also important to note that the planet we characterize today has evolved significantly since its formation. Planets hosted by relatively long-lived stars can provide insight into the end products of mass loss mechanisms such as photoe- vaporation (Owen & Wu 2013) and core-powered mass loss. Wolf 503 b, orbiting an 11 ± 2 Gyr old K-dwarf, likely experienced appreciable photoevaporation of its atmosphere. Neptune-class planets (M ∼ 20M ⊕ ) can have the majority of their atmosphere removed by its host star; in the most extreme cases the planet is left with H/He envelopes of fractions of a percent consistent with the H/He mass fractions found for Wolf 503 b. Much of the evaporation is thought to occur in the first 100 Myrs during a "saturation" phase early in the star's life when X-ray emission is at its peak and independent of rotation period (Owen & Wu 2017). Given the age of this system, Wolf 503 b could be an example of the end product of photoevaporation.

Potential for Atmospheric Characterization
Equipped with a precise mass measurement of Wolf 503 b, we are now able to more carefully consider the viability of this planet for transmission spectroscopy. A calculation of the transmission spectroscopy metric (TSM) of Kempton et al. (2018) places Wolf 503 b (TSM=63.9) in the top twenty best atmospheric follow-up targets in the size range 1R ⊕ < R P < 4R ⊕ (Guo et al. 2020). This immediately suggests Wolf 503 b as a potential target for atmospheric characterization.
A mass uncertainty below 20% decouples the similar effects that both high surface gravity and a high mean molecular weight composition have on atmospheric spectra allowing us to investigate the latter. However, with that degeneracy broken, we are potentially faced with another. The tentative correlation of the water absorption amplitude with T eq (Crossfield & Kreidberg 2017) suggest that hazes should be important considerations when modeling the atmospheres of planets like Wolf 503 b. A T eq of 790K places Wolf 503 b in a region where hazes might be commonplace for warm Neptunes, but its small size and low H/He mass fraction could indicate enhanced metallicity Fortney et al. (2013); Venturini et al. (2016). Both of these factors can have similar, flattening effects on transmission spectra.
We consider the ability of the James Webb Space Telescope (JWST) to distinguish between these effects by generating model spectra using ExoTransmit (Kempton et al. 2017) with varying degrees of aerosols ranging from a clear atmosphere to hazes with 100x and 1000x solar Rayleigh scattering or a cloud deck at 0.01 bar. For each of these aerosol compositions we simulated metallicities of 1x and 100x solar [M/H]. These spectra were then used to simulate JWST observations using NIRISS (Single object slitless spectroscopy covering 0.6−2.8µm) and NIRSpec (Bright object time series with G395H covering 2.87 − 5.27µm) instrument modes. Simulations were made with PandExo ) assuming a resolution of R = 35, a baseline equal in time to that of the transit and a zero noise floor. The model spectra were then smoothed to match the native resolution of the instrument and binned down to match the reso-lution of the simulated observations. Using a weighted least squares routine, each simulated JWST spectrum was then fit with both a linear model and the model spectra (including the model the simulation was generated from). The corresponding reduced χ 2 statistics and p-values (summarized in Tables 4 & 5) were calculated and used to compare the models.
The first question one would ask is whether these atmospheres are detectable (i.e is the linear model strongly rejected?). In this analysis we will consider a p-value > 0.05 to be a non-rejection of the model being fitted (or a non-detection in the linear case), 0.05 > p > 0.006 to be a weak rejection, and a p < 0.006 to be a strong rejection of the model. For NIRISS, each set of simulated spectra show an unambiguous detection with the exception of the cloud deck spectra for both metallicities ( Figure  8). This is not too surprising; the presence of clouds is expected to be a significant challenge when studying exoplanet atmospheres. Although, it is interesting to note that the cloudy 100x solar spectrum was a weak detection whereas its 1x solar counterpart was indistinguishable from the linear model. One would expect the combined effect of a cloudy, high mean molecular weight atmosphere would result in a stronger rejection than one of lower metallicity. For NIRSpec, the situation is slightly less optimistic. An aerosol-free composition was the only low metallicity atmosphere detectable but on the other hand, all high metallicity atmospheres were detectable.
Among the models that were detectable, we then ask whether these models are differentiable from one another. Both NIRISS and NIRSpec will be capable of distinguishing between atmospheres of different metallicity but NIRISS will be particularly useful for detecting possible hazes which is consistent with the results found by Batalha & Line (2017). For a solar metal content, NIRISS was able to resolve the differences in spectra due to various strengths of Rayleigh scattering, however for a metallicity 100 times solar the ability to detect these differences was lost. Cloud decks at lower pressures (higher altitude) would likely exacerbate this issue and, although not investigated here, we also have no reason to assume exoplanet atmospheres cannot contain both hazes and clouds potentially muting the effect of Rayleigh scattering.
From the TSM alone, Wolf 503 b proves to be a strong candidate for further atmospheric characterization. Our analysis shows that, at the very least, we could expect to differentiate a low mean molecular weight atmosphere from a higher one. Evidence of aerosols is also well within reach of JWST with a distinction between hazes and clouds being possible if the atmosphere has a close to solar metal content. Considering the increase in information to be gained from a low metallicity atmosphere, the relative metal-poorness of Wolf 503 b's host star only solidifies further this planet's potential as a follow-up target. Forming from a metal-poor disk may be helpful to keep the subsequent metallicity of the atmosphere low as well.

CONCLUSIONS
In this paper we characterized the sub-Neptune Wolf 503 b. Through radial velocity measurements we find that it is on an eccentric orbit (e = 0.41 ± 0.05) and determine its mass to be 6.26 +0.69 −0.70 M ⊕ . Employing stellar activity indicators, we find that the host star is indeed a well-behaved K-dwarf furthering this spectral class' reputation as the most amenable to radial velocity studies. We also compare Wolf 503 b to other K-dwarfs with recorded NUV measurements from the GALEX survey and find that it is a member of a large group with low NUV emission.
A joint analysis of previously unused short cadence K2 photometry and radial velocity data in combination with Gaia EDR3 data, provided us with a radius of 2.043±0.069 R ⊕ resulting in a bulk density of 2.92 +0.50 −0.44 gcm −3 . This low density helps confirm Wolf 503 b as a sub-Neptune with either a substantial H 2 O mass fraction of 45 +19 −16 % or an H-He mass fraction of 0.49±0.28%.
To enable future investigations of this planet, we utilized a Spitzer transit to further constrain ephemerides providing accurate transit predictions well into the JWST era. This analysis resulted in a 5-fold reduction in transit time uncertainty as compared to predictions made with values from Peterson et al. (2018).
We also explore the possibility of detecting a high metallicity atmosphere in addition to hazes finding that, in agreement with previous work by Batalha & Line (2017), that the NIRISS instrument will be an indispensable tool for atmospheric studies of Sub-Neptunes. The presence of clouds or the combination of strong haze effects with a high metallicity atmosphere understandably makes measurements less conclusive. We have found that Wolf 503 b offers itself as a good candidate for JWST follow-up observations and can act as a casestudy for planets orbiting old, metal-poor stars.