© The Authors 2025

1 Introduction

Since the first observation of a planet orbiting a solar-type star (Mayor & Queloz 1995), more than 5500 exoplanets have been confirmed, the majority detected by the Kepler satellite with the transit method (Borucki et al. 2010). In particular, this has opened the possibility of studying star-planet interactions shaping orbital architecture both at the scale of individual star-planet systems (e.g. Patra et al. 2017; Barker et al. 2024) and for larger samples (e.g. McQuillan et al. 2013; Messias et al. 2022; García et al. 2023).

Observing close-in planets, which exhibit orbital periods shorter than 1 day and experience strong tidal (e.g. Mathis 2018) and magnetic (e.g. Strugarek 2018) interactions with their host star, is crucial to calibrate dissipation processes related to these phenomena. Unfortunately, ultra-close-in planets, depending on the planetary mass, its orbital period, and the rotational period of the star, may be pushed away or experience orbital decay and disruption on very short timescales (e.g. Ahuir et al. 2021). This means that we are more likely to observe planets that are less dynamically perturbed and located further away from their host. Among all the planets detected using the transit method, approximately 97% have orbital periods exceeding one day¹.

This highlights the need for alternative detection methods to study ultra close-in planets. Statistically, non-transiting systems are likely far more numerous than transiting ones, suggesting that observing them could increase our chances of detecting planets very close to their star.

Having rather large modulations in their orbital phase curve, close-in planets can be detected with photometric observations, even if they are not transiting. Until now, only a few non-transiting companions have been found (Millholland & Laughlin 2017; Lillo-Box et al. 2021; Cullen & Bayliss 2024). These modulations are attributed to three main physical processes (Shporer 2017). The first is atmospheric, encompassing stellar light reflected by the companion’s atmosphere and thermal emission, due to heating of the companion’s atmosphere. The second is the tidal ellipsoidal distortion, induced by the tidal forces exerted by the companion on its host star. The third process is the Doppler boosting effect, caused by variations in the radial velocity between the object and the observer, leading to a modulation in the observed flux of the star. Fig. 1 shows an example of the three modulation components for a non-transiting planet, as well as the composite phase-curve, generated using Eq. (4), (7), (9), and (10) from Shporer (2017) and a black body emission law computed at the planetary equilibrium temperature, as a function of the star-planet-observer phase angle (Seager 2010). A realistic set of stellar and planetary parameters was chosen, analogous to the ones expected for a simplified bookcase-like system. A hot Jupiter with a radius of 1 R_jup, a mass of 1 M_jup, and an orbital period of 2 days was considered. The eccentricity was set to zero, while the inclination of the orbit on the line of sight was 60°, the maximum value taken for non-transiting systems (Winn 2014). We opted for a Bond albedo of 0.5 (e.g. Seager 2010), which is typical for a hot Jupiter-type planet, and picked an atmospheric recirculation rate of 50% (Fortney et al. 2021). We chose 1 M_⊙, 1 R_⊙, and T_eff = 6000 K as stellar parameters. In this configuration, the atmospheric processes outweigh the other two phenomena.

While the atmospheric processes are mostly driven by stellar irradiation, Doppler boosting and tidal ellipsoidal distortion are induced by the bodies’ gravitational interaction, and hence they depend on the stellar and planet masses. The greater the companion’s orbital velocity and mass, the stronger the Doppler boosting. However, depending on the mass ratio between the star and the planet, for short orbital periods, the ellipsoidal distortion effect might outweigh the other two phenomena in amplitude.

The layout of the paper is as follows. In Sect. 2, we present the sample of Kepler light curves we consider, and we explain how we look for a non-transiting companion signature. In Sect. 3, we perform a fit of the photometric modulations to extract additional constraints on the candidates’ global parameters. In Sect. 4, we discuss the perspectives that the characterisation of this sample may open.

Fig. 1 Simulated contributions of the three modulations to the phase-curve amplitude according to the phase angle. The modulation linked to atmospheric processes is shown in black and includes reflection and emission modulations. The ellipsoidal distortion is shown in purple, while the Doppler boosting is in orange. The overall contribution is shown by the composite phase-curve in yellow. Here, an eclipse would occur at the inferior conjunction, while the transit would happen at the superior conjunction with the observer.

2 Observational data and search methodology

In mono-band photometry and when blindly searched for in a large sample, non-transiting companion signatures may easily be confused with modulations arising from stable active longitudes of fast-rotating young stars. Therefore, prior knowledge of the rotation of the stellar host is necessary to undertake our analysis. To avoid confusing exoplanetary signals with the rotation period (e.g. García et al. 2014), we considered the catalogue by Santos et al. (2019, 2021), which consists of 55 232 main-sequence and subgiant FGKM stars.

To illustrate the procedure we describe below, we selected KIC 5697777, an F-type star with a rotation period of 3.90 days. The results obtained for this star are presented in Fig. 2 and are discussed in this section.

Fig. 2 Flux analysis of KIC 5697777. (a) Light curve calibrated according to the KEPSEISMIC method and filtered at 20 days. (b) Power spectral density of the light curve. The rotation period of the star is taken from Santos et al. (2021) and is located within the orange area. The green area shows the high-amplitude signal peak at around 0.5 days. The second high-amplitude signal peak located at around 0.25 days is the second harmonic of the orbital signal. (c) Phase-folding of the light curve with an additional filter to remove the stellar rotational contribution. The red curve is a smoothing of the black curve, obtained with a triangular filter. (d) Wavelet power spectrum of the light curve. A high-amplitude periodicity is highlighted by orange and red, whereas yellow and blue denote a low-power amplitude. The yellow bands with a constant period throughout time correspond to the high-amplitude peak and its second harmonic, identified in panel b. The reddish band corresponds to the rotation period and its harmonics, highlighted in orange in panel b. The shaded area represents the cone of influence of the reliable periodicities.

2.1 Unravelling signatures of non-transiting objects

We used Kepler light curves calibrated with the KEPSEISMIC method² (García et al. 2011; Pires et al. 2015), whereby signals with a period longer than 20 days were filtered out (see Santos et al. 2019) in order to be less affected by the instrumental noise and because we are interested in close-in objects. The light curve of KIC 5697777 is presented in panel a of Fig. 2.

The power spectral density (PSD) of a stellar light curve (Fig. 2, panel b) allows for the decomposition of the signal into its various frequency components. Identifying a high-amplitude peak, typically from a few tens to a few thousands of parts per million, is then akin to detecting a peak that corresponds to the passage of a non-transiting object. As is seen in Fig. 1, if tidal interactions dominate the flux modulation, the second harmonic of the orbital signal is expected to be the largest peak of the harmonic pattern (when still excluding the stellar rotation period signal). Nevertheless, it is the frequency of the first harmonic that corresponds to the orbital period of the non-transiting companion object.

To identify a non-transiting object signal, we computed the PSD of each of the 55 232 stars’ light curves. To remove false positives due to low-inclination stellar binary systems, we focused on systems with a signal at a period distinct from the rotation period of the targeted star. In a close binary stellar system, tidal effects typically tend to synchronise the rotational periods of the two stars (e.g. Zahn 1975). In contrast, a planetarymass companion generally lacks sufficient mass to induce such synchronisation. The aim here is to identify significant high-amplitude peaks in the PSD, different from the peak associated with the stellar rotation period. Given that only close-in objects should be detectable with our method, we focused on objects with short orbital periods and we searched for peaks with a frequency higher than 5 µHz (2.3 days). Below 5 µHz, the main signal in the PSD of a main-sequence solar-like star usually comes from the surface rotational modulation, granulation, and activity. In addition, because the object would be too distant from the star, the phase curve modulation amplitude would be too low to be detectable in visible photometry (e.g. Shporer 2017). A high-amplitude peak with a frequency higher than 5 µHz can then be associated with the presence of a non-transiting object. If we know the stellar rotation period and its harmonics, we do not expect any additional excess power of this nature associated with the dynamics of the star alone, when considering a main-sequence solar-type star (e.g. García & Ballot 2019).

A threshold limit was applied in the PSDs, to distinguish sharp high-amplitude peaks associated with non-transiting companions from background noise. The statistical distribution of the PSD function is a χ² distribution with two degrees of freedom (Woodard 1984). Following Appourchaux et al. (2000), we calculated the power level at which a peak due to noise has a probability, p_det, of appearing over a frequency range, Δ_det. Since the number of considered bins in the bandwidth will increase with the length of the time series, this threshold depends on the observation time,T , which is about 4 years for the Kepler observations. The ratio of the statistical χ² threshold to the average signal amplitude, $\tilde{S}$, is given by $s_{\det }/\tilde{s}\approx \ln (T)+\ln \left({\text{Δ}}_{\det }\right)-\ln \left(p_{\det }\right)$, with p_det the detection probability threshold, T in µs and Δ_det in microhertz. Here, Δ_det is within the range of 5 to 284 µHz, the Nyquist frequence of Kepler observations. After some empirical tests, we decided to consider a $s_{\det }/\tilde{s}$ value of approximately 260, which corresponds to an extremely low detection probability threshold, p_det, of 10⁻¹⁰⁰ (this corresponds to a sigma of 259). This choice allowed us to filter out most of the signal related to the harmonics of stellar rotation, while preserving the high-amplitude peaks we were looking for.

We computed $s_{\det }/\tilde{s}$ by dividing the PSD by a factor of α = (9/8)³ × median (PSD), corresponding to the corrected mean of the PSD. Nevertheless, the background signal of the PSD is non-uniform, predominantly influenced by the convective signal and the magnetic activity of the star, which results in higher amplitudes at low frequencies compared to high frequencies (García & Ballot 2019). Hence, adopting such a low probability threshold, p_det, reduces the likelihood of obtaining false positives and improves the accuracy of our signal analysis. With this method, we preselected 4788 stars.

Subsequently, we conducted a visual inspection of the PSDs displaying prominent peaks distinct from the stellar rotation period. This is illustrated in panel b of Fig. 2, in which the PSD exhibits a high-amplitude peak around 0.52 days (green area), while the stellar rotation period of KIC 5697777 corresponds to a period of 3.9 days (orange area). For our remaining stars, the PSDs presenting peaks that are part of the harmonic chain of the rotation period are excluded from the catalogue. However, we decided to keep in our sample targets for which the detected peak lies within the harmonic pattern of surface rotation but exhibits a clearly distinct morphology with respect to it; that is, for which we observe a sharp peak amidst the wide rotation harmonics. This analysis step reduced the considered sample to 362 stars.

2.2 Comparison with existing catalogues

To eliminate potential contaminants in our sample, we conducted a thorough cross-check with existing catalogues. The aim was to remove any stars previously identified as hosting binaries, brown dwarfs, transiting exoplanets, or identified false positives from the Kepler Object of Interest catalogue (KOI, Akeson et al. 2013). The first cross-analysis was performed using data from the National Aeronautics and Space Administration (NASA) Exoplanet Archive³, resulting in the removal of four stars confirmed to harbour transiting planets with an orbital period matching the ones detected in our signals (see Appendix A.1 for an example). We also removed 59 stars identified as false positives of transiting planets. A second comparison was carried out with the catalogue from Berger et al. (2018), allowing for the removal of 15 stars associated with binary systems. The detected signal in the PSD is indeed most likely attributable to the binary star system. Lastly, we compared our results to the Gaia DR3 non-single stars catalogue (Pourbaix et al. 2022), which allowed us to exclude another star from our sample, as it is part of a binary system. A total of 79 stars were removed from our sample following these cross-checks and upon completion of this step, our refined sample consists of 283 stars possibly hosting non-transiting exoplanets.

2.3 Phase-folding, wavelet analysis, and visual inspection

To further refine our study and pinpoint potential indicators of yet undetected transits, false positives, or binary systems, we phase-folded the light curves of the 283 candidates on the orbital period that we obtained from the threshold search conducted in Sect. 2.1. To remove the stellar contribution at lower frequencies, we applied a finite impulse response filter⁴ before phase-folding the signal. The phase-folding analysis for KIC 5697777 is presented in panel c of Fig. 2. The complex morphology of the phase curve suggests that different effects (see Fig. 1) combine to create a possible non-transiting modulation. A subsequent visual inspection was performed to exclude any remaining stars for which the phase-folding of the star signal fails to match the expected pattern of a non-transiting exoplanet passage (as in Fig. 2, panel c). As is shown in Appendix A.2, if the trend of the phase-folded light curve differs from a stable and coherent signal over time, the candidate was discarded. This step allowed for the refinement of our overall sample to 245 candidates.

We also computed the wavelet power spectrum (WPS, Torrence & Compo 1998; Liu et al. 2007) of the light curve, and we only kept stars that present stable modulations over time, as is shown in panel d of Fig. 2 for KIC 5697777. The aim is to determine whether the signal attributed to the non-transiting object exhibits a stable period throughout the different quarters of Kepler ‘s observations. If the signal disappears during at least one quarter, this suggests that the signal may result from contamination by a nearby star. An example of a wavelet power spectrum that does not match a non-transiting signal is provided in Appendix A.3. This step removed 92 stars from our sample.

To ensure that the Kepler signal of every star studied in this work was not contaminated by the light of a stellar neighbour with comparable brightness, we performed a search for potential contaminants using the Gaia catalogue (Gaia Collaboration 2023). We looked for Gaia objects within 40′′ of each target that have an apparent magnitude in the Gaia G band smaller than G_candidate + 1, where G_candidate is the apparent magnitude of our Kepler candidate in the Gaia G band. This results in the identification of 70 Gaia objects near 38 of our targets (see Table A.1 in Appendix A.4). For 11 of these 38 targets, the potential contaminants are located within a few arcseconds, which could result in the contamination of the Kepler light curve. Consequently, we chose to remove these 11 stars from our sample. For the remaining 27 targets, the potential contaminants are at least 10′′ away, making significant contamination issues unlikely. Similarly to what was done in Noraz et al. (2022), we also performed a visual inspection on the 142 stars remaining at this stage. We individually inspected the targets using the Two Micron All Sky Survey (2MASS) images provided on the Kepler Asteroseismic Science Operations Center (KASOC) database⁵ (Kjeldsen et al. 2010). If the star is within the halo of a brighter star, it was excluded from our sample because the short-period signal might come not from a non-transiting object but rather from the nearby star in question. An example of a star whose signal is likely contaminated by a neighbouring brighter star is shown in Appendix A.4. After this final step of selection, we obtained a list of 88 non-transiting exoplanet host candidates. At this point, it is important to underline that none of our candidates are in the list of 60 targets presented by Millholland & Laughlin (2017). We recall that they used an automated algorithm in the search for non-transiting short-period giant planets in Kepler, looking for targets with an orbital period between 1 and 5 days. While they modelled the non-transiting signal and stellar variability together, we used our prior knowledge of stellar rotation to filter out the stellar contribution in the light curve. More importantly, we extended the search towards shorter orbital periods.

3 Fit of the photometric data

As a final step, we attempted to extract additional information on the global parameters of our candidates by fitting their signature in the light curve. Two stars were excluded from this analysis, as their signatures could not be reliably fitted in the light curve. The signal of KIC 7379583 seems contaminated by residual instrumental drifts, leading to peaks in the phase-folded light curve that prevent the fit from converging. For KIC 11350907, the signal amplitude is 84 ppm, while the random noise is around 292 ppm, making it difficult for the fit to converge.

When only photometric data were available, the resulting signal presented in Fig. 1 depends on several free parameters, leading to degeneracies between different effects. To address this, we distinguish two cases. On the one hand, when the phase-folding of the light curve exhibits a modulation that is approximately sinusoidal, we assumed that the non-transiting signal is dominated by emission and reflection. We then fitted only this component in the data, using the following equation (see Shporer 2017; Lillo-Box et al. 2021): $F(\varphi )=A_0-A_{atm}\cos (2\pi [\varphi +\delta ]),$(1)

where A₀ is a flux offset parameter, and the phase, ϕ, is given by $$\varphi =\left(t-T_0\right)/P_{orb},$$(2)

where T₀ is taken to correspond to the time of inferior conjunction, and δ is a phase offset with respect to this conjunction.

The amplitude, A_atm, is given by $$A_{\text{atm}}=57{\alpha }_{\text{refl}}\sin i{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2/3}{\left(\frac{P}{\text{ day }}\right)}^{-4/3}{\left(\frac{R_{\text{p}}}{R_{\text{J}}}\right)}^2\text{ppm},$$(3)

where α_refl depends on the albedo and heat recirculation efficiency. We adopted this approach for the 77 targets listed in Table B.1.

On the other hand, when the modulation in the phase-folding was more complex (see panel c of Fig. 2), we attempted to perform the fit considering the three components of the signal, $$\begin{array}{cc}F(\varphi )& =A_0-A_{\text{atm}}\cos (2\pi [\varphi +\delta ])\\ & -A_{\text{ellip}}\cos (4\pi \varphi )+A_{\text{boost }}\sin (2\pi \varphi ),\end{array}$$(4)

where A_ellip is the amplitude of the tidal deformation, $$\begin{array}{cc}{A}_{\text{ellip }}=& 13{\alpha }_{\text{ellip }}\sin i{\left(\frac{R_{\star }}{R_{\odot }}\right)}^3\\ & \times {\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2}{\left(\frac{P}{\text{ day }}\right)}^{-2}\left(\frac{M_{\text{p}}\sin i}{M_{\text{J}}}\right)ppm.\end{array}$$(5)

with α_ellip a coefficient that depends on the limb and gravity darkening of the host star (Esteves et al. 2013). A_boost is the amplitude of the Doppler boosting contribution, $$A_{\text{boost }}=2.7{\alpha }_{\text{boost }}{\left(\frac{P}{\text{ day }}\right)}^{-1/3}{\left(\frac{M_{\star }}{M_{\odot }}\right)}^{-2/3}\left(\frac{M_{\text{p}}\sin i}{M_{\text{J}}}\right)\text{ppm},$$(6)

where α_boost is the photon-weighted bandpass-integrated boosting factor. After Millholland & Laughlin (2017), we used a parameterisation in terms of the effective temperature, T_eff, for α_ellip and α_boost, $$\begin{array}{c}{\alpha }_{\text{ellip }}=-\left(2.2\times {10}^{-4}{\text{K}}^{-1}\right)T_{\text{eff }}+2.6,\\ {\alpha }_{\text{boost }}=-\left(6\times {10}^{-4}{\text{K}}^{-1}\right)T_{\text{eff }}+\mathrm{7.2.}\end{array}$$(7)

We adopted this approach for nine targets presented in Table B.2 in Appendix B.

In practice, Eq. (1) and (4) were fitted to the data using Bayesian inference. The log-likelihood function is $$\ln {\mathcal{L}}_{\varphi }\left(S_{\varphi },\theta \right)=-\frac{\text{1}}{2}\times \sum _{i=1}^n\frac{{\left(S_{\varphi }\left(t_i\right)-{\mu }_{\varphi }\left({\varphi }_i,\theta \right)\right)}^2}{{\sigma }_{\varphi }{\left(t_i\right)}^2},$$(8)

where S _ϕ(t) is the light curve data and σ_ϕ(t) is the normally distributed noise, considering independent errors. µ_ϕ(ϕ, θ) is the model depending on the parameter, θ, as a function of the orbital phase, ϕ, as expressed in Eq. (1) when following the first approach, and Eq. (4) when following the second approach. We used the UltraNest sampler described in Buchner (2021), which is based on a nested sampling algorithm. To avoid degeneracies between T₀ and δ, we set δ to 0 and fitted for T₀, A_atm, and A₀ for the first approach, and T₀, A_atm, A₀, A_boost, and A_ellip for the second approach. We fixed the period, P_orb, at the value inferred from the PSD of the light curve. The priors are 𝒰 (−10, 10) for A₀ and 𝒰 (−10⁵, 10⁵) for A_atm, A_ellip, and A_boost, where 𝒰 (X, Y) denotes a uniform probability distribution between the bounds X and Y. We also considered uniform priors for T₀ and inferred them with the timing corresponding to the minimum amplitude of the phase-folded light curves. The posterior distributions of two datasets, chosen as examples of each approach, are shown in Fig. B.1 and B.2 of Appendix B.

At this stage, it is necessary to emphasise that the model used to perform the fit relies on a number of simplifying assumptions: the eccentricity is assumed null, the atmospheric contribution from the planet is modelled in the simplest possible way, and the stellar tidal response is considered to arise solely from equilibrium tides, with any potential contribution from dynamical tides being neglected. We therefore advise the reader to keep all these caveats in mind when they consider the results of the fit we present below. Some of the limitations of this approach are discussed in Shporer (2017).

An example of the fit result when the phase-folding of the light curve exhibits a modulation that is approximately sinusoidal is illustrated by Fig. 3. It shows the normalised light curve as a function of orbital phase for KIC 4930913, showing the observational data in grey, along with the best-fit model in orange. The model closely follows the general trend of the data despite significant scatter. The bottom panel displays the residuals, defined as the difference between the data and the model. The residuals are centred around zero and show no evidence of large-scale systematic trends, indicating a satisfactory fit.

An example of the fit result, based on the three components of the signal described by Eq. (4), is illustrated in Fig. 4 for KIC 4264634. The phase-folding of the light curve of this star shows a more complex pattern (in grey). As is seen in Fig. 3, the best-fit model (in orange) follows the general trend of the data and the residuals are again centred around zero.

The full results of the fits are displayed in Tables B.1 and B.2. For all 86 datasets, the posterior distribution of each parameter as well as the residuals of the fit follow a normal distribution. This indicates that the sampling successfully converged.

Using Eq. (3) and the fitted A_atm parameter, we derived the α_refl sin i (R_p /R_J )² parameter, obtaining a minimum value of 0.0 8 and a maximum value of 21.01 (see Table B.1). Provided that both α_refl and i are known, one can then infer the radius of the object considered. Through Eq. (5) and (6), and the fitted parameters, A_ellip and A_boost , we retrieved the (M_p sin² i/M_J) and (M_p sin i/M_J) parameters, respectively. For the 77 targets analysed in the first approach, the orbital inclination is unknown and cannot be derived from the fits. However, an upper limit can be estimated by considering the geometric condition required for a transit to occur. For a circular orbit, a transit event occurs for i ≥ i_max with $i\in \left[0;\frac{\pi }{2}\right]$, $$i_{\max }={\cos }^{-1}\left(\frac{R_{\star }}{a}\right),$$(9)

where a is the semi-major axis of the orbit. Assuming that all of our non-transiting companions follow a circular orbit (Mazeh 2008), they would then have an inclination inferior to i_max. To ensure that the inclination value is in compliance with the fit we inferred an inclination value using the ratio of (M_p sin² i/M_J) and (M_p sin i/M_J ), for the nine targets where A_ellip and A_boost are accessible. For three hosts (KIC 11668095, KIC 5440365, and (KIC 8844761), the uncertainties of the parameters (M_p sin² i/M_J) and (M_p sin i/M_J) are too large to infer sin i values. Regarding KIC 4373708, KIC 11702835, and KIC 5622796, the sin i values are inconsistent, suggesting that these stars are probably out of the range of validity of the model and would require more complex models (e.g. Penoyre & Sandford 2019). Nevertheless, we were able to infer the inclination angle for three systems. For KIC 11550689, we obtained an inclination of 6°, for KIC 5697777, 0.4°, and for KIC 4264634, 38°. All of these values are lower than the maximum inclination angles for these systems (see Table B.2). The corresponding masses that we can therefore derive in these three cases are 32, 1412, and 21 M_J, respectively. In the cases of KIC 11550689 and KIC 4264634, if the numerous assumptions made in this modelling approach are correct, this would suggest that the non-transiting objects might be brown dwarfs. In the case of KIC 5697777, the inferred mass is approximately 1.3 M_⊙, which is obviously inconsistent with the hypothesis of a non-transiting substellar companion. Again, this might be because the model we apply is not valid and the inclination we obtain from this inference is not correct.

Fig. 3 Phase-folded light curve of KIC 4930913, together with the best-fit model in orange. The red points correspond to data binned in 1-day intervals. The residuals show the difference between the model and the data points and are free of large systematic trends.

Fig. 4 Phase-folded light curve of KIC 4264634, together with the best-fit model in orange. The red points correspond to data binned in 1-day intervals. The scatter of the residuals is higher due to a larger dispersion of the data points but remains free of large systematic trends as well.

4 Discussion and conclusion

Among our 88 candidates, we have 13 F-type, 23 G-type, 37 K-type, and 15 M-type stars exhibiting a signature consistent with the presence of a non-transiting planetary companion in a close orbit. The properties of each candidate are presented in Table C.1. Among these candidates, we detected three systems for which the second harmonic is the largest, indicating that tidal ellipsoidal distortion outweighs the two other phenomena (see Appendix D).

About 62% of our sample (55 stars) exhibit orbital periods of less than a day, confirming the ability of this method to detect very short-period close-in substellar companions. In comparison, we recall again that approximately only 3% of planets detected with the transit method exhibit orbital periods shorter than 1 day.

Within our set of 88 candidates, two stars have been previously confirmed to harbour at least one planet, with an orbital period different from the one detected in this study (they would have been removed from our sample otherwise; see Sect. 2.2). Details of these two candidates are provided in the Appendix E. The method presented in this work could help increase the number of multiplanetary systems detected, which in turn can enable the scientific community to gain greater insights into the formation and evolution of exoplanetary systems.

McQuillan et al. (2013, hereafter M13) identified a depletion zone for planets close to fast rotators. García et al. (2023, hereafter G23) demonstrate that a physically motivated model that includes tidal and magnetic effects is able to reproduce the depleted region. Indeed, the amplitude of tidal effects increases dramatically as the semi-major axis gets smaller, meaning that planets in very short orbits decay very fast towards their host star when in an unstable configuration (Hut 1980). In particular, the main component of the tidal potential scales as 1/a³ (e.g. Mathis 2018), and under some conditions the orbital decay timescale is proportional to $a^{\frac{13}{2}}$ (Collier Cameron & Jardine 2018).

Fig. 5 displays the stellar rotation period P_rot as a function of the orbital period P_orb for our sample, linked to the findings of G23 in grey. Our 88 candidates are shown as a function of the effective temperature, T_eff, sourced from Berger et al. (2020). The dashed black line corresponds to the fit of the lower envelope of points obtained by M13. In this context, we demonstrate that our candidates, if confirmed, would be located well below this limit, thereby populating a previously deserted area in a region where we expect exoplanets to experience orbital decay on relatively short timescales. We remind the reader that the method used in this work strongly favours the detection of planets on very short orbits. In this study, we have focused on very short orbital periods below 2.3 days. If we had been able to detect non-transiting substellar companions with frequencies exceeding our limit of 2.3 days, the region above the dashed black line would contain a higher density of candidates.

Since the orbital evolution of close-in systems takes place over short timescales, studying non-transiting exoplanets would provide us with additional systems to refine our understanding of star-planet interactions. Given the scarcity of available data on planets in short orbits around fast rotators, such a discovery would be a major asset for improving the calibration of models computing the evolution of the architecture of star-planet systems (e.g. Benbakoura et al. 2019; Ahuir et al. 2021). In addition, by fitting phase curve models to the light curves, we have been able to extract additional constraints on the mass and radius for 86 of the detected objects. Subsequent spectroscopic ground-based follow-up is necessary to further characterise the sample identified in this work, in order to confirm the planetary nature of the candidates identified in our sample. Radial velocity observation campaigns could, for example, provide estimates of the masses of the non-transiting companions detected in this work. Such detections would be very useful to refine the set of assumptions that are used to model the signal of non-transiting planets in photometric light curves. In the future, a similar search for non-transiting objects could be done with the K2 (Howell et al. 2014), Transiting Exoplanet Survey Satellite (TESS, Ricker et al. 2015), and PLAnetary Transits and Oscillations of Stars (PLATO, Rauer et al. 2024) missions. Notably, the fact that PLATO will focus on bright target observations (Montalto et al. 2021) will greatly expand the number of systems accessible to ground-based follow-up, allowing for the detailed characterisation of a much larger sample of exoplanetary systems.

Fig. 5 Prot as a function of Porb for the 88 non-transiting-exoplanet candidate systems. The size of the dots is a function of the Kepler apparent magnitude of the host star, sourced from Brown et al. (2011), whereas the temperature Teff, sourced from Berger et al. (2020), is colour-coded. The grey sample consists of 576 confirmed single-planet-host main-sequence solar-like stars with a rotation period obtained by García et al. (2023). KIC 5697777 is represented by a diamond. The dashed black line corresponds to the fit of the lower envelope of points obtained by McQuillan et al. (2013).

Data availability

Full Tables B.1 and C.1 are available at the CDS via anonymous ftp to cdsarc.cds.unistra.fr (130.79.128.5) or via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/699/A226

Acknowledgements

C.G., S.N.B., and A.F.L. would like to thank Dr. I. Pagano for enabling C.G.’s research internship at INAF-OACT during Spring 2023. The authors acknowledge the critical reading and the constructive comments from an anonymous referee that significantly allowed improving the original version of this paper. A.S. and C.G. acknowledge funding from the European Research Council project ExoMagnets (grant agreement no. 101125367). A.S. acknowledges funding from the Origins project DynamEarths. A.S., R.G. and S.M. acknowledge funding from the Programme National de Planétologie (INSU/PNP). C.G, S.N.B., and A.F.L acknowledge support from PLATO ASI-INAF agreement no. 2022-28-HH.0 “PLATO Fase D”. R.A.G. acknowledges support from PLATO and SoHO/GOLF CNES grants. A.R.G.S acknowledges the support from the FCT through national funds and FEDER through COMPETE2020 (UIDB/04434/2020, UIDP/04434/2020 & 2022.03993.PTDC) and the support from the FCT through the work contract No. 2020.02480.CEECIND/CP1631/CT0001. S.M. acknowledges support by the Spanish Ministry of Science and Innovation with the grant no. PID2019-107061GB-C66 and through AEI under the Severo Ochoa Centres of Excellence Programme 2020–2023 (CEX2019-000920-S). This paper includes data collected by the Kepler mission, and obtained from the MAST data archive at the Space Telescope Science Institute (STScI). Funding for the Kepler mission is provided by the NASA Science Mission Directorate. STScI is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5–26555. This research has made use of the NASA Exoplanet Archive, which is operated by the California Institute of Technology, under contract with the National Aeronautics and Space Administration under the Exoplanet Exploration Program. Software: matplotlib (Hunter 2007), numpy (Harris et al. 2020), scipy (Virtanen et al. 2020), pandas (The pandas development team 2020), pymc (Wiecki et al. 2023), star-privateer (Breton et al. 2024). ultranest (Buchner 2021), corner (Foreman-Mackey 2016), astropy (Astropy Collaboration 2022), apollinaire (Breton et al. 2022).

Appendix A Contaminants

A.1 Transiting exoplanet

In Sect. 2.2, we were able to remove 4 stars confirmed to har-bour transiting exoplanets. For example, Morton et al. (2016) discovered Kepler-685b, a gas giant exoplanet orbiting around KIC 3351888 with an orbital period of 1.62 days.

Fig. A.1 (a) PSD of KIC 3351888. The Prot, shown in the orange area, was measured by Santos et al. (2021). The green area shows the high-amplitude peak signal detected in this work, at around 1.62 days. (b) Phase-folding analysis on the light curve of KIC 3351888. The red curve is a smoothing of the black curve, obtained with a triangular filter.

Panel (a) of Fig. A.1 shows the PSD of KIC 3351888 where a high-amplitude peak signal at 1.62 days can be observed. In panel b of Fig A.1, the phase-folding curve demonstrates a flux variation analogous to that observed during a transit. The signal detected in this work for KIC 3351888 therefore corresponds to the planet Kepler-685b. Because it has already been detected, this candidate is not part of our final sample.

A.2 Phase-folding

By phase-folding the light curve, we can examine the temporal stability of a signal corresponding to a given period. To refine our sample, we performed a visual inspection of the 283 stars obtained after Sects. 2.1 and 2.2. In these sections, we searched for high-amplitude peaks in the PSDs and compared our results with existing catalogues. In Sect. 2.3, we performed a phase-folding analysis followed by a visual inspection.

Fig. A.2 shows the phase-folded light curve of KIC 11154093, which exhibits a complex pattern with multiple distinct flux variations across the folding period. Visually, the observed pattern deviates significantly from the expected behaviour of a signal (see Fig. 1 and panel c of Fig. 2). This complex behaviour might be linked to the presence of multiple stars close to this target, such as seen when looking at the 2MASS image (see Appendix A.4) of the target in Fig. A.3. As a result, this target was excluded from our sample.

Fig. A.2 Phase-folding analysis on the light curve of KIC 11154093. The red curve is a smoothing of the black curve, obtained with a trian-gular filter.

Fig. A.3 Image from the KASOC website of KIC 11154093 using 2MASS (Skrutskie et al. 2006). The purple cross marks KIC 11154093.

A.3 Wavelet analysis

In Sect. 2.3, we computed the WPS of each of the 245 stellar light curves, and we only kept stars that present stable modulations along time, as shown in panel d of Fig. 2.

Figure A.4 presents the rotational analysis of KIC 2852669, a K-type star included in our sample of 245 stars. The top panel displays the PSD while the bottom panel shows the normalised time-averaged projection of the wavelet transform.

The different analyses all point to an orbital period of around 0.6 days (panel a, Fig. A.4). However, the signal shown on panel b of Fig. A.4 does not show stable modulations along time. We recall that the focal plane of Kepler is rotated after each quarter, so the photometric apertures used to monitor a given target have varying levels of background contamination. The periodic signal found in this study is therefore unlikely to be caused by a non-transiting object, and KIC 3233612 was excluded from our sample.

Fig. A.4 Rotational analysis of KIC 3233612. (a) PSD of the star. The stellar rotation period is shown in the orange area, while the green area shows the high-amplitude peak at around 0.6 days. (b) Period-time analysis of the stellar light curve. High-amplitude periodicity is highlighted by orange and red, whereas yellow and blue denote low-power amplitude. The yellow discontinuous structure appearing around 0.6 days results from background contamination. The shaded area represents the cone of influence of the reliable periodicities.

A.4 Field contaminants

In Sect. 2.3, targets with a nearby bright object in their aperture are removed from the sample. We first looked for potential con-taminants using the Gaia catalogue (Gaia Collaboration 2023) and found 70 Gaia objects near 38 of our targets. We excluded 11 targets for which the potential contaminants are located within a few arcseconds. The 27 remaining stars are located at least 10” away from the potential contaminant, making significant contaminations problems unlikely. The 10 closest contaminants, associated with these 27 retained targets, are listed in Table A.1.

Fig. A.5 Image from the KASOC website of KIC 12505054 using 2MASS (Skrutskie et al. 2006). The purple cross marks KIC 12505054.

We then carried out a visual inspection using the 2MASS images (Skrutskie et al. 2006), which can be accessed through the KASOC database (Kjeldsen et al. 2010). Similarly to what was done e.g. in Noraz et al. (2022), we excluded candidates if the star is within the halo of a brighter object.

Figure A.5 shows KIC 12505054, a star located within the halo of another bright object not listed in the KASOC database. As a result, this star was excluded from our final sample. Indeed, a nearby bright object is likely to contaminate the signal obtained for this star.

Appendix B Summary of photometric fit results

Fig. B.1 Corner plot for the fit of KIC 4930913 showing the posterior distributions of the fitted parameters A0, Aatm and T0. The values and error bars represent the median and 1σ uncertainty based on the posterior distributions.

In Sect. 3 we computed the fit of the photometric modulation data and retrieved the theoretical values of A_atm. By doing so, we were able to deduce the parameter α_refl sin i(R_p /R_J)², thanks to Eq. (3). These resulting parameters are presented in Table B.1 for the 15 first targets. Fig. B.1 shows the corner plot of the posterior probability distribution resulting from the fit of KIC 4930913, a star introduced in Sect. 3 to illustrate the fitting of light curves exhibiting approximately sinusoidal modulation. The diagonal panels display the marginalised one dimensional distribution for each parameter: the baseline amplitude A₀, the atmospherical amplitude A_atm and the reference time T₀. The off-diagonal panels present the two-dimensional joint posterior distributions, with contour levels corresponding to the 1σ, 2σ, and 3σ credible regions. The results indicate that all parameters are well constrained. In particular, T₀ shows a narrow distribution, implying high precision in the timing estimate. The distribution of A₀ is broader and slightly asymmetric, reflecting larger uncertainties. No significant correlations are observed between the parameters, suggesting they are independently constrained by the data. We also calculated with Eq. (9) the maximum inclination possible for each system, assuming they all follow a circular orbit. For two stars (KIC 3124224 and KIC 10396339), we do not provide i_max, because, using R_⋆ and M_⋆ from the reference catalogues adopted, we obtain a < R_⋆, which is obviously unphysical. Nevertheless, considering that R_⋆ and M_⋆ might be incorrect, we decided to keep these two stars in our sample.

Fig. B.2 Corner plot for the fit of KIC 4264634 showing the posterior distributions of the fitted parameters A0, Aatm, Aellip, Aboost and T0. The values and error bars represent the median and 1σ uncertainty based on the posterior distributions.

When the modulation of the phase-folded light curve was more complex, we performed a fit considering the three components of the signal. The corner plot in Fig. B.2 summarises the posterior probability distributions resulting for KIC 4264634, the example star chosen in Sect. 3 to illustrate the fit which accounts for the three components of the signal. As in Fig. B.1, the diagonal panels show the marginalised one-dimensional distributions for each fitted parameter: the baseline amplitude, A₀, the atmospheric amplitude, A_atm, the reference time, T₀, the ellipsoidal amplitude, A_ellip, and the Doppler boosting amplitude, A_boost. The off-diagonal panels illustrate the two-dimensional joint posterior distributions, with contours corresponding to the 1σ, 2σ, and 3σ credible intervals. All parameters are reasonably well constrained, with the amplitude parameters displaying broader distributions, especially A₀ and A_atm, which exhibit mild asymmetries. The joint distributions suggest weak correlations between some of the amplitude parameters, particularly between A_atm and A_ellip. Overall, the parameters remain largely uncorrelated. This supports the robustness of the fit and the ability of the model to independently constrain each physical contribution. The resulting parameters for all of the nine targets studied here are presented in Table B.2. In both Tables B.1 and B.2, the reported values of A_atm, A_beam and A_ellip correspond to their absolute values. Thanks to Eq. (5) and (6) we were also able to determine the parameters (M_p sin² i/M_J ) and (M_p sin i/M_J). With such calculations, we can in principle retrieve the value of inclination of the non-transiting object i. This was possible for three of our systems (KIC 11550689, 5697777, and 4264634). For three other stars, the uncertainties on our calculated parameters are too large to infer on the value i. For the last three stars, the uncertainties are not sufficiently high enough to explain the obtained inconsistent results. We are therefore led to believe that these stars are out of the validity range of the model used in this article and would probably require more complex models (see Penoyre & Sandford 2019).

Appendix C Final sample

We present in Table C.1 the overall properties of 16 out of our 88 selected candidates likely to host non-transiting companions. The orbital period and amplitude are derived from this work. The effective temperature T_eff, surface gravity log g, metallicity and mass measurements parameters are primarily sourced from Berger et al. (2020) while, for stars not included in this catalogue, the properties are taken from Mathur et al. (2017). The rotation periods are taken from Santos et al. (2019, 2021) and the Kepler apparent magnitudes Kpmag come from Brown et al. (2011). The semi-major axis is calculated with the Kepler’s third law, assuming the planetary mass is negligible compared to the stellar mass.

The flag column contains a sequence of three binary values (0 or 1) for each entry. Each digit represents a separate flag. If the first value is set to 1, it indicates that tidal forces outweighs atmospheric processes and Doppler boosting in the orbital signal. The second value, if set to 1, means that the target’s orbital period is a harmonic of the stellar rotation period (see Appendix D). Finally, the third value indicates whether an exoplanet has already been found orbiting around the targeted star, but at a different orbital period than the one detected in this work (see Appendix E). In such cases, the value is set to 1. The most common set of flags is 0, 0,0, indicating that the majority of the sample does not exhibit these characteristics.

We estimate the uncertainties on the orbital period P_orb and the periodogram amplitude A_periodogram of the modulation by modelling the main peak attributed in the PSD to the non-transiting planet with a Gaussian profile with central frequency ν₀ = 1/P_orb, height H, and standard deviation σ. Using the pyMC software (Wiecki et al. 2023), we perform a Hamiltonian Monte-Carlo exploration of the distribution of the set θ = {P_orb, H, σ}, assuming a χ² with two-degrees of freedom as likelihood function for the PSD (Woodard 1984). The periodogram amplitude A_periodogram of the modulation, in ppm, is computed as $A_{\text{periodogram }}^2=2\sqrt{\pi }\sigma H$. The parameter estimate is taken as the median of the distribution. The lower uncertainty is computed as the difference between the median value of the distribution and the 16th percentile, while the upper uncertainty is computed as the difference between the median value of the distribution and the 84th percentile.

Appendix D Systems with an orbital period corresponding to a harmonic of the stellar rotation period

In our 88 stars catalogue, the orbital period of 3 candidates corresponds to a harmonic of the stellar surface rotation period.

Fig. D.1 (a) PSD of KIC 7801559. The Prot, shown in the orange area, was measured by Santos et al. (2021). The green area shows the high-amplitude signal peak at around 1.34 days. (b) Phase-folding analysis on the light curve of KIC 7801559. The red curve is a smoothing of the black curve, obtained with a triangular filter.

Panel (a) of Fig. D.1 presents the PSD of KIC 7801559, an F-type star with a P_rot of 2.8 days. The companion’s orbital period detected in this work corresponds to the second harmonic of the stellar rotation period. On panel (a), the high-amplitude peak corresponding to the orbital signal (in the green area), seems embedded in the rotational harmonic pattern. Nevertheless, the signal shows a very stable periodic behaviour, distinct from the quasi-periodic signature observed in the case of spot modulation, both in the PSD and the light curve phase-folding in panel (b) of Fig. D.1. Therefore, we choose to keep this candidate in our catalogue and suggest that this could be the signature of systems in a spin-orbit resonance. This verification was also done for KIC 5560359 and KIC 8193786. The companion orbital period and the stellar rotation period of these candidates are shown in Table D.

Appendix E Systems with a confirmed exoplanet

In our 88 stars sample, we detected 2 stars with previously confirmed exoplanets, but at an orbital period different from the one detected in this study. The first star is KIC 10027247, an M-type star with a super-Earth exoplanet, named Kepler-1229 b (Morton et al. 2016). Its mass is 2.54 M_⊕, and it takes 86.8 days to complete one orbit of its star at a distance of 0.3006 Astronomical Unit (AU). Notably, its orbit is within the habitable zone (Torres et al. 2017). The second star in our catalogue that already hosts a confirmed exoplanet is KIC 9717943, an F-type star. Kepler-1169 b is a terrestrial exoplanet with a mass of 0.779 M_⊕ (Morton et al. 2016). Its orbital period is 6.11 days, and it is located at 0.0709 AU from its host star. The properties of these two systems are shown in Table E.

References

All Tables

All Figures

Fig. 1 Simulated contributions of the three modulations to the phase-curve amplitude according to the phase angle. The modulation linked to atmospheric processes is shown in black and includes reflection and emission modulations. The ellipsoidal distortion is shown in purple, while the Doppler boosting is in orange. The overall contribution is shown by the composite phase-curve in yellow. Here, an eclipse would occur at the inferior conjunction, while the transit would happen at the superior conjunction with the observer.

In the text

Fig. 2 Flux analysis of KIC 5697777. (a) Light curve calibrated according to the KEPSEISMIC method and filtered at 20 days. (b) Power spectral density of the light curve. The rotation period of the star is taken from Santos et al. (2021) and is located within the orange area. The green area shows the high-amplitude signal peak at around 0.5 days. The second high-amplitude signal peak located at around 0.25 days is the second harmonic of the orbital signal. (c) Phase-folding of the light curve with an additional filter to remove the stellar rotational contribution. The red curve is a smoothing of the black curve, obtained with a triangular filter. (d) Wavelet power spectrum of the light curve. A high-amplitude periodicity is highlighted by orange and red, whereas yellow and blue denote a low-power amplitude. The yellow bands with a constant period throughout time correspond to the high-amplitude peak and its second harmonic, identified in panel b. The reddish band corresponds to the rotation period and its harmonics, highlighted in orange in panel b. The shaded area represents the cone of influence of the reliable periodicities.

In the text

	Fig. 3 Phase-folded light curve of KIC 4930913, together with the best-fit model in orange. The red points correspond to data binned in 1-day intervals. The residuals show the difference between the model and the data points and are free of large systematic trends.
In the text

	Fig. 4 Phase-folded light curve of KIC 4264634, together with the best-fit model in orange. The red points correspond to data binned in 1-day intervals. The scatter of the residuals is higher due to a larger dispersion of the data points but remains free of large systematic trends as well.
In the text

Fig. 5 Prot as a function of Porb for the 88 non-transiting-exoplanet candidate systems. The size of the dots is a function of the Kepler apparent magnitude of the host star, sourced from Brown et al. (2011), whereas the temperature Teff, sourced from Berger et al. (2020), is colour-coded. The grey sample consists of 576 confirmed single-planet-host main-sequence solar-like stars with a rotation period obtained by García et al. (2023). KIC 5697777 is represented by a diamond. The dashed black line corresponds to the fit of the lower envelope of points obtained by McQuillan et al. (2013).

In the text

	Fig. A.1 (a) PSD of KIC 3351888. The Prot, shown in the orange area, was measured by Santos et al. (2021). The green area shows the high-amplitude peak signal detected in this work, at around 1.62 days. (b) Phase-folding analysis on the light curve of KIC 3351888. The red curve is a smoothing of the black curve, obtained with a triangular filter.
In the text

	Fig. A.2 Phase-folding analysis on the light curve of KIC 11154093. The red curve is a smoothing of the black curve, obtained with a trian-gular filter.
In the text

	Fig. A.3 Image from the KASOC website of KIC 11154093 using 2MASS (Skrutskie et al. 2006). The purple cross marks KIC 11154093.
In the text

Fig. A.4 Rotational analysis of KIC 3233612. (a) PSD of the star. The stellar rotation period is shown in the orange area, while the green area shows the high-amplitude peak at around 0.6 days. (b) Period-time analysis of the stellar light curve. High-amplitude periodicity is highlighted by orange and red, whereas yellow and blue denote low-power amplitude. The yellow discontinuous structure appearing around 0.6 days results from background contamination. The shaded area represents the cone of influence of the reliable periodicities.

In the text

	Fig. A.5 Image from the KASOC website of KIC 12505054 using 2MASS (Skrutskie et al. 2006). The purple cross marks KIC 12505054.
In the text

	Fig. B.1 Corner plot for the fit of KIC 4930913 showing the posterior distributions of the fitted parameters A0, Aatm and T0. The values and error bars represent the median and 1σ uncertainty based on the posterior distributions.
In the text

	Fig. B.2 Corner plot for the fit of KIC 4264634 showing the posterior distributions of the fitted parameters A0, Aatm, Aellip, Aboost and T0. The values and error bars represent the median and 1σ uncertainty based on the posterior distributions.
In the text

	Fig. D.1 (a) PSD of KIC 7801559. The Prot, shown in the orange area, was measured by Santos et al. (2021). The green area shows the high-amplitude signal peak at around 1.34 days. (b) Phase-folding analysis on the light curve of KIC 7801559. The red curve is a smoothing of the black curve, obtained with a triangular filter.
In the text