2.1 Data reduction of the raw spectra

Raw spectra were reduced using the custom pipeline of Krüger et al. (in preparation), which applies minimal preprocessing (i.e. no barycentric correction, normalisation or order merging). Bias was removed from each quadrant via overscan regions, and pixel-to-pixel sensitivity variations were corrected using a normalised master flat constructed from all fibre flats taken that night (0.1 s for Rosso, 1 s for Verde and 60 s for Azzurro). Scattered light subtraction was omitted, as it was negligible for our target. Spectral extraction proceeded by summing flux in the cross-dispersion direction around predefined traces modelled as fifth-order polynomials. Initial trace positions and dispersion-range definitions were taken from the Th traces in the Veloce Manual (Veloce Collaboration 2025). To mitigate edge effects, the summing aperture was extended by three pixels beyond the flux drop-off observed in the fibre flats.

Wavelength calibration – common to all exposures – was derived from ThXe lamp spectra extracted at the start of the night. Air wavelengths of Th lines from the NIST database (Kramida et al. Reference Kramida and Ralchenko2024) served as the absolute reference. An initial guess for line positions employed the converted-to-air solution provided by Chris Tinney in the Veloce Manual (Veloce Collaboration 2025). Pixel positions of Th lines were then determined by fitting generalised Gaussians to features nearest these guesses above a relative intensity threshold. We fitted a seventh-degree polynomial surface to a data cube comprising pixel positions, absolute order numbers and air wavelengths scaled by order number, with iterative outlier rejection via 3 $\sigma$ clipping of residuals. This solution was applied across the pixel range of each extracted order, yielding wavelength-calibrated spectra ready for CCF analysis. Figure A1 presents the resulting spectra for all three arms, plotted as wavelength against the square root of flux to emphasise variations in signal strength across the bandpass.

2.2 1D order stitching

We first cleaned each extracted 2D spectral order using a two-pass, sliding-window $\sigma$ –clip: a 51-pixel median filter was applied to estimate the local continuum, flagging and replacing pixels exceeding $3\,\sigma$ from this baseline with the median value. To maintain high signal-to-noise, we retained only the top 80 per cent of flux values in each order. The cleaned orders were then continuum-normalised by dividing by a Savitzky–Golay-filtered continuum (window length 601 pixels, 3rd-order polynomial). Each pixel was assigned an inverse-variance weight based on the local flux standard deviation (computed within a 51-pixel window). All orders were subsequently interpolated onto a common wavelength grid spanning the full observed wavelength range, with overlapping regions combined via inverse-variance weighting. Remaining outliers were identified and removed using Astropy’s sigma_clip (Astropy Collaboration et al. 2022). Finally, the stitched 1D spectrum was smoothed by convolution with a Gaussian kernel matched to the spectrograph’s resolution ( $R = 75\,000$ ).

2.3 Data post-processing

To correct for telluric absorption, we used the radiative transfer code Molecfit (Smette et al. Reference Smette2015; Kausch et al. Reference Kausch2015) to model and remove atmospheric features from each exposure individually. We aligned the one-dimensional spectrum using the telluric lines before performing a chi-squared fit to match the tellurics present in the data. A synthetic telluric model was then generated and divided out from each individual exposure, removing the imprint of Earth’s atmosphere. After telluric correction, we compute an out-of-transit master spectrum and divide each exposure by this baseline to isolate the planetary signal, producing a time series of residual spectra. These residuals are primarily shaped by variations in continuum levels between exposures, resulting in non-flat baselines. We flatten each spectrum using a high-pass filter with a width of 80 km s $^{-1}$ . A running median absolute deviation (MAD) filter with a width of 20 pixels is then applied: pixels deviating by more than 5 $\sigma$ from the local median are flagged as outliers and replaced via linear interpolation. Finally, the spectra are visually inspected for remaining contaminants, such as interstellar absorption features, which are manually masked with NaNs. This observational setup follows the procedure outlined in Borsato et al. (Reference Borsato2023).

3.1 CCF cleaning and detrending

Due to the sensitivity of cross-correlation, the process can amplify small residual artefacts in the data. Two dominant signals often appear: line distortions in the stellar spectrum caused by the Rossiter–McLaughlin (RM) effect, and residuals from imperfect telluric correction by Molecfit. To address these, we apply additional standard cleaning procedures (Hoeijmakers et al. Reference Hoeijmakers2019; Prinoth et al. Reference Prinoth2022; Borsato et al. Reference Borsato2023; Prinoth et al. Reference Prinoth2023). The RM effect appears in the CCF as a bright, time-dependent Gaussian feature. We remove this by forward-modelling a time-varying Gaussian using an MCMC sampler with informed priors on its amplitude and velocity offset. The best-fit Gaussian is then subtracted from the raw CCF map, effectively eliminating the RM signal. In principle, this feature can be modelled directly in the residual CCFs (as in tomographic studies such as Rainer et al. Reference Rainer2021) to extract the same parameters usually derived from RM analyses, such as the projected obliquity and $v \sin i_\star$ . In the present work, we remove it in order to recover a coherent planetary signal in the velocity–velocity diagrams. In future work, we plan to carry out a dedicated analysis of this dataset to investigate the RM effect in detail, alongside other observables that we speculate may be accessible (see Section 5).

To mitigate telluric residuals, we fit a fifth-degree polynomial along the orbital phase axis at each velocity step. This removes vertical patterns in the CCF. Prior to fitting, we mask the expected planetary trail using NaNs to ensure the fit is not biased. Subtracting this two-dimensional trend surface yields a fully cleaned CCF map. The full cleaning procedure is illustrated in Figure 1.

Figure 1. Step-by-step illustration of the cross-correlation function (CCF) cleaning pipeline applied to one night of WASP-189b observations using the combined template. (a) Raw CCF map showing the bright, time-varying Rossiter–McLaughlin (RM) effect. (b) Best-fit Gaussian model of the Rossiter–McLaughlin (RM) effect. (c) CCFs after subtracting the RM model. (d) Binary mask used to exclude the planetary signal from detrending. (e) Fifth-degree polynomial detrending surface fitted along the phase axis, masking the planetary trail. (f) Final cleaned CCF map. Colour represents signal strength in units of the continuum standard deviation, with blue indicating absorption and green indicating emission.

3.1.1 Constructing velocity–velocity diagrams

Velocity–velocity diagrams (often referred to as $K_p$ – $V_{\text{sys}}$ maps) enhance the detection fidelity of the cross-correlation function (Brogi et al. Reference Brogi2012). They represent co-added cross-correlation exposures, combining all observations into a single detection signal. Assuming a static and symmetrical atmosphere surrounding the planet, its radial velocity can be computed as:

(3) \begin{equation}v_r(\phi) = V_{\mathrm{rv}} + V_{\mathrm{grad}} \,\sin i \,\sin\bigl(2\pi\phi\bigr)\;;\quad\phi = \frac{t - T_0}{P}\end{equation}

where $v_r(\phi)$ is the planet’s radial velocity at orbital phase $\phi$ , t is the time stamp of observation, $T_0$ is the mid-transit time, and P is the orbital period. $V_{\mathrm{rv}}$ is a radial velocity offset, $V_{\mathrm{grad}}$ is a time-dependent velocity gradient, and i is the orbital inclination. By applying Equation (3) across a grid of trial $V_{\mathrm{grad}}$ values, we shift each CCF exposure into hypothetical planetary rest frames and average them. The resulting 2D map provides a clear visual confirmation of any atmospheric signal (see Figure 2). Note that $V_{\mathrm{rv}}$ and $V_{\mathrm{grad}}$ often correspond to the systemic velocity ( $V_{\mathrm{sys}}$ ) and orbital semi-amplitude ( $K_p$ ), respectively, but can differ due to atmospheric geometry and dynamical effects.

3.2 Model injection

To evaluate whether the absence of a signal reflects a true non-detection or stems from limitations in our analysis, we perform model injection. This technique has become standard in high-resolution exoplanet spectroscopy (Snellen et al. Reference Snellen, de Kok, de Mooij and Albrecht2010; Birkby et al. Reference Birkby2013; Hoeijmakers, Snellen, & van Terwisga Reference Hoeijmakers, Snellen and van Terwisga2018a).

For the injection–recovery tests, each of the 17 MANTIS templates at 4 000 K was injected individually into the observed spectra prior to CCF computation. This allows us to evaluate, species by species, whether a signal matching that template would be detectable given our S/N and observing conditions. Injected models were centred on the literature values of the systemic velocity ( $V_\textrm{sys} = -26.8$ km s $^{-1}$ ) and orbital semi-amplitude ( $K_p = 202$ km s $^{-1}$ ; Anderson et al. Reference Anderson2018).

Because the MANTIS templates assume a fiducial hot-Jupiter system ( $R_p = 1.5\,R_J$ , $R_\star = 1\,{\rm R}_\odot$ , Kitzmann et al. Reference Kitzmann2023), their absolute line depths must be rescaled to represent WASP-189b. We therefore apply a correction factor based on the squared radius ratio,

(4) \begin{equation} f = \left(\frac{(R_p/R_\star)_{\mathrm{WASP-189A}}}{(R_p/R_{\star})_{\mathrm{MANTIS}}}\right)^2 .\end{equation}

All templates are scaled by this factor before injection to ensure their line depths are consistent with WASP-189b’s true geometry. For WASP-189b ( $R_p = 1.619\,R_J$ , $R_\star = 2.362\,{\rm R}_\odot$ , Lendl et al. Reference Lendl2020), $f \approx 0.219.$

Injections allow us to verify whether our data quality and channel throughput are sufficient to recover the expected signal. In contrast, for newly discovered planets or species with no prior detections, theoretical templates generated via radiative transfer models can be injected to determine whether a non-detection reflects true atmospheric absence or insufficient sensitivity.

These models are broadened to account for instrumental resolution and planetary rotation, and resampled to match the Veloce wavelength grid. We then reprocess the injected spectra through our full pipeline, including telluric correction, continuum normalisation, and cross-correlation. By comparing the injected results to our original data, we assess whether a species should have been detectable given our signal-to-noise and observing conditions. This allows us to distinguish between non-detections driven by physical absence and those due to insufficient sensitivity or loss of signal from systematics.

3.3 Combining the channels of veloce

For each spectrograph arm (Azzurro, Verde and Rosso), we begin by loading the cross-correlation function (CCF) map. We measure the noise level by computing the standard deviation across the edges of the velocity grid, where no planetary signal is expected. We identify a small region in radial velocity–velocity gradient space around the expected planetary signal. Within this region, we measure the average strength of the normalised signal and, where available, compare it to the injected-model prediction. These values are used to assign a raw weight to each arm. Negative weights are discarded, and the remaining weights are normalised so that their sum equals one. This ensures that the contribution from each spectrograph arm is proportional to its signal content.

Figure 2. Detection of ionised calcium ( $\mathrm{Ca}^{+}$ ) in the atmosphere of WASP-189b. Top row: 2D cross-correlation function (CCF) maps for the data (left) and the expected signal from model injection (right). Bottom row: Corresponding velocity-velocity diagrams (data left; expected signal right). Dashed lines indicate the systemic velocity ( $V_\textrm{sys} = -26.8\,$ km s $^{-1}$ ) and semi-amplitude ( $K_p = 202\,$ km s $^{-1}$ ).

4.1 The detection of $\mathrm{Ca}^{+}$

We detect a clear atmospheric signal of ionised calcium ( $\mathrm{Ca}^{+}$ ) in WASP-189b, partially visible in the cross-correlation function (CCF) and robustly in the velocity–velocity diagram shown in Figure 2. The planetary trace appears slightly redshifted in $V_\textrm{rv}$ compared to the expected systemic velocity, while the peak $V_\textrm{grad}$ is offset to somewhat lower values than the expected orbital semi-amplitude. The signal in the velocity-velocity diagram exceeds $5\sigma$ in amplitude, based on the standard deviation of background regions away from the planetary trail. The injected model predicts a comparative detection significance as well, suggesting that with the current throughput of Veloce we would expect the signal to be this strong. Among known species in this planet’s atmosphere, $\mathrm{Ca}^{+}$ is the strongest contributor to optical absorption (Prinoth et al. Reference Prinoth2023), and our result represents a robust agreement to this fact using a single transit with Veloce.

4.2 Additional species detections

Beyond $\mathrm{Ca}^{+}$ , we also produce evidence for absorption at the expected velocities for H, Na, Mg, Ti, $\mathrm{Ti}^{+}$ and $\mathrm{Fe}^{+}$ (Figure 3). These panels exhibit coherent absorption peaks in the cross-correlation maps, indicative of possible planetary absorption. Other tested species–Ca, Fe and $\mathrm{Sr}^{+}$ –show broad depressions at the correct velocities but lack distinct signatures. While these signals warrant further investigation, robust confirmation will require stacking multiple transits to enhance signal-to-noise. Many of these signals have peaks at lower values than the expected $K_p$ of WASP-189b. This is a common outcome in cross-correlation studies and is attributed to the asymmetric geometry of the planet’s inflated, highly irradiated dayside (e.g. Prinoth et al. Reference Prinoth2022).

Figure 3. Mosaic of velocity-velocity diagrams for the species robustly detected in a single transit of WASP-189b with Veloce. Each panel shows the combined CCF signal for H, Na, Mg, $\mathrm{Ca}^{+}$ , Ca, Ti, $\mathrm{Ti}^{+}$ , $\mathrm{Fe}^{+}$ , Fe and $\mathrm{Sr}^{+}$ , arranged by increasing atomic number from top left to bottom right. White dashed lines mark the expected systemic velocity ( $V_\textrm{ sys}=-26.8\,$ km s $^{-1}$ ) and orbital semi-amplitude ( $K_p=202\,$ km s $^{-1}$ ). The colour scale indicates cross-correlation strength in units of the continuum standard deviation.

We compare the relative performance of Veloce’s three arms using the $\mathrm{Fe}^{+}$ signal as a reference (see Figure B1). The green (Verde) arm contributes the strongest $\mathrm{Fe}^{+}$ absorption, and the blue (Azzurro) arm also shows a significant contribution, reflecting deeper $\mathrm{Fe}^{+}$ line strengths in covered wavelength region. In contrast, the red (Rosso) arm contributes the least, as the $\mathrm{Fe}^{+}$ lines it covers are relatively shallow, despite the higher signal fidelity. This comparison highlights the importance of tailoring future observational strategies to optimise sensitivity to ionised metals.

For completeness and transparency, we present the mosaic of non-detections in Figure B2. These panels demonstrate the absence of clear signals for V, Cr, Mn, Ni and $\mathrm{Ba}^{+}$ , underscoring the need for higher S/N or multiple-transit datasets to explore these species further. While TiO appears to show a feature, we do not interpret this as a detection. Our injected models do not predict a measurable TiO signal, and prior studies (e.g. Prinoth et al. Reference Prinoth2022) have shown that robust TiO detections typically require many stacked transits.