EXPRES. IV. Two Additional Planets Orbiting ρ Coronae Borealis Reveal Uncommon System Architecture

The detection of 51 Peg b (Mayor & Queloz 1995) was enabled by the ground-breaking radial velocity (RV) precision of ∼15 m s⁻¹ (Campbell et al. 1988; Walker et al. 1995). The RV technique continued to improve and soon reached a new state-of-the-art precision of ∼3 m s⁻¹ (Butler et al. 1996). Several additional gas giant planets in short-period orbits were detected, and the seventh of these was a Jupiter-mass planet in a 39.6 day orbit around ρ CrB b (Noyes et al. 1997).

In the ensuing years, lower-mass planets in wider orbits were detected. From the earliest days of exoplanet surveys, it was understood that obscuring RV scatter was correlated with strong chromospheric activity (Campbell et al. 1991; Queloz et al. 2001). However, it it was not clear whether the ultimate precision for radial velocity measurements would be set by instrumental stability or velocity noise from the stellar photosphere.

The High Angular Resolution Planetary Spectrograph (HARPS) achieved the next breakthrough in RV precision when it reached a single-measurement precision of 1 m s⁻¹ (Queloz et al. 2001; Mayor et al. 2003; Pepe et al. 2003). This precision improvement and a relatively high-cadence program made a big difference in the ability to detect lower-amplitude and multiplanet signals, epitomized by the addition of a sub-Neptune planet with a velocity amplitude of only 4 m s⁻¹to the μ Arae system (Butler et al. 2001; McCarthy et al. 2004; Santos et al. 2004; Pepe et al. 2007).

The latest generation of spectrographs, including the Echelle Spectrograph for Rocky Exoplanets and Stable Spectroscopic Observations (Pepe et al. 2021), the EXtreme PREcision Spectrograph (EXPRES; Jurgenson et al. 2016), and NEID (Schwab et al. 2016) have established a new state-of-the-art RV precision of 0.3 m s⁻¹. The story of ρ CrB illustrates the significance of this increased precision and high observing cadence.

Noyes et al. (1997) identified a 39.6 day, 67 m s⁻¹ sinusoidal signal around ρ CrB in only 41 nights of observations over one observing season. This planet was well within the capabilities of all of the planet-hunting spectrographs of the time. The Hamilton spectrograph at Lick Observatory had been upgraded to 3 m s⁻¹ precision in 1994, and ρ CrB was added to the program in 1997. Over the next 10 yr, an additional 26 of these higher-precision measurements had been obtained, refining the orbital parameters slightly but identifying no new signals (Butler et al. 2006). The rms scatter to the one-planet fit was 6.9 m s⁻¹, and the "stellar jitter" was estimated to be 3.9 m s⁻¹.

By the time of the Butler et al. (2006) follow-up on ρ CrB, there were multiple spectrographs with ∼1–3 m s⁻¹ measurement precision (Fischer et al. 2016, and references therein). ρ CrB was added to the target list at the Keck High Resolution Echelle Spectrometer (HIRES) and then the Automated Planet Finder (APF) as part of the Eta-Earth RV Survey (Howard et al.2009). Fulton et al. (2016) used 519 Keck/HIRES measurements and 157 APF velocities, taken over 8 yr, to identify a second planet. Despite their long time baseline, the Whipple and Lick velocities were omitted; their comparatively low cadence and low precision did not contribute to the significance of the detection. This second planet, ρ CrB c, was found to be a 25 M_⊕ planet on a 102 day orbit with an RV semi-amplitude of 3.74 m s⁻¹ (Fulton et al. 2016), a signal smaller than the "stellar jitter" reported a decade earlier. The revised estimate of the stellar jitter by Fulton et al. (2016) was 2.57 m s⁻¹. They saw no evidence of additional planets but mentioned possible structure to the residuals with a period ≳10 yr.

With every improvement in instrumental precision, our estimate of the "stellar noise floor" has changed. The residual velocity scatter attributed to stellar jitter has proven to be a combination of unknown instrumental errors, stellar variability, and multiple planetary signals (Isaacson & Fischer 2010; Luhn et al. 2020). With the new generation of extreme precision radial velocity (EPRV) spectrographs (Jurgenson et al. 2016; Schwab et al. 2016; Pepe et al. 2021), the instrumental contribution is often smaller than the stellar component (Isaacson & Fischer 2010; Luhn et al. 2020). We now know that nearly all stars have planets (Howard et al. 2010b; Lovis et al. 2011; Mayor et al. 2011; Burke et al. 2015; Hsu et al. 2019). However, our knowledge about the detailed planetary architectures are limited by our lack of sensitivity to a large portion of the mass–period (or radius–period) parameter space. Current statistics suggest that there are an average of three planets per stellar host and that the average planet size is between that of Earth and Neptune (Borucki et al. 2011; Foreman-Mackey et al. 2014; Winn & Fabrycky 2015; Christiansen et al. 2016; Foreman-Mackey et al. 2016; Winn 2018; Hsu et al. 2019, 2019; Christiansen et al. 2020; Zhu & Dong 2021, and references therein).

Planet multiplicity greatly increases the number of observations needed to disentangle the Keplerian signals (He et al. 2021). If those are obtained over long timescales, then different phases in the activity cycle of the star are sampled, making it more difficult to detect low-amplitude signals. Combining high cadence with high instrumental precision can help us identify the small signals that may be lurking in our data.

Radial velocity data has been collected for hundreds of stars for more than 25 yr (e.g., Cochran & Hatzes 1993; Butler et al. 2001; Queloz et al. 2003; Pepe et al. 2004; Howard et al. 2010a; Fischer et al. 2014). For the last 10 yr, the precision of these surveys has been ∼1 m s⁻¹ (Fischer et al. 2016). This precision, combined with relatively low observing cadence in most surveys has made it challenging to detect planets with M < 30M_⊕ on periods beyond ∼100 days. For systems with planets on orbits shorter than 40 days, we have still been limited to planets of a few Earth masses for G and K dwarfs, missing many of the smaller planets that have been identified by transit surveys. It is more difficult to disentangle multiplanet signals, especially when one or more are close to the instrumental precision (He et al. 2021), preventing detection of architectures similar to the solar system.

The 100 Earths Survey is a high-cadence EPRV survey (Brewer et al. 2020) designed to locate the small and intermediate mass planets in orbits that have so far eluded both transit and RV surveys. It specifically aims to detect planets in the habitable zones of Sun-like (G and K dwarf) stars. The program combines the extreme stability of the EXPRES spectrograph (Jurgenson et al. 2016; Blackman et al. 2020; Petersburg et al. 2020) with very high cadence to identify low-amplitude signals in multiplanet systems on relatively short timescales. The rapid switching capabilities of the Lowell Discovery Telescope (LDT; Levine et al. 2012; DeGroff et al. 2014) allows observing in quarter night increments, and has even been used for very-high-cadence observations with allocations of just 30 minutes. Brewer et al. (2020) showed that for HD 3651 it was possible to recover the known planet parameters with residuals of 58 cm s⁻¹ using only 60 observations taken over ∼5 months, roughly two orbital periods of HD 3651b. Here, we will present the first new detections by this program: two additional planets around ρ CrB (HD 143761). The architecture of the four-planet system differs from most previously detected systems in arrangement and variety of masses.

Since August 2019 the 100 Earths Survey has collected 163 observations of HD 143761 on 89 separate nights (Table 1). We use an exposure meter to stop all observations at the same signal-to-noise ratio (S/N) in order to reduce the effects of charge transfer inefficiency (CTI) on the radial velocities; at S/N = 300, CTI contributes about 10 cm s⁻¹ offset in the measured RV (Blackman et al. 2020). When the program began, we obtained three consecutive observations for ρ CrB , each with S/N = 250 at 500 nm. This strategy allowed us to evaluate the short-term RV scatter in our data for our target stars. This cadence was updated in 2020 with two consecutive observations at S/N = 310, improving the duty cycle of our observations while maintaining the same nightly RV precision. Clouds and bad seeing sometimes limit the number of exposures we can obtain, and we stop all exposures at 1200 s to ensure a reasonable correction for the barycentric velocity; integrations not reaching the required S/N in this time are discarded. Five observations of HD 143761 during 2020 were obtained using the lower per observation S/N, while the remainder were at the higher value as well as higher cadence. We increased the observing priority for this star in 2021. All observations that met the lower S/N requirement in less than 20 minutes were included in this analysis: 153 observations on 82 nights. The spectra were extracted using a flat relative optimal extraction, and radial velocities were derived using forward modeling (Petersburg et al. 2020).

During the night, we obtain laser frequency comb (LFC) observations every 20–30 minutes and Thorium–Argon emission lamp (ThAr) observations every hour. Our pipeline (Petersburg et al. 2020) uses the ThAr observations for a broad wavelength solution over the entire spectrum and uses the cross-correlation technique to derive an absolute radial velocity. We also derive a separate, more precise, wavelength solution using the LFC frames between ∼5000 and 7200 Å to construct a hierarchical, nonparametric model (Zhao et al. 2021). With these wavelengths, we use forward modeling of individual ∼2 Å chunks to derive precise relative velocities with single-measurement precision of 30 cm s⁻¹ at S/N = 250 (Brewer et al. 2020; Petersburg et al. 2020).

2.1. Stellar Properties

Following the procedure of Brewer et al. (2016), we use Spectroscopy Made Easy (SME) to forward model ∼350 Å of the spectrum to determine effective temperature, surface gravity, projected rotational velocity, and overall metallicity along with precise abundances for 15 elements (C, N, O, Na, Mg, Al, Si, Ca, Ti, V, Cr, Mn, Fe, Ni, and Y; Table 2). Some RV surveys use an iodine cell for wavelength calibration on each science frame, prohibiting use of those observations for abundance analysis. Because the wavelength calibration for EXPRES is performed in frames adjacent to our stellar observations, we are able to analyze every observation to constrain the uncertainties in the derived parameters.

Table 2. Stellar Properties for HD 143761

Parameter	Value
Identifier	HD 143761
Teff [K]	5817±24
$\mathrm{log}g$ [m s−2]	4.25±0.05
$v\sin i$ [km s−1]	0.8±0.3
V [mag]	5.39 a
B − V	0.61 a
BP [mag]	5.55±0.05 b
RP [mag]	4.75±0.05 b
Distance [pc]	17.497±0.015 b
[C/H]	−0.14±0.03
[N/H]	−0.29±0.04
[O/H]	+0.07±0.04
[Na/H]	−0.24±0.02
[Mg/H]	−0.11±0.01
[Al/H]	−0.07±0.03
[Si/H]	−0.16±0.01
[Ca/H]	−0.12±0.02
[Ti/H]	−0.08±0.01
[V/H]	−0.15±0.03
[Cr/H]	−0.24±0.02
[Mn/H]	−0.40±0.02
[Fe/H]	−0.20±0.01
[Ni/H]	−0.22±0.01
[Y/H]	−0.25±0.03
Mass [M⊙]	0.95±0.01
Luminosity [L⊙]	1.82±0.08
Radius [R⊙]	1.34±0.04
Age [Gyr]	10.2±0.5

Notes.

^a Leeuwen (2007). ^b Lindegren et al. (2021).

Download table as:  ASCIITypeset image

After deriving the stellar parameters, we combine the T_eff, $\mathrm{log}g$, [Fe/H], and [α/Fe] measurements with the V magnitude, B − V color, and Gaia DR3 distance measurement (Fabricius et al. 2021) to derive the stellar luminosity and radius (Table 2). We use those same inputs as prior constraints in our isochrone analysis to fit for the mass and age. In addition to the Yonsei-Yale isochrone (Y2 isochrones; Demarque et al. 2004) analysis as described in Brewer et al. (2016), we performed a more detailed analysis for the sake of comparison with the stellar properties from the previous planet detection around ρ CrB (Fulton et al. 2016).

Although most of the stellar properties are consistent with those reported in earlier studies, there are a couple of items of interest in the context of its planets. As an early system with a hot Jupiter, ρ CrB stands out with its relatively low metallicity of [Fe/H] = −0.2. This makes it a rarity in the context of the giant planet metallicity correlation (Fischer & Valenti 2005). The abundance pattern itself is also interesting, with oxygen significantly elevated compared to scaled solar values. However, an investigation into possible causes for this lies outside the scope of this paper.

In addition to having subsolar metallicity, the star is evolved. It is an early subgiant with a radius of 1.34 R_⊙ and near solar mass. However, there is slight discrepency in the T_eff, [Fe/H], and mass between our analysis and those used in Fulton et al. (2016) based on the interferometric measurements of Braun et al. (2014). In that analysis, they find the T_eff is ∼200 K cooler, it is 0.1 dex more metal poor, and the mass of 0.89 M_⊙ 6%–10% smaller. That temperature falls outside the distribution of literature values, which have a mean of 5798 ± 58 K that is consistent with our value.

2.2. Isochrone Analysis

Over its main-sequence lifetime, elemental diffusion causes an apparent depletion in surface abundances (Dotter et al. 2017). Determining the early main-sequence properties of a turnoff star, where this depletion is at its maximum, requires a more careful analysis. The isochrones ⁷ (Montet et al. 2015) Python package using the Modules for Experiments in Stellar Astrophysics (MESA) Isochrones and Stellar Tracks (MIST; Choi et al. 2016) allows us to determine the initial parameters of the star and provides consistent evolutionary tracks to look at the main-sequence behavior. Using our spectroscopic parameters combined with Gaia BP and RP colors and parallax, we found that the best fit was a 0.95 ± 0.01 M_⊙ at 10.2 ± 0.5 Gyr with initial [Fe/H] = −0.11 (Figure 1). The mass and age are close to that derived from the Y2 isochrones (0.98 ± 0.01M_⊙ and 9.5 ± 0.4 Gyr). Both isochrone analyses derived a slightly lower $\mathrm{log}g$ (4.16) than that derived via spectroscopy.

Zoom In Zoom Out Reset image size

We then performed the same analysis using the parameters from Fulton et al. (2016) and find a mass of M = 0.90 ± 0.02M_⊙ at an age of ${12.4}_{-0.52}^{+0.73}$ Gyr and an initial metallicity of [Fe/H] = −0.26. At this age, that would imply this star was in one of the first generations of stars in the Milky Way, and yet it is relatively metal rich. Although this is possible, we find that our parameters cover a more reasonable age range and so adopt the parameters in Table 2 for the rest of this work. The adopted distance from the prior work (closer by 0.25 pc) may have influenced their adoption of parameters for a smaller, cooler star. The ∼5% higher stellar mass will result in a 4% increase in the planet minimum masses.

We use a Lomb–Scargle periodogram to identify the strongest period in the velocities with a false-alarm probability (FAP) below 10⁻³. We fit and subtract a Keplerian signal with that period then repeat the procedure until there are no periods with FAP <10⁻³. Using the resulting parameters as starting values, we perform a simultaneous fit for the four Keplerian signals. For the two new planets, we placed a Gaussian prior on the eccentricity centered at 0.05 with σ = 0.05 in keeping with the other two planets. Adding a linear trend with a flat prior gave an improved fit with a trend of about −0.6 m s⁻¹ yr⁻¹. To test whether or not a stellar jitter term was warranted, we added jitter as a free parameter to our model with an initial value of 0.35 m s⁻¹, corresponding to the intranight scatter for this star. The resulting model has nearly identical parameters to the previous best-fit model with the exception of planet d, which has a period two days shorter. The reduced χ² of this model was 1.3. Model comparison analysis showed that a model with eccentricity fixed to zero for the longest-period planet (d) was preferred over allowing it to vary. This is likely due to the incomplete phase coverage in our data for this planet; the period is longer than an observing season.

We use this model to set the priors on a Markov Chain Monte Carlo (MCMC) analysis using the python software package RadVel (Fulton et al. 2018) to obtain uncertainties on the parameters. To keep the scales of the parameters comparable, we fit with some parameters combined and scaled: $\sqrt{e}\cos \omega$, $\sqrt{e}\sin \omega$ and $\mathrm{log}K$. The MCMC analysis used 50 walkers and the first 300,000 steps of burn-in were discarded; the chains were well mixed after 1,920,000 steps. The posterior distributions were relatively tight and generally Gaussian in shape (Appendix B).

The highest peak in the residuals to a three-planet fit is at 12.9 days and has an FAP = 1.0 × 10⁻⁶ (Figure 2). After adding a fourth Keplerian at 12.9 days, there are no significant peaks in the periodogram of the residuals.

Zoom In Zoom Out Reset image size

We used ΔBIC and ΔAIC to compare two, three, and four-planet models. In each case, the model with the additional planet was overwhelmingly preferred over the simpler system. Our final four-planet fit (Table 3) has a ΔBIC_3p−4p = 13 compared to the best three-planet fit. The derived planetary masses and semimajor axes can be found in Table 4.

Table 3. MCMC Posteriors

Parameter	Credible Interval	Maximum Likelihood	Units
Orbital Parameters
Pb	39.8438 ± 0.0027	39.844	days
Tconjb	55479.62 ± 0.29	55479.6	JD
Tperib	${55498.7}_{-39.0}^{+0.75}$	55499	JD
eb	0.038 ± 0.0025	0.0379
ωb	$-{1.577}_{-0.058}^{+0.06}$	−1.577	radians
Kb	67.28 ± 0.19	67.28	m s−1
Pc	${102.19}_{-0.22}^{+0.27}$	102.18	days
Tconjc	${55629.3}_{-11.0}^{+8.6}$	55629.5	JD
Tperic	${55609}_{-14}^{+13}$	55609	JD
ec	${0.096}_{-0.054}^{+0.053}$	0.09
ωc	${0.17}_{-0.54}^{+0.44}$	0.16	radians
Kc	${4.0}_{-0.17}^{+0.18}$	4.02	m s−1
Pd	${282.2}_{-3.7}^{+2.2}$	281.4	days
Tconjd	${55583}_{-31}^{+54}$	55596	JD
Tperid	${55512}_{-32}^{+55}$	55525	JD
ed	≡0.0	≡0.0
ωd	≡0.0	≡0.0	radians
Kd	2.2 ± 0.14	2.21	m s−1
Pe	12.949 ± 0.014	12.949	days
Tconje	${55498.5}_{-4.7}^{+4.6}$	55498.3	JD
Tperie	${55496.2}_{-5.0}^{+5.7}$	55495.4	JD
ee	${0.126}_{-0.078}^{+0.054}$	0.073
ωe	$-{0.01}_{-1.0}^{+0.86}$	−0.01	radians
Ke	1.14 ± 0.015	1.142	m s−1
Other Parameters
γ	$-{0.74}_{-0.33}^{+0.32}$	−0.75	m s−1
$\dot{\gamma }$	$-{0.00146}_{-0.00047}^{+0.00048}$	−0.00144	m s−1 d−1
$\ddot{\gamma }$	≡0.0	≡0.0	m s−1 d−2
jitter	${1.111}_{-0.063}^{+0.067}$	1.054

Note. 1,620,000 links saved. Reference epoch for γ, $\dot{\gamma }$, $\ddot{\gamma }$: 58982.73654927462. Dates are Barycentric JD—2,400,000.

Download table as:  ASCIITypeset image

Table 4. Derived Posteriors

Parameter	Credible Interval	Maximum Likelihood	Units
${M}_{b}\sin i$	1.093 ± 0.023	1.066	MJup
ab	${0.2245}_{-0.0024}^{+0.0023}$	0.2219	AU
${M}_{c}\sin i$	28.2 ± 1.5	28.5	M⊕
ac	${0.4206}_{-0.0045}^{+0.0044}$	0.416	AU
${M}_{d}\sin i$	21.6 ± 2.5	19.5	M⊕
ad	0.827 ± 0.011	0.81	AU
${M}_{e}\sin i$	${3.79}_{-0.54}^{+0.53}$	3.71	M⊕
ae	0.1061 ± 0.0011	0.1049	AU

Download table as:  ASCIITypeset image

The resulting rms scatter in the residuals to median MCMC model is 1.20 m s⁻¹ (Figure 3), which is less than twice that of the quietest star in the 100 Earths Survey but 3 times our mean empirically derived single-measurement uncertainties for this star of 0.37 m s⁻¹ (Figure 4). The extra scatter is likely due to stellar activity and possibly additional unresolved planets. We also tested a model with the eccentricity of the fourth planet fixed to zero. All of the parameters were consistent with the 1σ uncertainties; the most notable change was a slightly lower eccentricity for the planet with the 102 day period. The ΔBIC = 2 shows little support for this model, and the resultant rms scatter was higher at 1.25 m s⁻¹. We chose to keep the model with the free eccentricity for planet (e), but note that additional data will help constrain the eccentricities of all planets.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.1. Dynamical Stability

The system is stable out to at least 10 million years for system inclinations of 90° (i.e., edge-on orbits) and 50° but becomes unstable with a system inclination of 20°. We ran a dynamical simulation of the system, with the parameters given in Table 3 and different inclinations, using REBOUND's implementation of the Wisdom–Holman Fast integrator (Wisdom & Holman 1991; Rein & Liu 2012; Rein & Tamayo 2015). Each time step was of 0.1 days, and we integrated out to 10 million years. Though there was some change in the orbital parameters for the four different planets, all orbits remained stable with no planets ejected, crossing orbits, or being sent into the host star for 90° and 50° inclination. Figure 5 shows the evolution of each planet's semimajor axis and eccentricity throughout the 10 million year simulation.

Zoom In Zoom Out Reset image size

With a system inclination of 90° (i.e., edge-on orbits for all planets), the eccentricities of the two lower-eccentricity planets, b and d, stay below 0.04 and 0.1 respectively throughout the simulation. Planet c, with an initial eccentricity of 0.1 ± 0.009, oscillates between 0 and 0.10 in eccentricity. Planet e's eccentricity ranges between 0.05 and 0.08 over the 10 million year simulation. The semimajor axes of planets b, c, and d change by much less than 1% throughout the simulation, while planet e's semimajor axis varies by less than 2%. With an inclination of 50°, the situation is nearly identical and the system remains stable over the 10 million years simulated. At an inclination of 20° the system is unstable; all planets experience increased variability in both eccentricity and position. Planet e undergoes changes in semimajor axis of up to 10% and eccentricity varying between 0 and 0.6.

Stellar rotation and activity cycles can often masquerade as planetary signals, although high-cadence observations mitigate this issue to some extent. Short-period planetary signals will necessarily have constant phase; limited spot lifetimes will induce phase variations on the rotational signal. Although this adds noise to high-cadence observations, the phase variations weaken the periodic signal. In particular, the ∼13 day signal has a low amplitude and potentially a harmonic of a longer rotation period. We examined activity indicators in the spectroscopic data and analyzed two independent photometric data sets to evaluate the robustness of the identified RV signals.

4.1. Spectroscopic Activity Indicators

As part of the spectroscopic analysis, we measured the H-α emission, H-α equivalent widths, and derived S-values, and during the radial velocity analysis measured the FWHM and bisector span of the cross-correlation function. All of these activity indicators are obtained for every observation. Periodogram analysis shows no significant periodicity in any of the indicators with the exception of the FWHM of the cross-correlation function (CCF; Figure 6).

Zoom In Zoom Out Reset image size

The CCF FWHM is a measure of line shape changes typically associated with stellar activity. We identified three significant periods at 28.9, 84.9, and 385 days. The 85 day signal is a possible alias of the 29 day signal. None of the periods we identified coincides with any of the periods in the radial velocities and the period of planet e is not a low-order harmonic of 29 days. Either of the slowest periods at 28.9 and 84.9 days are reasonable rotation periods for a subgiant of this mass (Nascimento et al. 2012; Saders & Pinsonneault 2013). As the 28.9 day period has a narrower and stronger peak (with an FAP ∼10 times lower), we assume that the longer period is an alias.

4.2. Photometric Variability

4.3. Stellar Rotation Period Search

The high-cadence and increased precision of the EXPRES 100 Earths Survey allows us to identify low-mass and long-period planetary signals. We have demonstrated that, with less than three years of data, we were able to find the two known planets around ρ CrB as well as two previously undetected planets. The lower-mass interior planet (e) fits well with those discovered in transit surveys but rarely found in prior RV surveys. The outer Neptune-mass planet (d) lies in a region of parameter space sparsely populated by either transit or RV surveys due to the prior technical limitations and biases of both. This opens a window into system architectures that have been missed in prior studies. In addition, this evolved star is showing conflicting indicators on rotation that could inform our understanding of late stage evolution of planetary systems.

5.1. Not Peas-in-a-Pod

Most of the planets in the solar-system orbit within a couple of degrees of orbital inclination of each other, with Mercury being an outlier at 6° away from the orbital plane. Despite this, the large orbital distances between the planets mean that a transit survey that saw one of our planets would likely miss several others, and the survey would have to run for decades to catch the outer planets. However, a commonly found system architecture in the Kepler survey (Borucki et al. 2011) is a compact system of multiple planets (Winn & Fabrycky 2015). Subsequent studies have found that, within these systems, the planets tend to have very similar radii (Weiss et al. 2018) and masses (Millholland et al. 2017); the so-called "peas-in-a-pod."

Recently, Millholland et al. (2022) have found that at some point, either the alignment or the tight period spacing breaks down. The next planet in the pattern is missing, even though that planet should have been detectable given its predicted period and radius. Their analysis was based on all systems with four or more planets with measured radii. Only four of their 59 systems had detected planets beyond 100 days.

To place ρ CrB in the context of those systems, we retrieved data from the NASA Exoplanet Archive for comparison. We have measured minimum planet masses, not radii, so we selected only confirmed systems with four or more planets that have measured masses. In Figure 7 we plot all those with no detected planets with orbital periods beyond 100 days. All but three of these systems were discovered via transits, and we see the same types of systems as in Millholland et al. (2022), which focused on Kepler compact multiplanet systems. Additional planets in systems like Kepler-80, Kepler-107, and Kepler-223 should have been detectable if the regular period spacing continued.

Zoom In Zoom Out Reset image size

We also selected all systems with four or more planets with measured masses from the Exoplanet Archive that have at least one planet with a period longer than 100 days; these would have been more challenging to detect via transits. The distributions of planets in these systems look very different (Figure 8) than the "peas-in-a-pod" systems. There are still some very uniform systems that, inside 100 days, look like the compact multiplanet systems. However, many have a wider dispersion in the masses, and the exterior planets are generally more massive. The outer planets also break the period spacing pattern of the transit system.

Zoom In Zoom Out Reset image size

ρ CrB has a warm Jupiter with an interior super-Earth and two exterior Neptune-mass planets. A similar scrambling can be seen in GJ 876 and 55 Cnc (Figure 8). As even 20 M_⊕planets such as ρ CrB d have eluded detection in previous RV surveys, it is possible that there are many more such systems. These previously overlooked architectures appear more like random assortments of planets than peas-in-a-pod. EPRV surveys such as the 100 Earths Survey run with EXPRES will be able to resolve outstanding questions about the compact multiplanet systems, as well as identify these new more widely spaced architectures.

5.2. Rotation and v sin i

5.3. Formerly Habitable Zone Neptune

Planet d resides in a 281 day orbit around this 0.95 M_⊙ star. Although this nine-month period raises hopes of a habitable zone planet, ρ CrB is currently a subgiant with a luminosity more than 75% greater than the Sun (Table 2). We can calculate the equilibrium temperature for planet (d):

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={\left(\displaystyle \frac{1-{A}_{B}}{1-{A}_{B,\oplus }}\right)}^{\displaystyle \frac{1}{4}}{\left(\displaystyle \frac{{R}_{* }}{{R}_{\odot }}\right)}^{\displaystyle \frac{1}{2}}{\left(\displaystyle \frac{{M}_{* }}{{M}_{\odot }}\right)}^{\displaystyle \frac{1}{3}}{\left(\displaystyle \frac{P}{{P}_{\oplus }}\right)}^{-\tfrac{2}{3}}\left(\displaystyle \frac{{T}_{\mathrm{eff}}}{{T}_{\mathrm{eff},\odot }}\right){T}_{{eq},\oplus }.\end{eqnarray} \tag{ 1 }$$

The giant planets in the solar system and Earth all have similar albedos (∼0.3). Absent additional information we will use Earth's albedo; for G-type stars, this is a reasonable assumption (Genio et al. 2019). We find T_eq = 345 K, almost 90 K hotter than the Earth. However, while the star was on the main sequence, this region would have been much more temperate. Using the 0.95M_⊙ MIST evolutionary track found using our best-fit isochrone, we calculated the stellar flux at the location of the planet over its main-sequence lifetime (Figure 9). Early in its life, planet d would have resided in the conservative habitable zone, and even out to 7 Gyr would have been in the optimistic habitable zone (Kane et al. 2016). More massive planets can accommodate a higher flux without entering a runaway greenhouse (Kopparapu et al. 2014), as illustrated by the 5 M_⊕ inner edge in Figure 9. Although we might not expect liquid water on the surface of giant planets, Neptune-mass planets in this region are currently rare.

Zoom In Zoom Out Reset image size

5.4. An Ancient Metal-poor Exoplanetary System

We have yet to detect planetary systems that look like the solar system. These systems have been hidden by lack of RV precision, low cadence, and bias in transit detections. A primary goal of the 100 Earths Survey is to determine the frequency of systems like ours. Only by fully exploring the parameter space will we be able to reveal the diversity and frequency of planetary system architectures. There are already several EPRV spectrographs with high-cadence exoplanet programs taking science data and several more in commissioning or planned (e.g., Jurgenson et al. 2016; Schwab et al. 2016; Thompson et al. 2016; Gibson et al. 2020; Pepe et al. 2021; Seifahrt et al. 2022). This work shows that a high-cadence survey with an extreme precision spectrograph like EXPRES can begin to give us those answers. In only a few observing seasons, we recovered two known planets around ρ CrB and identified two new planets: a hot super-Earth and a temperature Neptune-mass planet. The system has an uncommon architecture, hosting a series of close-in planets that have a wide range of masses.

We have had great success over the last 27 yr in finding thousands of planetary systems. Most of those discoveries look nothing like what we expected. Instead of solar-system analogs, we have found hot Jupiters, sub-Neptunes, super-Earths, and systems with tightly packed inner planets; none of which exist in the solar system. It might seem that systems like ours are rare, but despite the long search, limits in our precision and biases in our detection methods have left large regions of the mass and period parameter space unexamined. We are now beginning to probe that space. In these early days, we are not yet finding systems similar our own, but systems like ρ CrB are not like those that are so prevalent in transit surveys either. Revealing what lies hidden at longer periods can help us better understand planet formation and migration.

These results made use of data provided by the EXPRES team using the EXtreme PREcision Spectrograph at the Lowell Discovery telescope, Lowell Observatory. Lowell is a private, nonprofit institution dedicated to astrophysical research and public appreciation of astronomy and operates the LDT in partnership with Boston University, the University of Maryland, the University of Toledo, Northern Arizona University and Yale University. EXPRES was designed and built at Yale with financial support from NSF MRI-1429365, NSF ATI-1509436, and Yale University. Research with EXPRES is possible thanks to the generous support from NSF AST-2009528, NSF 1616086, NASA 80NSSC18K0443, the Heising-Simons Foundation, and an anonymous donor in the Yale alumni community.

We gratefully acknowledge the support for this research from NASA grant 80NSSC21K0009, NASA XRP-80NSSC21K0571, and NSF 2009528. R.M.R. acknowledges support from the Heising-Simons 51 Pegasi b Postdoctoral Fellowship. We also acknowledge the generous support for telescope time from the Heising-Simons foundation. G.W.H. acknowledges long-term support from NASA, NSF, Tennessee State University, and the State of Tennessee through its Centers of Excellence Program.

This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program.

This publication makes use of The Data & Analysis Center for Exoplanets (DACE), which is a facility based at the University of Geneva (CH) dedicated to extrasolar planets data visualization, exchange and analysis. DACE is a platform of the Swiss National Centre of Competence in Research (NCCR) PlanetS, federating the Swiss expertize in Exoplanet research. The DACE platform is available at https:dace.unige.ch.

A.1. TESS Photometry

ρ CrB was observed by TESS (Ricker et al. 2014) during Sectors 24, 25, and 51 at 2 minute cadence as well as during Sector 51 at a 20 s cadence. For the analysis discussed in Section 4.3, 2 minute observations are sufficient, so the 20 s cadence observations in Sector 51 are not considered. The TESS observing periods are only 27.4 days per sector, limiting our activity analysis with these data to periods shorter than ∼14 days.

We removed cotrending basis vectors (CBVs) from the simple aperture photometry (SAP) light curves. This process preserves signals from stellar astrophysics while accounting for systematic variations affecting TESS observations. For the 2 minute cadance observations the first four CBVs were removed. Data points with nonzero quality flags were also removed. From the 2 minute Sector 51 data, the data points between 2459712.8167376 and 2459712.9347927, inclusive, were removed due to likely systematic issues that were not addressed by the CBVs or the quality flags. The data with CBVs and flagged data points removed are shown in Figure 10. The data and CBVs were retrieved from the Barbara A. Mikulski Archive for Space Telescopes (MAST).

Zoom In Zoom Out Reset image size

A periodogram analysis of the TESS light curve was performed combining Sectors 24 and 25. Sector 51 was not included in this analysis due to the large gap between the sectors and the large data gap in the sector itself. The strongest signal in the periodogram of the TESS data occurs at 14.8 days. However, this corresponds to a strong peak in the window function and is where we would therefore expect to see a signal given the sampling of the data.

A.2. Automated Photoelectric Telescope Photometry

Ground-based photometry of ρ CrB was obtained at the Fairborn Observatory, AZ with the T4 0.75 m APT (Henry 1999). Observations were obtained of both ρ CrB and comparison star HD 139389 in Strömgren b and y pass bands. The data were combined into a single pass band of (b + y)/2, and differential magnitudes for ρ CrB are included in Table 5. Observations were obtained between 2000 January 13 and 2020 June 26, with an average of 54 data points in each of 21 observing seasons.

Table 5. APT Photometry

Heliocentric Julian Date	(b + y)/2 Differential	Trend
(HJD—24000000)	Magnitude	Removed
51557.0280	−1.0087	−1.0073
51572.0254	−1.0073	−1.0072
51579.0347	−1.0071	−1.0072
51585.9417	−1.0066	−1.0072
51586.9829	−1.0078	−1.0072
⋯	⋯	⋯

Only a portion of this table is shown here to demonstrate its form and content. A machine-readable version of the full table is available.

Download table as:  DataTypeset image

A trend was removed from the photometry to account for any signatures that might be long term, leaving those attributable to the rotation of surface features. The data were smoothed with a Gaussian kernel with a FWHM of 100 days (following Cabot et al. 2021; Roettenbacher et al. 2022). The data are included in Figure 11 and Table 5.

Zoom In Zoom Out Reset image size

We used an MCMC analysis to determine the uncertainties in the parameters. We used 50 walkers and the first 300,000 steps of burn-in were discarded; the chains were well mixed after 1,920,000 steps. The eccentricity of planet d was fixed to zero, and the remaining eccentricities were constrained to be less than 0.2. Gaussian priors were placed on the periods and RV semi-amplitudes of all planets with μ and σ set to the parameters and uncertainties returned from the initial maximum likelihood fit. The jitter used a Gaussian prior with μ = 0.35 and σ = 0.2.

In Figure 12 we plot the posterior distributions for the parameters of planets d and e along with the offset (γ), linear drift ($\dot{{gamma}}$), and jitter. The linear drift has been converted to m s⁻¹ yr⁻¹ to aid legibility. The distributions of the parameters for planets b and c are small and nearly Gaussian shaped and were also excluded to aid legibility.

Zoom In Zoom Out Reset image size

In addition, we have included the posteriors for the derived parameters (Figure 13).

Zoom In Zoom Out Reset image size