© The Authors 2026

1 Introduction

Exoplanetary reflection signatures in high-resolution spectroscopic datasets continue to elude detection efforts. Despite advances in instrumentation, the faint nature of these signals poses a significant challenge. The first attempt at detecting an exoplanet’s reflection signature in high-resolution spectroscopy was conducted as a series of case studies of the τ Bootis system (Charbonneau et al. 1999; Collier Cameron et al. 1999, 2000). While resulting in non-detections, these efforts yielded gradually improving constraints on the radius and albedo of τ Bootis b (Leigh et al. 2003; Rodler et al. 2010, 2013; Hoeijmakers et al. 2018; notable examples targeting other systems include Collier Cameron et al. 2002, Rodler et al. 2008, and Martins et al. 2015). One of the major complications in those early studies was that the planet’s orbit was not precisely known, making it impossible at the time to shift the obtained spectra to the correct planetary rest frame. Incidentally, the detection of water vapour in the planet’s atmosphere by Lockwood et al. (2014) has helped to better constrain the orbital solution. Whether or not the water vapour detection is robust is still a matter of debate, however (e.g. Pelletier et al. 2021; Webb et al. 2022).

Despite these difficulties, pursuing reflection signals remains a worthwhile endeavour and has recently begun to re-attract some attention (e.g. Scandariato et al. 2021; Spring et al. 2022; Vaughan et al. 2026). Once successfully realised, the observation of reflected light would represent a direct planetary detection that can constrain the orbital solution of a given system without the need for a planetary transit. Accordingly, such a detection can help break the mass degeneracy associated with an unknown inclination. As the investigations into τ Bootis have demonstrated, even a non-detection can place upper limits on a planet’s albedo, a parameter with limited observational constraints in exoplanetary research. Beyond hot Jupiters, reflection spectroscopy could become particularly valuable for studying temperate exoplanets, which have low transit probabilities around Sun-like stars. Additionally, crucial biomarkers such as oxygen imprint spectral signatures in reflected light (Meadows et al. 2018). Ultimately, advanced reflection studies can complement both transmission and emission spectroscopy, much like reflected-light observations of Solar System planets have been key to our understanding of atmospheres closer to home (see e.g. Walker 2017; Irwin et al. 2022).

In this paper we introduce a new approach for assessing a system’s suitability for reflection spectroscopy, the reflection spectroscopy metric (RSM). We then demonstrate its use by applying it to real data in an attempt to detect or constrain the reflection signal of a specific exoplanet, KELT-9 b. This planet orbits an early-type star on a highly misaligned orbit, introducing the challenge of having to account for variable line broadening. This stems from the fact that throughout the course of an orbit around its rapidly rotating host, the planet will perceive different degrees of stellar rotation. We discuss how neglecting such rotational effects can strongly impact the detectability of reflected light.

The rest of this manuscript is organised as follows. Section 2 describes the methodology, and Sect. 3 explains how the RSM guides target selection, culminating in the choice of KELT-9 b. In Sect. 4 we introduce our observations with CARMENES and describe the employed data reduction procedure. Section 5 deals with the analysis and results, including the injection tests we carried out to assess detectability limits. Finally, we discuss the results and offer our conclusions in Sects. 6 and 7, respectively.

2 Methodology

2.1 Reflected flux and albedo

The planet-to-star contrast in reflected light, can be written as $$\epsilon :=\frac{F_{\mathrm{r}}}{F_{\mathrm{s}}}=A_{\mathrm{g}}\mathrm{\Phi }{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2,$$(1)

where F_r and F_s are the reflected and stellar flux, respectively, while A_g is the geometric albedo, Φ is the phase function, R_p is the planetary radius, and a is the orbital semi-major axis (Charbonneau et al. 1999; Kreidberg 2018). In principle, Φ depends on the orbital phase angle, reaches its maximum when the planet is fully illuminated from our line of sight (just before and after secondary eclipse), and vanishes at the primary transit midpoint. Inferring A_g from a successful detection or a stringent upper limit on ε thus provides insights into atmospheric composition, condensate clouds, or reflection from gas-phase scattering (e.g. Collier Cameron et al. 2002; Leigh et al. 2003).

2.2 Doppler motion and rotational broadening

High-resolution Doppler spectroscopy leverages the planet’s radial motion to separate the faint planetary signal from the much stronger stellar lines and telluric contamination (Rodler et al. 2010; Hoeijmakers et al. 2018). For close-in planets, the orbital velocity can be hundreds of kms⁻¹, significantly shifting the planetary lines over the course of an orbit.

Beyond the orbital motion, in general an exoplanet also sees its host star rotate. For a planet orbiting in the host’s equatorial plane, the perceived stellar rotational velocity amounts to $$v_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{s},\mathrm{p}}=2\pi R_{\mathrm{s}}(\frac{1}{P_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{s}}}-\frac{1}{P_{\mathrm{o}\mathrm{r}\mathrm{b}}}),$$(2)

where R_s is the stellar radius, P_rot,s is the star’s rotation period, and P_orb is the planetary orbital period. For a planet itself rotating with period P_rot,p, an additional velocity component, $$v_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{p}}=2\pi \sin i\frac{R_{\mathrm{p}}}{P_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{p}}},$$(3)

further broadens the lines (Rodler et al. 2010). Here, i is the inclination of the planet’s rotational axis. The net line-broadening velocity relevant to the reflected spectrum is then $$v_{\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{l}}=\sqrt{v_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{s},\mathrm{p}}^2+v_{\mathrm{r}\mathrm{o}\mathrm{t},\mathrm{p}}^2}.$$(4)

As we demonstrate later and as was shown by Spring et al. (2022), neglecting this broadening can decrease the sensitivity of cross-correlation procedures to a planet’s reflected signal, especially in strongly misaligned systems.

3 Target selection

Central to guiding this and future studies in terms of their target selection is the development of a suitable metric that quantifies the relative favourability of different systems for reflection studies. To this end, we introduce the reflection spectroscopy metric (RSM), a figure of merit mirroring the definitions of the transmission spectroscopy metric and emission spectroscopy metric in Kempton et al. (2018). In contrast to these, the RSM is conceived to inform and facilitate searches using high-resolution spectroscopy, however. The favourability of a given system for such studies is governed by

• (i)

  the apparent brightness of the host star in the visual wavelength domain, m_V, since it is directly related to the achievable signal-to-noise (S/N) of the reflection signature;

• (ii)

  the solid angle of the planet, Ω_p ∝ (R_p/a)², as seen from the host star since it defines the size of the reflecting planetary disc;

• (iii)

  the planetary radial velocity semi-amplitude, K_p, since the spectral separation of the direct and reflected stellar lines increases with higher orbital radial velocities;

• (iv)

  the number of spectral lines, N, in the host star’s spectrum since cross-correlation techniques, whose performance improves with N, will be employed in the detection attempt;

• (v)

  the geometric albedo of the planet, A_g, since it defines the reflected fraction of the total light incident on the planetary disc.

In an attempt to construct a working definition of the RSM, we combined the different factors listed above into $$\mathrm{R}\mathrm{S}\mathrm{M}={10}^{-m_V/5}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2{K}_{\mathrm{p}}\sqrt{N}{A}_{\mathrm{g}}.$$(5)

Note that the number of spectral lines is accounted for via the square root of N, which reflects a weighting by its Poisson uncertainty. While in theory the usage of Eq. (5) would facilitate the ideal target selection, in practice the number of lines, N, as well as the planetary geometric albedo are difficult to estimate or unknown for many target systems. We therefore commenced our target selection by applying the simplified metric $${\mathrm{R}\mathrm{S}\mathrm{M}}^{\mathrm{\prime }}={10}^{-m_V/5}{\left(\frac{R_{\mathrm{p}}}{a}\right)}^2{K}_{\mathrm{p}}$$(6)

to all planets within the NASA exoplanet archive (Akeson et al. 2013). Figure 1 shows the declination of the most favourable targets as a function of their respective metric score. As such, the candidates most favourable for a reflection signature search as suggested by the simplified metric can be picked depending on the chosen observation site. To illustrate this point, the latitudes of two high-resolution spectrographs, CARMENES (Quirrenbach et al. 2014, 2016) at the Calar Alto Observatory and ESPRESSO (Pepe et al. 2021) at the Very Large Telescope on Cerro Paranal, are indicated by dashed horizontal lines in Fig. 1. Since a given target culminates in the zenith, and thus under optimal viewing conditions, when its declination coincides with the observer’s latitude, picking targets close to the horizontal dashed line corresponding to the instrument one is planning to use will maximise the achievable S/N.

From Fig. 1 it is evident that v And b, τ Bootis b, HIP 65 A b, WASP-33 b, and KELT-9 b represent the most suitable target systems for a reflection signature search. Of these, apart from achieving the lowest metric score, v And b and τ Bootis b do not show transits. Restricting the target selection to transiting planets only ensures precise knowledge of the planet’s orbital solution, which will make the reflection signature search that is to follow considerably simpler. HIP 65 A b on the other hand, is far removed from the optimal declination range for CARMENES, the instrument picked for conducting this case study. While WASP-33 b achieves the highest metric score, its host star exhibits intense pulsations (Collier Cameron et al. 2010; Smith et al. 2011), which can severely complicate the data reduction and interpretation. Cont et al. (2021) comprehensively discuss the effects of this star’s pulsation on cross-correlation procedures. Caution therefore suggests picking KELT-9 b as the case study’s target. The planet’s and its host star’s characteristics are presented in Table 1. KELT-9 b is a prominent specimen among the ultra-hot Jupiter class of exoplanets in that it exhibits the hottest known day-side temperature of ~4600 K and possesses an extended hydrogen envelope of a size comparable to the Roche lobe (Yan & Henning 2018). Hooton et al. (2018) set and upper limit of 0.14 on its geometric albedo in the U band.

Fig. 1 Declination of the 50 targets presenting the largest RSM′ values as a function of their respective score. The ten most favourable targets are annotated with their names. The areas of the data points scale with the planetary radii. For reference, an exemplary data point corresponding to 1 RJup is shown in the bottom right. The colour of a given data point indicates the respective planetary equilibrium temperature. Since non-transiting planets do not have any measured equilibrium temperatures listed in the NASA exoplanet archive (Akeson et al. 2013), they are shown as empty data points. Planets that had previously been targeted for reflection signature searches are marked by dotted crosshairs (see text for references). The latitudes of two exemplary high-resolution spectrographs, CARMENES and ESPRESSO, are visualised as horizontal dashed lines, the declination ranges within which targets culminate at airmasses lower than 1.5 are shaded in the respective colour.

4 Observations and data reduction

This study used observational data obtained by the CARMENES spectrograph at the Calar Alto Observatory during the night of 30 to 31 August 2021 (PI: Fei Yan). They consist of 43 exposures, each lasting 200 s and spanning orbital phases between φ = 0.58 and 0.66, where we define the orbital phase to assume the value 0 at the primary transit midpoint increasing to 1 during the course of a full orbit. The orbital coverage is visualised in Fig. 2.

While CARMENES spectra are simultaneously obtained in a visible channel (VIS; 0.52–0.96 μm) and a near-infrared (NIR) channel (0.96–1.71 μm), in this study we focussed on the VIS channel only since, compared to the NIR channel, is not as troubled by telluric contamination and covers a wavelength domain in which the stellar emission and, therefore, the expected reflection signal strength is higher. Due to the cross-dispersive instrument design, the VIS spectrum is made up of 61 distinct but overlapping one-dimensional spectral orders. These were extracted from the raw frames using the reduction pipeline caracal v2.20 (Zechmeister et al. 2014; Caballero et al. 2016). The exposures were corrected for telluric absorption using molecfit (Smette et al. 2015; Kausch et al. 2015). During the initial assessment of the dataset it was noticed that the S/N ratio of exposure 30 far exceeded the rest of the time series. A more thorough investigation revealed an additional broad peak in the spectrum centred around 6620 Å. This was shown to have been caused by a red LED that briefly lights up the pinhole entry of the science fibre to enable the re-acquisition of the signal in cases when the link between the light beam and pinhole is lost. Evidently, this automated procedure was triggered during the exposure in question and contaminated the recorded science spectrum. We decided to discard the exposure, thereby incurring a loss of 2.3% of the total data.

Fig. 2 Orbital coverage of the CARMENES dataset used in this study. The spectra were taken on 30 August 2021 between orbital phases of 0.58 and 0.66, here indicated in orange. The secondary eclipse range around φ = 0.5 is indicated in grey.

5 Analysis and results

5.1 Search for reflection signatures

To assess whether reflection signatures are present in the CARMENES dataset we used high-resolution Doppler spectroscopy, a method that exploits the time-varying Doppler shift experienced by an exoplanet’s signature to disentangle it from the stationary imprint of the Earth’s atmosphere and the signal of the stellar host (e.g. Snellen et al. 2010; Rodler et al. 2012; Brogi et al. 2012; Birkby 2018). While commonly used to search for planetary signatures in transmission spectra or emission signals originating from the planet’s dayside (see Appendix B for a detection of iron emission in the CARMENES dataset), it is also readily applicable to probe for reflected light.

Processing the CARMENES dataset into a suitable format to apply these techniques to, required several reduction steps. To begin with, the relative motion between the observatory and the target system, consisting of the motion of the observatory around the Solar System barycentre, the star’s radial velocity variation induced by the orbiting planet, and the target’s systemic radial velocity of (−20.567 ± 0.100) km/s (Gaudi et al. 2017), was then corrected for. To streamline the process and avoid unnecessary re-binning steps, we corrected for both sources of relative motion simultaneously. To ensure compatibility, the shifted spectra were sampled onto the original wavelength solution. This aligned the different exposures throughout the night to the same frame of reference, the stellar frame.

Next, we median-combined the exposures to produce a high S/N stellar template. This approach ensures that any faint planetary signal (if present at detectable levels) is effectively blurred out due to the planet’s orbital radial velocity, which shifts its signal to different wavelengths in each exposure. As a result, the stellar template provides an accurate representation of the star’s direct light, surpassing the fidelity of a model spectrum or that of an analogue star by capturing the unique characteristics and potential idiosyncrasies of the specific star being observed.

Subtraction of the stellar template from the individual exposures throughout the night yielded residual exposures, effectively lifting the veil of stellar glare from each spectrum. Three arbitrary residual orders are visualised in Fig. A.1. Crucially, these residuals may also harbour the faint planetary signal, awaiting detection. To remove any lingering low-frequency continuum signal that might result from an inadequate removal of the stellar component, we applied a high-pass filter to the residual exposures.

These were subsequently transformed from the stellar frame into the planetary frame by applying the orbital radial velocity at each mid-exposure time to reverse the Doppler shift affecting the planetary signal. It is important to note that this transformation is only possible when the planet’s orbital parameters are known, which was not the case for the initial studies of the τ Bootis system mentioned in Sect. 1. After applying the required shift to the wavelength axis of each exposure, we resampled the resulting planetary frame residual exposures onto the original wavelength solution.

Having thus processed the data into a suitable format, we still required a template that describes the reflection spectrum to compare with the residual exposures. As described in Sect. 2, depending on the rotational velocity derived in Eq. (4), we expect the reflected component to present a different degree of rotational broadening as compared to the direct stellar component. The KELT-9 system offers an intriguing case regarding the stellar rotational velocity as observed from the planet. This is due to the fact that the orbit of KELT-9 b is severely misaligned with respect to the rotational plane of its host. The angle between the planetary orbit’s normal and the star’s rotational axis is measured to be (−84.8 ± 0.3) deg (Stephan et al. 2022). The consequence of this peculiar configuration is that the rotational velocity of KELT-9 as seen from KELT-9 b varies in the course of its orbit. Correspondingly, the rotational broadening experienced by the reflection spectrum will vary throughout the planetary orbit. When traversing the rotational plane, the lines in the reflection spectrum will be broadened according to the full equatorial rotational velocity of the star, v_rot = (111.4 ± 1.3) km/s (Gaudi et al. 2017). As the planet nears the rotational poles, the intensity of line broadening is expected to diminish. This unique configuration makes the KELT-9 system an ideal testbed for investigating the effects of rotational line broadening on the reflection spectrum of an exoplanet.

Due to the varying stellar rotational velocity perceived by the planet, the foundation of a reflection template adequate for the KELT-9 system must lie in the an un-broadened stellar spectrum. Consequently, the stellar template used earlier to subtract the direct stellar component and obtain the residual exposures is unsuitable for this purpose. Instead, we resorted to spectral models generated without incorporating any rotational broadening. To mimic KELT-9, we used a PHOENIX model spectrum (Husser et al. 2013), assuming an effective temperature of T_eff = 10 200 K, surface gravity of log(g) = 4.0, metallicity of [Fe/H] = 0, and an alpha element abundance of [α/Fe] = 0. To adapt the PHOENIX spectrum into a reflection template customised for our analysis, we chunked it into the CARMENES orders and scaled it by the planet-star contrast, ε, defined in Eq. (1). Since KELT-9 b’s phase function, Φ(φ), varies over time, ε is expected to vary for each exposure. For a first-order approximation, both A_g and Φ(φ) were set to 1, representing a perfectly reflective, fully illuminated planetary dayside as viewed from Earth. Under this assumption, ε simplifies to the planet’s solid angle seen from the star, that is, ε ≈ (R_p/a)² ≈ 7 × 10⁻⁴. Finally, the template was convolved with the CARMENES instrument profile. These steps transformed the generic PHOENIX spectrum into our un-broadened reflection template shown in Fig. 3. In addition, we constructed a binary reflection template by identifying the absorption line positions within the PHOENIX spectrum. For each wavelength channel containing an absorption line, a relative flux value of −1 was assigned in the binary template spectrum, while all other wavelength channels were given a relative flux value of 0.

At this stage, we did not account for the rotational broadening, and the templates were constructed assuming the most favourable geometric albedo as well as a phase function fixed at one. As such, attempting to detect these templates in the residual exposures constitutes the simplest possible implementation of a reflection signature search.

To this end we employed a normalised and uncertainty-weighted cross-correlation technique, which allowed us to combine the numerous spectral lines into – in the ideal case – a single, prominent peak, thereby enhancing the detectability of any signal present in the data. To achieve this, we computed the cross-correlation function (CCF) between the residual exposures and the reflection template that was Doppler shifted along a grid of radial velocity deviations, Δv, ranging from −300 to 300 km/s with a step size of 1 km/s, according to $$\mathrm{C}\mathrm{C}\mathrm{F}=\frac{\sum _i{\alpha }_i{r}_i{t}_i(\mathrm{\Delta }v)}{\sum _i{({\alpha }_i{r}_i{t}_i(\mathrm{\Delta }v))}^2}.$$(7)

Here, r denotes the residual exposures in the planetary frame, σ_r is the uncertainty array associated with the exposure, and t represents the reflection template spectrum shifted by Δv in radial velocity space. Here, r_i and t_i denote the residual and template flux in wavelength channel i, respectively, while α_i represents the applied weights. Defining these as ${\alpha }_i={\sigma }_{r_i}^{-1}$, where σ_ri is the residual flux error in a given wavelength channel, we ensured that maximising the CCF corresponds to maximising the likelihood (Lockwood et al. 2014).

Collapsing the CCFs resulting from the individual residual exposures yielded a one-dimensional CCF representing signal content in the entire dataset. Recall that prior to the cross-correlation we converted the residual spectra from the stellar to the planetary frame by correcting for the orbital motion of the planet. This operation was based on the planet’s known K_p value. Repeating this frame conversion and the subsequent cross-correlation procedure for a series K_p values covering an interval of 600 km/s around the true K_p value at a step size of 1 km/s we built up a two-dimensional CCF map that is spanned by the Δv and K_p grids. By dividing each K_p row by its standard deviation under the exclusion of a 100 × 100 km/s region around the position at which the true signal was expected, we normalised the CCF map to obtain a S/N map. The two panels of Fig. 4 show these S/N maps when using the PHOENIX and binary template, respectively. There is no peak discernible at the expected semi-amplitude and radial velocity deviation in either map. The planet’s reflection signature could therefore not be detected within the CARMENES dataset, using these templates. Both panels in Fig. 4 exhibit an anti-correlation feature at K_p = 0 km/s. This artefact is most likely caused by the imperfect removal of the direct stellar light component.

Fig. 3 Reflection template of KELT-9 b. The inset in the lower right shows a close-up of the template compared to two rotationally broadened versions of the model.

Fig. 4 S/N map produced by cross-correlating the spectral residuals with the reflection templates as a function of the planetary radial velocity semi-amplitude, Kp, and the radial velocity deviation, Δv. The top and bottom maps present the results when using the PHOENIX-based and binary reflection template, respectively. The dashed black lines indicate where a genuine reflection signal would be expected to manifest.

5.2 The effect of rotational line broadening

In the above search for reflection signatures we assumed a high geometric albedo of 1, as well as a constant phase function, Φ = 1, for all exposures. The resulting non-detection prompted us to re-evaluate these first order assumptions. An injection study can yield upper limits on the albedo of the planet as well as shed light on the effect of the rotational line broadening and the sensitivity of the search method to varying planetary reflectivity and orbital dynamics, refining our understanding of the detectability of such signals given our data.

To this end, we computed a momentary radial velocity resulting from the fictitious semi-amplitude K_inj = K_p/2 = 126.5 km/s during every exposure. Here, K_p denotes the true semi-amplitude of the planet that results from K_p = (M_s/M_p)K_s, where M_s and M_p are the stellar and planetary masses, respectively, and K_s = (276 ± 79) m/s is the semi-amplitude of KELT-9 (Gaudi et al. 2017). The use of a fictitious semi-amplitude avoids confusion between the injected signal and any residual genuine signal present at the true semi-amplitude of the planet. The reflection template was shifted according to the velocity corresponding to each exposure, scaled by an injection factor, I, and added to the telluric-corrected but not yet normalised dataset. The full injection routine is thus described by $$f_{k,i}^{\mathrm{i}\mathrm{n}\mathrm{j}}=f_{k,i}+(I\cdot t_{k,i}).$$(8)

Here, f_k,i is the relative flux in wavelength channel k and exposure i while t_k,i denotes the same for the reflection template exposure i. Accordingly, $f_{k,i}^{\text{inj}}$ is the relative flux that results from the injection procedure. The injection factor, I, is given by the product of the geometric albedo assumed in the injection run and the phase function, Φ, which changes for each exposure, i. We used $$\mathrm{\Phi }(\phi )=\frac{1+\cos (\pi +2\pi \phi )}{2}$$(9)

to parametrise Φ as a function of the orbital phase, φ. This formulation assumes a simplified phase curve, approximating the reflected light intensity as a cosine function of the orbital phase. Unlike the Lambertian phase function (e.g. Madhusudhan & Burrows 2012), which exhibits a steeper increase in brightness near full phase due to isotropic scattering, this approach provides a less elaborate but nonetheless effective first-order approximation of the planet’s phase curve. The choice of this parametrisation simplifies the analysis by avoiding the additional complexities of integrating the Lambertian phase function, which would require precise modelling of the planet’s scattering properties and reflectivity.

To explore the detectability of planetary signals, we gradually decreased the albedo, varying A_g from 1 to lower values. This approach mimics the injection of increasingly darker planets, allowing us to investigate the limits of the data in constraining the planet’s albedo. As was shown in Eq. (4), the reflected spectrum is subject to rotational line broadening resulting from the combination of the perceived stellar rotation and the planetary rotation velocity. To include this effect in our analysis we repeated the injection procedure for different net rotational velocities with which we broadened the spectral lines of the injected reflection signal as well as the template used for cross-correlation. While this assumes the same rotational broadening for the two, such an approach simplifies the analysis by circumventing the need for a complex model that accounts for the variation of rotational broadening during the planet’s orbital motion. Such modelling would require precise knowledge of the relative latitudinal stellar rotation and its impact on the reflected spectrum, which is currently unfeasible without introducing artefacts. Instead, this approach serves as a preliminary investigation to explore the impact of rotational broadening and the possibility to derive an upper limit for the geometric albedo, A_g, given the limited signal-to-noise of our data. The S/N maps resulting from each injection run are shown in Fig. 5.

A note on the computation of the S/N at the injection position: the underlying S/N map resulting from cross-correlating the original (non-injection-treated) dataset with the reflection template exhibits autocorrelation features (see Fig. 4). Thus, the base S/N level at the injection position must not be zero. To obtain the corrected signal-to-noise ratio, S/N_c, with which the injected signal can be recovered, we corrected for the base S/N value by subtracting it from the value encountered at the injection position. These are the values reported in the corner of the panels that make up Fig. 5. Furthermore, it should also be noted that the simple CCF to S/N conversion employed here can break down for strong rotational broadening (Parker et al. 2024). These caveats imply that the approximate S/N maps depicted in Fig. 5 should not be read as a formal statistical significance measure.

Running such injection recoveries on a more extensive grid spanned by a set of geometric albedo values (steps of 0.1 between 0 and 1) as well as by several rotational velocities (0, 1, 5, 10 and 50 km/s) before linearly interpolating between the recovered S/N values, we obtained the two-dimensional S/N plane visualised in Fig. 6. Here, we have also indicated where the planet’s genuine reflection signature can be expected. This results from the upper limit on the geometric albedo, A_g < 0.14, reported in Hooton et al. (2018), and the rotational velocity of the star as seen by the planet during the CARMENES exposures. These velocities can be approximated by assuming the projected rotational velocity to vary according to $$v_{\mathrm{r}\mathrm{o}\mathrm{t},i}={\hat{v}}_{\mathrm{r}\mathrm{o}\mathrm{t}}\left(\frac{1+\cos (4\pi \cdot {\phi }_i)}{2}\right).$$(10)

Here, v_rot,i is the momentary projected rotational velocity, φ_i is the orbital phase during the exposure i and ${\hat{v}}_{\text{rot}}$ = (111.4 ± 1.3) km/s is the maximum rotational velocity (Gaudi et al. 2017). This parametrisation implies that the projected rotational velocity vanishes when the planet is located above the stellar rotational pole at φ = 0.25 and 0.75. Consequently, it assumes its maximum value, ${\hat{v}}_{\text{rot}}$, when traversing the stellar rotational plane at φ = 0 and 0.5. Strictly speaking, the system’s spin-orbit-misalignment measures (−84.8 ± 0.3) deg (Stephan et al. 2022) instead of a clean 90 deg, however.

6 Discussion

The application of the simplified RSM as defined in Eq. (6) suggested selecting KELT-9 as the target system for a reflection study. Its significant spin-orbit-misalignment combined with the rapid relative rotation of the host negatively impact the system’s favourability, however. Apart from these factors, the inclusion of a spectral type weighting in the RSM seems advisable, since early-type stars, such as KELT-9, exhibit fewer spectral lines and consequently, their reflection signatures will manifest at lower S/N when applying cross-correlation routines compared to later-type stars. Overall, to prevent future investigations from falling victim to the same pitfalls, these experiences motivate a future refinement of the RSM. A modified version of the metric, including the stellar rotational velocity, spin orbit alignment as well as a weighing factor that incorporates the star’s spectral type, is left for a follow-up study. In summary, KELT-9 b is an unfavourable target for reflection studies.

Our simple search for an un-broadened reflection signature resulted in a non-detection. The primary reasons for this outcome are the low albedo of KELT-9 b and the significant rotational broadening of its spectrum. The former means the reflection signal is too weak relative to the noise of the CARMENES dataset to be confidently detected. The latter severely reduces the effectiveness of the cross-correlation technique at detecting faint signals. A secondary reason for the non-detection is the looseness and consequent inadequacy of the PHOENIX model used as a basis for the template building.

The residual orders shown in Fig. A.1 provide insight into the adequacy of the stellar template that was used to remove the direct stellar light component. While they are noise-dominated, they also retain some low-frequency structure that likely stems from imperfect stellar or telluric removal. These artefacts can contribute to the autocorrelation visible in Figs. 4 and 5 and further reduce our sensitivity to genuine reflection signatures. The imperfect removal of the direct stellar light likely also gives rise to the anti-correlation feature around K_p = 0 km/s in Fig. 4.

The effect altering the reflection spectrum most strongly is the rotational broadening of its lines. To assess its significance on reflection studies, we performed a series of injection recoveries. Figure 5 demonstrates visually how the detection of the injected signal depends on the planet’s geometric albedo and the rotational velocity. Both dependences are to be expected intuitively. The local peaks visible in the S/N maps in Figs. 4 and 5 are confined to amplitudes below 2 and are thus not significant. Figure 6 helps make sense of why the planet’s genuine reflection signal could not be recovered. When comparing the parameter space where a tentative detection can be realised to the region where the planet’s signal is expected to reside, it comes as no surprise that the search ultimately resulted in a non-detection. These findings are of consequence for reflection studies targeting systems in which the planet perceives a non-negligible stellar rotational velocity.

The favourability of KELT-9 b as the target for a reflection study is further reduced by the fact that its high temperature results in significant thermal emission within the CARMENES VIS range. As shown in Appendix B, the search for emission signatures relies on different spectral templates and probes distinct signal morphologies as compared to the reflection study presented here.

For the sake of simplicity, it is advisable to restrict oneself to target systems in which the planetary orbital period is equal or close to the stellar rotation period. Misaligned systems such as KELT-9 represent a special case, however: on the one hand the variable rotational broadening experienced by the reflected spectrum constitutes a challenge to model, on the other hand – in the case of near polar orbits – they can offer periods in which the reflected spectrum is not or only barely affected by rotational line broadening. In this sense, misaligned systems offer a chance at increased detection possibilities during polar passage.

While rotational broadening significantly affects the detectability of reflected light, another relevant factor is the wavelength-dependent reflectivity of the planet. The geometric albedo of an exoplanet varies due to molecular absorption, Rayleigh scattering, and the presence or absence of clouds or hazes. These factors shape the planet-star contrast in high-resolution spectroscopy and play a crucial role in determining whether a reflection signal can be detected (e.g. Hu 2019; Damiano & Hu 2020, 2022). For hot Jupiters, strong alkali absorption (e.g. Na and K) tends to suppress reflectivity in the optical, leading to low albedos (often <0.1; e.g. Demory 2014). In contrast, cooler planets with high-altitude clouds can exhibit much higher albedos (~0.3–0.6) in the visible and near-infrared due to enhanced scattering from condensates (Damiano & Hu 2022; Sudarsky et al. 2000). The spectral region where reflection spectroscopy is most effective depends not only on the planet’s albedo spectrum but also on the host star’s spectral energy distribution; for Sun-like stars, the visible band is optimal, whereas for M-dwarf hosts, the infrared is more favourable (e.g. Kitzmann et al. 2011). In light of these dependences, future studies may benefit from prioritising planets with significant cloud coverage, where higher albedos could improve detectability, or focussing on spectral regions where Rayleigh scattering or cloud reflection enhances the contrast with the host star.

Assuming the reflection template is a valid and sufficiently precise description of the reflected spectral component, the parameter space where tentative detections are achievable (the top-left corner in Fig. 6) can be enlarged in two ways: by collecting more data that cover orbital phases where the reflected component is well distinguishable from the direct stellar light, or by basing the search on spectra possessing a higher S/N ideally obtained through larger telescopes. In any case, the large separation between the sensitivity region of our data and the expected signal position implies that it will be exceedingly difficult to achieve overlap between the two for KELT-9 b.

Fig. 5 Injection-recovery tests of KELT-9 b for different geometric albedo, Ag, and net rotational velocity, vrot, values. The panels show S/N maps resulting from the cross-correlation of the reflection template with the CARMENES dataset, into which the respective artificial reflection signals were injected beforehand. The rotational velocity used to broaden the injected signal as well as the template increases towards the right. The geometric albedo with which the artificial signal and the template were scaled drops off towards the bottom. The positions at which the artificial signal was injected are indicated by the dashed black lines. The S/N value encountered at this position is shown in the white box in the top-left corner of each respective panel. Note that this value corresponds to the corrected S/N (hence the c subscript), that is, the signal strength normalised by the encountered signal strength at the same position in the absence of an injected signal.

Fig. 6 Sensitivity region of a reflection study applied to the CARMENES dataset. The blue colour map shows the corrected S/N as a function of the planet’s geometric albedo as well as the velocity used for the rotational line broadening of the artificial signal. The dashed white lines indicate different corrected S/N contours, with the S/Nc = 5 contour delineating the region where a robust detection based on the CARMENES data at hand could be realised. The black dots mark grid points at which injection-recovery tests were performed. The hatched regions show where the reflection signal is expected to manifest and are defined by the upper limit on the planet’s geometric albedo (Ag < 0.14; Hooton et al. 2018) as well as the stellar rotational velocity as perceived by the planet during the CARMENES exposures (note that this vrot region extends to higher values than shown here, but the corresponding S/N achievable using the CARMENES dataset vanishes completely). The genuine reflection signal of the planet can be assumed to be located where the two hatched regions overlap.

7 Conclusion

This study introduces the RSM, which can be used to evaluate the suitability of planetary systems for high-resolution reflection studies, offering a valuable tool for prioritising targets in future research. While KELT-9 was identified as a promising target, the spin-orbit misalignment and significant rotational broadening of its reflection spectrum presented unforeseen challenges, leading to a non-detection for KELT-9 b. These findings underscore the critical importance of considering rotational broadening effects in target selection and highlight the need for refined metrics that account for spectral complexities.

Looking ahead, advancements in instrumentation, such as the advent of ESO’s Extremely Large Telescope and next-generation high-resolution spectrographs, will enhance sensitivity and resolution, potentially enabling the first robust detections of reflected light in high-resolution spectra. Coupled with improved metrics and methodologies, these tools will open new frontiers in exoplanet research, including the study of temperate planets and their atmospheric properties.

Reflection spectroscopy remains a promising avenue for probing exoplanetary atmospheres, and continued investment in data quality, target selection strategies, and analytical frameworks will ensure its growth into a cornerstone technique in the coming decades.

Acknowledgements

This work constitutes a Master’s thesis project completed in 2022 at Landessternwarte Königstuhl. E.N. acknowledges the support from the Deutsches Zentrum für Luft- und Raumfahrt (DLR, German Aerospace Center) - project number 50OP2502. CARMENES is an instrument at the Centro Astronómico Hispano en Andalucía (CAHA) at Calar Alto (Almería, Spain), operated jointly by the Junta de Andalucía and the Instituto de Astrofísica de Andalucía (CSIC). The authors wish to express their sincere thanks to all members of the Calar Alto staff for their expert support of the instrument and telescope operation. CARMENES was funded by the Max-Planck-Gesellschaft (MPG), the Consejo Superior de Investigaciones Científicas (CSIC), the Ministerio de Economía y Competitividad (MINECO) and the European Regional Development Fund (ERDF) through projects FICTS-2011-02, ICTS-2017-07-CAHA-4, and CAHA16-CE-3978, and the members of the CARMENES Consortium (Max-Planck-Institut für Astronomie, Instituto de Astrofísica de Andalucía, Landessternwarte Königstuhl, Institut de Cirències de l’Espai, Institut für Astrophysik Göttingen, Universidad Complutense de Madrid, Thüringer Landessternwarte Tautenburg, Instituto de Astrofísica de Canarias, Hamburger Sternwarte, Centro de Astrobiología and Centro Astronómico Hispano-Alemán), with additional contributions by the MINECO, the Deutsche Forschungsgemeinschaft (DFG) through the Major Research Instrumentation Programme and Research Unit FOR2544 “Blue Planets around Red Stars”, the Klaus Tschira Stiftung, the states of Baden-Württemberg and Niedersachsen, and by the Junta de Andalucía. We acknowledge financial support from the Agencia Estatal de Investigación (AEI/10.13039/501100011033) of the Ministerio de Ciencia e Innovación and the ERDF “A way of making Europe” through projects PID2022-137241NB-C4[1:4], PID2021-125627OB-C31, RYC2022-037854-I, and the Centre of Excellence “Severo Ochoa” and “María de Maeztu” awards to the Instituto de Astrofísica de Canarias (CEX2019-000920-S), Instituto de Astrofísica de Andalucía (CEX2021-001131-S) and Institut de Ciències de l’Espai (CEX2020-001058-M). This work was also funded by the LMU München Fraunhofer-Schwarzschild Fellowship and DFG under Germany’s Excellence Strategy (EXC 2094–390783311), Generalitat de Catalunya/CERCA programme, and the National Natural Science Foundation of China (grant no. 42375118). We thank the Calar Alto Observatory for allocation of director’s discretionary time to this programme (DDT.A21.280). We thank the anonymous referee for their thorough and insightful feedback and suggestions.

References

Appendix A Residual orders

Fig. A.1 Main panels: Arbitrary residual orders resulting from the subtraction of the stellar template from CARMENES exposure 10. Right panels: Distributions of the residual flux marginalised over all wavelength channels in the respective orders, with the dotted grey lines indicating a residual relative flux of zero.

Appendix B Iron emission detection

Before investigating the planetary reflection signal, we intended to assess the feasibility of detecting the spectral signature of KELT-9 b in our CARMENES data. Prior work has shown that the thermal emission spectrum of KELT-9 b is dominated by a dense forest of atomic iron (FeI) lines (e.g. Pino et al. 2020, 2022; Kasper et al. 2021). Therefore, we targeted the emission lines from atmospheric FeI applying the same procedures as described in previous work. This procedure differs somewhat from the one performed in Sect. 5.1, where we employed a simpler template based on a copy of the stellar spectrum. For an in-depth description of the analysis methods used here, we refer the reader to the literature (e.g. Cont et al. 2022a,b; Yan et al. 2023).

For each spectral order, the spectra were arranged chronologically to form a two-dimensional spectral matrix. The individual spectra were then fitted with a second-order polynomial and normalised using the resulting fit function. This process also included the removal of outliers and the masking of saturated telluric absorption lines. To address contamination from telluric and stellar lines, the detrending algorithm SYSREM (Tamuz et al. 2005) was applied. Ten consecutive iterations of the algorithm were performed, resulting in a residual spectral matrix for each iteration. The planet’s atmospheric signature is expected to be hidden within the noise of the residual spectra. Using the radiative transfer code petitRADTRANS (Mollière et al. 2019) we computed a template of the emission spectrum emerging from KELT-9 b’s dayside atmosphere. The atmospheric structure assumed for these calculations covers a pressure range between 10⁻¹⁰ and 10² bar and consists of 130 layers. We approximated the temperature-pressure (T-p) profile and the FeI volume mixing ratio with the values used in the work on KELT-20 by Johnson et al. (2023). This T-p profile is defined by a low-pressure point (T₁, p₁) and a high-pressure point (T₂, p₂). Between the two pressure points the temperature changes linearly with log p; outside this range the atmosphere is assumed to be isothermal. The two points were set as (T₁, p₁) = (5000 K, 10⁻³ bar) and (T₂, p₂) = (1800 K, 10^−1.3 bar), implying an atmospheric thermal inversion layer. We assumed a FeI volume mixing ratio of 5.4 × 10⁻⁵. The FeI opacities needed for the radiative transfer calculation were obtained from the Kurucz line database (Kurucz 2018). The resulting petitRADTRANS template spectrum was divided by the stellar flux, normalised to the continuum level, and convolved with the instrument profile. The left and middle panels of Fig. B.1 show the used T-p profile and the resulting planet-to-star flux ratio (F_p/F_s).

By means of the cross-correlation procedure briefly presented in Sect. 5.1, we successfully extracted the faint FeI emission lines originating from KELT-9 b’s atmosphere. To this end, we shifted the template spectrum over a radial velocity deviation, Δv, interval ranging from −520 km/s to 520 km/s with steps of 1 km/s. We arranged the resulting CCFs of each echelle order in an individual two-dimensional array. In the following, the arrays from the different echelle orders were co-added, yielding the final CCF map. This CCF map was aligned over a range of different K_p values. The tested K_p values covered a range between 0 km/s and 400 km/s with steps of 0.5 km/s. For each individual K_p the mean value of the shifted two-dimensional CCF map was computed along its time evolution axis, which resulted in a one-dimensional CCF. Subsequently, the one-dimensional CCFs that originated from alignments with different K_p values were arranged in a two-dimensional array. Normalising by the maps standard deviation under the exclusion of the region around the strongest signal peak yields the S/N map of the FeI emission lines shown in the right hand panel of Fig. B.1.

The maximum FeI detection significance is reached after three SYSREM iterations with a S/N of 5.2. The detected orbital parameters are K_p = ${251.5}_{-69.5}^{+27.5}$ km/s and Δv = $8_{-59}^{+21}$ km/s, which is in line with the expected values. We conclude that the quality of the CARMENES observational data used in this work is sufficient to investigate the planetary reflection signal.

Fig. B.1 Left: Normalised FeI emission spectrum emerging from the dayside atmosphere of KELT-9 b as modelled with petitRADTRANS. Middle: Adopted T-p profile used for generating the emission template spectrum. Right: S/N map of FeI thermal emission from KELT-9 b’s dayside atmosphere. The signal peaks with a S/N of 5.2 after three consecutive SYSREM iterations. We indicate the coordinates of the strongest signal with the dashed white lines. Cross-sections of the S/N peak are shown in the horizontal and vertical panels.

All Tables

All Figures

Fig. 1 Declination of the 50 targets presenting the largest RSM′ values as a function of their respective score. The ten most favourable targets are annotated with their names. The areas of the data points scale with the planetary radii. For reference, an exemplary data point corresponding to 1 RJup is shown in the bottom right. The colour of a given data point indicates the respective planetary equilibrium temperature. Since non-transiting planets do not have any measured equilibrium temperatures listed in the NASA exoplanet archive (Akeson et al. 2013), they are shown as empty data points. Planets that had previously been targeted for reflection signature searches are marked by dotted crosshairs (see text for references). The latitudes of two exemplary high-resolution spectrographs, CARMENES and ESPRESSO, are visualised as horizontal dashed lines, the declination ranges within which targets culminate at airmasses lower than 1.5 are shaded in the respective colour.

In the text

	Fig. 2 Orbital coverage of the CARMENES dataset used in this study. The spectra were taken on 30 August 2021 between orbital phases of 0.58 and 0.66, here indicated in orange. The secondary eclipse range around φ = 0.5 is indicated in grey.
In the text

	Fig. 3 Reflection template of KELT-9 b. The inset in the lower right shows a close-up of the template compared to two rotationally broadened versions of the model.
In the text

Fig. 4 S/N map produced by cross-correlating the spectral residuals with the reflection templates as a function of the planetary radial velocity semi-amplitude, Kp, and the radial velocity deviation, Δv. The top and bottom maps present the results when using the PHOENIX-based and binary reflection template, respectively. The dashed black lines indicate where a genuine reflection signal would be expected to manifest.

In the text

Fig. 5 Injection-recovery tests of KELT-9 b for different geometric albedo, Ag, and net rotational velocity, vrot, values. The panels show S/N maps resulting from the cross-correlation of the reflection template with the CARMENES dataset, into which the respective artificial reflection signals were injected beforehand. The rotational velocity used to broaden the injected signal as well as the template increases towards the right. The geometric albedo with which the artificial signal and the template were scaled drops off towards the bottom. The positions at which the artificial signal was injected are indicated by the dashed black lines. The S/N value encountered at this position is shown in the white box in the top-left corner of each respective panel. Note that this value corresponds to the corrected S/N (hence the c subscript), that is, the signal strength normalised by the encountered signal strength at the same position in the absence of an injected signal.

In the text

Fig. 6 Sensitivity region of a reflection study applied to the CARMENES dataset. The blue colour map shows the corrected S/N as a function of the planet’s geometric albedo as well as the velocity used for the rotational line broadening of the artificial signal. The dashed white lines indicate different corrected S/N contours, with the S/Nc = 5 contour delineating the region where a robust detection based on the CARMENES data at hand could be realised. The black dots mark grid points at which injection-recovery tests were performed. The hatched regions show where the reflection signal is expected to manifest and are defined by the upper limit on the planet’s geometric albedo (Ag < 0.14; Hooton et al. 2018) as well as the stellar rotational velocity as perceived by the planet during the CARMENES exposures (note that this vrot region extends to higher values than shown here, but the corresponding S/N achievable using the CARMENES dataset vanishes completely). The genuine reflection signal of the planet can be assumed to be located where the two hatched regions overlap.

In the text

	Fig. A.1 Main panels: Arbitrary residual orders resulting from the subtraction of the stellar template from CARMENES exposure 10. Right panels: Distributions of the residual flux marginalised over all wavelength channels in the respective orders, with the dotted grey lines indicating a residual relative flux of zero.
In the text

Fig. B.1 Left: Normalised FeI emission spectrum emerging from the dayside atmosphere of KELT-9 b as modelled with petitRADTRANS. Middle: Adopted T-p profile used for generating the emission template spectrum. Right: S/N map of FeI thermal emission from KELT-9 b’s dayside atmosphere. The signal peaks with a S/N of 5.2 after three consecutive SYSREM iterations. We indicate the coordinates of the strongest signal with the dashed white lines. Cross-sections of the S/N peak are shown in the horizontal and vertical panels.

In the text