© The Authors 2025

1 Introduction

Following the discovery of the first confirmed extrasolar planet orbiting a main sequence star (Mayor & Queloz 1995), close to 6000 exoplanets have been identified¹. One of the key points of these past decades is the fact that the Solar system is not unique in many aspects: it is simply one of the ever-growing number of known planetary systems in our Galaxy. Based on our understanding of the Solar System, we expect that there are perhaps orders of magnitude more satellites orbiting exoplanets compared to the number of exoplanets. Yet at the time of writing, there are no confirmed exomoons, but there are a number of possible explanations for this (see Szabó et al. 2024, and references therein). With the introduction of space telescopes such as Kepler (Borucki et al. 2010) and CHaracterising ExOPlanets Satellite (CHEOPS; Benz et al. 2021), we have the theoretical capability of finding transiting exomoons below the size of Earth.

The Kepler space telescope was the first instrument recognised to be able to detect an exomoon (Szabó et al. 2006; Simon et al. 2007). Kipping et al. (2015) and Kipping et al. (2022) performed extensive analyses of Kepler light curves with the aim of finding transiting exomoons, yielding Kepler-1708b-i as a candidate. Another method thought to be useful for detection of satellites is based on the transit timing tariations (TTVs) of the planets (Simon et al. 2007; Kipping 2009) caused by exomoons. Fox & Wiegert (2021) identified six candidates based on their Kepler datasets, which were found to be false positives by Kipping (2020). Kipping & Yahalomi (2023) found no strong TTV signals that could be linked to satellites in a large set of suitable planets observed by Kepler. Simon et al. (2015) predicted that Earth-sized exomoons could likely be detectable with CHEOPS. For that reason, a considerable effort has been devoted within the CHEOPS Science Team to observing systems that are likely to sustain exomoons for long periods of time (see e.g. Dobos et al. 2022) in pursuit of detecting these objects. Ehrenreich et al. (2023) found that there are likely no transiting bodies above the size of Mars orbiting ν² Lupi d, based on 1.5 transits observed with CHEOPS.

In this paper, we present the necessary steps in the search for moons orbiting the bright (V = 8.62 mag) star HD 95338, including refining the stellar and planetary parameters, analysing synthetic data, and placing upper limits on the hypothesised satellite. Díaz et al. (2020) reported the discovery of a dense Neptune-sized planet (${3.89}_{-0.20}^{+0.19}$ R_⊕, ${42.44}_{-2.08}^{+2.22}$ M_⊕, yielding a density of ${3.98}_{-0.64}^{+0.62}$ g cm⁻³), HD 95338b. Based on a single transit observed by the Transiting Exoplanet Survey Satellite (TESS; Ricker et al. 2015) and a large number of radial velocity measurements, Díaz et al. (2020) found an orbital period of 55.087 ± 0.020 days. This meant that the planet could be far enough from its host star that an exomoon could survive for several gigayears (Dobos et al. 2022), yet close enough that it can be observed multiple times during the life cycle of CHEOPS, leading to tighter constraints on its parameters, as well as on those of its hypothetical companion. However, due to the Sun-synchronous, dusk-dawn low Earth orbit of the satellite, long-period planets might only be observable once a year.

This paper is structured as follows. In Sect. 2, we describe in detail the photometry of CHEOPS and TESS observations, with a special emphasis on the various detrending steps needed to account for the known instrumental effects. We also provide details of the modelling of the transit of HD 95338b. New and improved stellar parameters are provided in Sect. 3. In Sect. 4, we present the refined planetary parameters, based on five transits observed with CHEOPS and TESS. In Sect. 5, we present a novel method for chasing transits linked to exomoons. We argue for the need of a thorough noise treatment to avoid false positive exomoon detections. We exclude the presence of Mars-sized transiting exomoons at 1σ. In Sect. 5.3, we present the analyses of synthetic data, which are used to justify the algorithm utilised in the search for additional transits. We show that even from a single transit, the 3σ detection of a less than 1 R_⊕ transit would be feasible in this system.

2 Methods

2.1 CHEOPS photometry

CHEOPS observed four transits of HD 95338b (during four separate visits) in March 2021, April 2022, May 2023, and April 2025 (Table A.1). The spacecraft is on a low Earth orbit, with an ≈100 min period. Consequently, targets cannot be observed continuously throughout an orbit due to Earth occultations, crossing over the South Atlantic Anomaly, and too high stray-light contributions, as is often discussed with regards to CHEOPS data (e.g. Lacedelli et al. 2022; Krenn et al. 2024). This is reflected in the so-called observational efficiency (Table A.1). We extracted light curves using the CHEOPS Data Reduction Pipeline (DRP; version 14.1.2; Hoyer et al. 2020). The DRP performs the bias, read-out noise, flatfield, and smearing-corrections, while also estimating the telescope jitter and background contamination. The pipeline also carries out aperture photometry using a circular apertures, with radii between 15 and 40 pixels (in 1-pixel steps). In this work, we rely on light curves extracted using the R = 25 px radius aperture (the so-called default mode).

Due to the spacecraft rolling around its line of sight, the well-known triangular point spread function (as discussed in e.g. Szabó et al. 2021; Kiefer et al. 2023) of any given star (Fig. 1) falls onto different pixels throughout each orbit (causing a part of the flux to fall in and out of the circular aperture as the telescope rolls). Several factors, including the inhomogeneous nature of the CCD, stray light, smearing trails, and thermal variability of the instrument result in the appearance of periodic or quasi-periodic systematic effects that are correlated with the angular orientation (and possibly the jitter) of the spacecraft. Additionally, the flux contribution of nearby stars also introduces a noise source that changes throughout each orbit and is therefore correlated in nature. The background variation thus also has to be estimated for every exposure.

We emphasise, that in the context of nominal science observations, a subset of the full frame images, known as ‘window images’, are downlinked. In order to minimise the noise caused by bad pixels, it is necessary to move the position of the window in the CCD to different regions of the CCD when required. In May 2024, the position of the window was relocated to the central location of the CCD². This resulted in the fourth visit being observed in a distinct region of the CCD compared to the initial three observations. This modification has resulted in a shift in the configuration of the star pattern within the full frame. This alteration has had an impact on the correlation between the flux and the roll angle of the telescopes (and the number of used harmonics of the roll angle in the analysis). This is due to the fact that the smearing trails of the background stars within the full frame intersect the window as the field rotates around the target star.

Fig. 1 The Subarray image depicts the point spread function (PSF) of HD 95338b and the neighbouring stars, with their respective magnitudes of brightness indicated. The target is centred on the coordinates (180, 178) on the CCD.

2.2 TESS photometry

In addition to the CHEOPS photometry described above, we also obtained the short-cadence SAP (Simple Aperture Photometry) light curves from sectors 10 and 63 of TESS (Transiting Exoplanet Survey Satellite; Ricker et al. 2015). These include one transit of HD 95338b each. The transit from sector 10 was also analysed by Díaz et al. (2020). We corrected for contamination in the aperture using the CROWDSAP keyword. We downloaded the light curves from the two sectors from MAST. Additionally, we corrected for the pointing instability of the telescope using the pipeline-estimated position of the photocenter. As a final preparatory step, we cut out a 1.5 and a 3-day window (in sectors 10 and 63, respectively) of the light curves – that are centered on the two transits – in order to minimise the computational time needed for the analysis described below.

2.3 Modelling the transit of HD 95338b from CHEOPS

There are a number of algorithms and software packages capable of modelling a transiting exoplanet and their satellites, including planetplanet (Luger et al. 2017), the phase-folding method by Kipping (2021), the analytical framework of Saha & Sengupta (2022), and Pandora (Hippke & Heller 2022). There are, however, none that incorporate the advanced noise filtering algorithms such as those based on Gaussian processes (GPs; e.g. Rasmussen & Williams 2006; Roberts et al. 2012; Ebden 2015) or a wavelet-based technique (e.g. Carter & Winn 2009). These tools are used to deal with the time-correlated noise arising from both astrophysical and instrumental noise sources. Their benefits in avoiding distortions of the light curve parameters of transiting exoplanets are clearly demonstrated in Barros et al. (2020), Csizmadia et al. (2023), and Kálmán et al. (2023b). Regarding the detection of exomoon transits specifically, Kálmán et al. (2024a) and Szabó et al. (2024) discuss the implications of insufficient noise treatment.

The two perhaps most well-known exomoon candidates, Kepler-1625b-i (Teachey & Kipping 2018) and Kepler-1708b-i (Kipping et al. 2022), identified by a so-called photodynamical method, where the parameters of the planet and its satellite are fitted jointly, are not confirmed through independent investigations. In fact, Heller et al. (2019) and Kreidberg et al. (2019) provide alternative explanation for the exomoon signal observed in Kepler-1625b-i, involving the presence of red noise. More recently, Heller & Hippke (2024) presented evidence suggesting that ‘large exomoons are unlikely’ around both Kepler-1625b and Kepler-1708b, although this is again refuted (Kipping et al. 2025). It is reasonable to expect that a simultaneous fitting of the planetary and lunar signals with the time-correlated noise would be beneficial in discerning false-positives from bona fide exomoons.

Lacking such a tool, we developed a sequential scheme to search for exomoon-like signals in the CHEOPS light curves of HD 95338b. In the first step, we model the transit of HD 95338b, then in the second step, we search the residuals for signs of additional transits. The transit signal that could be attributed to an exomoon has an amplitude on the order of ≈100 ppm (as is discussed in Section 5). Thus a thorough detrending of the known systematic effects is essential to avoid false-positive detections. We therefore constructed a strict framework to handle the known (quasi-periodic) systematic effects as follows:

• 0.

  Remove light curve points that are heavily affected by cosmic rays and remove four points before and after each gap (only for the fourth visit).

• 1.

  Conduct preliminary transit modelling of HD 95338b.

• 2.

  Find the best-fit model of the instrumental effects in the residuals of the previous step; find the >3σ outliers in the residuals of this step.

• 3.

  Subtract the best-fit systematics model from the raw light curve and mask the outliers identified in the previous step.

• 4.

  Model the transit of HD 95338b in the pre-whitened light curve.

• 5.

  Check the residuals for transits of an exomoon.

In step 2, we constructed a basis vector (F_BV) for detrending (following the common recipe, e.g. Szabó et al. 2021; Smith & Csizmadia 2022) in the form of $$F_{\mathrm{B}\mathrm{V},N}=\sum _{i=1}^N(a_i\sin (i\varphi )+b_i\cos (i\varphi )),$$(1)

where ϕ is the roll angle of the spacecraft, a_i and b_i are amplitudes of the sinusoidal terms. We then used the Bayesian Information Criterion (BIC Schwarz 1978) to find the optimal number of harmonics N, where N ∈ {1 . . . 10}, to be included in the detrending process. Defining the optimal set of detrending vectors is also part of the optimisation problem, and besides the roll angle harmonics, all telemetry data and by-products of the photometry reduction pipeline have initially been considered. We concluded that the background contamination is among the most important detrending parameters, and therefore, we included this vector into the noise model as well.

We analysed the transits of HD 95338b using the Transit and Light Curve Modeller (TLCM; Csizmadia 2020). TLCM uses a two-dimensional Gauss-Legendre integration approach to calculate the transit curves. These are described in terms of the scaled semi-major axis of the planet, a/R_⋆, the planet-to-star radius ratio, R_p/R_⋆, the time of mid-transit, T₀, and the orbital period, P. Díaz et al. (2020) found that the planet has an eccentric orbit, specifically, an eccentricity of e = 0.127 ± 0.045 and an argument of the periastron ω = 39.4° ± 18.7°. In TLCM, it is possible to account for the orbital eccentricity, where the ingress and egress durations can differ. In that case, instead of the commonly used impact parameter, the so-called (inferior) conjunction parameter contains the information of the orbital inclination of the planet (i). It is given by Csizmadia (2020) as: $$b^{\mathrm{\prime }}=\frac{a}{R_{\star }}\frac{(1-e^2)\cos i}{(1+e\sin \omega )}.$$(2)

We fixed e and ω at their reported values during the light curve modelling.

Stellar limb darkening is taken into account via the quadratic law (Wade & Rucinski 1985), described by the linear (u₁) and quadratic (u₂) coefficients using the re-parametrisation introduced in Kálmán et al. (2024b) and Harre et al. (2024): $$A=\frac{1}{4}(u_1(\frac{1}{\alpha }-\frac{1}{\beta })+u_2(\frac{1}{\alpha }+\frac{1}{\beta })),$$(3) $$B=\frac{1}{4}(u_1(\frac{1}{\alpha }+\frac{1}{\beta })+u_2(\frac{1}{\alpha }-\frac{1}{\beta })),$$(4)

where $\alpha =\frac{1}{2}\cos \left({77}^{\circ }\right)$, and $\beta =\frac{1}{2}\sin \left({77}^{\circ }\right)$. In step 1 of the algorithm, the A and B parameters were fitted using the [1.06, 1.56] and [1.40, 1.80] uniform priors, respectively. These priors are set according to the theoretical limb-darkening of the star (with stellar parameters from Table 1), as calculated by LDCU³. This software uses a modified version of the routines of Espinoza & Jordán (2015) to calculated the limb darkening coefficients, based on stellar intensity profiles calculated via the ATLAS (Kurucz 1979) and PHOENIX (Husser et al. 2013) models. In step 4, when the definitive planetary parameters are computed, the A and B coefficients were treated as free parameters of the fit.

Furthermore, TLCM uses the wavelet-based formalism of Carter & Winn (2009) to handle time-correlated noise (also commonly known as ‘red noise’) that comprises of any remaining instrumental effects or astrophysical variations (Csizmadia et al. 2023). This technique introduces two additional fitting parameters, σ_w for the standard deviation of the white noise component and σ_r for characterising the red component. These are fitted simultaneously with the parameters describing the transits. To avoid overcorrecting the light curves, the red noise fit is constrained by prescribing that the standard deviation of the white noise (i.e. the residuals) has to be equal to the mean photometric uncertainty of the input light curve (Csizmadia et al. 2023). In step 1, this technique is used to handle the instrumental noise sources, as the wavelet-based noise model is fitted jointly with the transits. Theoretical limb darkening coefficients are needed at this phase to retrieve plausible transit shapes. This allows for a clear detection of the instrument-related effects in step 2. The wavelet-based noise handling is used again in step 4 to obtain the final set of planetary parameters. At this stage, it is expected that astrophysical noise sources will make up the majority of the autocorrelated effects. We emphasise that, from the perspective of the planetary transit fitting, any additional transits – including those from an exomoon – are treated as noise by the wavelet technique or other non-parametric methods, and were therefore removed. For this reason, the search for moons needs to occur on a residual light curve where the time-correlated effects (that are not clearly attributed to instrumental effects) are not removed. Therefore, the complex detrending process outlined above is warranted.

A much less computationally efficient method would be to fit the F_BV jointly with the (planetary) transits and the wavelet-based noise model⁴. This would, however, involve the introduction of up to 66 additional fitting parameters, since the a_i and b_i coefficients of Eq. (1), as well as the amplitudes of the background contamination are independent variables in the four visits. In the steps outlined above, the systematic effects are fitted using a simple least squares method, which is several orders of magnitude faster than the joint F_BV plus red noise plus transit fit in our experience. This way, we are able to drastically narrow down the parameter space of the transit modelling itself, thus reducing the possibilities for the MCMC chains to get stuck in so-called local minima, and increasing the reliability of the solutions.

2.4 Joint TESS+CHEOPS modelling

In order to increase the precision of the recovered planetary parameters, we also performed a joint TESS+CHEOPS LC analysis. We assume that the passband-independent orbital parameters, including the eccentricity, the argument of periastron, the semi-major axis and the orbital period remain constant on the timescale relevant in these observations. We also assume that although the TESS passband is redder compared to CHEOPS (see e.g. Scandariato et al. 2024), the possible atmospheric features on HD 95338b are not detectable with the present noise levels, meaning that there is no wavelength-dependence of the transit depth. Therefore, we assumed that R_p/R_⋆ is the same in the two passbands. We further assumed that the transit chord does not change between the different observations, meaning that the scaled semi-major axis, the conjunction parameter, and the orbital period are the same for all five modelled transits. The limb-darkening coefficients are passband-dependent, which neccesitates the inclusion of A^CHEOPS, B^CHEOPS, A^TESS, and B^TESS into the analysis. It is understood that red noise from an astrophysical origin decreases with increasing wavelength. This manifests in an expectedly lower σ_r (see e.g. Kálmán et al. 2023a) in the TESS observations compared to the CHEOPS light curves. On the other hand, the white noise level depends primarily on the aperture size and the exposure time. We therefore include ${\sigma }_r^{\text{CHEOPS}},{\sigma }_w^{\text{CHEOPS}},{\sigma }_r^{\text{TESS}}$, and ${\sigma }_w^{\text{TESS}}$ in the global LC modelling. The likelihood of the joint TESS+CHEOPS modelling is then the sum of the log-likelihoods of the light curve fits of the data from the two satellites.

In addition to increasing the number of transits that can be analysed, the two different passbands also allowed us to break the degeneracy between the limb-darkening coefficients and the transit parameters, yielding more precise estimates on both.

During the revision of this paper, Saha (2025) published an independent analysis based on the first three transits from CHEOPS and the two TESS transits.

2.5 Sampling the posterior distributions

In every TLCM-based light curve analysis described in this paper, we followed the same recipe. We used wide, uninformed priors on every parameter. A genetic algorithm is then able to find the (global) minimum in the likelihood space, with N_init = 2000 randomly selected values for each parameter (see Csizmadia 2020, for more details). Having found starting points that are suspected to be good, we applied a Differential Evolution Markov-chain Monte Carlo (DE-MCMC; Ter Braak 2006) to get robust estimates of the uncertainty ranges of the fitting parameters. In every case, we used ten chains (also known as walkers) each consisting of 20 000 possible steps. The convergence was monitored via the Gelman–Rubin statistic (Gelman & Rubin 1992) and the effective sample size (e.g. Ripley 1987). If convergence is not reached within these 20 000 steps, the chain length is automatically extended to 10 · 20 000.

Fig. 2 Demonstration of the fitting and detrending steps. The raw CHEOPS light curves of HD 95338b are shown in the top row, along with the initially fitted transit curves (step 1). The residuals of the initial transit fit, along with the best-fit systematics model is shown in the second row, with N number of harmonics of the roll angle included. The final transit fit (step 4) and its residuals (later used in the search for additional transits) are shown in the bottom two rows, where the full white circles represent 200-point bins (with an average uncertainty of ≈20 ppm in each bin).

3 Stellar parameters

We used the ARES+MOOG spectroscopic analysis methodology to derive the stellar spectroscopic parameters (T_eff, log g, micro-turbulence, [Fe/H]). This methodology is described in detail in Sousa et al. (2021); Sousa (2014); Santos et al. (2013). The equivalent widths (EW) were consistently measured using the ARES code⁵ (Sousa et al. 2007, 2015). We used the list of iron lines appropriate for stars cooler that 5200 K which was presented in Tsantaki et al. (2013). For this spectral analysis we used a combined HARPS spectrum of HD 95338. To converge for the best set of spectroscopic parameters for each spectrum we use a minimisation process to find the ionisation and excitation equilibrium. This process makes use of a grid of Kurucz model atmospheres (Kurucz 1993) and the latest version of the radiative transfer code MOOG (Sneden 1973). We also derived a more accurate trigonometric surface gravity using recent Gaia data following the same procedure as described in Sousa et al. (2021) which provided a consistent value (4.44 ± 0.04 dex) when compared with the spectroscopic surface gravity.

We computed the stellar radius of HD 95338 utilising a MCMC modified infrared flux method (IRFM – Blackwell & Shallis 1977; Schanche et al. 2020) that allows us to produce synthetic photometry from a constructed spectral energy distribution (SED) using stellar atmospheric models (Castelli & Kurucz 2003) and our stellar spectroscopic parameters as priors. These simulated data are compared to broadband fluxes in the following bandpasses: Gaia G, G_BP, and G_RP, 2MASS J, H, and K, and WISE W1 and W2 (Skrutskie et al. 2006; Wright et al. 2010; Gaia Collaboration 2023) that are listed in Table 4 of Díaz et al. (2020) to derive the stellar bolometric flux. The last step in the MCMC is to convert this flux to the stellar effective temperature and angular diameter, that is combined with the offset-corrected Gaia parallax (Lindegren et al. 2021) to produce the stellar radius.

To accurately characterise the stellar properties of HD 95338 and account for some modelling systematics, we derived two pairs of {M_⋆, t_⋆} from two different approaches. The first mass and age pair were estimated by applying the isochrone placement algorithm (Bonfanti et al. 2015, 2016) that interpolates the input stellar parameters within pre-computed grids of isochrones and tracks of PARSEC⁶ v1.2S (Marigo et al. 2017). The second pair of mass and age estimates were derived from a Levenberg-Marquardt minimisation scheme (as detailed in Salmon et al. 2021) computing stellar models on the fly with the CLÉS code (Code Liégeois d’Évolution Stellaire – Scuflaire et al. 2008). This minimisation scheme optimises over the age and mass to produce the best agreement between model and observed radii and effective temperatures. As detailed in Bonfanti et al. (2021), we checked the mutual consistency between the two respective pairs of outcomes via a χ²-based criterion and then merged the results obtaining M_⋆ = 0.848 ± 0.043 M_⊙ and t_⋆ = ${9.5}_{-5.7}^{+4.3}$ Gyr.

The derived stellar parameters are listed in Table 1. In general, they are in good agreements with the values reported by Díaz et al. (2020), although we present more precise estimates on R_⋆ and log g, which are used in the transit modelling.

In order to increase the precision and accuracy of the planetary model, we incorporated the stellar parameters of Table 1 as inputs to the light curve analysis. By a rearrangement of Kepler’s third law of planetary motion, the scaled semi-major axis is connected to the stellar density as well, since (under the assumption of a spherical star) $${\left(\frac{a}{R_{\star }}\right)}^3=P^{2}G(1+q)\frac{{\rho }_{\star ,\text{transit}}}{3\pi },$$(5)

where q = M_p/M_⋆, and ρ_⋆,transit is a geometrically constrained density, since a/R_⋆ is determined from the transit duration. Using the estimated T_eff, [Fe/H] and ρ_⋆,ransit, we are able to provide estimate M_⋆,transit and R_⋆,transit using the empirical formulae⁷ of Southworth (2011). Then, according to Csizmadia (2020), the goodnes-of-fit metric (in this case, the logarithmic likelihood) can be modified as $$\begin{array}{rl}-\log {\mathcal{L}}_{\text{mod\hspace{0.17em}}}=& -\log \mathcal{L}+\frac{1}{2}{\left(\frac{R_{\star }-R_{\star ,\text{transit\hspace{0.17em}}}}{\sqrt{{\left(\mathrm{\Delta }{R}_{\star }\right)}^2+{\left(\mathrm{\Delta }{R}_{\star ,\text{transit}}\right)}^2}}\right)}^2\\ & +\frac{1}{2}{\left(\frac{\log g-{\log g|}_{\text{transit\hspace{0.17em}}}}{\sqrt{(\mathrm{\Delta }\log g{)}^2+{\left({\mathrm{\Delta }\log g|}_{\text{transit\hspace{0.17em}}}\right)}^2}}\right)}^2,\end{array}$$(6)

where log g|_transit = 4.43–2 log R_⋆,transit + log M_⋆,transit, where log g_⊙ = 4.43. Given that both the second and third terms contain information about ρ_transit, a/R_⋆ is constrained via the input stellar parameters as well.

4 Results

4.1 Planetary parameters

The best-fit CHEOPS and TESS light curves (from the joint fit) are shown in Figs. 3 and 4. The corresponding parameters are shown in Table 2. The parameters describing the planetary transit agree well with each other between the CHEOPS only and the joint TESS plus CHEOPS modelling. We note that the uncertainty ranges do not decrease with the inclusion of the two additional TESS transits, which implies that they are dominated by the uncertainty in the input stellar parameters. The best-fit conjunction parameter is in a <1σ agreement with the impact parameter presented by Díaz et al. (2020), and the relative planetary radius shows a discrepancy of only ≈1.2σ. There is a 4.8σ disagreement between the estimated a/R_⋆ values, and a 2.8σ difference between the T₀ from our joint TESS+CHEOPS analysis and the respective values shown by Díaz et al. (2020). Using the S/N criteria established by Csizmadia et al. (2023) for a/R_⋆, we know that from the CHEOPS analysis alone, we are able to constrain the scaled semi-major axis well within 2% of the truth (since for a 2% accuracy, a S/N = 3 is needed, and in our case, S/N = 22). Consequently, we argue that the uncertainty range presented in Díaz et al. (2020) is likely underestimated. Díaz et al. (2020) find that the first TESS transits occurs 166 ± 52 seconds earlier than our estimates, which is likely a combination of inadequate noise treatment (Kiefer et al. 2023) and the fact that we analyse five transits together. However, the discrepancy in T₀ is <3σ. A brief TTV (transit timing variation) analysis is shown in Sect. 4.2.

The smaller aperture of the TESS cameras implies that the TESS observations are more noisy – this is confirmed by the higher ${\sigma }_w^{\text{TESS}}$ compared to ${\sigma }_w^{\text{CHEOPS}}$. The red noise parameters are not directly comparable between the two instruments, although we expect the this noise type to decrease at longer wavelengths (e.g. Kálmán et al. 2023a). The retrieved A^CHEOPS and B^CHEOPS are within 3σ of the theoretical values. Both limb darkening coefficients in the TESS passband are in agreement with those presented by Díaz et al. (2020) (after conversion to the formalism used in this paper).

Our results are also in agreement with the parameters presented by Saha (2025). Due to the additional transit included in this work and the refined stellar parameters, we estimate slightly lower uncertainties in R_p/R_⋆, b′, P, and T₀. A direct comparison of the confidence intervals of a/R_⋆ and the limb-darkening coefficients is not applicable, since Saha (2025) did not use these as fitting parameters.

We find that the planetary orbit has an inclination of 89.654° ± 0.0155°, the planet has a radius of R_p = 0.3811 ± 0.0035R_Jup, with a transit duration of ${5.587}_{-0.016}^{+0.018}$ hours.

Fig. 3 Phase-folded detrended CHEOPS light curve of HD 95338 (top panel, orange; same data as the third row of Fig. 2.) The rednoise-corrected light curve is shown with blue, while the best-fit transit light curves of HD 95338b is represented by the solid red curve.

Fig. 4 Phase-folded detrended TESS light curve of HD 95338 (top panel, orange, the same data as the third row of Fig. 2.) The rednoise-corrected light curve is shown with blue, while the best-fit transit light curves of HD 95338b is represented by the solid red curve.

Fig. 5 Transit timing data of HD 95338b. The dashed purple lines represent the uncertainty range from the global fit.

4.2 Transit timing

We extracted the timing of each individual transit (two from TESS, four from CHEOPS) by fixing all parameters to the best-fit values from Table 2 with the exception of σ_w, σ_r and T₀, assuming a constant orbital period. We also included a possible quadratic trend in the modellings. The resultant transit timing data are shown in Fig. 5 and are compiled in Table 3. There is no apparent TTV, meaning that the transit timings are explained well by a constant period. The T₀ value from sector 63 is in a 2.5σ disagreement with the global modelling, which is not statistically significant. A deeper exploration of the TTV signals is beyond the scope of this paper.

5 The search for moons

We chose to restrict the search for exomoons on the CHEOPS dataset alone. This is done because the smaller (mirror) aperture of TESS in comparison to CHEOPS manifests as higher white noise levels, which are detrimental to the quest for moons with the methodology presented here (191 ppm for CHEOPS compared to 289 ppm in TESS despite the exposure time that is more than 4 times shorter, Table 2). For that reason, and to keep the analysis self-consistent, we use the planetary orbital parameters (P, b and a/R_⋆) that are derived from the CHEOPS photometry alone in the processes described below.

5.1 Sensitivity maps

We examined the CHEOPS residual light curves for shallow transit signals (that could be linked to an exomoon) following the approach detailed in Ehrenreich et al. (2023). This method explores a grid of moon transit epochs and durations {T_0,Moon, W}. For each point of the grid, we fit for a box-shaped transit depth (allowed to be negative) together with a trend (linear or quadratic) and a multiplicative noise term (factor applied to the data error bars). The fit is done using the Markov chain Monte Carlo algorithm (MCMC) implemented in the Python package emcee (Foreman-Mackey et al. 2013). From the posterior distribution on the transit depth, we derive the best-fit value of the moon radius and its 3-σ lower limit. We also compute the Bayesian Information Criterion and the Akaike Information Criteria (AIC; Akaike 1998) defined as: $$\mathrm{B}\mathrm{I}\mathrm{C}=k\ln (n)-2\ln (\mathcal{L})$$(7) $$\mathrm{A}\mathrm{I}\mathrm{C}=2k-2\ln (\mathcal{L}),$$(8)

where n is the number of data points, k is the number of parameters included in the fit, ℒ is the best-fit (maximum) likelihood value. The BIC and AIC values are compared to their respective values in the case where the moon transit depth is fixed to 0, i.e. a ‘no-moon’ reference fit.

The parameter space scanned by the {T_0,Moon, W} grid must cover mid-transit times across the full Hill-sphere of HD 95338b⁸ while allowing for either shorter or longer transit duration (e.g. if the moon impact parameter is larger or smaller) than that of its planet. With a transit duration of the Hill sphere of 17.6 ± 0.9 hours (i.e. <20.3 h at 3 σ), we considered T_0,Moon values that satisfy −10.5 h ≤ ΔT₀ ≤ + 10.5 h, where ΔT₀ = T_0,Moon − T_0,Planet. The range of moon transit durations W covers values from 30 min (grazing moon) to 7 h, with a planetary transit duration of ~5.6 h. The grid resolution (increment) is set to 15 min in both dimensions.

In order to provide conservative estimates of a possible exomoon signal in the T₀−W parameter space, we propagate the uncertainties from the planetary transit fit. To that end, we calculate the standard deviation of 100 red noise plus planetary transit (see Sect. 4.1) models (with parameters selected randomly from the posterior distribution of the fit) at every time stamp of the observations, and add this to the flux uncertainties estimated by DRP in quadrature. In order to increase computational efficiency, we fit F_BV via a least-squares algorithm. Consequently, we do not propagate the uncertainties of F_BV into the residuals.

5.2 Box-shift method

We performed a more detailed transit modelling after the planetary transit in the third visit to try to explore the possibility of an additional transit being there. To do that, we used narrow box priors on T_0,Moon, which allowed us to perform an in-depth scan of the selected areas, by shifting them through the light curve. The priors are chosen to be 3 hours wide (≈1/2 transit duration), and at each new step, their borders are increased by 1.5 hours (≈1/4 transit duration), yielding an overlapping configuration for the transit modelling. In total, we performed 8 different transit modellings in every light curve.

We assumed that the inclination of the orbit of the moon is such that the impact parameter is the same as for HD 95338b. We further assume that the separation of the planet and its companion is so small that the star–moon distance can be described sufficiently well with the semi-major axis of the planetary orbit, and that the relative velocity of the moon is negligible compared to the orbital velocity of the system. The difference in the transit duration, although difficult to estimate as it depends on the mass of the moon and its semi-major axis in the orbit around the planet. Based on the specific configuration, it is possible that a moon within the Hill sphere would not be transiting or would show a grazing transit (Ehrenreich et al. 2023), or that the duration of the lunar transit would exceed that of the planet (Kálmán et al. 2024a). It is also established that the limb-darkening coefficients can more reliably be measured in case of deeper transits (e.g. Eq. (34) of Csizmadia et al. 2023). For these reasons, we only fit the depth and epoch of the lunar transit, while we adopt all other transit parameters from the analysis of the planetary transits (Table 2).

5.3 Injection-and-retrieval tests

In order to verify the performance of the sequential exomoon-detecting algorithm presented above, we performed a number of injection-and-retrieval tests. We injected transits with the a/R_⋆ and b values of Díaz et al. (2020). We tested a wide range of moon radii, with a grid consisting of R_{Moon,injected} ∈ {0.8, 0.9, 1.05, 1.2, 1.35} R_⊕. The transits were injected at two epochs (BJD 2459301.25 in the first visit and BJD 2459687.05 in the second visit) so that they overlap with the transit of HD 95338b. We followed the same steps as described in Sect. 2.1, then proceeded to generate the sensitivity maps, and, having identified the transits on these we also performed the box-shifting approach as well. There is an infinite number of possible T_0,Moon − R_Moon combinations that could be used in such tests. These specific cases are useful to demonstrate that the technique described above works well for finding lunar transits. We only use data from the first three CHEOPS visits for the injection-and-retrieval tests.

In Figs. B.2 and B.3, we show the sensitivity maps for the R_{Moon,injected} = 0.8 R_⊕ and R_{Moon,injected} = 1.35 R_⊕ cases. The second rows, showing the 3σ lower limit on the moon radii show on clear strip for about 4 hours before and 5 hours after the conjunction of HD 95338b in the first and second visits, respectively. The presence of the signal is also confirmed by the ΔBIC and ΔAIC values (third and fourth rows of Figs. B.2 and B.3, respectively). The ΔBIC and ΔAIC metrics also show that at certain times (≈2 hours after and ≈1 hour before the midtransit of the planet) positive bumps (i.e. negative moon radii) are preferred to the combination of GPs and the quadratic trend. Such signals do not correspond to the transit of hypothetical exomoon. Comparing with Fig. B.1, we may conclude that they appear as a consequence of the injected lunar transits. It is likely that these two additional signals act as weights on the planetary transit, causing minor distortions that are then detected during the sensitivity mapping.

After identifying these regions of interest (between ≈2301.04 and ≈2301.26 for the first visit, ≈2686.93 and ≈2697.15 for the second visit), we also conducted the box-shifting parameter extraction similarly to the study of the real dataset (Sec. 5.2). We include a quadratic trend (Eq. (6)) in the analyses. For the cases of R_{Moon,injected} = 0.8 R_⊕ and R_{Moon,injected} = 1.35 R_⊕, the light curves from this method are shown on Figs. B.4 and B.5. We can observe the transit appearing and then disappearing again as the searchbox for T_0,Moon is shifted through the light curves.

The extracted T_0,Moon and R_Moon values are listed in Tables B.1 and B.2 for all tested moon sizes, along with the respective deviations from the nominal values. The transit depths are recovered even for R_{Moon,injected} = 0.8 R_⊕ case in both visits. We find that the significance of this parameter (and consequently the transit itself) is ≈5.5 R_⊕ in the first visit and ≈4.2 in the second visit for R_{Moon,injected} = 0.8 R_⊕. The significances of the detected transits naturally increase with R_{Moon,injected}. Table 3 suggests that in the second visit, the time-correlated noise level is higher at σ_r = 1356.0 ± 259.8 ppm compared to just σ_r = 589.8 ± 255.2 ppm. We emphasise that anything not explicitly included in the modelling (via a fitting parameter) is described by the wavelets. For this reason, σ_r depends heavily on the underlying light curve (whether it contains injected transits or not) and the actual fitted quantities. Consequently, comparison of the red noise levels between the various R_{Moon,injected} cases is not possible. The ≈900 ppm difference in the red noise levels may explain the larger uncertainties on R_Moon seen in the second visit. Similarly, the in the third visit, we find σ_r = 2443.8 ± 282.3 ppm (Table 3), which may shroud any evidence for an exomoon. A comprehensive analysis of the underlying causes of these discrepancies is beyond the scope of this paper, but it should be noted that the instrument displays a certain degree of sensitivity loss in relation to its ageing process (Fortier et al. 2024, Sect. 4).

We note that the new window on the CCD yielded considerably lower red noise (Table 3).

Although one of the injected moon transits is on the ingress of the planetary transit and the other one is on the egress, they are not symmetric with respect to T₀. Consequently, the shape of the modelled planetary transit gets distorted in comparison with the original dataset, and this discrepancy is more evident the larger the injected exomoon radius is. In this sequential approach of the chase for moons, we are in fact fitting the photocenter (Simon et al. 2007) as the ‘planet’. The retrieved ‘planetary’ parameters are therefore more distorted if we increase R_{Moon,injected}. This is shown in Table 4. We observe a shift in T₀ to earlier epochs, and increase in R_p/R_⋆, and a decrease in a/R_⋆. Although the recovered parameters agree in a 1σ with each other and even the original dataset (Table 5), these minor inconsistencies may lead to differences in the residuals that are then converted into the parameter biases seen in Tables B.1 and B.2. Furthermore, the additional signals (i.e. transits) can also cause and increased variance in the 100 models chosen from the posteriors for the uncertainty propagation. As a result, the hint of an exomoon seen on the left column Fig. B.1 disappears (left columns of Figs. B.2, B.3). A more detailed investigation of the distortion of the planetary parameters is beyond the scope of this paper.

Given that the distortion of the transit parameters and the flux uncertainties increases with R_{Moon,injected}, the recovered R_Moon and T_0,Moon are also further from their input values at the higher moon radii (Tables B.1, B.2). At the same time, the lower moon sizes yield {T_0,Moon, R_Moon} combinations that are within 1σ of the injections.

6 The search for moons orbiting HD 95338b

Having verified that we are able to achieve >3σ detections of additional transiting bodies with 0.8 R_⊕ based on single transits (Sect. 5.3), we explored the residual light curves of HD 95338 from the four CHEOPS visits. The sensitivity map for each visit of HD 95338 is shown in Fig. 6 (with an underlying linear trend), which hints and the presence of a signal that occurs slightly later than the planetary transit during the second visit (middle column), and a strong signal in the third visit (right column). The signals in these two visits can be described by a box with a depth corresponding to 0.8 R_⊕ and ≈1.2 R_⊕, respectively. We also estimate 3-σ lower limits at ≈0.6 R_⊕ (in the second visit) and ≈0.8 R_⊕ (in the third visit). We find that the cases where a box-shaped signal is included in the fit is favoured by ΔBIC ≥ 10 and ΔAIC ≥ 10. In the first visit (Fig. 6, left column), the lower limit on the hypothetical transit depth is negative at every (ΔT₀, W) grid point. In the case of the second transit, the detection statistics are more sporadically distributed in the parameter space and not in a manner that we expect for a transit based on the injection–and–retrieval tests (see Sect. 5.3).

The spread-out detections (non-zero 3σ lower limits) before the transit in the second visit are likely noise. In the first and fourth visits, the best-fit box does not have a 3σ detection significance (i.e. the 3σ lower limit on the radius is negative), and the inclusion of the box is disfavoured by the BIC. For these reasons, we conclude that a moon signal is not detectable in the first two visits or the last visit, while the signal seen in the third visit warrants further investigations. We note that at the bottom edge of the maps (Fig. 6) the width of the boxes (i.e. the duration of the supposed lunar transit) is commensurable with the duration of the gaps in observations, as discussed in Sect. 2.1. For this reason, any signal seen there can also safely be disregarded.

We also constructed the sensitivity maps by using a quadratic trend (Fig. B.1). Firstly, we note that no signals appear in the first or fourth visits. Secondly, we observe that the features seen in the second and third visits (Fig. 6) disappear, as the 3σ lower limit (second row of Fig. B.1) no longer implies a clear detection. The BIC and AIC values still suggest that the inclusion of a dip is preferred, although lower with probabilities (ΔBIC ≈8 instead of ≥10 as seen in the case when a linear trend is included).

After narrowing down the possible locations of the exomoon transit via the sensitivity maps (Fig. 6), we pinned down the exact position and duration of the signal. We identified the relevant light curve section to be between ≈246072.427 BJD and ≈246072.777 BJD. We constructed eight overlapping segments (each spanning 3 hours, or roughly half of the transit duration of HD 95338b). These segments are used as uniform priors (or searchboxes) for mid-time of the lunar transit. To scan the above-mentioned portion of the light curve, we shift these searchboxes by 45 minutes (or one-eighth of the planetary transit duration) for each individual fit. We thus performed eight individual fits in search of the lunar transit, as seen on Fig. 7. We also included a quadratic trend (for the flux–time detrending) in the form of $$F_{\text{trend\hspace{0.17em}}}=p_0+p_1\cdot (t-T_{\text{ref\hspace{0.17em}}})+p_2\cdot {(t-T_{\text{ref\hspace{0.17em}}})}^2,$$(9)

where p₀, p₁, and p₂ are the constant, linear, and quadratic coefficients respectively, and t_ref = T_0,Moon is the reference time of the trend. We also allow for the two parameters of the wavelet-based noise filtering to vary in an unconstrained manner. The results of the box-shifting method are shown on Fig. 7, and the best-fit parameters are listed in Table 5. The noise-related parameters (σ_r, σ_w, p₀, p₁ and p₂) are consistent with each other in all eight cases, while R_Moon/R_⋆ is consistent with 0 within 1σ. We therefore do not see a significant detection of a transit-like dip.

Fig. 6 Sensitivity maps for the four CHEOPS visits of HD 95338. The first row represents the best-fit box depth (i.e. moon radius) in the ΔT0−W space. The 3σ lower limit of the moon radius is shown in the second row, while the third and fourth rows represent the ΔBIC and ΔAIC values between a solution with and without a moon. The empty (blank) areas are values that are below a threshold of 0 R⊕ for the moon radius (first 2 rows) and a threshold of 2 for the ΔBIC and ΔAIC (last 2 rows). The grey squares (and corresponding dotted lines) show the position of the transit of HD 95338b in this parameter space for each visit. An underlying linear trend is assumed for each of the four visits.

Fig. 7 Results of the box-shifting method. The shaded gray areas represent the searchboxes within which the mid-time of the lunar transit is constrained. The best-fit moon radius and its 1σ uncertainty is shown on the bottom of each plot.

7 Discussion and concluding remarks

In this work, we have carried out a search for possible satellites orbiting the eccentric Neptune-sized planet HD 95338b in a sequential approach. Using four transits from CHEOPS and two from TESS, we present a detailed transit modelling, followed by a novel approach in the quest for exomoons. We re-derived the stellar parameters of the host star HD 95338, and present considerable improvements on some of them in comparison with the discovery paper (Díaz et al. 2020). We present improvements on the basic transit parameters of HD 95338b. Most notably, we reduce the uncertainty on the planetary radius by a factor of 10 compared to the results presented by Díaz et al. (2020). Given that the precision of these parameters does not change with the inclusion/exclusion of the two TESS transits, we argue that the stellar parameters and their uncertainties are properly taken into account in TLCM (as described in Csizmadia 2020) and that the uncertainties of the modelling parameters are limited by our knowledge of the star. The most remarkable improvement is achieved regarding the orbital period, which we estimate with a precision of ≈0.3 s compared to ≈1700 s in the discovery paper.

The quest for exomoons is difficult as we are trying to find transits whose timing and duration are not known, and that are shallow (especially in comparison with the instrumental noise effects of CHEOPS). A thorough detrending, incorporating all known noise sources is therefore a pre-requisite for such a project, as outlined in Sect. 2. We scanned the residual light curves after subtracting the planetary signal for signs of additional transit-like features that could be attributed to an exomoon orbiting HD 95338b. For a thorough characterisation of an exomoon which would incorporate the combination of two Keplerian orbits (Hippke & Heller 2022), at least three transit detections would be needed. We are able to rule out the presence of additional signals below ≈0.6 R_⊕ at 1σ, thus placing an upper limit on the size of a hypothetical satellite in the system.

There are obvious long-term trends in the residuals seen in Fig. 2. With the assumption that these are linear (and letting the complex noise handling tools, either GP-based or wavelet-based, handle all other possible non-white noise realisations), we observed the detection of a transit-like feature (Fig. 6). However, with the addition of a quadratic trend (Figs. B.1 and 7), the transit detection is no longer significant. The deviation between the fitted linear and quadratic trends is the greatest around the position of the detected transit. The fact that the significance of the detection decreases after including a higher-order time-dependent term – which is commonly associated with astrophysical noise sources, such as spots – implies that the transit-like feature is likely not originated by a transit of an exomoon. Increasing the order of polynomials included in the modelling introduces the risk that the shorter timescale noise components are enhanced in a way that may introduce additional transit-like features. Thus we do not test for these possibilities, relying instead on the noise handling algorithms that are thought to be more efficient at higher frequencies. In any case, the retrieved moon size is consistent with a no-moon case at 3σ. Additionally, the lack of obvious TTVs (Fig. 5) also serves as counterargument against the presence of a satellite (Simon et al. 2007; Kipping 2009).

Based on the analysis of the synthetic LCs, we conclude that the method presented in Sect. 5 is viable for searching for exomoon transits in a non-photodynamical way. The injection-and -retrieval tests suggest that the combination of the sensitivity maps and the box-shifting method does not induce false positive detections, however, it does recover the size and position of the injected exomoon transits (within 2σ of their respective input values), and leads to >3σ detections even for the smallest tested moon sizes.

HD 95338b is an interesting target and an important milestone in the quest for exomoons. Future observations with CHEOPS, PLATO (PLAnetary Transits and Oscillations of stars Rauer et al. 2014, 2025), Ariel (Atmospheric Remote-sensing Infrared Exoplanet Large-survey Tinetti et al. 2018, 2022), or JWST (James Webb Space Telescope Gardner et al. 2006) may further reduce the current upper limit of 0.6 R_⊕ on the radius of a possible moon in the system. The methodology presented in Sections 2 and 5 can also be regarded a stepping stone in the hunt for exomoons using the telescopes of the future. In the event of detections of additional signals that are similar to the one seen on Fig. 7, a fully photodynamic characterisation may be able to provide deeper insights into the nature of this system.

The discovery of the first confirmed exomoon will open a new chapter in exoplanet research. The above-told tale about the exomoon hunt around HD 95338b, together with the earlier examples mentioned in Section 1, sends a strong message to the research community regarding this subject. Detecting exomoons is currently situated at the frontiers of contemporary research, and success in confirming such discoveries can only be achieved through a comprehensive understanding and in-depth exploration of the intricacies involved in signal design, measurement techniques, and data analysis. Our discussion also outlines potential pathways to be explored during this investigation, and we are convinced that our considerations in the above analysis will provide significant methodological tools to this demanding task.

Data availability

The raw and detrended photometric time series data are available at the CDS via https://cdsarc.cds.unistra.fr/viz-bin/cat/J/A+A/701/A240

Acknowledgements

CHEOPS is an ESA mission in partnership with Switzerland with important contributions to the payload and the ground segment from Austria, Belgium, France, Germany, Hungary, Italy, Portugal, Spain, Sweden, and the United Kingdom. The CHEOPS Consortium would like to gratefully acknowledge the support received by all the agencies, offices, universities, and industries involved. Their flexibility and willingness to explore new approaches were essential to the success of this mission. CHEOPS data analysed in this article will be made available in the CHEOPS mission archive (https://cheops.unige.ch/archive_browser/). This work has made use of data from the European Space Agency (ESA) mission Gaia (https://www.cosmos.esa.int/gaia), processed by the Gaia Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium). Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement. On behalf of the ‘Searching for transiting exomoons in CHEOPS data’ project Sz.K. is grateful for the possibility to use HUN-REN Cloud (see Héder et al. 2022, https://science-cloud.hu/) which helped us achieve the results published in this paper. ASi and C.Br. acknowledge support from the Swiss Space Office through the ESA PRODEX program. A.De. Gy.M.Sz. acknowledges the support of the Hungarian National Research, Development and Innovation Office (NKFIH) grant K-125015, a a PRODEX Experiment Agreement No. 4000137122, the Lendúlet LP2018-7/2021 grant of the Hungarian Academy of Science and the support of the city of Szombathely. This project has received funding from the Swiss National Science Foundation for project 200021_200726. It has also been carried out within the framework of the National Centre of Competence in Research PlanetS supported by the Swiss National Science Foundation under grant 51NF40_205606. The authors acknowledge the financial support of the SNSF. T.Wi. acknowledges support from the UKSA and the University of Warwick. S.G.S. acknowledge support from FCT through FCT contract no. CEECIND/00826/2018 and POPH/FSE (EC). The Portuguese team thanks the Portuguese Space Agency for the provision of financial support in the framework of the PRODEX Programme of the European Space Agency (ESA) under contract number 4000142255. The Belgian participation to CHEOPS has been supported by the Belgian Federal Science Policy Office (BELSPO) in the framework of the PRODEX Program, and by the University of Liège through an ARC grant for Concerted Research Actions financed by the Wallonia-Brussels Federation. G.Sc., L.Bo., V.Na., I.Pa., G.Pi., R.Ra., and T.Zi. acknowledge support from CHEOPS ASI-INAF agreement no. 2019-29-HH.0. ABr was supported by the SNSA. YAl acknowledges support from the Swiss National Science Foundation (SNSF) under grant 200020_192038. We acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2019-107061GB-C61, PID2019-107061GB-C66, PID2021-125627OB-C31, and PID2021-125627OB-C32, from the Centre of Excellence “Severo Ochoa” award to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), from the Centre of Excellence “María de Maeztu” award to the Institut de Ciències de l’Espai (CEX2020-001058-M), and from the Generalitat de Catalunya/CERCA programme. DBa, EPa, and IRi acknowledge financial support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación MCIN/AEI/10.13039/501100011033 and the ERDF “A way of making Europe” through projects PID2019-107061GB-C61, PID2019-107061GB-C66, PID2021-125627OB-C31, and PID2021-125627OB-C32, from the Centre of Excellence “Severo Ochoa’ award to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), from the Centre of Excellence “María de Maeztu” award to the Institut de Ciències de l’Espai (CEX2020-001058-M), and from the Generalitat de Catalunya/CERCA programme. SCCB acknowledges the support from Fundação para a Ciência e Tecnologia (FCT) in the form of work contract through the Scientific Employment Incentive program with reference 2023.06687.CEECIND. ACC acknowledges support from STFC consolidated grant number ST/V000861/1, and UKSA grant number ST/X002217/1. ACMC acknowledges support from the FCT, Portugal, through the CFisUC projects UIDB/04564/2020 and UIDP/04564/2020, with DOI identifiers https://doi.org/10.54499/UIDB/04564/2020 and https://doi.org/10.54499/UIDP/04564/2020, respectively. A.C., A.D., B.E., K.G., and J.K. acknowledge their role as ESA-appointed CHEOPS Science Team Members. P.E.C. is funded by the Austrian Science Fund (FWF) Erwin Schroedinger Fellowship, program J4595-N. This project was supported by the CNES. This work was supported by FCT – Fundação para a Ciência e a Tecnologia through national funds and by FEDER through COMPETE2020 through the research grants UIDB/04434/2020, UIDP/04434/2020, 2022.06962.PTDC. O.D.S.D. is supported in the form of work contract (DL 57/2016/CP1364/CT0004) funded by national funds through FCT. B.-O.D. acknowledges support from the Swiss State Secretariat for Education, Research and Innovation (SERI) under contract number MB22.00046. A.De., B.Ed., K.Ga., and J.Ko. acknowledge their role as ESA-appointed CHEOPS Science Team Members. C.B. acknowledges support from the Swiss NCCR PlanetS. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40182901 and 51NF40205606. J.K. acknowledges support from the Swiss National Science Foundation under grant number TMSGI2_211697. MF and CMP gratefully acknowledge the support of the Swedish National Space Agency (DNR 65/19, 174/18). DG gratefully acknowledges financial support from the CRT foundation under Grant No. 2018.2323 “Gaseousor rocky? Unveiling the nature of small worlds”. M.G. is an F.R.S.-FNRS Senior Research Associate. CHe acknowledges the European Union H2020-MSCA-ITN-2019 under GrantAgreement no. 860470 (CHAMELEON), and the HPC facilities at the Vienna Science Cluster (VSC). MNG is the ESA CHEOPS Project Scientist and Mission Representative. BMM is the ESA CHEOPS Project Scientist. KGI was the ESA CHEOPS Project Scientist until the end of December 2022 and Mission Representative until the end of January 2023. All of them are/were responsible for the Guest Observers (GO) Programme. None of them relay/relayed proprietary information between the GO and Guaranteed Time Observation (GTO) Programmes, nor do/did they decide on the definition and target selection of the GTO Programme. K.W.F.L. was supported by Deutsche Forschungsgemeinschaft grants RA714/14-1 within the DFG Schwerpunkt SPP 1992, Exploring the Diversity of Extrasolar Planets. This work was granted access to the HPC resources of MesoPSL financed by the Region Ile de France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir supervised by the Agence Nationale pour la Recherche. This work has been carried out within the framework of the NCCR PlanetS supported by the Swiss National Science Foundation under grants 51NF40_182901 and 51NF40_205606. AL acknowledges support of the Swiss National Science Foundation under grant number TMSGI2_211697. ML acknowledges support of the Swiss National Science Foundation under grant number PCEFP2_194576. PM acknowledges support from STFC research grant number ST/R000638/1. This work was also partially supported by a grant from the Simons Foundation (PI Queloz, grant number 327127). NCSa acknowledges funding by the European Union (ERC, FIERCE, 101052347). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. V.V.G. is an F.R.S-FNRS Research Associate. JV acknowledges support from the Swiss National Science Foundation (SNSF) under grant PZ00P2_208945. EV acknowledges support from the ‘DISCOBOLO’ project funded by the Spanish Ministerio de Ciencia, Innovación y Universidades undergrant PID2021-127289NB-I00. NAW acknowledges UKSA grant ST/R004838/1. Project no. C1746651 has been implemented with the support provided by the Ministry of Culture and Innovation of Hungary from the National Research, Development and Innovation Fund, financed under the NVKDP-2021 funding scheme.

Appendix A Information about the CHEOPS observations

The metadata of the four CHEOPS observations are shown in Table A.1. The list of nearby stars, that may interfere with the aperture photometry, is shown in Table A.2.

Appendix B Sensitivity maps

B.1 Quadratic trend

The sensitivity maps in case of a quadratic trend (for the real observations) are shown on Fig. B.1.

Fig. B.1 Same as Fig. 6 but with a quadratic trend instead of a linear trend.

B.2 Synthetic transits

The results of injection-and-retrieval tests are shown in Figs. B.2 – B.5, and Tables B.1 – B.1.

Fig. B.2 Same as Fig. 6 but with a quadratic trend instead of a linear trend and with an RMoon = 0.8 R⊕ moon injected into the first two visits. The injection-and retrieval tests were carried out only on the first three CHEOPS visits.

Fig. B.3 Same as Fig. B.2 but with a quadratic trend instead of a linear trend and with an RMoon = 1.35 R⊕ moon injected into the first two visits.

Fig. B.4 Performance of the box-shifting approach on the injected exomoon transits. Top rows: Recovery of an RMoon = 0.8 R⊕ in the first visit. Bottom row: Recovery of an RMoon = 0.8 R⊕ in the second visit.

Fig. B.5 Performance of the box-shifting approach on the injected exomoon transits. Top row: Recovery of an RMoon = 1.35 R⊕ in the first visit. Bottom row: Recovery of an RMoon = 1.35 R⊕ in the second visit.

References

All Tables

All Figures

	Fig. 1 The Subarray image depicts the point spread function (PSF) of HD 95338b and the neighbouring stars, with their respective magnitudes of brightness indicated. The target is centred on the coordinates (180, 178) on the CCD.
In the text

Fig. 2 Demonstration of the fitting and detrending steps. The raw CHEOPS light curves of HD 95338b are shown in the top row, along with the initially fitted transit curves (step 1). The residuals of the initial transit fit, along with the best-fit systematics model is shown in the second row, with N number of harmonics of the roll angle included. The final transit fit (step 4) and its residuals (later used in the search for additional transits) are shown in the bottom two rows, where the full white circles represent 200-point bins (with an average uncertainty of ≈20 ppm in each bin).

In the text

	Fig. 3 Phase-folded detrended CHEOPS light curve of HD 95338 (top panel, orange; same data as the third row of Fig. 2.) The rednoise-corrected light curve is shown with blue, while the best-fit transit light curves of HD 95338b is represented by the solid red curve.
In the text

	Fig. 4 Phase-folded detrended TESS light curve of HD 95338 (top panel, orange, the same data as the third row of Fig. 2.) The rednoise-corrected light curve is shown with blue, while the best-fit transit light curves of HD 95338b is represented by the solid red curve.
In the text

	Fig. 5 Transit timing data of HD 95338b. The dashed purple lines represent the uncertainty range from the global fit.
In the text

Fig. 6 Sensitivity maps for the four CHEOPS visits of HD 95338. The first row represents the best-fit box depth (i.e. moon radius) in the ΔT0−W space. The 3σ lower limit of the moon radius is shown in the second row, while the third and fourth rows represent the ΔBIC and ΔAIC values between a solution with and without a moon. The empty (blank) areas are values that are below a threshold of 0 R⊕ for the moon radius (first 2 rows) and a threshold of 2 for the ΔBIC and ΔAIC (last 2 rows). The grey squares (and corresponding dotted lines) show the position of the transit of HD 95338b in this parameter space for each visit. An underlying linear trend is assumed for each of the four visits.

In the text

	Fig. 7 Results of the box-shifting method. The shaded gray areas represent the searchboxes within which the mid-time of the lunar transit is constrained. The best-fit moon radius and its 1σ uncertainty is shown on the bottom of each plot.
In the text

	Fig. B.1 Same as Fig. 6 but with a quadratic trend instead of a linear trend.
In the text

	Fig. B.2 Same as Fig. 6 but with a quadratic trend instead of a linear trend and with an RMoon = 0.8 R⊕ moon injected into the first two visits. The injection-and retrieval tests were carried out only on the first three CHEOPS visits.
In the text

	Fig. B.3 Same as Fig. B.2 but with a quadratic trend instead of a linear trend and with an RMoon = 1.35 R⊕ moon injected into the first two visits.
In the text

	Fig. B.4 Performance of the box-shifting approach on the injected exomoon transits. Top rows: Recovery of an RMoon = 0.8 R⊕ in the first visit. Bottom row: Recovery of an RMoon = 0.8 R⊕ in the second visit.
In the text

	Fig. B.5 Performance of the box-shifting approach on the injected exomoon transits. Top row: Recovery of an RMoon = 1.35 R⊕ in the first visit. Bottom row: Recovery of an RMoon = 1.35 R⊕ in the second visit.
In the text