DMPP-4: candidate sub-Neptune mass planets orbiting a naked-eye star Open Access

ABSTRACT

We present radial velocity measurements of the very bright (V ∼ 5.7) nearby F star, DMPP-4 (HD 184960). The anomalously low Ca ii H&K emission suggests mass-loss from planets orbiting a low activity host star. Periodic radial velocity variability with ∼10 m s⁻¹ amplitude is found to persist over a >4 yr time-scale. Although the non-simultaneous photometric variability in four TESS sectors supports the view of an inactive star, we identify periodic photometric signals and also find spectroscopic evidence for stellar activity. We used a posterior sampling algorithm that includes the number of Keplerian signals, N_p, as a free parameter to test and compare (1) purely Keplerian models (2) a Keplerian model with linear activity correlation and (3) Keplerian models with Gaussian processes. A preferred model, with one Keplerian and quasi-periodic Gaussian process indicates a planet with a period of |$P_\textrm {b} = 3.4982^{+0.0015}_{-0.0027}$| d and corresponding minimum mass of |$m_\textrm {b}\, \textrm {sin}\, i = 12.2^{+1.8}_{-1.9}$| M_⊕. Without further high-time resolution observations over a longer time-scale, we cannot definitively rule out the purely Keplerian model with two candidates planets with |$P_\textrm {b} = 2.4570^{+0.0026}_{-0.0462}$| d, minimum mass |$m_\textrm {b}\, \textrm {sin}\, i = 8.0^{+1.1}_{-1.5}$| M_⊕ and |$P_\textrm {c} = 5.4196^{+0.6766}_{-0.0030}$| d and corresponding minimum mass of |$m_\textrm {b}\, \textrm {sin}\, i = 12.2^{+1.4}_{-1.6}$| M_⊕. The candidate planets lie in the region below the lower-envelope of the Neptune Desert. Continued mass-loss may originate from the highly irradiated planets or from an as yet undetected body in the system.

1 INTRODUCTION

The dearth of planets, now commonly dubbed the ‘Neptune Desert,’ is a very clear region in the planetary mass versus orbital period plane (Mazeh, Holczer & Faigler 2016). It is delineated by a sharp upper edge with a planet mass inversely proportional to the orbital period, the trend first noted by Mazeh, Zucker & Pont (2005), and a lower boundary with masses that are roughly linearly proportional to the period. Several mechanisms to explain both boundaries were discussed by Mazeh et al. (2016) and more recently by Vissapragada et al. (2022). The upper boundary can be explained as a death line, where inwards migrating planets lose much of their mass due to insolation from the host star. They consequently move down below the lower boundary of the desert where a large number of transiting Kepler planets were found in a ridge in the planet radius versus orbital period diagram (Mazeh et al. 2016). The lack of intermediate-mass planets at short periods exists because the time-scale for this mass-loss is thought to be short. A number of planets have nevertheless been found in the Neptune Desert (West et al. 2019; Armstrong et al. 2020; Burt et al. 2020; Díaz et al. 2020; Dreizler et al. 2020; Jenkins et al. 2020; Jordán et al. 2020; Kanodia et al. 2021; Murgas et al. 2021; Smith et al. 2021; Mori et al. 2022). These planets orbit stars with 11.1 < V < 17.0, except for a single example with V = 9.8 (Jenkins et al. 2020). Identification of planetary systems with brighter stellar hosts that are either in, or have potentially transitioned the Neptune Desert, would be advantageous for follow-up characterization.

Chromospheric emission has been found to be depressed in stars harbouring mass-losing planets (Vidal-Madjar et al. 2003; Haswell et al. 2012). The Ca ii H&K emission in these stars is below the basal level for inactive stars owing to absorbing gas that is lost from close-orbiting planets (Haswell et al. 2020). Hence, a priori identification of stars that harbour hot ablating planets is possible. The Dispersed Matter Planet Project (DMPP) targets bright stars within 100 pc, with typical apparent magnitudes of m_v < 10, making follow-up characterization easier (Haswell et al. 2020). A total of 39 host stars with anomalously low Ca ii H&K emission have been carefully selected by DMPP. Because this gas is concentrated in the orbital plane of a planet (Debrecht et al. 2018), the method preferentially picks out edge-on systems. Consequently, a significantly higher than average fraction of these planets are expected to transit. Most of the DMPP planets reported to date lie just below the lower Neptune desert boundary (Barnes et al. 2020; Haswell et al. 2020; Staab et al. 2020). One explanation is that they are the stripped cores of more massive planets that have crossed the Neptune Desert. It is possible that the DMPP planets identified so far are part of the same population as the dense concentration or ‘ridge’ of Kepler planets in the period–radius diagram (Mazeh et al. 2016). The DMPP planets thus offer valuable opportunities to understand the mechanisms sculpting exoplanet demographics, and thus the evolution of individual hot planets.

The compact multiplanet system, DMPP-1, comprises some of the most irradiated rocky planets yet known (Staab et al. 2020). The DMPP planets are not typically mass-losing hot Jupiters; rather, it is likely that the quenching of Ca ii H&K emission comes from the low-mass, potentially rocky planets, or from additional planets that are below the radial velocity (RV) detection threshold. Jones et al. (2020) identified a candidate transiting signature in TESS observations of DMPP-1 that would correspond to such a low-mass planet with an extended mass-losing atmosphere. This planet was not identified in the RVs reported by Staab et al. (2020). DMPP target stars could thus host analogues and precursors of catastrophically disintegrating planets (Rappaport et al. 2012, 2014; Sanchis-Ojeda et al. 2015).

DMPP has proved very successful, with an effectively 100 per cent success rate (Haswell et al. 2020) on targets with sufficient observations. While most targets are in the Southern hemisphere, the DMPP sample comprises a number of northern targets. Here, we present SOPHIE and HARPS-N RV observations of our brightest target, the late-F star, DMPP-4, which is visible to the naked eye. We present system parameters for DMPP-4 in Section 2. In Section 3, we discuss ground based and TESS photometric observations and search for evidence of the stellar rotation period. Spectroscopic activity indicators are analysed in Section 4 before searching for Keplerian signals in the RVs in Section 5. A summary and further discussion in Section 6 is followed by a brief conclusion in Section 7.

2 SYSTEM PARAMETERS OF DMPP-4

At V ∼ 5.7, DMPP-4 is visible to the naked eye. It is a Northern hemisphere star (J2000 coordinates 19:34:19.8,  +51:14:11.8). Stellar observed and derived parameters are listed in Table 1. The Gaia Data Release 3 parallax is π = 39.323 mas, from which a directly estimated distance (d = 1/π) of 25.43 pc is implied.¹

A range of spectral types between F5 (Roman 1949) and F8V (Eggen 1960) have been estimated.² Table 1 includes estimates of V and also B – V = 0.4344 ± 0.0076 derived from photometric conversions³ documented in Gaia Data Release 3 (Gaia Collaboration et al. 2016, 2023), which is more consistent with an earlier spectral type of F5V (Pecaut, Mamajek & Bubar 2012; Pecaut & Mamajek 2013). Our own estimates of T_eff, M_*, and T_* using our HARPS-N observations and species (Soto & Jenkins 2018) are also more consistent with an earlier spectral-type classification of F5V. We also find an age of 2.15 Gyr using species. With careful consideration of macroturbulence, v_mac, the equatorial rotation velocity, v sin i and v_mac have been obtained following the method described in Murphy et al. (2016). Macroturbulence contributes significantly to the broadening, with v_mac = 7.0 ± 2.0 km s⁻¹. The derived v sin i|$=7.5^{+0.5}_{-2.0}$|  km s⁻¹ is consequently significantly lower than would be measured without including v_mac.

3 PHOTOMETRIC OBSERVATIONS

3.1 LCOGT photometry

Photometric observations of DMPP-4 were collected using the Las Cumbres Observatory global telescope network (LCOGT) robotic telescopes (Brown et al. 2013). A total of 761 exposures of ∼1 s duration each were collected over 46 observing nights. Between 3 and 63 observations were made each night, with 1 h cadence between batches of 3−9 observations. The upper panels of Fig. 1 shows the data as mean subtracted flux in parts per thousand (ppt) collected between 2016 June 17 and August 28. The most significant peak in the corresponding periodogram (Fig. 2, upper panel) is found at 5.119 d, with an additional short period peak at 1.979 d and a longer 10.95 d peak. The peak at 1.245 d is the 1-d alias of the 5.119 d peak. The peak at 10.95 d corresponds to a sinusoidal signal with amplitude 5640 ppm. When phase-folded on this periodicity, the data are found to be clustered into 11 distinct groups. The 5.119 d peak corresponds to a sinusoidal signal with semi-amplitude of 7.56 ppt (7560 ppm). We modelled a B-band light curve using the image code DoTS (Collier Cameron 2001) for a single spot with T_phot−T_spot = 2000 K (Berdyugina 2005), finding that a r_spot = 5.04° is required to reproduce this semi-amplitude. A small spot group with the same effective area as a single spot could also potentially yield the same photometric amplitude, though a more widely distributed collection of spots would likely need to cover a larger total effective area to yield the same sinusoidal or near sinusoidal photometric amplitude. The implied effective spot area appears to be at odds with the low log (⁠|$R^{\prime }_{\rm HK}$|⁠) activity of DMPP-4, especially given that the largest sunspots are never this large when the solar maximum log (⁠|$R^{\prime }_{\rm HK}$|⁠) is still considerably greater than is seen on DMPP-4. Moreover, we find that a 5.04° spot would induce a K = 82 m s⁻¹ RV signal, which is more than an order of magnitude greater than we find in our RV data (Section 4). We therefore conclude that the periodicities recovered from the LCOGT B-band data are not genuine activity-related signals.

3.2 TESS photometry

The Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) monitored DMPP-4 in Sectors 14, 15, 40, and 41. The Mikulski Archive for Space Telescopes⁴ provides detrended Pre-search Data Conditioning Simple Aperture Photometry (PDCSAP) light curves. Release notes for each TESS sector⁵ give the Combined Differential Photometric Precision (CDPP) for 1 h time-scales on ∼20 000 targets. The TESS magnitude for DMPP-4 is 5.22 (Table 1), indicating an expected CDPP of ∼15 ppm. This scales to ∼10 ppm for the 2.4 h (0.1 d) bins we used, and is in close agreement with the formal errors on our binned data of 11.4−11.6 ppm.

Photometric semi-amplitude variability on a scale of <100 ppm (Fig. 1) can be seen at times in the TESS observations; considerably less than the 7560 ppm seen in the LCOGT data. However, periodic variability is only easily discerned during part of Sectors 15 and 41. Log-likelihood periodogram searches using sinusoidal signals, are shown in Fig. 2 for Sectors 14 and 15, Sector 40 and 41 and all sectors combined. Highly significant peaks appear in the periodograms with the most significant periods at, respectively, 6.33, 6.47, and 6.40 d. The respective half width at half-maxima (HWHMs) of these peaks of 0.33, 0.44, and 0.33 d give an indication of the spread of power, but should not be treated as period uncertainties, as discussed by VanderPlas (2018). We also note that the longer period at 11.6 d (HWHM  = 1.4 d) is also significant in the Sectors 40 and 41 periodogram, along with a shorter 4.770 d peak. As with the 6.40 d strongest peak when using all sectors, the longer peak at 12.1 d is strongly aliased due to sampling (VanderPlas 2018). The envelope of the aliased peak possesses HWHM  = 1.7 d. The 12.1 d peak is marginally significant with ΔlogL = 13.0 (FAP =0.003).

To investigate the periodicities further, we performed a wavelet analysis using the scaleogram package,⁶ which is based on the PyWavelets library (Lee et al. 2019). The resulting continuous wavelet transform in Fig. 3 confirms the visual inspection of periodicities in Fig. 1. The periodicities identified in the log likelihood periodograms are not present throughout the extent of the time series and the periodicities are often localized. The global spectrum reveals peak power in Sectors 14 and 15 at 5.8 d and in Sectors 40 and 41 at 6.0 d and 10.7 d (peaking at 11.1 d for the maximum power).

We again used DoTS to model the photometric variability expected from TESS observations and the absorption-line profiles for a cool spot on DMPP-4 observed with 550 nm central wavelength. The TESS flux half-amplitude of ∼29 ppm found by phasing Sectors 15 and 41 is commensurate with an equatorial cool starspot of radius 0.34°. A spot of this size is expected to induce RV variability of |$K \langle 16\,{\rm cms}^{-1}$| at 550 nm and is thus below the level that could easily be detected by an RV instrument achieving ∼1 m s⁻¹ precision. If a Solar-like facular/spot area ratio of 11.6 is present on DMPP-4 (Shapiro et al. 2014), a 0.082° region of plage surrounding a 0.024° spot could induce a flux half-amplitude of ∼29 ppm. In this case, the expected RV amplitude is K ∼ 38 cm s⁻¹ at 550 nm. The RV amplitude rises by an order of magnitude for a cool spot with radius 0.2° and surrounding plage of 0.68°. Although the TESS observations and our spot estimates provide further evidence in support of the low activity of DMPP-4, we cannot rule out the possibility of higher activity levels at other epochs when our RV observations were made.

3.3 Spot evolution and the orbital period of DMPP-4

A study of sunspot group data spanning over a century by Berdyugina (2005) found that two active longitudes persist and are separated by ∼180°. One region usually dominates, periodically alternating between the two. More recent work by Basri & Nguyen (2018) has looked at the phenomenon of ‘double-dipping,’ where a stellar photometric light curve spends some of time exhibiting a single sinusoid behaviour at the rotation period and some time in a double-dip mode with two dips or sinusoids within one rotation. The phenomenon was found in the analysis of over 34 000 Kepler light curves by McQuillan, Mazeh & Aigrain (2014), where periodogram peaks at the rotation period and half the rotation period were identified. While two distinct active groups may be present, any spot distribution will typically yield a photometric light curve with either one to two dips. For more active stars, spots may thus be distributed over a range of longitudes rather than in two distinct longitude regions. The exact light-curve morphology is determined by either the relative strength, the distribution of the spots, or by both effects. Basri & Nguyen (2018) used the Kepler light curves to examine how much time stars spend in single-dip mode and how much in double-dip mode. For their hottest sample corresponding to stars with T_eff = 6000−6200K, the time spent in double-dip mode is 1.8 times that spent in single-dip mode. This bimodal behaviour is seen in the TESS observations of DMPP-4. At the TESS Sectors 14 and 15 epochs, two spot groups of similar intensity could yield the periodic signatures seen at ∼5.5−6.5 d, while in Sectors 40 and 41, there appears to be evidence for one spot group for some of the time and two spot groups for some of the time. Evolution of spot groups appears to be relatively rapid and on the time-scale of a single rotation, resulting in a rapidly changing light curve. The light curve in Sectors 40 and 41 can be explained by a single weak/small spot that grows in strength with the simultaneous appearance of a second spot group resulting in a switch from single-dip to double-dip mode.

In summary, we believe that DMPP-4 exhibits small spots, as evidenced by the low log(⁠|$R^{\prime }_{\rm HK}$|⁠), photometric amplitude and spot size estimate. It possibly follows the solar paradigm of two persistent and distinct active regions rather than that of potentially more distributed activity on an active star. However, interpreting periodicities at P/2 from this scenario may be simplistic since other effects such as differential rotation and spot evolution can result in both double dipping and multiple periodicities. A rotation period of ∼11 d rather than ∼6 d thus seems more likely given the estimated stellar rotation velocity and radius.

3.4 Period and inclination distributions

We expect DMPP to find systems with high orbital inclinations where the stellar Ca ii H&K emission would be more readily obscured (Haswell et al. 2012, 2020). So, under the assumption that the stellar spin axis and planetary orbital planes are aligned, we expect a stellar rotation axis which is highly inclined to our line of sight. The upper plot in figure  shows the distribution of expected periods from the R_* = 1.38 ± 0.01 R_⊙ (Gaussian distribution) and v sini i  = |$7.5^{+0.5}_{-2.5}$| m s⁻¹ (asymmetric Gaussian distribution) estimates listed in Table 1. A uniform distribution of axial inclinations of 0° < i < 90° was simulated. We find a resulting modal period of P  = 8.89 d with a range between 2.74 d and 11.10 d at the 16 and 84 per cent intervals. A simulation with 60° < i < 90° only affects the most likely period slightly (P  = 9.02 d), but skews the posterior distribution to higher values of P (8.69 d to 13.50 d). The lower panel in figure shows the axial inclination distributions, i, assuming period distributions with P = 6.0 d and 11.0 d from the TESS periodicities identified in Section 3.2 (assuming Gaussian distributions of σ = ±1 d and ±2 d). For the P = 6.0 d case, a low axial inclination with modal value, i = 32.8° (23.2° to 43.8° at 16 per cent and 84 per cent intervals) is found. For P = 11.0 d, inclinations >49.7° (16 per cent) are found with a modal sin i > 1. The large uncertainty in our v sin i estimate, particularly in the lower limit is a likely explanation for the difference between the 11 d photometric rotation period and Monte Carlo modal period of 8.89 d. If i = 90°, v sin i must be ≤6.35 km s⁻¹ when P = 11 d.

4 ACTIVITY INDICATORS

Spectroscopic observations of DMPP-4 were made with SOPHIE (Perruchot et al. 2008) at the 1.9 m Observatoire de Haute Provence Telescope during four semesters in 2015 and 2016. Further data were obtained with HARPS-N (Cosentino et al. 2014) at the Telescopio Nazionale Galileo telescope. Observations were made in 2015 April–May on three nights over a time span of four nights; and again later in 2015 Nov on five nights. Similarly, observing runs in 2016 May and 2016 November–December each resulted in four nights of observations over a five night time span. We refer, respectively, to the individual SOPHIE runs as SO-15A, SO-15B, SO-16A, and SO-16B. In 2016, HARPS-N observations on four nights (HN-16A) were made in 2016 July between the SO-16A and SO-16B runs. The last set of HARPS-N observations were made in 2019 Aug on a single night (HN-19A). Exposures were monitored to ensure relatively stable counts with SOPHIE. Because DMPP-4 is very bright, observation times of 200–900 s were used, depending on observing conditions. When observations were cycled with other targets, exposures on DMPP-4 were made for at least ∼900 s in a single block to ensure that possible 5 min stellar oscillations were minimized.

The RVs were derived using the matched-template code, harps-terra (Anglada-Escudé & Butler 2012). We also used the activity indices provided by the standard data reduction software pipelines for both instruments. A total of 111 observations were made with SOPHIE, while 87 observations were made with HARPS-N. After optimally weighted binning into 900 s observations, the respective data sets comprise 45 and 27 observations. The binned data are tabulated in Tables A1 and A2.

The low altitude of the Observatoire de Haute Provence, small telescope aperture and instrumental effects such as reference/target fibre drift, potential chromatic κ-correlation (see below) and charge transfer inefficiency effects limit the precision that is achieved with SOPHIE. The formal uncertainties for our 2015 and 2016 observations, before 900 s binning, are 2.25 and 2.26 m s⁻¹. The binned data uncertainties are, respectively, 1.69  and 1.60 m s⁻¹. A measurement precision of ∼1−2 m s⁻¹ has been demonstrated with SOPHIE on RV standard stars on time-scales of a few tens of days (Bouchy et al. 2013). For HARPS-N, the formal uncertainties for single 200 s observations is 1.17 and 0.73 m s⁻¹ for 900 s binned data.

4.1 Kappa correlation with RVs

We investigated correlations between possible diurnal chromatic systematics and the RVs. Bourrier, Lecavelier des Etangs & Vidal-Madjar (2014) and Berdinas et al. (2016) identified a significant intra-night systematic effect on observations made with the HARPS-N spectrograph. The primary cause appears to be an incomplete correction of differential atmospheric refraction by the atmospheric dispersion corrector (ADC), causing colour-dependent (chromatic) flux-losses. The effect is thus potentially important for targets observed throughout the night at different airmasses. We measured the ‘chromatic index,’ κ, defined as the slope of a linear portion of the pseudo spectral energy distribution about a specified wavelength. In practice, we measured the slope in mean flux per échelle order across several orders in a similar manner to Berdinas et al. (2016); see Haswell et al. (2020) for full details. Fig. 5 (left-hand panel) shows the chromatic index, κ, plotted against the RVs for all SOPHIE and HARPS-N observations. The degree of correlation is measured via Pearson’s r coefficient and listed in Table 2; only weak or very weak correlations are seen during most observing runs. There is a moderate correlation in the HN-16A data set, but the significance is also moderate, as indicated by the student’s t test p-value. This suggests that there is no clear evidence of a correlation. The correlation is driven by the two outlying RVs (greatest positive values); removing them yields r = 0.05 and p = 0.86.

4.2 RV correlations with activity measures

The correlation of RVs with the full width at half-maximum (FWHM), line Bisector Inverse Span (BIS) and Ca ii H&K S-index are also plotted in Fig. 5. Correlation measurements are listed in Table 2. We include the HN-19A data in Table 2 for completeness; since this run comprised a single night, activity correlations are likely to be less meaningful.

The BIS shows significant moderate anticorrelation in the SO-16A and SO-16B data sets, but is less pronounced in SO-15B. We note that the SO-15A data is of significantly lower precision than the subsequent observations. For BIS in SO-15A, removing the data point with largest positive RV = 8.49 m s⁻¹ results in r = −0.64 and p = 0.17. For individual observing runs, the BIS anticorrelation is largest and most significant in the SO-16A and SO-16B epochs. In SO-16B, the RVs also show a significant (p = 0.02) moderate anticorrelation with FWHM.

HARPS-N shows moderate correlations, but the significance is low. The HN-16A run shows a moderate S-index variability, although the probability of no correlation is still nearly 10 per cent (i.e. p = 0.097). Of the three activity indicators, only BIS shows a consistently moderate-good correlation.

4.3 Periodicities in activity measures

The right-hand panels in Fig. 5 show periodograms for the κ-correlation and each activity index. Aliases and integer fractions of a day in the Kappa and FWHM periodograms are seen. The κ-correlation and FWHM show no other significant peaks. For BIS, we find a 2.54 d peak with FAP  = 2.7 × 10⁻⁶ and Δ LogL = 27.2 relative to the scenario with no signal. A neighbouring 2.68 d peak with logL_2.54−2.68 < 5 (Δ LogL = 22.9) relative to the main peak is also present. This periodicity, taken with the moderate anticorrelation between BIS and RV, suggests that the RVs may be affected by stellar activity. It is nevertheless difficult to reconcile this finding with our analysis that reveals very low photometric variability and the low log(⁠|$R^{\prime }_{\rm HK}$|⁠), suggesting that DMPP-4is an inactive star. The variable degree of correlation between BIS and RV at different observing epochs suggests that the star shows variable activity. This should however be treated with some caution since each observing epoch comprises only a few nights of observations and so may not sample a complete stellar rotation.

Further, the most significant BIS periodicity is likely to be a fraction of the true period. Period analysis of apparent RV shifts due to line-shape changes as described by BIS, are a proxy measure of the third central moment or line skewness. We find that for stars with realistic spot and facular distributions and solar activity levels, RV periodicities may appear with dominant power at harmonics of the rotation period, including P_rot/2 and P_rot/4. This behaviour was first noted by Boisse et al. (2011). For instance, a single high latitude spot can easily induce RV and BIS variability at P_rot/2, while a low latitude spot could simultaneously induce RV variability at P_rot/2 and BIS variability at P_rot/4. A pair of low-latitude spots located at active longitudes 180° apart, as has been observed for the sun (Berdyugina & Usoskin 2003), could then also be expected to lead to predominant periodicities in RV and BIS at P_rot/4. The exact periodicities are likely to vary from star to star, but depend on the exact spot patterns, presence of absence of faculae and the combination of spot latitude and stellar axial inclination and |$v\, \textrm {sin}\, i$|⁠.

The BIS periodicities of 2.54 and 2.68 d may thus be manifestations of active regions that have true periodicities at twice or four times these values. In other words, P_rot inferred from BIS may be either 5.08 d and 5.36 d (2×) or 10.2 d and 10.7 d (4×). The longer periods are close to photometric periodicities identified from the TESS Sector 40 and 41 observations in Section 3.2. There is tentative evidence for low significance peaks (1– 10 per cent FAP) at ≥10 d in the S-index periodogram, which shows a most significant peak at 11.1 d (Δ LogL = 15.9), again in close agreement with the TESS observations.

5 RV ANALYSIS

We began by assessing the coherence of periodic RV signals over time as more data points were added. Fig. 6 shows the stacked Bayesian general Lomb–Scargle periodograms of the RVs and the residual RVs after accounting for the first Keplerian (Mortier et al. 2015; Mortier & Collier Cameron 2017). Increasing the number of observations over the six observing campaigns, spanning 2015 and 2016 with the additional night in 2019 steadily increases the significance of each signal. There is some aliasing of peaks, which split into sub-peaks with addition of the final data points taken at HARPS-N in 2019. The long-term signal coherence suggests that the periodicities are Keplerian in origin and can be attributed to planets. Nevertheless, since the periodicities found in the BIS data are similar to those of the second RV signal, we investigate a number of scenarios below.

5.1 Maximum likelihood searches

We carried out likelihood period analyses on the combined SOPHIE and HARPS-N RVs, assuming circular Keplerian orbits with e = 0. The Keplerian model includes offset terms for each instrument. In addition to a period search on the RVs only (RV model), we also performed a period search that incorporated the BIS values into the likelihood model via a linear correlation coefficient (RV + linBIScorr model). Keplerians are added recursively into the model when searching for additional signals. The likelihood model is described in detail in Anglada-Escudé et al. (2013) and Anglada-Escudé et al. (2016) and includes white noise terms to account for noise in addition to the formal uncertainties and systematic RV offset terms for each data set.

Fig. 7 (top left) shows that periodicity at P = 3.498 d is found for the 1-Keplerian models (N_p = 1) for both RV and RV + linBIScorr models. Table 3 lists the likelihood values, logL, and change in likelihood, ΔlogL, of recursively added periodicities. The corresponding estimates of the Bayesian Information Criteria, ΔBIC = BIC_k − BIC_{k − 1} were used to estimate Bayes Factor values, |$\textrm {BF} = e^{(\textrm {BIC}_{k-1}-\textrm {BIC}_{k})}/2$|⁠. The ratio of Bayesian evidence of competing models, the Bayes Factor (BF), provides a measure of support in favour of one model over another and is obtained directly from the samples obtained. We use the scale defined by Trotta (2008) and discussed in Standing et al. (2022), where 3 ≤ BF < 12 implies weak evidence, 12 ≤ BF < 150 implies moderate evidence and BF ≥ 150 indicates strong evidence. We note that these BF values are approximations and should not be directly compared with the values reported from posterior sampling in the following sections. For N_p = 1, ΔBIC for the RV + linBIScorr model is smaller than for the RV model. Nevertheless, both ΔBIC and BF indicate very highly significant periodicity at P = 3.498 d for both models.

For the 2-Keplerian solutions, with N_p = 2, significant peaks at 2.459 d are found (Table 3). For the RV model, Δ LogL = 17.3 (ΔBIC = −21.9). While this periodicity also contains the highest power in the RV + linBIScorr model, the significance is decreased to ΔLogL = 14.6. The corresponding ΔBIC = −14.3 and BF = 1274 still indicates very strong evidence for this model (Raftery 1995), though the RV model is preferred over the RV + linBIScorr model. Although the linBIScorr correlation coefficient for HARPS-N is −1.1 ± 0.5, the coefficient for the SOPHIE data is consistent with zero at −0.2 ± 0.7. This suggests that the linear correlation coefficient is not adequate for modelling long baseline time series.

The changing correlation between the RVs and simultaneous BIS values at different epochs, as evidenced by Pearson’s r statistic in Table 2, might be indicative of an active star with continually evolving activity. Since DMPP-4 shows evidence of low level photometric modulation and evidence for BIS variability at a period that matches the signal at ∼2.5 d, we investigated whether Gaussian process (GP) modelling can be used to characterize the effects of activity on the RVs.

5.2 Signal recovery and model comparison using posterior sampling

We used kima (Faria et al. 2018) to search for Keplerian signals and investigate various model scenarios. We include the maximum likelihood values in our model solutions in Table 4, but note that they should not be used to distinguish models since they do not take into consideration the number of model parameters or data points. Kima uses Diffusive Nested Sampling (Brewer, Pártay & Csányi 2011) to sample from the joint posterior distribution and enables the marginal likelihood (Bayesian evidence) of a model to be directly estimated (in contrast to the BF estimates derived from the BIC values in Section 5.1). This has the benefit of correctly enabling intercomparison of models (Feroz, Balan & Hobson 2011; Brewer & Donovan 2015). With kima, it is possible to either fix the number of Keplerians, N_p, or include N_p as a free parameter. The likelihood model includes the systemic velocity, γ, of a reference data set (in our case, the HARPS-N RV data, γ_RV,HN) and relative offsets for other data sets. Since our RV data are derived via matched spectrum template derived from the observations, these parameters are expected to be close to |$0\, \textrm {ms}^{-1}$|⁠. As with our likelihood periodogram searches, additive white noise, σ, are included as free parameters for each data set. Further details can be found in Standing et al. (2022).

5.2.1 1-Keplerian (N_p = 1) solution

Using kima to obtain planetary parameters as described in Standing et al. (2022), we find |$P_\textrm {b}= 3.49733^{+0.00093}_{-0.23513}$| when we restrict the search to 1 Keplerian (N_p = 1), in agreement with the likelihood period search. The skewed uncertainties (i.e. a relatively large negative uncertainty) arise from posterior sampling of the alias peaks, which can be seen at P = 3.19 and P = 3.26 d in Fig. 7 (a). The peak aliasing is a consequence of the sampling, since observations span ≤5 nights at any given observing epoch with SOPHIE or HARPS-N. The N_p = 1 solution parameters are shown in Table 4, column 1; the RVs and solution curve are shown in Fig. 8. We find BF > 7676 in favour of a single-planet model compared with a zero-planet model. [N.B.: Since no samples were obtained during the run for the N_p = 0 model, we estimate a lower limit on the BF by setting the number of samples with N_p = 0 equal to 1 (Standing et al. 2022; Triaud et al. 2022; Baycroft et al. 2023; Standing et al. 2023).]

5.2.2 2-Keplerian (N_p = 2) solution

We allowed kima to sample posterior space with N_p as a free parameter. In this way, it is possible to recover the preferred number of Keplerians. We place a maximum N_p = 3, and find the solution with N_p = 2 has the highest posterior evidence, with BF = 2397, compared with the N_p = 1 model, and is shown in Fig. 9. There is not sufficient evidence for further Keplerians. However, the two preferred periods differ from the those found by the likelihood model, which adds Keplerian signals recursively. Instead of the P_b = 3.498 d and P_c = 2.459 d identified in Section 5.1, the N_p = 2 model recovers |$P_\textrm {b} = 2.45703^{+0.00293}_{-0.00082}$| and |$P_\textrm {c} = 5.4193^{+0.0060}_{-0.0032}$| d (Table 4, column 2). With the additive white noise terms, σ_RV,SO and σ_RV,HN, the effective uncertainties in the N_p = 2 model are 1.68 and |$1.36\, \textrm {ms}^{-1}$|⁠, suggesting that the formal uncertainties, particularly for HARPS-N (see Section 4), may be underestimated.

We investigated the difference between the periods identified by the N_p = 2 recursive likelihood method in Section 5.1 and the posterior sampling method of kima. We performed MCMC posterior sampling on the two periods we identified in Section 5.1 and, then performed the same analysis on the periods identified by kima. Moderate evidence in favour of the two periodicities recovered by kima was found over the recursively added second periodicity, with ΔBIC + = −5.04 and BF = 12.2. Similarly, we tried restricting the first signal to the N_p = 1 period identified by kima and used kima to detect any further signals. The posterior sampling yields a second period at 2.23 d with a moderate BF = 29.2, somewhat lower than for the preferred N_p = 2 solution.

The discrepancy between the periods recovered by the two approaches arises because of the data sampling, which is not optimal for recovering multiple periodicities with short observing epochs spanning only a few days. In this scenario, the multimodal posterior sampling used by kima may perform better than a method that adopts recursive addition of signals.

5.3 Quasi-periodic GP

kima has also been adapted to enable the use of a GP. The hyperparameters of the GP are inferred together with the orbital parameters (Faria et al. 2023). The GP quasi-periodic kernel is defined in the form given by Rasmussen & Williams (2006, see also Haywood et al. 2014) as

which describes the correlation between RVs at times t_i and t_j. Four hyperparameters define the covariance matrix. The amplitude, η₁, describes the deviation of the GP models from the mean function. The parameter, η₂ is the long-term evolution time-scale, while η₃, the characteristic period of the GP, can be related directly to the stellar rotation period. As noted by Rajpaul et al. (2015) and Barragán et al. (2022), it is desirable that η₃ < η₂, to ensure the validity of using the quasi-periodic GP with a periodic component. The factor, η₄, is proportional to the inverse of the square root of the harmonic complexity, Γ (e.g. see Rajpaul et al. 2015, Barragán et al. 2022 and Nicholson & Aigrain 2022). From equation (1), |$\Gamma = 2/\eta _4^2$|⁠, and describes the roughness or complexity of variations within each characteristic period, η₃.

5.3.1 Keplerian  + quasi-periodic GP

Because our data comprise only a few nights per observing run, a careful choice of priors for the GP hyperparameters, is required. We tested solutions informed by the photometry and activity signatures investigated in Sections 3 and 4. Following our finding of photometric periodicity (Section 3), we first trialled a GP prior for η₃ of |$\mathcal {G}[11.6\pm 1.5]$| d (i.e. a Gaussian distribution). Given the likely activity signatures at P_rot/2 and P_rot/4 discussed above in Section 4.3, we expect characteristic activity signatures on timescales shorter than P_rot. This suggests some degree of harmonic complexity is required to account for any structure related to the characteristic period, η₃. Nicholson & Aigrain (2022) note, for instance, that analysis of stellar light curves typically yield harmonic complexity of 0.5 < Γ < 2. To investigate the posterior distribution of the harmonic complexity, we used a fairly wide prior of 0.1 < η₄ < 10, corresponding to harmonic complexity of 0.02 < Γ < 200.

With this configuration, we do not find evidence for any planetary signals. Although the posterior distribution returns significant periodicities at ∼2.5, 3.5, and 5.4 d, the posterior samples for N_p = 1, 2, and 3 return low respective Bayes Factors of BF  = 1.2, 1.3, and 1.4. We found that the prior distributions are recovered in the posteriors for η₂ and η₃. This not surprising since the few-day span of each observing epoch is much less than the observationally informed priors. The posterior distribution for |$\eta _4 = 0.15^{+0.17}_{-0.04}$| corresponds to harmonic complexity of |$\Gamma = 89^{+76}_{-69}$|⁠. This high harmonic complexity thus accounts for most of the RV variability. However, the sampling rate of the observations leads us to reject solutions that fit RV variability on very short time-scales. There are too many fitting parameters and not enough data points to constrain a credible solution. Essentially, the model is too flexible.

5.3.2 A simplified quasi-periodic GP

The intensive sampling strategy we adopted to enable searches for shorter period planets means that several observations were made over a few hours on each night. By contrast, the significant periodicities in the data are typically a few days, with the effect that the 71 data points essentially behave as 21 high S/N data points. Since the N_p = 1  + GP model in Section 5.3.1 contains 13 fitting parameters, the flexibility of the GP model and thus the degeneracy with planetary signals is high (Nicholson & Aigrain 2022). Thus, it is not surprising that the hyperparameter priors either dominate the posteriors in Section 5.3.1, or lead to overfitting when relatively uninformative and wide.

In simulations, Nicholson & Aigrain (2022) noted that there is no single starspot parameter in their models that relates directly to Γ. They adopted initial guesses of Γ = 1 when fitting simulated light curves, as this value represents the transition between low and high harmonic complexity (Barragán et al. 2022). Informed by the simultaneously derived BIS periodicity of 2.54 d, we obtained posterior samples with a prior on η₃ of |$\mathcal {G}[2.5\pm 0.5]$| d. In light of the posterior distribution for η₄ that we found in Section 5.3.1, we followed Nicholson & Aigrain (2022), but fixed Γ = 1 (i.e. we fixed |$\eta _4 = \sqrt{2}$|⁠). With these priors, posterior sampling leads to a preference for N_p = 1, with |$P = 3.4982^{+0.0015}_{-0.0027}$| d and evidence of BF = 111.4. As in Section 5.3.1, the posterior distribution is dominated by the priors for η₂, though the data appears to have had slight influence on the posterior GP periodic component |$\eta _3 = 2.42^{+0.25}_{-0.24}$|⁠. The RV solution curves and folded RVS are plotted in Fig. 10 (a, b). Hence, there is moderate evidence for a single Keplerian with the same period as the purely Keplerian N_p = 1 model. The posteriors show that inclusion of the GP resulted in reduced significance of the periodicities at 2.457 and 5.419 d found for the purely Keplerian N_p = 2 model.

5.4 Alternative models using activity indicators

With the evidence for moderate activity correlations between the RV and BIS data, we used kima to obtain simultaneous posterior samples using both data sets and a simultaneous single GP model. Here, an additional GP hyperparameter, η_1,2 is required to model the BIS amplitude. The number of fitting parameters for the model is 17. This is similar to the approach adopted by Barragán et al. (2022), although we do not consider the time derivative of the GP amplitude terms here. Further, adding the BIS time series doubles the size of the data set to 142, or if we consider the effective data set size, to 42 as in Section 5.3.2. The flexibility of the model is thus reduced with inclusion of the activity time series.

5.4.1 Model with simultaneous BIS activity indicator

We used the same model priors as in Section 5.3.1. We again found that the η₂ and η₃ posterior distributions matched the priors while a high harmonic complexity is preferred, with |$\eta _4 = 0.15^{+0.05}_{-0.04}$| (⁠|$\Gamma = 89^{+76}_{-39}$|⁠). The posterior samples demonstrate an overdensity at P ∼ 3.5 d, though for N_p = 1, 2, and 3, respective Bayes Factors of BF  = 1.0, 0.8, and 0.8 indicate very weak evidence for any Keplerian signals. Even with the inclusion of the BIS time series, the model is too flexible.

5.4.2 Model with simultaneous BIS activity and fixed harmonic complexity

Finally, with fixed Γ = 1 (⁠|$\eta _4 = \sqrt{2}$|⁠), as in Section 5.3.2, we used a prior on η₃ of |$\mathcal {G}[2.5\pm 0.5]$| d. The solution is shown in Table 4, column 4 and in Fig. 10 (c, d, and e). A BF = 1038 indicates strong evidence in favour of a N_p = 1 solution with with P = 3.4982 d. The inclusion of BIS thus appears to have provided further evidence for the 2.5 d activity periodicity, yielding a higher BF compared with the same model and priors on the RV-only time series in Section 5.3.2.

Fig. 10 (e) demonstrates that there is considerable scatter in the BIS at some epochs, particularly in the poorer quality data in the SO-15A run. The white noise terms for the BIS data (Table 4) reflect this. Modelling the RV and BIS data simultaneously has also likely resulted in the higher RV white noise terms compared with the N_p = 2 solution; accounting for σ_RV,SO and σ_RV,HN yields respective effective uncertainties of |$2.22$| and |$1.48\, \rm {m\,s}^{-1}$|⁠.

6 SUMMARY AND DISCUSSION

Although there is clear evidence for a single candidate Keplerian (N_p = 1) with |$P_\textrm {b}~= 3.49791^{+0.00046}_{-0.00143}$| d and m_b sin i  = |$11.6_{-1.5 }^{+1.6 }$| M_⊕, the purely Keplerian model indicates strong evidence favouring a 2-Keplerian (N_p = 2) solution, with a relative Bayes Factor, BF = 2397. The preferred N_p = 2 solution does not include the N_p = 1 period, but rather recovers periods of |$P_\textrm {b}~=~2.4570^{+0.0026}_{-0.0462}$| and |$P_\textrm {c}~=~5.4196^{+0.6766}_{-0.0030}$| d with respective derived minimum masses of |$m_\textrm {b}\, \textrm {sin}\, i~= 8.0_{-1.5 }^{+1.1}$| M_⊕ and |$m_\textrm {c}\, \textrm {sin}\, i~= 12.2_{-1.6}^{+1.4 }$| M_⊕. The N_p = 2 solution derived from likelihood periodogram searches, with P_b =  3.498 d and P_c =  2.459 d, suggests that there may be conditions where the common approach of recursively adding signals does not lead to the best solution; here, it is likely that the length of the candidate periods and relative length of the observing runs has resulted in ambiguity. Our finding may also in part be a limitation of the period search algorithm, which does not explore maximum likelihood space for a global maximum, but rather, iteratively adjusts the Keplerians already identified when searching for additional signals. In this respect, an algorithm capable of finding a global marginalized maximum likelihood in a potentially multimodal posterior space is to be preferred, particularly when considering multiple Keplerian signals. It is important to note that the posterior samples from kima also contained a significant number of samples around the 3.498 d period (around 40 per cent compared with the 2.4570 d period and 50 per cent compared with the 5.4196 d period).

Despite the low log (⁠|$R^{\prime }_{\rm HK}) = -5.24$| that we used to select DMPP-4 for follow-up RV observations, there is some evidence for activity-induced RV variability. This is consistent with the intrinsic stellar activity corresponding to a higher value of log(⁠|$R^{\prime }_{\rm HK}$|⁠) than observed. The hypothesis underlying DMPP is that circumstellar absorption depresses the log (⁠|$R^{\prime }_{\rm HK}$|⁠), so this is as expected. Because the N_p = 2 solution recovers periods that are close to one half and one quarter integer factors of the inferred rotation period (P_rot/2 and P_rot/4), as discussed in detail in Section 4.3, we were prompted to investigate RV models that contained an activity component.

The GP models we investigated for DMPP-4 were necessarily restrictive because of the nature of the data sampling, the relatively short observing runs and the effective size of the data sets. When a quasi-periodic GP is included in the model, we find moderate evidence (BF = 111.4) for a single planet in the RVs with |$P_\textrm {b}~= 3.4982^{+0.0015}_{-0.0027}$| d and |$m_\textrm {b}\, \textrm {sin}\, i~= 12.6^{+1.6}_{-1.8} \textrm {M}_{\oplus }$|⁠. Including the simultaneous BIS measurements augments the Bayesian evidence, yielding |$P_\textrm {b}~= 3.4982^{+0.0022}_{-0.2327}$| d and |$m_\textrm {b}\, \textrm {sin}\, i~= 12.2^{+1.8}_{-1.9} \textrm {M}_{\oplus }$| with a BF = 1038.

Thus, although the purely Keplerian model prefers N_p = 2, the model restricted to N_p = 1 and the models with a GP and no N_p restriction all provide moderate to strong evidence for a single planet. The periods and minimum masses are consistent within the uncertainties for these models. For the Keplerian models with a GP, with and without simultaneous BIS data, either greater model flexibility, the adopted tighter priors, or both, result in lower overall Bayesian evidence. Without further data it is thus difficult to distinguish between the preferred N_p = 2 model and models that return single Keplerians in the presence of variable activity. Nevertheless, there is (i) evidence for activity periodicities that coincide with the identified RV periodicities in the purely Keplerian N_p = 2 model. Further, there is (ii) a not-insignificant appearance of the 3.498 d periodicity in the model posterior samples for the N_p = 2 model. And finally, (iii) a preferred periodicity at 3.4982 d is found in the Keplerian  + GP models. These three observations might be taken as reasonable evidence to prefer the solutions with a single 3.498 d Keplerian.

We based our choices of GP hyperparameter priors on starspot distribution hypotheses and photometrically inferred periodicities. Recently, however, Nicholson & Aigrain (2022) have suggested that harmonic complexity priors derived from photometry are not necessarily appropriate. Simplifying the model by removing the harmonic complexity, Γ, and relying only on simultaneously derived and modelled BIS periodicities makes intuitive sense, especially in light of the relatively few available effective data points. Nicholson & Aigrain (2022) also investigated the effect of reducing the number of data points and hence increasing the flexibility of the model. The resulting increased degeneracy highlights the need for sufficient data points.

The following sections further investigate the nature of the planet candidates. Despite a preference for the models with a single planet, we also considered the N_p = 2 solution.

6.1 Orbital stability of planets in the N_p = 2 model

We checked the stability of the orbital solution for the N_p = 2 model using the IAS15 integrator within the N-body orbital integrator, rebound (Rein & Liu 2012; Rein & Spiegel 2015). Simulations were started with e = 0 and the semimajor axes listed in Table 4 for the N_p = 2 solution. For both planet candidates, only very small fluctuations in a and e are seen. For DMPP-4b, we find semimajor axis r.m.s. variability over 25 000 years of 5.5 × 10⁻⁷ au (i.e. 0.0014 per cent of the semimajor axis). A very low degree of eccentricity of e_b = 0.00020 ± 0.00009 is established. For the putative DMPP-4c, we find corresponding rms semimajor axis variability of 4.7 × 10⁻⁷ au (i.e. 0.0007  per cent of the semimajor axis) and eccentricity e_c = 0.00007 ± 0.00003. For planets that formed in circular orbits, or evolved into non-eccentric orbits via steady tidal circularization, we expect a stable configuration. Guided by the upper limit to the eccentricities we found in Table 4, we conducted simulations that started with e = 0.065. The respective semimajor axis rms variabilities for DMPP-4b and DMPP-4c are 0.0023 per cent and 0.0019  per cent and the corresponding eccentricities were e_b = 0.071 ± 0.030 with a full range of 0.017 < e_b < 0.107 and e_c = 0.055 ± 0.019) with a range of 0.024 < e_c < 0.079. Apart from the 680 yr periodic eccentricity variability, the orbits did not show any further evidence of evolution on the 25 000 yr time-scale of the simulations. The stability of the orbits is not surprising given that the closest approach of the planets is 0.021 au, while the respective Hill’s radii of DMPP-4b and c are r_H = 0.0007 and 0.0013 au.

6.2 Fundamental properties of the exoplanet candidates

Fig. 11 shows the period-mass diagram for planets with known mass determined with ≤ 20 per cent uncertainty. Only archival planets⁷ with true dynamical mass estimates and tabulated T_eff and R_* are plotted (e.g. planets with lower mass estimates, |$m_\textrm {p}\, \textrm {sin}\, i$|⁠, are excluded). The N_p = 1, DMPP-4b, solution using RV and BIS data (Table 4, column 4) and the N_p = 2 solution (Table 4, column 2) are shown. The planet candidates are closer to the Neptune desert than the planets orbiting DMPP-1, −2, and −3 (Barnes et al. 2020; Haswell et al. 2020; Staab et al. 2020). All DMPP planets possess lower mass estimates to better than 20 per cent with the exception of DMPP-1e, which has an uncertainty of 22 per cent. The points in Fig. 11 are colour coded according to the planet equilibrium temperatures, T_eq. The value obtained depends on the assumed albedo. In the Solar system Bond albedos, A_B vary between ∼0.088 for Mercury (Mallama 2017), through a mean value of 0.36 for gas giant planets (Hanel et al. 1983; Pearl et al. 1990; Pearl & Conrath 1991; Li et al. 2018) up to 0.76 for Venus (Haus et al. 2016). Using these values in the equation |$T_{\rm eq}~= T_{\rm eff}\sqrt{R_*/2a}(1-A_{\rm B})^{0.25}$| therefore gives a range of plausible equilibrium temperatures. For the N_p = 1, respective T_eq  ∼ 1608, 1472, and 1152 K are implied for DMPP-4b. For N_p = 2, the corresponding temperatures for DMPP-4b are T_eq  ∼ 1809, 1655, and 1295 K and for DMPP-4c are T_eq  ∼ 1388, 1270, and 994 K. The equilibrium temperatures for the Mercury-like, low albedo cases are thus around 400−500 K higher.

Using the probabilistic Forecaster (Chen & Kipping 2017) and the minimum mass for our preferred N_p = 1 solution (Table 4), we estimate the minimum radius of DMPP-4b as |$R_{\rm b} = 3.49^{+1.47}_{-1.04}\, \rm {R}_{\hbox{$\oplus $}}$|⁠. For the N_p = 1 solution, minimum radii of |$R_{\rm b} = 2.72^{+1.17}_{-0.81}\, \rm {R}_{\hbox{$\oplus $}}$| and |$R_{\rm c} = 3.49^{+1.46}_{-1.02}\, \rm {R}_{\oplus }$| are predicted. DMPP-4b and c would thus be likely to possess few per cent H₂ dominated atmospheres (Zeng et al. 2019). The planets may have lost much of their atmospheres, through radiation-driven mass-loss, evolving downwards through the Neptune Desert (Fig. 11). Given this probable evolutionary history, and the ongoing high irradiation from the F7V host star, DMPP-4b and c may deviate from expectations based on the general planet population. If DMPP-4b and DMPP-4c are rocky ‘Chthonian’ planets stripped of their atmospheres, they are likely hot enough that their surfaces comprise liquid magma oceans. Temperatures of 1100–1500 K are typical for molten magma on the Earth (Sigurdsson et al. 2015), as noted by Staab et al. (2020). This more extreme scenario could mean that DMPP-4b and c lie below the exoplanet radius valley, with radii of less than |$1.7\, -\, 2$| R_⊕ (Fulton et al. 2017; Zeng et al. 2021), in turn implying densities of between 1.6 × ρ_Earth and 2.5 × ρ_Earth if R ≤ 1.7 R_⊕. This could potentially affect mass loss since the surface gravities would be 2.8−4.2 greater, but needs detailed modelling to quantify. If substantial mass-loss does not occur, and the planets are indeed atmosphereless, to satisfy the DMPP hypothesis, we might expect that there are additional interior and as yet undetected mass-losing planets in the DMPP-4 system.

Direct size determinations from transits would reveal bulk properties and potential atmospheric compositions. If the rotation period of DMPP-4 is 11.6 d and the orbital angular momentum of the planets is aligned with the stellar rotational angular momentum, there is a high transit probability (Haswell et al. 2020). We used the Lightcurve Analysis Tool for Transiting Exoplanets (LATTE) to search for transits in the 2 min cadence TESS photometry (Eisner, Lintott & Aigrain 2020). LATTE corrects the data for residual systematics using an iterative non-linear filter (Aigrain & Irwin 2004). A Box-Least-Squares period search is performed on the flattened light curves. We do not find evidence for transits in any of the TESS sectors (see Section 3) observed to date. Of course, these planets may not be aligned exactly edge-on, so may not transit. In this case, information from phase curves, e.g. using JWST could be revealing.

7 CONCLUSION

Because it is so bright, DMPP-4 is an important target for follow-up characterization. Although the current data do not enable us to unambiguously distinguish between a system with one or two planets, there is significant evidence to favour models with a single planet. Our currently preferred solution indicates a 12.2 M_⊕ planet in a 3.498 d orbit. The potential activity contributions, amplitude modulations and periods of the DMPP-4RV signals we observe indicate that an intensive RV monitoring campaign is needed. This would better enable us to distinguish between multiplanet solutions and activity signals on the time-scales of the stellar and planet candidate periods. It may also reveal additional planets below our current mass detection threshold.

ACKNOWLEDGEMENTS

JRB, MRS, and CAH were funded by STFC under consolidated grants ST/T000295/1 and ST/X001164/1. The RV data were obtained at through The OPTICON Trans-National Access Programme (FP7II, 2013–2016). This work also used observations from the Las Cumbres Observatory global telescope network. This paper includes data collected by the TESS mission. Funding for the TESS mission is provided by the NASA’s Science Mission Directorate. Data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia) have been used. The Gaia data were processed by the Gaia Data Processing and Analysis Consortium (DPAC; https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. This research 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 (see https://exoplanetarchive.ipac.caltech.edu for further details of data sources).

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author. The reduced radial velocities and activity indices are tabulated in Appendix A.

Footnotes

References

APPENDIX A: DATA

See Tables A1 and A2.