A Perfect Tidal Storm: HD 104067 Planetary Architecture Creating an Incandescent World

The original discovery of planet b orbiting HD 104067 by Ségransan et al. (2011) employed 88 RV observations obtained using the HARPS spectrograph. Since then, HARPS ceased monitoring of the target and the HARPS temporal coverage of this target remains around 2271 days from 2004 February to 2010 April. In 2020, a rereduction of all HARPS data carried out by Trifonov et al. (2020) corrected systematics, such as nightly zero-point RVs and intranight drifts, that slightly improved the overall HARPS precision. These revised HARPS reductions were utilized for our subsequent analysis, and the HD 104067 data yield a mean RV uncertainty of ∼0.89 m s⁻¹. In addition, the Keck RV program has been carrying out observations of HD 104067 for over two decades using the HIRES spectrograph, and were used in an independent analysis of the system by Rosenthal et al. (2021). In total, 72 HIRES RVs were collected, spanning from 1997 January to 2023 May (∼9600 days) with an average RV uncertainty of ∼1.28 m s⁻¹.

HD 104067 has been observed during TESS sectors 10, 36, and 63. To search for signs of stellar activity, we utilized the 2 minute cadence Simple Aperture Photometry (SAP) lightcurves that have been processed by the Science Processing Operations Center (Jenkins et al. 2016) and are available on the Mikulski Archive for Space Telescopes (MAST):10.17909/t9-nmc8-f686. Data points flagged as poor quality were removed, in addition to 5σ outliers. A Lomb–Scargle periodogram was used to estimate a stellar variability period of 18.3 ± 4.9 days from each single-sector lightcurve. We further discuss a possible interpretation of the observed variability in Section 3.1. To constrain the planet candidate's orbital and physical parameters, we used the Pre-search Data Conditioning (PDC) lightcurves, since the flux measurements are adjusted for dilution from nearby stars. We describe the preparation of the lightcurve and fitting process in Section 3.3.

Similar to the Kepler spacecraft, TESS is subject to known observing and instrumental systematics during observations (Jenkins et al. 2016). These include fluctuations in pauses in data collection after each orbit, changes in flux due to temperature changes, and periodic pointing corrections. While the pipeline that produces the PDCSAP lightcurves corrects these known systematics in the majority of TESS targets, some astrophysical changes in flux can also be removed. In particular, changes in flux that are on the order of the length of an individual sector tend to be flattened in the PDCSAP lightcurves. We find this to be the case for HD 104067, and thus have selected to alternatively analyze the SAP lightcurves. Some known systematics persist, such as changes in flux due to temperature fluctuations near gaps in the observations, but the long-period variations are sufficiently large in amplitude and persist over a multiyear baseline such that we consider these variations to be physically associated with the star.

Given the brightness of HD 104067, the star is also known by numerous aliases, such as GJ 1153, HIP 58451, TIC 428673146, and Gaia DR3 3494677900774838144. Ségransan et al. (2011) refer to HD 104067 as a moderately active K2V star with a spectroscopically determined rotation period of 34.7 days. We derived stellar properties for the star by applying the SpecMatch (Petigura et al. 2015) and Isoclassify (Huber et al. 2017) software packages to the obtained template Keck-HIRES spectrum. These parameters are provided in Table 1. In addition, we include in Table 1 the parallax and distances extracted from the third data release of the Gaia mission (Gaia Collaboration et al. 2021). In addition to the discussion of stellar properties in the context of the known exoplanet detection, provided by Ségransan et al. (2011) and Rosenthal et al. (2021), the star has also been the subject of other surveys. For example, Suárez Mascareño et al. (2015) spectroscopically determined a rotation period of 29.8 ± 3.1 days, roughly consistent with that found by Ségransan et al. (2011).

Table 1. HD 104067 Derived Stellar Parameters

Parameter	Units	Values	Source
V	V-band magnitude	7.93	This work, Specmatch
M⋆	Mass (M⊙)	0.82 ± 0.03	This work, Specmatch
R⋆	Radius (R⊙)	0.78 ± 0.01	This work, Specmatch
$\mathrm{log}g$	Surface gravity (cgs)	4.56 ± 0.10	This work, Specmatch
Teff	Effective Temperature (K)	4952 ± 100	This work, Specmatch
[Fe/H]	Metallicity (dex)	0.11 ± 0.06	This work, Specmatch
$v\sin i$	Projected rotational velocity (km s−1)	2.32 ± 1.0	This work, Specmatch
Prot	Stellar rotation period (days)	18.3 ± 4.9	This work, TESS
ϖ	Parallax (mas)	49.1470 ± 0.0235	Gaia DR3
d	Distance (pc)	20.35 ± 0.01	Gaia DR3

Download table as:  ASCIITypeset image

As mentioned in Section 2.2, we used a Lomb–Scargle periodogram to examine stellar variability that may be associated with the rotation period of the star from each individual sector of TESS photometry. Shown in Figure 1 are the photometry, periodogram analysis, and folded lightcurves for TESS sectors 10 (top row), 36 (middle row), and 63 (bottom row). These analyses reveal a consistent periodic variability signature of 18.3 ± 4.9 days for the three sectors, despite the observations each being separated by ∼1 yr. The second period identified in sector 36 is approximately half of the primary variability period, and thus is likely an alias. If the 18.3 days signal is indeed the signature of the stellar rotation period, then it is quite different from that derived spectroscopically by Ségransan et al. (2011) and Suárez Mascareño et al. (2015). A possible cause of this discrepancy in rotation period measurements is differential rotation, similar to that observed for the known exoplanet host HD 219134 (Johnson et al. 2016). Furthermore, there is growing evidence that spectroscopic and photometric rotation period measurements do not always agree, possibly the result of activity differences in the chromosphere and photosphere (Lafarga et al. 2021; Schöfer et al. 2022; Isaacson et al. 2024). However, there are several other pieces of information to note here. First, the R_⋆ and $v\sin i$ values shown in Table 1 result in a calculated rotation period of ∼17 days, which is more consistent with the rotation period derived from TESS photometry than those previously determined spectroscopically. Ségransan et al. (2011) measured a $v\sin i=1.61$ km s⁻¹; slower than the 2.32 km s⁻¹ value shown in Table 1, and enough to account for their larger rotation period of ∼30 days. Second, the TESS photometry is substantially separated in time from the HARPS RV data, such that the star may be located at a significantly different place within its magnetic activity cycle, as has been photometrically observed for numerous other stars (Henry et al. 2000; Dragomir et al. 2012). Third, it is worth acknowledging that aliases of the true rotation period can be difficult to disentangle, depending on the sampling rate of the observations, and so some of the reported rotation periods may be examples of such aliases.

Zoom In Zoom Out Reset image size

3.2.1. Radial Velocity Solutions

We utilized both HARPS and HIRES data in our RV analysis. In our modeling process, we divided the HIRES data set into two separate data sets (HIRES1 and HIRES2) due to a major upgrade carried out on HIRES in 2004 August (Butler et al. 2017) that could introduce a velocity offset between the two data sets. We first used an iterative planet search tool, RVSearch, to look for potential planetary signals within the combined data. We refer readers to Rosenthal et al. (2021) for details regarding the working process of RVSearch. We set the algorithm to search within a period space between 2 days and five times the combined observing baseline of HARPS and HIRES, which is around 48,000 days. The search process includes instrument offsets as well as linear and quadratic trends, and we only considered periodogram signals that have a false-alarm probability (FAP) less than 0.1%, as described by Howard & Fulton (2016) and Fulton et al. (2018). RVSearch returned two significant signals from the search: one belonging to the known b planet at ∼58 days with an FAP that is consistent with zero, and the other yielding a new period signature of ∼14 days with an FAP of 3.29 × 10⁻⁷. It is worth noting that the 14 days signal is primarily detectable within the HARPS data, which is unsurprising given the high cadence nature of the HARPS observations as well as their slightly higher RV precision compared to the HIRES data.

The RVSearch results were then used as input for the RV modeling toolkit RadVel (Fulton et al. 2018) to sample posteriors of Keplerian orbital parameters of the returned signals, as well as to estimate the associated uncertainties via a Markov chain Monte Carlo (MCMC) analysis. All parameters, including instrument offset and trends, were allowed to vary as free parameters. Chain convergences were evaluated under four criteria: Gelman–Rubin statistic (<0.01), minimum autocorrelation time factor (>40), maximum relative change in autocorrelation time (<0.03), and minimum independent draws (>1000). The MCMC chain rapidly converged, and we present the full results of our two planet model in Table 2. The planetary masses and semimajor axes were derived using the stellar mass shown in Table 1. The model includes both a linear and a quadratic trend, with significances of 3σ and 2σ, respectively. The RV data are shown in Figure 2, where the blue lines indicate the best fit to the data. The top panels show the full observational baseline and the residuals from the fit to the data. The bottom panels show the individual fits to the two planetary signatures detected in the RV data: the known 55.8 days period planet (planet b), and the newly detected 13.9 days period planet (planet c).

Zoom In Zoom Out Reset image size

Table 2. Radial Velocity Derived Planetary Parameters

Parameter	Credible Interval	Maximum Likelihood	Units
Orbital Parameters
Pb	55.851 ± 0.017	55.851	days
Tconjb	${2454167.3}_{-0.87}^{+0.94}$	2454167.0	BJD
Tperib	${2454159.5}_{-3.4}^{+3.5}$	2454159.7	BJD
eb	${0.123}_{-0.051}^{+0.048}$	0.136
ωb	${0.5}_{-0.41}^{+0.4}$	0.5	radians
Kb	12.0 ± 0.57	11.98	m s−1
Pc	${13.8992}_{-0.0037}^{+0.0047}$	13.8985	days
Tconjc	${2455192.9}_{-0.56}^{+0.59}$	2455192.64	BJD
Tperic	${2455191.6}_{-1.1}^{+1.2}$	2455191.6	BJD
ec	${0.29}_{-0.13}^{+0.12}$	0.32
ωc	${0.64}_{-0.56}^{+0.57}$	0.66	radians
Kc	${4.25}_{-0.63}^{+0.62}$	4.49	m s−1
Other Parameters
${\gamma }_{\mathrm{HIRES}2}$	$-{1.81}_{-0.74}^{+0.72}$	−1.59	m s−1
${\gamma }_{\mathrm{HIRES}1}$	${0.0}_{-2.8}^{+2.6}$	0.4	m s−1
γHARPS	${1.02}_{-0.82}^{+0.84}$	1.14	m s−1
$\dot{\gamma }$	${0.0013}_{-0.00043}^{+0.00042}$	0.00142	m s−1 day−1
$\ddot{\gamma }$	2.7e − 07 ± 1.3e − 07	2.7e − 07	m s−1 day−2
${\sigma }_{\mathrm{HIRES}2}$	${4.34}_{-0.52}^{+0.6}$	3.85	m s−1
${\sigma }_{\mathrm{HIRES}1}$	${8.2}_{-1.2}^{+1.1}$	8.1	m s−1
σHARPS	${4.05}_{-0.33}^{+0.37}$	3.83	m s−1
Derived Posteriors
${M}_{b}\sin i$	${62.1}_{-3.2}^{+3.3}$	61.5	M⊕
ab	${0.2674}_{-0.0033}^{+0.0032}$	0.2671	au
${M}_{c}\sin i$	13.2 ± 1.9	13.2	M⊕
ac	0.1058 ± 0.0013	0.1057	au

Note. 80,000 links saved. Reference epoch for γ, $\dot{\gamma }$, $\ddot{\gamma }$: 2455275.946932.

Download table as:  ASCIITypeset image

3.2.2. Planet Search Completeness

We carried out an injection and recovery test on the combined HARPS and HIRES RV data to examine the sensitivity of RV data to potential additional companions in the system. In total, 3000 synthetic planets were injected into the data, with orbital eccentricities drawn from a beta distribution using parameters from Kipping (2013), and orbital periods and minimum masses drawn from log-uniform distributions. We then used the results from the injection and recovery test to compute a search completeness contour, shown in Figure 3 as a function of mass (${M}_{p}\sin i$) and semimajor axis (a). The two recovered signals described in Section 3.2.1 are marked as solid black dots with the ∼14 days signal lying on the 50% completeness curve, indicated by the black line. The individual blue and red dots indicate whether the recovery was successful at each injected location. The results suggest that additional giant planets more massive than the already discovered planet b are clearly recoverable if present in the inner part of the system, but were not detected within our data, suggesting that such an additional planetary presence is extremely unlikely within 1–2 au of the host star. However, our data do not rule out their presence beyond 10 au, which corresponds roughly to the observing baseline of the combined RV data set. On the other hand, smaller planets less massive than 10 Earth masses lie below the 50% completeness curve for the majority of the parameter space, hinting the possibility that small planets may yet be present within the system. These small planets sit below the sensitivity of our RV data and their discoveries require further observations and/or more precise measurements. Overall, our injection-recovery analysis is able to exclude Saturn-mass and Jupiter-mass planets out to 3 au and 8 au, respectively. It is worth noting that the completeness analysis described here refers to ${M}_{p}\sin i$, and so depends on the inclination of the synthetic planetary orbits relative to the plane of the sky.

Zoom In Zoom Out Reset image size

3.2.3. Stellar Activity

As described in Section 2.2, the TESS photometry consists of three sectors: 10, 36, and 63. For the analysis of the transit signal, designated as TOI-6713.01, we used the LightKurve package (Lightkurve Collaboration et al. 2018) to download the relevant sectors of the PDC and SAP time series measurements from MAST. We prepared the data by first using a Savitzky–Golay filter (Savitzky & Golay 1964) to temporarily flatten the lightcurve after masking out the transits (with one transit duration on either side). An outlier mask was created by performing five iterations of 3σ clipping on the flattened lightcurve. We then also removed segments of the lightcurve that were heavily affected by momentum dumps. This, in addition to the outlier mask, removed ∼8% of the lightcurve data.

We fit the prepared lightcurve using the exoplanet package (Foreman-Mackey et al. 2021), which uses a Hamiltonian Monte Carlo (HMC) routine to explore posterior probability distributions. To model the transit, we assumed a circular orbit with orbital period (P), time of inferior conjunction (T_conj), scaled planet radius (R_p /R_⋆), impact parameter (b), transit duration (T₁₄), and mean flux offset (μ) as free parameters. Quadratic limb darkening coefficients were held constant at u₀ = 0.45 and u₁ = 0.18 (Claret 2017). Gaussian priors were set on P and T_conj based on values reported by the Exoplanet Follow-up Observing Program.

The stellar variability was modeled using the RotationTerm kernel included in celerite2 (Foreman-Mackey et al. 2017), which is a mixture of two simple harmonic oscillators. For the secondary oscillation quality factor (Q), the difference in quality factors between the primary and secondary modes (dQ), fractional amplitude of secondary mode relative to primary mode (f), and the standard deviation of the GP (σ), we sampled the logarithmic value of the parameter. A Gaussian prior was placed on the rotation period (P_rot) informed by the stellar activity analysis in Section 3.1. However we note the 18 days signal was not well-preserved in the PDC lightcurve.

We used the values obtained from minimizing a negative log-likelihood function as initial positions for two parallel chains and ran the HMC for 2000 tuning steps and 4000 sampling steps per chain. In Table 3 we provide the median values with their 1σ uncertainties for each fitted and derived parameter. The model fit to binned TESS data is shown in Figure 4. The transit depth is 231 ppm and the out-of-transit PDCSAP photometry has a standard deviation of ∼300 ppm. As there are ∼1300 measurements acquired during transit, the transit signal-to-noise ratio (S/N) for this detection is S/N ∼ 27, according to the methodology described by Pont et al. (2006). Note that we examined both the PDC and SAP photometry from the individual TESS sectors and found that the transit signal was present in all cases, although the signal becomes less significant in the SAP data. More follow-up observations are needed to validate and confirm the planetary nature of the photometric signal.

Zoom In Zoom Out Reset image size

Table 3. TOI-6713.01 TESS Transit Parameters

Parameter (Units)	Credible Interval	Source
Transit Parameters
P (days)	2.1538197 ± 0.0000041	This work, fit
Tconj (BJDTDB–2457000.0)	1571.52211 ± 0.00081	This work, fit
Rp /R⋆ (–)	0.0152 ± 0.0014	This work, fit
T14 (hr)	1.667 ± 0.065	This work, fit
b (–)	0.52 ± 0.29	This work, fit
e (–)	≡0	This work
ω (–)	≡0	This work
Rp (R⊕)	1.30 ± 0.12	This work, derived
a/R⋆ (–)	8.42 ± 0.15	This work, derived
a (au)	0.03054 ± 0.00037	This work, derived
i (deg)	86.5 ± 2.0	This work, derived
Gaussian Process Parameters
$\mathrm{ln}(Q)$ (–)	−6.59 ± 0.53	This work, fit
$\mathrm{ln}({dQ})$ (–)	−1.43 ± 0.75	This work, fit
Prot (days)	8.9 ± 2.6	This work, fit
$\mathrm{ln}(\sigma )$ (–)	−8.30 ± 0.12	This work, fit
f (–)	0.093 ± 0.025	This work, fit

Download table as:  ASCIITypeset image

Given the eccentric orbits of the RV-detected giant planets, and the prospect of a further inner terrestrial planet (see Section 3), the architecture of the system presents an interesting opportunity to explore planetary dynamical interactions within a relatively compact system. The dynamical simulations conducted for this work make use of the Mercury Integrator Package (Chambers 1999), using similar methodology to that conducted by Kane (2019) and Kane et al. (2021c). We also adopt the hybrid symplectic/Bulirsch–Stoer integrator with a Jacobi coordinate system, which generally provides more accurate results for multiplanet systems (Wisdom & Holman 1991; Wisdom 2006). We used a time step of 0.05 days (1.2 hr) to ensure adequate sampling of the inner planetary orbit and ran the simulation for a 10⁷ yr. Although neither of the two RV planets are known to transit, we assumed that the system is roughly coplanar. Adopting the parameters for the inner planet described in Section 3.3, we assumed a circular orbit for that planet. Given the estimated radius for the inner planet of 1.30 ± 0.12 R_⊕, we adopted a planetary mass of 2 M_⊕, based on predictions from exoplanet mass–radius relationships (Weiss & Marcy 2014; Chen & Kipping 2017).

Our simulations found that the system is able to retain long-term stability for the full 10⁷ yr integration. The results of our simulations are represented in Figure 5, showing the eccentricity evolution of the transiting planet candidate, TOI-6713.01 (top panel), the newly detected planet c (middle panel), and the previously known planet b (bottom panel). In order for the structure of the eccentricity evolution to be easily visible, we have restricted the data shown in the plot to the first 10⁵ yr of the full 10⁷ yr simulation. Note also that we use identical eccentricity scales for each panel to emphasize the relative amplitudes of the eccentricity variations. Figure 5 shows that the two giant planets experience a regular exchange of angular momentum that results in eccentricity ranges of 0.12–0.16 and 0.10–0.29 for planets b and c, respectively. The eccentricity range for TOI-6713.01 is 0.00–0.16, with a clear beating pattern, indicating that the presence of the inner planet causes a slight difference in frequency between the angular moment exchange between planets b and c. Since TOI-6713.01 lies so close to the host star, even a relatively small eccentricity excitation may have significant consequences for tidal forces experienced by the planet.

Zoom In Zoom Out Reset image size

Based on the eccentricity evolution resulting from the dynamical simulations from Section 4.1, we calculated the tidal heating of the inner planet candidate and the consequences for the planetary surface temperature. The power contributed to a rigid body via tidal heating is given by

$$\begin{eqnarray}&&{P}_{\mathrm{tides}}=\dot{E}=-\mathrm{Im}({k}_{2})\displaystyle \frac{21}{2}\displaystyle \frac{{G}^{3/2}{M}_{\star }^{5/2}{R}_{p}^{5}{e}^{2}}{{a}^{15/2}}\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, M_⋆ is the mass of the host star, R_p is the radius of the planet, e is the orbital eccentricity, and k₂ is the second-order Love number (Squyres et al. 1983; Meyer & Wisdom 2007). For our calculations, we adopt a Love number of k₂ = 0.35, based on the terrestrial bodies (Zhang 1992). Equation (1) clearly demonstrates the extreme sensitivity of tidal heating and power output to the semimajor axis of the body, since it scales with a^−15/2. For example, the tidal heating of Io results in a total tidal dissipation of 10¹⁴ W (Bierson & Steinbrügge 2021).

The proximity of TOI-6713.01 to the host star will also result in significant incident flux of stellar radiation that contributes to the overall energy budget of the planet. The power absorbed by the planet due to the star is given by

$$\begin{eqnarray}&&{P}_{\mathrm{star}}=\displaystyle \frac{{L}_{\star }}{4\pi {a}^{2}}(1-A)\pi {R}_{p}^{2}\end{eqnarray} \tag{ 2 }$$

where L_⋆ is the stellar luminosity and A is the Bond albedo of the planet. In order to set an upper limit on the stellar contribution to the planetary heating, we adopt a Bond albedo of A = 0.0. For epochs of the simulated dynamical evolution where the planet has a nonzero eccentricity, we consider only the average flux received and so evaluate Equation (2) at the semimajor axis.

To calculate the total power emitted by the planet, we consider the planet using a blackbody approximation such that the power radiated by the planet may be expressed as

$$\begin{eqnarray}&&{P}_{\mathrm{rad}}=4\pi {R}_{p}^{2}\epsilon \sigma {T}_{p}^{4}\end{eqnarray} \tag{ 3 }$$

where is the emissivity of the planet, σ is the Stefan–Boltzmann constant, and T_p is the blackbody equilibrium temperature of the planet. We assume an emissivity of = 1.0 which treats the planet as a perfect emitter. As a blackbody, the power radiated equates to the total power absorbed (P_star + P_tides). The equilibrium temperature of the planet is therefore provided by

$$\begin{eqnarray}&&{T}_{p}={\left(\displaystyle \frac{{P}_{\mathrm{star}}+{P}_{\mathrm{tides}}}{4\pi {R}_{p}^{2}\epsilon \sigma }\right)}^{\tfrac{1}{4}}\end{eqnarray} \tag{ 4 }$$

where the power calculations of Equations (1) and (2) have been incorporated.

Figure 6 shows the results of our tidal effects calculations for TOI-6713.01 for the first 10⁵ yr of the eccentricity evolution (see Section 4.1). The top panel shows an expanded view of the eccentricity evolution of the planet candidate. The middle panel shows the tidal energy dissipation rate (power) emitted by the planet (Equation (1)). The bottom panel shows the equilibrium temperature of the planet, combining the effects of tidal energy and incident stellar flux (Equation (4)). Clearly, the proximity of the planet to the host star results in a highly sensitive response of the planet's tidal energy to the eccentricity variations. The mean tidal power produced by the planet is 8.3 × 10²⁰ W, or almost 7 orders of magnitude larger than Io. Amazingly, using the above mentioned assumptions regarding Love number, Bond albedo, and emissivity, these calculations lead to a very high equilibrium temperature for the planet, where the temperature varies in the range 1202–2646 K. For an emissivity of = 0.5, the equilibrium temperature raises even higher, reaching peak values of 3130 K. In either scenario, the peak value of the temperature evolution is well above that needed for the surface to be in an entirely molten state (Boukaré et al. 2022).

Zoom In Zoom Out Reset image size