DMPP-4: Candidate sub-Neptune mass planets orbiting a naked-eye star

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 − 1 amplitude is found to persist over a > 4 year timescale. Although the non-simultaneous photometric variability in four tess sectors supports the view of an inactive star, we identify periodic photometric signals and also ﬁnd 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 b = 3 . 4982 +0 . 0015 − 0 . 0027 d and corresponding minimum mass of m b sin i = 12 . 2 +1 . 8 − 1 . 9 M ⊕ . Without further high time resolution observations over a longer timescale, we cannot deﬁnitively rule out the purely Keplerian model with 2 candidates planets with P b = 2 . 4570 +0 . 0026 − 0 . 0462 d, minimum mass m b sin i = 8 . 0 +1 . 1 − 1 . 5 M ⊕ and P c = 5 . 4196 +0 . 6766 − 0 . 0030 d and corresponding minimum mass of m b 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.


INTRODUCTION
The dearth of planets, now commonly dubbed the "Neptune Desert", is a very clear region in the planetary mass vs orbital period plane (Mazeh et al. 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 et al. (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 inward 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 vs orbital period diagram (Mazeh et al. 2016).The lack of in-termediate mass planets at short periods exists because the timescale 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;Jenkins et al. 2020;Díaz et al. 2020;Armstrong et al. 2020;Burt et al. 2020;Jordán et al. 2020;Dreizler et al. 2020;Smith et al. 2021;Murgas et al. 2021;Kanodia 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 characterisation.
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 mv < 10, making followup characterisation 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 (Haswell et al. 2020;Barnes 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 masslosing 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(Rappaport et al. , 2014;;Sanchis-Ojeda et al. 2015).
DMPP has proved very successful, with an effectively 100% 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 §2.In §3, we discuss ground-based and tess photometric observations and search for evidence of the stellar rotation period.Spectroscopic activity indicators are analysed in §4 before searching for Keplerian signals in the RVs in §5.A summary and further discussion in §6 is followed by a brief conclusion in §7.

Parameter value Reference
Spectral Type F5 -F8V 1,2,3 Mean G [mag] 5.5974 ± 0.0005 4 Mean BP [mag] 5.8437 ± 0.0009 4 Mean RP [mag]  A range of spectral types between F5 (Roman 1949) and F8V (Eggen 1960) have been estimated2 .Table 1 includes estimates of V and also B-V = 0.4344 ± 0.0076 derived from photometric conversions3 documented in Gaia Data Release 3 (Gaia Collaboration et al. 2016Collaboration et al. , 2023)), which is more consistent with an earlier spectral type of F5V (Pecaut et al. 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, vmac, the equatorial rotation velocity, v sin i and vmac have been obtained following the method described in Murphy et al. (2016).Macroturbulence contributes significantly to the broadening, with vmac = 7.0 ± 2.0 kms −1 .The derived v sin i = 7.5 +0.5  −2.0 kms −1 is consequently significantly lower than would be measured without including vmac.

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 hour 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 Aug 28.The most significant peak 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 lightcurve using the image code DoTS (Collier Cameron 2001) for a single spot with T phot − Tspot = 2000 K (Berdyugina 2005), finding that a rspot = 5.04 o 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 HK ) activity of DMPP-4, especially given that the largest sunspots are never this large when the solar maximum log(R HK ) is still considerably greater than is seen on DMPP-4.Moreover, we find that a 5.04 o spot would induce a K = 82 ms −1 RV signal, which is more than an order of magnitude greater than we find in our RV data ( §4).We therefore conclude that the periodicities recovered from the LCOGT B band data are not genuine activity-related signals.

TESS photometry
The Transiting Exoplanet Survey Satellite (tess, 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 d, 6.47 d, 6.40 d.The respective half width at half maxima (HWHM) of these peaks of 0.33 d, 0.44 d and 0.33 d give an indication of the spread of power, but should not be treated as period uncertainties, as discussed by Vander-Plas (2018).We also note that the longer period at 11.6 d (HWHM = 1.4 d) is also significant in the Sector 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 (Van-derPlas 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 package6 , 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 timeseries and the periodicities are often localised.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 o .A spot of this size is expected to induce RV variability of K < 16 cms −1 at 550 nm and is thus below the level that could easily be detected by an RV instrument achieving ∼ 1 ms −1 precision.If a Solar-like facular/spot area ratio of 11.6 is present on DMPP-4 (Shapiro et al. 2014; see also Barnes et al. 2023), a 0.082 o region of plage surrounding a 0.024 o spot could induce a flux half-amplitude of ∼ 29 ppm.In this case, the expected RV amplitude is K ∼ 38 cms −1 at 550 nm.The RV amplitude rises by an order of magnitude for a cool spot with radius 0.2 o and surrounding plage of 0.68 o .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.1, assuming uniformly distributed periods with the indicated ranged.Bottom: Distribution of stellar axial inclination assuming tess periodicities identified in Fig. 2.

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 o .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 lightcurve 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 lightcurves by McQuillan et al. (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 lightcurve 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 lightcurve morphology is determined by either the relative strength, the distribution of the spots, or by both effects.Basri & Nguyen (2018) used the Kepler lightcurves 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 Sector 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 timescale of a single rotation, resulting in a rapidly changing lightcurve.The lightcurve 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 doubledip mode.
In summary, we believe that DMPP-4 exhibits small spots, as evidenced by the low log(R 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.

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(Haswell et al. , 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 Fig. 4 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 ms −1 (asymmetric Gaussian distribution) estimates listed in Table 1.A uniform distribution of axial inclinations of 0 o < i < 90 o was simulated.We find a resulting modal period of P = 8.89 d with a range between 2.74 d to 11.10 d at the 16 per cent and 84 percent intervals.A simulation with 60 o < i < 90 o 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 Fig. 4 shows the axial inclination distributions, i, assuming period distributions with P = 6.0 d and 11.0 d from the TESS periodicities identified in §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.8o (23.2 o to 43.8 o at 16 per cent and 84 percent intervals) is found.For P = 11.0 d, inclinations > 49.7 o (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 o , v sin i must be ≤ 6.35 kms −1 when P = 11 d.

ACTIVITY INDICATORS
Spectroscopic observations of DMPP-4 were made with 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 Apr-May on 3 nights over a time span of 4 nights; and again later in 2015 Nov on 5 nights.Similarly, observing runs in 2016 May and 2016 Nov-Dec each resulted in 4 nights of observations over a 5 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 4 nights (HN-16A) were made in 2016 Jul 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 s to 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 minute stellar oscillations were minimised.
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  & 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 ms −1 and 2.26 ms −1 .The binned data uncertainties are respectively 1.69 ms −1 and 1.60 ms −1 .A measurement precision of ∼ 1 − 2 ms −1 has been demonstrated with SOPHIE on RV standard stars on timescales of a few tens of days (Bouchy et al. 2013).For HARPS-N, the formal uncertainties for single 200 s observations is 1.17 ms −1 , and 0.73 ms −1 for 900 s binned data.

Kappa correlation with RVs
We investigated correlations between possible diurnal chromatic systematics and the RVs.Bourrier et al. (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 colourdependent (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 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.

RV correlations with activity measures
The correlation of RVs with the full width at half-maxiumum (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 anti-correlation 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 ms −1 results in r = −0.64 and p = 0.17.For individual observing runs, the BIS anti-correlation 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 anti-correlation 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% (i.e.p = 0.097).Of the three activity indicators, only BIS shows a consistently moderate-good correlation.

Periodicities in activity measures
The right 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 −6 and ∆ LogL = 27.2 relative to the scenario with no signal.A neighbouring 2.68 d peak with logL2.54−2.68< 5 (∆ LogL = 22.9) relative to the main peak is also present.This periodicity, taken with the moderate anti-correlation 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 HK ), suggesting that DMPP-4 is 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.Barnes et al. (2023) finds 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 Prot/2 and Prot/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 Prot/2, while a low latitude spot could simultaneously induce RV variability at Prot/2 and BIS variability at Prot/4.A pair of low-latitude spots located at active longitudes 180 o 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 Prot/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 sin i.
The BIS periodicities of 2.54 d 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, Prot 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 §3.2.There is tentative evidence for low significance peaks (1 -10 per cent FAP) at ≥ 10 d in the Sindex periodogram, which shows a most significant peak at 11.1 d (∆ LogL = 15.9),again in close agreement with the tess observations.

RV ANALYSIS
We began by assessing the coherence of periodic RV signals over time as more data points were added.Figure 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 6 ob-   Parameter Model fits to RV data only RV + BIS Systemic RVs and white noise parameters Table 4. Posterior RV Model parameters and derived planet masses with 68.3% confidence uncertainties.The log L statistics in the final row are given for the best fitting (maximum posterior) models.Models with 1-Keplerian (Np = 1) and 2-Keplerians (Np = 2) are tabulated in columns 1 and 2. Column 3 is the Np = 1 + GP solution using only the RVs.Column 4 tabulates parameters for the Np = 1 + GP solution which uses simultaneous RVs and BIS.Eccentricities are all consistent with 0 and therefore upper limits are quoted at 2σ.

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  shows that periodicity at P = 3.498 d is found for the 1-Keplerian models (Np = 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, BF = e (BIC k−1 −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 Np = 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 Np = 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.3and 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 lin-BIScorr 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 timeseries.
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 characterise the effects of activity on the RVs.

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 use to distinguish models since they do not take into consideration the number of model parameters or data points.Kima uses Diffusive Nested Sampling (Brewer et al. 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 §5.1).This has the benefit of correctly enabling inter-comparison of models (Brewer & Donovan 2015;Feroz et al. 2011).With kima, it is possible to either fix the number of Keplerians, Np, or include Np 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 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).

1-Keplerian (Np = 1) solution
Using kima to obtain planetary parameters as described in Standing et al. (2022), we find P b = 3.49733 +0.00093 −0.23513 when we restrict the search to 1 Keplerian (Np = 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  we estimate a lower limit on the BF by setting the number of samples with Np = 0 equal to 1 (Standing et al. 2022;Triaud et al. 2022;Baycroft et al. 2023;Standing et al. 2023).

2-Keplerian (Np = 2) solution
We allowed kima to sample posterior space with Np as a free parameter.In this way it is possible to recover the preferred number of Keplerians.We place a maximum Np = 3, and find the solution with Np = 2 has the highest posterior evidence, with BF = 2397, compared with the Np = 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  4, column 2).With the additive white noise terms, σRV,SO and σRV,HN, the effective uncertainties in the Np = 2 model are 1.68 ms −1 and 1.36 ms −1 , suggesting that the formal uncertainties, particularly for HARPS-N (see §4), may be underestimated.
We investigated the difference between the periods identified by the Np = 2 recursive likelihood method in §5.1 and the posterior sampling method of kima.We performed MCMC posterior sampling on the two periods we identified in §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 Np = 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 Np = 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.

Quasi-periodic Gaussian Process
kima has also been adapted to enable the use of a Gaussian Process (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 ti and tj.Four hyperparameters define the covariance matrix.The amplitude, η1, describes the deviation of the GP models from the mean function.The parameter, η2 is the long term evolution timescale, while η3, 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 η3 < η2, to ensure the validity of using the quasi-periodic GP with a periodic component.The factor, η4, is proportional to the inverse of the square root of the harmonic complexity, Γ (e.g.see Rajpaul et al. 2015, Barragán et al. 2022and Nicholson & Aigrain 2022).From equation 1, Γ = 2/η 2 4 , and describes the roughness or complexity of variations within each characteristic period, η3.

Keplerian + quasi-periodic Gaussian Process
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 §3 and §4.Following our finding of photometric periodicity (Section 3), we first trialled a GP prior for η3 of G[11.6 ± 1.5] d (i.e. a Gaussian distribution).Given the likely activity signatures at Prot/2 and Prot/4 discussed above in §4.3, we expect characteristic activity signatures on timescales shorter than than Prot.This suggests some degree of harmonic complexity is required to account for any structure related to the characteristic period, η3.Nicholson & Aigrain (2022) note, for instance, that analysis of stellar lightcurves 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 < η4 < 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 d, 3.5 d and 5.4 d, the posterior samples for Np = 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 η2 and η3.This not surprising since the few-day span of each observing epoch is much less than the observationally informed priors.The posterior distribution for η4 = 0.15 +0.17 −0.04 corresponds to harmonic complexity of Γ = 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 timescales.There are too many fitting parameters and not enough data points to constrain a credible solution.Essentially, the model is too flexible.

A simplified quasi-periodic Gaussian Process
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 Np = 1 + GP model in §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 §5.3.1, or lead to over-fitting 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 lightcurves, 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 η3 of G[2.5 ± 0.5] d.In light of the posterior distribution for η4 that we found in §5.3.1, we followed Nicholson & Aigrain (2022), but fixed Γ = 1 (i.e.we fixed η4 = √ 2).With these priors, posterior sampling leads to a preference for Np = 1, with P = 3.4982 +0.0015 −0.0027 d and evidence of BF = 111.4.As in §5.3.1, the posterior distribution is dominated by the priors for η2, though the data appears to have had slight influence on the posterior GP periodic component η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 Np = 1 model.The posteriors show that inclusion of the GP resulted in reduced significance of the periodicities at 2.457 d and 5.419 d found for the purely Keplerian Np = 2 model.

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 timeseries doubles the size of the data set to 142, or if we consider the effective data set size, to 42 as in §5.3.2.The flexibility of the model is thus reduced with inclusion of the activity timeseries.

Model with simultaneous BIS activity indicator
We used the same model priors as in §5.3.1.We again found that the η2 and η3 posterior distributions matched the priors while a high harmonic complexity is preferred, with η4 = 0.15 +0.05 −0.04 (Γ = 89 +76 −39 ).The posterior samples demonstrate an over-density at P ∼ 3.5 d, though for Np = 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 timeseries, the model is too flexible.

Model with simultaneous BIS activity and fixed harmonic complexity
Finally, with fixed Γ = 1 (η4 = √ 2), as in §5.3.2, we used a prior on η3 of G[2.5 ± 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 Np = 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 timeseries in §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 Np = 2 solution; accounting for σRV,SO and σRV,HN yields respective effective uncertainties of 2.22 ms −1 and 1.48 ms −1 .

SUMMARY AND DISCUSSION
Although there is clear evidence for a single candidate Keplerian (Np = 1) with P b = 3.49791 +0.00046  −0.00143 d and m b sin i = 11.6 +1.6  −1.5 M⊕, the purely Keplerian model indicates strong evidence favouring a 2-Keplerian (Np = 2) solution, with a relative Bayes Factor, BF = 2397.The preferred Np = 2 solution does not include the Np = 1 period, but rather recovers periods of P b = 2.4570 +0.0026 −0.0462 d and Pc = 5.4196 +0.6766 −0.0030 d with respective derived minimum masses of m b sin i = 8.0 +1.1 −1.5 M⊕ and mc sin i = 12.2 +1.4  −1.6 M⊕.The Np = 2 solution derived from likelihood periodogram searches, with P b = 3.498 d and Pc = 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 marginalised 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% compared with the 2.4570 d period and 50% compared with the 5.4196d period).
Despite the low log(R 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 HK ) than observed.The hypothesis underlying DMPP is that circumstellar absorption depresses the log(R HK ), so this is as expected.Because the Np = 2 solution recovers periods that are close to one half and one quarter integer factors of the inferred rotation period (Prot/2 and Prot/4), as discussed in detail in §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 b = 3.4982 +0.0015 −0.0027 d and m b sin i = 12.6 +1.6  −1.8 M⊕.Including the simultaneous BIS measurements augments the Bayesian evidence, yielding P b = 3.4982 +0.0022 −0.2327 d and m b sin i = 12.2 +1.8 −1.9 M⊕ with a BF = 1038.
Thus, although the purely Keplerian model prefers Np = 2, the model restricted to Np = 1 and the models with a GP and no Np 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 Np = 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 Np = 2 model.Further, there is (ii) a not-insignificant appearance of the the 3.498 d periodicity in the model posterior samples for the Np = 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 the Np = 2 solution.
6.1 Orbital stability of planets in the Np = 2 model We checked the stability of the orbital solution for the Np = 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 Np = 2 solution.For both planet candidates, only very small fluctuations in a and e are seen.For DMPP-4 b, we find semi-major axis r.m.s.variability over 25,000 years of 5.5 × 10 −7 AU (i.e.0.0014% of the semi-major axis).A very low degree of eccentricity of e b = 0.00020 ± 0.00009 is established.For the putative DMPP-4 c, we find corresponding r.m.s.semi-major axis variability of 4.7 × 10 −7 AU (i.e.0.0007% of the semi-major axis) and eccentricity ec = 0.00007 ± 0.00003.For planets that formed in circular orbits, or evolved into non-eccentric orbits via steady tidal circularisation, 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 semi-major axis r.m.s.variabilities for DMPP-4 b and DMPP-4 c are 0.0023% and 0.0019% and the corresponding eccentricities were e b = 0.071 ± 0.030 with a full range of 0.017 < e b < 0.107 and ec = 0.055±0.019)with a range of 0.024 < ec < 0.079.Apart from the 680 year periodic eccentricity variability, the orbits did not show any further evidence of evolution on the 25,000 year timescale 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 radii of DMPP-4 b and c are rH = 0.0007 AU and 0.0013 AU.

Fundamental properties of the exoplanet candidates
Fig. 11 shows the period-mass diagram for planets with known mass determined with ≤ 20% uncertainty.Only archival planets 7 with true dynamical mass estimates and tabulated T eff and R * are plotted (e.g.planets with lower mass estimates, mp sin i, are excluded).The Np = 1, DMPP-4 b, solution using RV and BIS data (Table 4, column 4) and the Np = 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 (Staab et al. 2020;Haswell et al. 2020;Barnes et al. 2020).All DMPP planets possess lower mass estimates to better than 20% with the exception of DMPP-1 e, which has an uncertainty of 22%.The points in Fig. 11 are colour coded according to the planet equilibrium temperatures, Teq.The value obtained depends on the assumed albedo.In the Solar System Bond albedos, AB vary between ∼ 0.088 for Mercury (Mallama 2017), through a mean value of 0.36 for gas giant planets (Li et al. 2018;Hanel et al. 1983;Pearl et al. 1990;Pearl & Conrath 1991) up to 0.76 for Venus (Haus et al. 2016).Using these values in the equation Teq = T eff R * /2a(1 − AB) 0.25 therefore gives a range of plausible equilibrium temperatures.For the Np = 1, respective Teq ∼ 1608, 1472 and 1152 K are implied for DMPP-4 b.For Np = 2, the corresponding temperatures for DMPP-4 b are Teq ∼ 1809, 1655 and 1295 K and for DMPP-4 c are Teq ∼ 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 Np = 1 solution (Table 4), we estimate the minimum radius of DMPP-4 b as R b = 3.49 +1.47  −1.04 R⊕.For the Np = 1 solution, minimum radii of R b = 2.72 +1.17 −0.81 R⊕ and Rc = 3.49 +1.46 −1.02 R⊕ are predicted.DMPP-4 b and c would thus be likely to possess few per cent H2 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-4 b and c may deviate from expectations based on the general planet population.If DMPP-4 b and DMPP-4 c 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-4 b 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 minute cadence tess photometry (Eisner et al. 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 lightcurves.We do not find evidence for transits in any of the tess sectors (see §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.

CONCLUSION
Because it is so bright, DMPP-4 is an important target for follow-up characterisation.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-4 RV signals we observe indicate that an intensive RV monitoring campaign is needed.This would better enable us to distinguish between multi-planet solutions and activity signals on the timescales of the stellar and planet candidate periods.It may also reveal additional planets below our current mass detection threshold.

Figure 1 .
Figure 1.LCOGT B band lightcurve and tess PDCSAP lightcurves for sectors 14, 15, 40 and 41.The data are binned with a 0.1 d interval and are mean-subtracted.The change in flux, ∆F , is plotted in parts per thousand for the LCOGT data and parts per million for the tess data.

Figure 3 .
Figure 3. Continuous wavelet transform (CWT) for the tess photometry.The most significant continuous wavelet transform (CWT) powers are shown in red.The shaded area denotes the wavelet cone of influence region.The right hand panel shows the global spectrum for Sectors 14 and 15 (dotted line), Sectors 40 and 41 (dashed line) and all sectors (solid line).

Figure 4 .
Figure 4. Monte Carlo simulations for stellar rotation period, P , and axial inclination, i Top: The distribution of stellar rotation periods based on the parameters in Table1, assuming uniformly distributed periods with the indicated ranged.Bottom: Distribution of stellar axial inclination assuming tess periodicities identified in Fig.2.

Figure 5 .
Figure5.Left: Activity index correlations with RV for SOPHIE (panels a-c) and HARPS-N (panels d-f) observations.The points are colour coded according to the observation semester.Right: Activity index periodograms for line full width at half maximum (FWHM), bisector inverse span (BIS) and Calcium S-index.FAP at 0.1%, 1% and 10% are (top to bottom) are indicated by the dashed horizontal lines.

Figure 6 .
Figure 6.Stacked periodograms of the RVs (top) and residual RVs after accounting for the first signal (bottom).The dotted lines indicate the end of each observing run.

Figure 7 .
Figure 7. Radial velocity log likelihood periodograms using (a) RVs only and (b) with a BIS linear correlation (RV+linBIScorr).The periodograms for single (Np = 1) Keplerian signals are shown in the top panels.The period searches with a second recursively added signal (Np = 2) are shown in the bottom panels.The vertical dashed line indicates the BIS period of 2.54 d.Horizontal dashed lines are the 0.1%, 1% and 10% FAPs as in Fig. 5.

Figure 8 .
Figure 8. Single Keplerian RV solution (Np = 1) for combined SOPHIE (SO) and HARPS-N (HN) data sets taken in semesters 15A, 15B, 16A, 16B and 19A.(a) RVs vs observation time and (b) phased RVs showing fit with 68% and 95% model uncertainties.The light blue and pink data uncertainties indicate the combined formal uncertainties and additive white noise.

Fig. 7 (
Fig. 7 (top left)  shows that periodicity at P = 3.498 d is found for the 1-Keplerian models (Np = 1) for both RV and RV+linBIScorr models.Table3 lists the likelihood values,

Figure 9 .
Figure 9.As in Figure 8 for the 2-Keplerian (Np = 2) solution.(a) The RVs with solution curve and (b) the phased RVs plotted for the 5.419 d candidate signal (top) and the 2.457 d signal (bottom).

Figure 10 .
Figure 10.Single Keplerian solution (Np = 1) with Gaussian Process (GP).(a) Solution using only the RVs with a GP.The RV component is represented by a short dashed line and the GP by a long dashed line.As in Figure 8, the RV + GP curve shaded regions show 68% and 95% model uncertainties.(b) The phased RVs corresponding to panel (a).(c) RVs for the RV + BIS solution with a GP and (d) corresponding phased RVs.(e) The BIS timeseries for the RV + BIS solution showing the GP as a long dashed line.

Figure 11 .
Figure 11.The Planet mass vs orbital period diagram for planets with known mass.Each planet is colour coded according to its estimated equilibrium temperature, Teq (assuming a Bond Albedo A B = 0.36).The preferred Np = 1 solution in column 4 of Table 4 for DMPP-4 b is shown by the large square symbol.The Np = 2 candidate planets, DMPP-4 b and c, are shown by the large upright triangle and circle.DMPP-1 b, c, d & e, DMPP-2 b and DMPP-3A b are also highlighted.

Table 2 .
correlations between RV and chromatic index, κ, full width at half maximum (FWHM), RV vs bisector span (BIS) and RV vs Ca H&K S-index.Pearson's r, and student's p statistic are shown for each correlation for all SOPHIE observations combined and all HARPS-N observations combined (top rows).The remaining rows show the same statistics for individual observing runs.

Table 3 .
Candidate maximum likelihood periodogram search summary for (a) RV period search and (b) RV period search with a linear BIS correlation (RV+linBIScorr model).The corresponding Bayesian Information Criteria (BIC) are also given along with the ∆BIC values for BIC(Np − (Np − 1)) in parentheses and directly derived estimates of the corresponding Bayes Factors (BF).

Table A1 .
SOPHIE RVs and activity indices with their respective uncertainties.