Issue 
A&A
Volume 673, May 2023



Article Number  A158  
Number of page(s)  31  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202245121  
Published online  24 May 2023 
The Mantis Network
III. Expanding the limits of chemical searches within ultrahot Jupiters: New detections of Ca I, V I, Ti I, Cr I, Ni I, Sr II, Ba II, and Tb II in KELT9 b
^{1}
Lund Observatory, Department of Astronomy and Theoretical Physics, Lund University,
Box 43,
221 00
Lund,
Sweden
email: n.borstato@astro.lu.se
^{2}
Department of Astronomy, School of Science, The University of Tokyo,
731 Hongo,
Bunkyoku, Tokyo
1130033,
Japan
^{3}
University of Bern, Center for Space and Habitability,
Gesellschaftsstrasse 6,
3012,
Bern,
Switzerland
^{4}
UniversitätsSternwarte, Fakultät für Physik, LudwigMaximiliansUniversität München,
Scheinerstr. 1,
81679
München,
Germany
Received:
3
October
2022
Accepted:
1
April
2023
Crosscorrelation spectroscopy is an invaluable tool in the study of exoplanets. However, aliasing between spectral lines makes it vulnerable to systematic biases. This work strives to constrain the aliases of the crosscorrelation function to provide increased confidence in the detections of elements in the atmospheres of ultrahot Jupiters (UHJs) observed with highresolution spectrographs. We use a combination of archival transit observations of the UHJ KELT9 b obtained with the HARPSN and CARMENES spectrographs and show that it is possible to leverage each instrument’s strengths to produce robust detections at a substantially reduced signaltonoise. Aliases that become present at low signaltonoise regimes are constrained through a linear regression model. We confirm previous detections of H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, and Fe II, and detect eight new species, Ca I, Cr I, Ni I, Sr II, and Tb II, at the 5σ level, and Ti I, V I, and Ba II above the 3σ level. Ionised terbium (Tb II) has never before been seen in an exoplanet atmosphere. We further conclude that a 5σ threshold may not provide a reliable measure of confidence when used to claim detections, unless the systematics in the crosscorrelation function caused by aliases are taken into account.
Key words: planets and satellites: atmospheres / planets and satellites: gaseous planets / techniques: spectroscopic
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Exoplanets and their inherent diversity have dramatically transformed the field of astronomy. We now understand that most stars host exoplanets. Beyond just detecting these celestial bodies, we can now observe them in intricate detail, revealing their atmospheric and orbital characteristics, formation histories, and prevalence. Examining exoplanet atmospheres is of significant interest to both observers and theorists, as it allows us to further characterise these entities. Nevertheless, discerning the atmospheric signals of an exoplanet from those of its host star remains a challenge. Since the initial successful detection of the ultrahot Jupiter HD 209458 b’s atmosphere (Charbonneau et al. 2002), instrumental advancements have made it possible to identify the absorption signatures of individual ionised, atomic, and molecular species.
KELT9 b, which was discovered by Gaudi et al. (2017), is an ultrahot Jupiter (UHJ) with the highest recorded equilibrium temperature of 4000 K. The planet has been observed with many instruments, including the HARPSNorth (Hoeijmakers et al. 2019; Pino et al. 2022), CARMENES (Yan & Henning 2018), PEPSI (Cauley et al. 2019; Pai Asnodkar et al. 2022b), FIES (BelloArufe et al. 2022), and TRES (Lowson et al. 2023) spectrographs. These observations have led to the detection of the following species: H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II.
Detecting heavier elements, such as those produced by neutroncapture processes (Burbidge et al. 1957) in exoplanets, is desirable due to the information they can reveal about planet formation and the host stars they orbit (Ek et al. 2020). Studies of meteorites in our Solar System have shown that the inner Solar System (e.g. the Earth) is relatively enriched in isotopes produced by the slower neutroncapture process (sprocess) compared to the much more rapid neutroncapture process (rprocess), which contradicts current established theories on planet formation (Burkhardt et al. 2011; Budde et al. 2016, 2019; Saji et al. 2020). The inhomogeneity of these heavy isotopes in Solar System meteorites may be driven from where the protoplanetary material, which makes up the meteorites, originates (Dauphas et al. 2010; Lugaro et al. 2016; Ek et al. 2020). In conjunction with this observation, massive planets are commonly known to migrate inwards towards their host stars, and form on faster timescales than terrestrial planets (Pollack et al. 1996; Chambers & Wetherill 1998). Therefore, giant planets should sample the protoplanetary material used in the formation of terrestrial planets as they migrate, leading to the postulation that there could be a coupling of the inhomogeneity trend with the Jupitermass planet’s atmosphere.
Furthermore, Meléndez et al. (2009) noted that the Sun, as compared to solar twins, is depleted in several refractory elements – elements with a high condensation temperature – such as yttrium and barium. The possibility that planet formation imprints chemical signatures on the host star is an active area of research, where differences in observed stellar abundance trends with condensation temperature are likely associated with the histories of planet formation (e.g. Booth & Owen 2020; Liu et al. 2020, 2021). Additionally, Kruijer et al. (2017) show that measurements of heavy elements have proven to be a useful tool in determining the age of Jupiter, which is directly linked to the age of the Sun. Age determination is an ongoing challenge in stellar astrophysics (Pont & Eyer 2004; Sahlholdt et al. 2019), and an independent approach through planet age determination could help to validate current age estimate approaches. Thus, the study of heavy transiron elements in exoplanets has the potential to help explain how planets form, and to infer more information about the host stars they inhabit. Detection studies and further analysis of heavier elements in UHJs would present the first steps in constructing a sample size to further test the ideas of planet formation in a context outside of the Solar System.
In order to study the chemical composition of exoplanets, it is necessary to remove auxiliary signals, such as telluric and stellar signals. The removal of auxiliaries results in a noisedominated spectrum with the planet’s absorption spectrum buried inside but all for the deepest transmission lines. These deep lines, however, are routinely analysed (e.g. Ehrenreich et al. 2008; Wyttenbach et al. 2017; Salz et al. 2018; Cauley et al. 2019; AlonsoFloriano et al. 2019b; Žák et al. 2019; Seidel et al. 2020, 2022; Turner et al. 2020; Lampón et al. 2021; Pai Asnodkar et al. 2022b), and they can directly measure the massloss rate, exosphere, temperature profile, and atmospheric dynamics, and constrain where in the atmosphere a particular species is located. The crosscorrelation technique (Snellen et al. 2010) can be used to leverage all the known lines of a given species in a spectrum to produce an average line profile. This method increases the signaltonoise ratio (S/N) enough so that the species becomes detectable, although it removes information about the individual lines and prevents studies of individual line profiles. Nevertheless, crosscorrelation enables the acquisition of chemical inventories and the analysis of their atmospheric characteristics, such as temperaturepressure profiles, chemical abundances, and surface gravity (e.g. Brogi & Line 2019; Fisher et al. 2020; Pelletier 2021; Gibson et al. 2022).
The primary use of the crosscorrelation technique is to detect chemical species in an exoplanet’s atmosphere. With it, a diverse range of ionised, atomic, and molecular species have been identified, including metals such as iron and titanium (e.g. Hoeijmakers et al. 2018a, 2019, 2020a,b; CasasayasBarris et al. 2019; Ehrenreich et al. 2020; Pino et al. 2020; Cabot et al. 2021; Merritt et al. 2021; BelloArufe et al. 2022), and molecules such as water, carbon monoxide, hydroxide radicals (e.g. Brogi et al. 2012, 2013, 2014; de Kok et al. 2013; Birkby et al. 2017; AlonsoFloriano et al. 2019a; Landman et al. 2021; Giacobbe et al. 2021; Nugroho et al. 2021; van Sluijs et al. 2023; Guilluy et al. 2022), and titaniumoxide (Nugroho et al. 2017; Cont et al. 2021; Prinoth et al. 2022).
A substantial proportion of chemical species detections have come from the same class of planets, UHJs, which are shortperiod giant planets with equilibrium temperatures comparable to the effective temperatures of M dwarfs. Their extreme equilibrium temperature creates an extended, cloudfree atmosphere that maximises the potential flux contribution that can come from the optically thin region of an exoplanet’s atmosphere (Heng 2016; Stevenson et al. 2016; Kitzmann et al. 2018). Thus, UHJs have effectively become laboratories for experimenting with exoplanet spectra, enabling the field to develop and refine techniques. Provided that an adequate line list exists, it is possible to search for any element, although not all are equally detectable (Kesseli et al. 2022). For heavier elements, abundances generally decrease, entailing weaker absorption from the given species. Furthermore, absorption is dictated by the combined line opacity of a species (Heng & Kitzmann 2017), so species with fewer lines may also be difficult to detect. If the scientific goal is to detect a particular species, as much spectroscopic information as possible should be leveraged to obtain the highest chance of detection.
Crosscorrelation results rely heavily on the chosen line list, with different lists containing slight variations in the calculated crosssections of a given species (GharibNezhad et al. 2021). Incorrect crosssections can lead to inefficient or even null detections if the mismatch is large enough (Hoeijmakers et al. 2015). Linelist inaccuracies are likely the cause of why the molecular species TiO has been challenging to detect until recently (e.g. Hoeijmakers et al. 2015; Nugroho et al. 2017; Sedaghati et al. 2021; Prinoth et al. 2022). Therefore, accurate line lists that are representative of the atmosphere are required (Kitzmann et al. 2023). Combining lines that would not be present in a spectrum, or assigning too much weight to a specific line reduces the efficiency of the crosscorrelation process. In all cases, the result is a noisier signal with a high risk of nondetection for faint signals.
Another inhibiting factor that has a tangible effect on the crosscorrelation signal is the wavelength range of the spectrograph. Typically, studies only analyse data from one instrument (e.g. Yan & Henning 2018; Hoeijmakers et al. 2019; Ehrenreich et al. 2020; Merritt et al. 2021; Prinoth et al. 2022; Stangret et al. 2022; BelloArufe et al. 2022), as it is relatively easy to combine observations from the same instrument. Techniques that combine spectrographic observations from different highresolution instruments are still in their infancy, but they have been used in the past (e.g. Turner et al. 2020; Pai Asnodkar et al. 2022b; Pino et al. 2022). As archives collect multiple transit observations of the same target, developing robust combination techniques becomes important. Combining the information from multiple spectrographs enhances the quality of observations beyond the capabilities of a single instrument, as it allows for a more extensive coverage of wavelength ranges that may be more conducive to detecting certain species over others.
An important topic discussed less in crosscorrelation studies is aliases. Aliases are systematic artefacts of crosscorrelation functions, either from autocorrelation with the search species (Hoeijmakers et al. 2018b) or correlations with alternative species in the spectrum (Hoeijmakers et al. 2018a, 2019). The systematic variation distorts the noise in the baseline of the crosscorrelation function and can, in turn, severely undermine the uncertainty analysis. This problem becomes more apparent when searching for less abundant species or species with fewer lines, as interference from species with stronger lines becomes more likely.
The growing number of observations of KELT9 b leads to a relative decrease in what is achievable in a single night of observations, while increasing the motivation to combine multiple spectroscopic observations. This study combines archival data from the HARPSN and CARMENES spectrographs to detect exotic species in the atmosphere of KELT9 b. The defined approach combines crosscorrelation functions from multiple observations, and constrains the predicted alias contributions to yield more statistically robust detections. First searches for the same elements found in the publications of Yan & Henning (2018), Yan et al. (2019), and Hoeijmakers et al. (2019) serve as verification before searching for new undetected species in the combined data sets. The aim is to show that signals increase in their S/Ns when multiple spectrographs are combined, and also begin to resolve alias structures, while constraining these structures can lead to new detections of atoms in KELT9 b’s atmosphere.
This paper is structured as follows: in Sect. 2 we summarise the data and processing steps. In Sect. 3, we introduce the crosscorrelation technique and our new method for extracting the planetary signals from the data; additional supplementary information is provided in Appendix A. We present the results of our attempts to detect new species in KELT9 b’s atmosphere in Sect. 4, with additional results discussed in more detail in Appendix B, and validating these detections in Appendix C. Finally, we discuss the implications of the results in Sect. 5, and draw conclusions in Sect. 6.
Summary information of the transit observations used in this study.
2 Data collection and processing
The first two of the four analysed data sets originate from observations taken with the HARPSNorth spectrograph (HARPSN, Mayor et al. 2003) located on the Telescopio Nazionale Galileo. The original studies that used these data were Hoeijmakers et al. (2018a, 2019). HARPSN is an optical spectrograph spanning a wavelength range between 387.4 and 690.0nm, at a spectral resolution of 115 000. The two transit observations were conducted on 31 July 2017 and 20 July 2018, respectively, with 600 s exposures. Both observations cover the 3.9 h transit and additional baseline measurements totalling 49 and 46 exposures, with median S/N values over the echelle orders of 74 and 71, respectively. The primary science products obtained from these observations are the individually extracted echelle orders from the 2D echellogram and the blazecorrected, stitched and resampled 1D spectra.
The two additional data sets originate from the CARMENES spectrograph (Quirrenbach et al. 2014) from the CalarAlto Observatory. CARMENES is an optical and nearinfrared spectrograph with the optical arm covering a wavelength range between 520 and 960 nm and an infrared wavelength between 960 and 1710 nm at a resolution of 94 600 and 80 400. The first publications to analyse this data were Yan & Henning (2018) and Yan et al. (2019). The observations were taken on 6 August 2017 and 16 June 2018, covering the entire transit. Both observations cover the entire transit with 300 s exposures for the first night and 111 s for the second. The present study uses the optical arm of the instrument with median S/N values of 35 and 34. The primary science products drawn from the archive are the individually extracted and blazecorrected echelle orders, estimated errors and blaze profiles. The 2D orders were stitched together using a weighted average to construct a 1D spectrum of each exposure. Crosscorrelation is more efficient when the blaze profile has not been removed, as this eliminates the need to reweigh lines in the template by the effective throughput of the spectrograph (i.e. near order edges). Therefore, the profiles are divided back into spectral orders. This information is summarised in Table 1.
In order to optimise the crosscorrelation search, we removed any systematic signals that do not belong to the planet itself. The adopted process follows the same methodology applied in Hoeijmakers et al. (2020b). This reference provides an indepth walkthrough of the data preparation. The systematics primarily consists of atmospheric telluric contamination, velocity offsets, the stellar spectrum, and any additional effects such as interstellar medium (ISM) lines or residual continuum offsets.
The telluric lines produced by the H_{2}O and O_{2} absorption bands in the atmosphere are removed using molecfit (Smette et al. 2015). molecfit creates a model for the telluric transmission spectrum over the entire wavelength of the spectrograph range using the 1D spectra of each observation. The models are then interpolated onto the wavelength solutions of the 2D echelle orders and divided out. This process applies to spectra without blaze correction since the operation is multiplicative. The correction is enough for HARPSN and the CARMENES visible arm to bring the telluric corrections to the noise level of the spectra.
In addition to the telluric removal, the stellar spectra require two Doppler corrections to correctly place each spectral order into the host star’s rest frame. The first is the centre of mass motion of the EarthSun system, the barycentric velocity, while the second corrects the star’s radial velocity induced by the planet’s gravitational pull. The correction leaves the spectra at a constant radial velocity shift set by the systemic velocity of approximately −18 km s^{−1} (Gaudi et al. 2017; Borsa et al. 2019; Hoeijmakers et al. 2019; Pai Asnodkar et al. 2022b), all other system parameters were taken from Gaudi et al. (2017). We removed the stellar signal from the time series spectra by calculating an outoftransit baseline of the star and dividing each spectrum by it. The spectral orders were then flattened using a median filter in order to remove broadband continuum variations, and with a running standard deviation applied in suit. Values greater than 5σ were flagged as NaNs. The process was then repeated instead using a Gaussian highpass filter of 80 km s^{−1}, with any remaining 5σ outliers flagged as NaNs. The resultant NaNs created an outlier mask which avoided problematic pixels in the crosscorrelation process (Hoeijmakers et al. 2018a).
3 Crosscorrelation analysis for transit spectroscopy
Crosscorrelation combines the expected contributions of the absorption lines for a particular species in spectrum. The procedure requires using a spectral template T of a specific species over a defined wavelength axis i (T_{i}), which predicts the transmission spectrum of a specific species. The template is converted to a set of weights by normalisation: (1)
where T_{i}(υ) represents the template species Doppler shifted to a specific velocity, and υ represents the Doppler shift of the template. Crosscorrelation is implemented using the following equation: (2)
where c(υ) represents the weighted average of line positions at a given Doppler shift, N represents the total number of wavelength bins, and x_{i} is the spectral data which shares the same wavelength range i as the template. The template and crosscorrelation functions are velocitydependent. During transit, a planet’s signal Doppler shifts to its instantaneous radial velocity and a crosscorrelation will form a maximum at this point, assuming the atmosphere’s velocity profile is comparable to the planet’s. With a radial velocity step of 1 km s^{−1}, we sample a radial velocity range of −1000 km s^{−1} to 1000 km s^{−1}. This option covers more aliasing lines in the crosscorrelation function, which allows for tighter constraints on the aliasregression models without compromising computational efficiency. We introduce this in Sect. 3.3.
During transit, the stellar disk becomes partially obscured and this creates a Doppler shift due to the rotation of the star (Rossiter 1924; McLaughlin 1924). In a crosscorrelation map, this results in a bright apparent emission feature called the Doppler shadow that requires removal. We remove the Doppler shadow by masking where the planet’s signal is expected to lie on the crosscorrelation plane, assuming a Keplerian circular orbit: (3)
where υ_{rv} is the planetary radial velocity during transit, υ_{orb} is the orbital velocity of the system, ϕ is the orbital phase of the planet, i is the system inclination, and υ_{(sys)} is the star’s systemic velocity. We use this expression to mask the planet’s signal along the trace by masking pixels in the crosscorrelation function that lie within 20 km s^{−1} of the expected planet signal. Removal of the shadow follows the same procedure as Hoeijmakers et al. (2020b), which fits an empirical model consisting of two Gaussians as a profile that has a positive central core, negative sidelobes and a free scaling factor.
Furthermore, correlated structures caused by the incomplete removal of strong stellar and telluric lines can remain in the crosscorrelation function. These stellar aliases create nearvertical structures. Removing these artefacts follows a modified approach to Prinoth et al. (2022), which fits a 1D polynomial from each velocity value and subtracts the fit. Instead, rather than fitting for a single velocity value, a fit is applied to a sliding average on a column of ten velocity values, and subtraction is applied to the central column. The column window is then moved across the crosscorrelation map, subtracting over the whole radial velocity range. The modified approach utilises the partial correlation between neighbouring pixels to better sample vertical trends. Figure 1 illustrates the applied cleaning steps.
The exoplanet’s signal can be strengthened further by coadding all pixels in the planet’s trace. To do this, the individual exposures are shifted to the timeindependent rest frame of the planet. Since the planet’s orbital phase and inclination are well known, it is possible to apply Eq. (3) by assuming a range of orbital velocities, in this case spanning 0 to 400 km s^{−1} and shifting the pixels in the crosscorrelation map accordingly. Coadding intransit crosscorrelation functions creates a 1D crosscorrelation function for every potential orbital velocity (Nugroho et al. 2017; Prinoth et al. 2022). Stacking these 1D crosscorrelation functions produces a K_{p} – V_{sys} map (Brogi et al. 2012). If an absorption spectrum of the search species exists then the resulting map should show a peak at the orbital velocity of the planet and the stellar systemic velocity.
Fig. 1 Overview of the crosscorrelation cleaning process. Panel A: the raw crosscorrelation function timeseries applied after the outoftransit spectra have been divided and the residual broadband features removed. Both the Doppler shadow and the alias structures are present. Panel B: the planet mask covers the expected planet signal, the halfwidth of the planet mask is 20 km s^{−1}. Panel C: the fitted Dopplershadow profile is a twocomponent Gaussian. Panel D: the crosscorrelation function after removal of the Doppler shadow, residual stellar aliases are still present. Panel E: the resulting crosscorrelation map after removing the stellar aliases. The signature of the planet appears as the dark slanted trace. In this plot, the sign is negative; however, it is flipped in preceding figures to denote absorption. 
Fig. 2 Comparison plot demonstrating the improved signal from applying the combination method for the species Fe I. Left panel: K_{p} – V_{sys} map of a single transit observation, this specifically is the first HARPSN transit. The red border marks the window where the average signal strength is computed, and the region outside where the standard deviation is computed. Right panel: K_{p} – V_{sys} map when combining all four transit observations. The result is a much higher S/N. 
3.1 Crosscorrelation templates
Crosscorrelation requires a template mask to define the weights with which pixels are coadded together, as specified in Eq. (2). We constructed masks through forward modelling KELT9 b’s transmission spectrum, assuming that the planet contained the absorption lines of the search species. The calculation of these templates primarily follows that of Hoeijmakers et al. (2019). We used HELIOSK^{1} (Grimm & Heng 2015) to compute the opacity functions for each species, drawing from the transmission tables of the VALD (Ryabchikova et al. 2015) and Kurucz (Kurucz 2018) databases. The continuum opacity function for the mask comprised HHe, H_{2}–H_{2} and H_{2}–He collision induced absorption obtained from HITRAN (Richard et al. 2012), while the H^{−} boundfree and freefree absorption contributions were taken from John (1988). Crosscorrelation templates are computed following the method described in Kitzmann et al. (2023). The chemical composition is determined with the equilibrium chemistry code FASTCHEM^{2} (Stock et al. 2018, 2022) assuming solar element abundances and an isothermal temperature profile with 4000 K in the pressure range 10–10^{−15} bar. These high temperatures create partially ionised atoms and a low abundance of molecules. The templates for the detected species can be found in Fig. A.2.
3.2 Stacking observations
Combining spectroscopic observations from multiple nights with different spectrographs is nontrivial. When using a single spectrograph, it is possible to shift the observed spectra to the planetary rest frame for each night, and compute an average spectrum before applying the crosscorrelation function (Hoeijmakers et al. 2020b). Alternatively, observations can be averaged in K_{p} – V_{sys} space as the process projects the crosscorrelation functions on the same grid for each night. In either case, an average dampens the noise leading to clearer signals. However, a standard mean does not consider variations in data quality caused by weather and airmass variations. Additionally, when using different instruments, each spectrograph’s unique wavelength coverage also creates differences in results. Therefore, additional steps are required to take these aspects into consideration.
To deal with the variations in data quality from the same instrument, we were able to mask out the planet’s signal in the K_{p} – V_{sys} map. This was done by constructing a window of values surrounding the expected area of the signal. In our case, the window ranged from −50 km s^{−1} to 10 km s^{−1} for the systemic velocity, and 100 km s^{−1} to 300 km s^{−1} for the orbital velocity (see Fig. 2). The window selected deliberately encompassed a wide range of velocity points to ensure any signal which may be offset, due to different atmospheric dynamics (Prinoth et al. 2022), was captured (see Fig. 6). We then maskedout the signal, and took the standard deviations of the remaining pixels. Taking the standard deviation of the masked K_{p} – V_{sys} map measured the noise and any systematics that remained in the data without contribution from the signal of the planet. It was then possible to invert the mask and apply it to the K_{p} – V_{sys} map such that only the expected region where one would expect the signal to lie was extracted. If a signal existed in this region taking the mean would lead to a value which was larger on average than the background, if no signal existed the mean would average the noise and any systematics that lay in that region. A pseudoS/N value could be constructed by taking the ratio of the mean and the standard deviation, whose value was dependent on the night’s quality.
The pseudoS/N calculation to K_{p} – V_{sys} maps is not sufficient to compare the quality of observations between different spectrographs. Different instruments cover different wavelength ranges and thus cover different lines. Since crosscorrelation uses the line position of all lines across the wavelength range, this will create a difference in the result even if the data quality is identical. It is possible to account for this difference through the templates applied in crosscorrelation. Applying a simple sum of the total absorbed flux of the template over the wavelength range achieves an approximate measurement for this value. This number will be the same for the same spectrograph but will differ for other instruments. Instruments that contain deeper more readily observed lines obtain a higher value.
With the sum of the line depths of the instrument and the pseudoS/N of the masked K_{p} – V_{sys} map, it is possible to compute a statistic that enables comparison between nights and normalise the values to create a set of weights, (4)
where i refers to the night of observation, Σ_{T} is the sum of the template flux, and SNR_{pseudo,i} is the pseudoS/N of the K_{p} – V_{sys} map. Figure A.1 outlines the weight extraction process for Fe I and Na I. For Fe I, higher preference is given to the HARPSN nights. This is likely due to its higher quality produced by the telescope, its high efficiency in resolving blue wavelengths and greater abundance of deeper lines. For Na I, however, the weights the CARMENES share a larger proportion of the total sum, as the combined flux from both spectrographs for this species is comparable. However, the data quality is still of better quality for the HARPSN nights leading to stronger weight values.
The last step required before the K_{p} – V_{sys} maps can be combined is to normalise each map by subtracting the mean of the masked map and dividing by the standard deviation for each K_{p} – V_{sys} map. The transformation creates S/N maps commonly used in the literature (Giacobbe et al. 2021; Guilluy et al. 2022). The act of scaling removes the implicit dependence of the line depth, which vary with different instruments. It is then possible to apply a weighted average using the precalculated weights to all nights to create one master S/N map. Finally, to reconvert the master map to one with a line strength: the master S/N map is normalised to a maximum value of one. This array of values is multiplied with the K_{p} – V_{sys} map which produces the most significant weight. At this point, all maps have been combined to produce a signal that mimics the profile of the best observation but at a larger S/N. Figure 2 illustrates the improved signal of Fe I, moving from the best HARPSN night to the combined observation.
3.3 Aliasing in crosscorrelation
A potential vulnerability with the crosscorrelation technique is the presence of aliases, which occur when lines in a template correlate with the lines of other species. In a 1D crosscorrelation function, this appears as an apparent signal located away from the systemic velocity, which is easy to spot by eye (Hoeijmakers et al. 2019, 2020b). However, peaks that lie close to the systemic velocity are more difficult to untangle, given that the dynamic effects of an UHJ’s atmosphere can also shift signals of the order of a few km s^{−1} (Ehrenreich et al. 2020; Prinoth et al. 2022; Pai Asnodkar et al. 2022a).
A simple toy model of a spectrum illustrates how aliases may cause false detections. In the following example, presume that a transmission spectrum consists of three different species over 400–420 nm. One species has many absorption lines (species A), for example, Fe I, which has approximately 500 lines in this region. The other two species (species B & C) have 100 and 20 lines, respectively. A summation of a uniform distribution of Gaussian peaks, with random height, width, and centre positions creates the transmission spectrum: (5)
where i is the total number of species, l is the total number of lines that each species has, and A, λ and σ represent the wavelength amplitudes, centroid position and standard deviation, respectively.
The amplitude of each line varies from 0.001 to 0.005 of a transit depth, with a standard deviation of 0.003–0.005, with the exception of species C which has twenty lines all with the same line depth of 0.005. The produced data sets represent a simplified transmission spectrum with the planet centred at 0 km s^{−1}. The different species can be extracted with the crosscorrelation technique. Figure 3 displays the crosscorrelation results of searching for species C in the transmission spectrum.
In the K_{p} – V_{sys} map, we see a central peak where we expect to see the planetary signal. However, two additional peaks lie at 50 km s^{−1} and 100 km s^{−1}. Converting to a 1D crosscorrelation function, we see that the central peak barely sits above the other two. Peaks of the crosscorrelation function no longer represent random noise but are artificial aliases caused by crosscorrelating with the other species that reside within the spectrum. Strong aliases thus carry the potential to bury signals of fainter species. This situation is most potent in two cases: for species with few lines and species with weak lines. In the case of only a few lines, each line position is given more weight in the crosscorrelation function, making the misscorrelation with other lines in the spectrum more significant. For faint species, the alias signals become more substantial since the more dominant species add a stronger contribution to the overall alias structure. For more dominant species with many lines, aliases are expected to be less pronounced. Since many lines contribute to the crosscorrelation function, the lines do not correlate as well with other absorption species in the spectrum, leading to a noise that much better approximates the theoretical photon noise. While all species have alias profiles, the contribution is more negligible for species with strong lines.
Observations have routinely demonstrated that signals from different species detected through crosscorrelation can have offsets in systematic velocity (Pai Asnodkar et al. 2022b). Peaks are known to deviate of the order of a few km s^{−1}, typically explained through atmospheric dynamics (Ehrenreich et al. 2020; Prinoth et al. 2022). This effect is slightly noticeable in Fig. 3, as species A has a peak close to the planetary rest frame at 2.8 km s^{−1}, which has increased the signalpeak height and produced a positive offset of 0.3 km s^{−1}. This result implies that other species can potentially contribute to the crosscorrelation signal and, if a signal is not present, appear to be one. This scenario outlines the need for an analysis method when applying crosscorrelation, especially for line lists with relatively few lines.
The root of the alias problem also betrays its solution. It is possible to predict their positions and relative strengths, and leverage those predictions to more accurately determine if an expected peak is, in fact, a signal. This can be done by crosscorrelating the model spectrum of the species suspected of producing aliases with the searchspecies templates to produce models of alias profiles. Additionally, correlating the search species template with itself creates an expected autocorrelation function. The sum of the aliases and the expected signal comprises the entire crosscorrelation function. Figure 3’s final panel contains the alias profiles of species A and B.
Assuming that the crosscorrelation function is a linear sum of the aliases and the true signal, we propose the following numerical model: (6)
where Y_{i} is the response variable and represents each data value i, k is the total number of predictor variables, X_{ki} are the predictor variables, β_{0} is a numerical intercept of the model, and β_{k} is a scaling factor, N is the total number of data values allowed into the model, and e_{i} is the error on each measurement. If e_{i} measurements are independent and identically distributed, with a zero mean E(e_{i}) = 0, and variance VAR(e_{i}) = σ^{2}, then regression can be used to estimate the β_{i} parameters. The aliasregression model is estimated by multiple least squares (MLS; Olive 2017).
Applying multiple least squares regression to this data provides estimates for the coefficient values whose standard errors can be measured through jackknife resampling (Quenouille 1949), and creates an empirical model to compare with the final results. The division of any coefficient value by its standard deviation provides a measure of the detection significance. Under idealised conditions, the error distribution follows a Student’s tdistribution and can therefore be subjected to a Wald ttest (Olive 2017). The use of tstatistics is routine in the verification of crosscorrelation signals (Brogi et al. 2012; Nugroho et al. 2017). Additionally, a statistic which describes the overall constraining power of the model is the number of degrees of freedom: (7)
where df is the degrees of freedom of the model, n is the total number of data points, and p is the total number of predictor variables. If df ≥ 30, a tdistribution can be accurately approximated to follow a Gaussian distribution N(0, 1) (Olive 2017). As is currently conducted in the literature, claiming a detection requires the tstatistic to be ≥ 5σ, weaker signals lower than this value are usually classed as ‘tentative’ detections.
It is essential to be aware that 1D crosscorrelation functions can strongly correlate with neighbouring data points depending on one’s choice of radial velocity step, and this can overconstrain any fit applied. Under the statistical assumptions applied in MLS, it is expected that the data points are independently distributed (Olive 2017). Therefore, we mitigate this issue by only fitting every fourth sample of the 1D crosscorrelation function (Hoeijmakers et al. 2019, 2020b), corresponding to a stepsize of 4 km s^{−1}.
Fig. 3 Crosscorrelation results from correlating species C with Species A, Species B, and the whole toy spectrum. Top left panel: the generated K_{p} – V_{sys} map using the toy spectrum. In the absence of noise, clear signals occur at the same orbital velocity but various systemic velocities, which cannot be related to the search species signal. Top right panel: the extracted 1D CCF, while a peak has formed at the expected location of the signal, the resulting crosscorrelation function is highly structured with two large peaks forming close to 50 km s^{−1} and 100 km s^{−1}. Bottom panel: the decomposition of the 1D crosscorrelation into its signal and alias components, we see that species A is largely responsible for the signal structure away from the planet’s signal, including the two offset peaks. The regression model specified in Eq. (8) has successfully reconstructed the 1D crosscorrelation function. 
Summarised tscores for the constructed toy model.
Applying aliasregression to the theoretical case
Following the outlined methodology, the following MLS regression model is constructed for the hypothetical sample spectrum for species A, B, and C: (8)
In this ideal case, the estimated alias profiles match the output of the toy model exactly. Implying that the coefficient values have a value of unity and do not need to be scaled, note the values in the parenthesis denote the standard error measurements of the coefficients. Figure 3 illustrates this model as a dashed line and demonstrates that the regression approach recreates the crosscorrelation function, and the significance of each coefficient in the model is summarised in Table 2. Template C has a tstatistic of t = 9.8, which passes the detection threshold (≥ 5σ). Another result of this fit is the detection significance of the A and B line profiles, both of which are detected, with species A obtaining the highest significance value. This suggests that aliases peaks carry enough information to enable detections by themselves, and an example of such aliases has previously been identified in Hoeijmakers et al. (2019).
Since the model reproduces the aliases; noise, systematics, or omitted alias contributors are likely responsible for deviations. The toy model presumes no stellar or other systematic contributions exist in the transmission spectrum, and it is common for residual stellar signals to be present within crosscorrelation maps (Hoeijmakers et al. 2020b). Finally, aliasregression also operates under the assumption that all alias contributors are known. As it is challenging to know every species that may contribute to the alias structure ahead of time, deviations from the model should be expected. All detections found in this paper were included in the aliasregression model. Detections follow the same detection criteria (σ) as for the toy model.
3.4 Validating detections
Candidate signals must undergo scrutiny to test their robustness. A good primary step is to confirm that the signals are not produced by systematic artefacts and originate exclusively from the intransit exposures, which is achieved through bootstrapping the in and outoftransit exposures (see Wyttenbach et al. 2015, and references therein). Following this confirmation, the aliasregression models can be statistically examined through an Analysis of Variance (ANOVA) Ftest and visually examined through response and residual plots.
3.4.1 Bootstrapping
We tested the robustness of detections through a bootstrapping approach applied to the 2D crosscorrelation functions as described in Hoeijmakers et al. (2020b). The method tested the distribution of candidate signals caused by correlated noise in the crosscorrelation functions to quantify the probability that the signal was caused by systematic noise, by selecting intransit values of the crosscorrelation function away from the expected planet signal and shifting these values to a random radial velocity. It was possible to test for this. Using a random Gaussian probability distribution, we sampled the shifts with 10 km s^{−1} and 20 km s^{−1} fixed widths. Without systematics, the profiles were expected to be Gaussiancentred at an average line depth of zero, with the same widths specified by the fixed widths used to sample the random shifts. The signal was validated by plotting the measured line strength of the signal as a vertical line, and the line should be distinctly offset from the two Gaussians if a signal was present. Any systematics present could cause an offset in the Gaussian histograms leading to overlapping with the signal, which would suggest that the signal came from this source as opposed to the planet signal.
We conducted an analysis to assess whether the signal originated from the intransit exposures. Beginning with the cleaned crosscorrelation function, the intransit and outoftransit exposures were separated into two groups and averaged. The average of the outoftransit group should have been zero, while the intransit one should have experienced an offset due to absorption by the planet’s atmosphere. Fifty percent of in and outoftransit exposures were randomly selected, averaged and divided by the master in and outoftransit crosscorrelation function. Repeating this process 20 000 times created three distributions of inin, inout, and outout residuals. Plotting these as histograms enabled one to identify if the signal indeed did come from the intransit exposures, as the inout transit distribution should have been offset from the inin and outout distributions, which should have been centered around zero. An offset confirmed that the signal came from the intransit exposures.
3.4.2 Testing aliasregression through Fstatistics, response, and residual plots
A MLS model comprised of the predicted signals and alias functions should describe the 1D crosscorrelation function. How well this is achieved can be evaluated through a response plot, residual plot, and an ANOVA Fstatistic. A response plot is a scatter plot of the predicted values of a model against the measured data values. In this plot, all values should be distributed around the identity line y = x, and measure how well the model represents the data, with the gradient representing the overall predictive power of the model (Olive 2017). The residual plot demonstrates the difference between the two models against the x_{i} data points. Ideally, these plots should show a random distribution of points along a y_{i} value of zero (Olive 2017).
The ANOVA Fstatistic examines if the predictors are needed in the MLS model, that is, whether or not the model Y_{i} offers a better prediction than the sample mean Ȳ. In other words, ANOVA tests the coefficients against the null hypothesis that β_{1} = … = β_{k} = 0 (Olive 2017). The Fstatistic is calculated as the ratio of the explained variance and the unexplained variance of the model (9)
where Ȳ_{i} denotes the sample mean of a specific predictor variable group, Y_{ij} is the jth observation out of the predictor variable groups K, n_{i} is the number of observations in the ith group, Ȳ denotes the overall mean of the data, and N is the overall sample size. The Fstatistic follows an Fdistribution, which does not approximate a normal distribution with higher degrees of freedom. The statistic will be significant if the explained variance is large compared to the unexplained variance, which is unlikely if the null hypothesis is true (Lowry 2014).
A onetailed hypothesis test with a 5σ cutoff under a Gaussian assumption has a probability of P(z ≥ 5σ) ≤ 5.724 × 10^{−7}. Software packages such as scipy (Virtanen et al. 2020), can calculate an equivalent probability value for the Fstatistic. This value will be used as the cutoff to evaluate its significance.
In addition to the Fstatistic, it is also possible to test how statistically significant some parameters of the model are in combination through the use of a reduced regression model (Olive 2017). In this case, we want to evaluate if the alias coefficients, in combination with each other, constrain the crosscorrelation function. Computing a second regression model omitting the searchspecies coefficient enables this. It is then possible to compute an Fstatistic which evaluates the constraining power of the reduced model using the following formula (10)
where RSS refers to the residual sum of squares for the reduced and full model, p is the number of predictors removed, n is the total number of observations in the data set, and k is the number of coefficients, including the intercept in the full model. The Fstatistic follows an F(d_{1}, d_{2}) with the degrees of freedom defined as (Olive 2017) (11)
As the alias profiles lie away from the planetary signal, they are better constrained through crosscorrelating a wide radialvelocity range. The range −1000 km s^{−1} to 1000 km s^{−1} was enough to achieve this.
4 Results
The four transit observations of KELT9 b were subject to a chemical survey using the procedures outlined in Sect. 3. Detections fell into one of two groups: confirmed detections or new detections. Confirmed detections are species previously found by Hoeijmakers et al. (2019), Yan et al. (2019) or Turner et al. (2020). The purpose was to show that the method reproduces these results while enhancing the S/N. We explicitly go into more detail with the Mg I detection, and the tentative Ba II detection to give the reader a clearer view of how the model fits the data. All detections are subject to the statistical verification techniques outlined in Sect. 3.4 with the figures located in Appendix C.
4.1 Detection confirmations
The original detections from Hoeijmakers et al. (2019), Yan et al. (2019) and Turner et al. (2020) were recovered, with the noted exception of Y II which is recovered at a significance of 2.7σ. The K_{p} – V_{sys} maps and 1D crosscorrelations for the species H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II are displayed in Fig. 4. The centroid velocities for this subset are close to the accepted systemic velocity values quoted in the literature (Pai Asnodkar et al. 2022b) with the notable exception of Ca II and H I, which have strongly blueshifted velocities. These signals likely originate from a higher region of the planet’s atmosphere. An additional point of interest for these plots is that there appears to be a secondary trough located around −50 km s^{−1}; this is likely a residual of the RossiterMcLaughlin effect (Rossiter 1924; McLaughlin 1924) caused by incomplete removal of the stellar signal.
Additionally, Sc II and Cr II detections are verified. While a signal appears present in the K_{p} – V_{sys} plots, their structures are irregular. However, the signals’ locations are at the expected velocities. The peak of Cr II appears to have some substructure, with a secondary peak at a slightly higher velocity, suggesting aliasing.
4.2 New detections
In addition to reconfirming the species already detected in the archived data, five new detections above the 5σ detection threshold are provided. These are Ca I, Cr I, Ni I, Sr II, and Tb II, and three tentative detections (above 3σ) for Ti I, V I, Ba II. There is what appears to be a secondary signal in the Sr II detection located at a K_{p} velocity of 130 km s^{−1}. This signal only occurs in the second observation using the HARPSN spectrograph, while the expected signal appears in both. Suggesting that the feature is specific to that observation and not the planet in general. The detection results are displayed in Fig. 5. The line lists for these species are plotted in Appendix A. All these species display significant aliasing in their K_{p} – V_{sys} and 1D crosscorrelation function, although the aliasing is not severe for Ca I. All the signals are located close to the accepted planetary signals with overlapping peak widths. Ionised Tb is the first detection of these species in any exoplanet atmosphere. Other observations have reported detections of the remaining seven, but never in KELT9 b. Figure 6 provides an outline of the measured peak positions and widths for each signal, at each K_{p} value, all values fall well within the averaging window selected in K_{p} – V_{sys} space. There appears to be significant variability in atmospheric signals, with signals clustering into three groups. This may be due to the signals occurring at different atmospheric layers.
4.3 Verification tests
Figure 7 provides a summary of the statistical significance of each search species, and the largest tscore of the alias coefficients for each model, while Table B.1 provides a summary of the statistical significance of all coefficient values for the aliasregression models constructed for each species. Most of these alias coefficient values lie below the significance threshold of 5σ, with the exception of the predicted aliases of Fe II with searching for Mg I and Ba II which are detected at significance’s of 8.4σ and 8.0σ, respectively. However, there are multiple occurrences of alias coefficients obtaining significance values above 3σ, which suggests these artefacts are contributing to the overall signal of the crosscorrelation function. The alias modelling appears to have aided in constraining the significance of the detections of each species.
4.3.1 The aliases of magnesium and ionised barium
Figure 8 illustrates each predicted aliascontribution for the aliasregression model for Mg I and Ba II, overplotted with the 1D crosscorrelation function and the model itself. The models capture the Mg I and Ba II signals and multiple Fe II aliases through the crosscorrelation function. In the panels for Mg I and Ba II in Fig. C.3 we show that the response plots have gradients of 0.406 and 0.278 respectively. This suggests that the model explains a significant proportion of the underlying crosscorrelation function, and the residual plot shows that most features have been removed. The other coefficients of the functions do not reach a significance above 3σ; however, in combination, they constrain the data to a statistically significant degree (see Table 3). The results show that regressing alias profiles adds constraining power to crosscorrelation detections.
4.3.2 Bootstrapping to confirm the location of the signal
Bootstrapping results for the eight new species are displayed in Figs. C.1 and C.2. All detections show that the calculated line strength of the signal is visually offset from the two sampled Gaussians. The offset suggests that any systematics sampled from this bootstrapping approach are not able to explain the signals seen in the K_{p} – V_{sys} maps of the newly detected signals. It is also the case for the confirmed detections seen in Fig. 4.
4.3.3 Verifying the aliasregression models
The ANOVA Fstatistic of the aliasregression models are presented in Table 3. All models have Fstatistics which generate pvalues that are more significant than above a 5σ threshold of the normal distribution. The result suggests that all coefficient values in combination add constraining power to the 1D crosscorrelation function. The results presented by the Fstatistics agree with the results presented in Figs. C.3 and C.6, all models obtain positive gradients, with lines of worst fit regions away from zero. The strongest gradient of the species is H I, which almost has a onetoone modeltodata relationship, while the weakest is Tb II which, while obtaining a positive gradient, leaves much of the residual data left unexplained. Additionally, the residual plots show that the models remove the signals partially (see Fig. C.3). Based on the results of this section and the bootstrapping, it can be concluded that the new detections are possible detections, but the models could be improved to better explain the signals. Table 3 also contains reduced model Fstatistics; in all cases, the alias coefficients in combination remain statistically significant above the equivalent >5σ value. Moreover, the statistics are much more significant for the confirmed species than for the newly detected species, which suggests that the crosscorrelation functions for these results contain additional variation that the model does not explain.
Fig. 4 Confirmed species detections in KELT9 b. Each panel presents each species’ 2D and 1D crosscorrelation signatures in the following order: H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. In the K_{p} – V_{sys} maps, the white dashed lines locate the maximum peak of the signal. Onedimensional crosscorrelation functions are extracted from the K_{p} – V_{sys} map by selecting the row of pixels at the maximal location of the signal. The dashed dot represents a Gaussian fit to the crosscorrelation function, and the vertical dotted lines mark the location of the peak for each example. The blueshaded region represents the standard deviation of all points away from the peak of each function. 
Fig. 5 Newly detected species detections in KELT9 b. Each panel presents each species’ 2D and 1D crosscorrelation signatures in the following order: Ca I, V I, Cr I, Ni I, Sr II, Ba II, and Tb II. The plots are constructed in the same manner as Fig. 4. 
5 Discussion
The purpose of aliasregression is to explain the signals in crosscorrelation maps. The methods and validation techniques in this study illuminate a path forwards that can potentially yield detections of such species. Based on the results of this study, new observations should be combined with archival data to increase the S/N. Atmospheric observations of exoplanets are photon starved, motivating the development of techniques that combine existing and new observations to increase the contrast of known signals. This method also has the potential to produce additional detections and reduce the photon noise in the data to the point where the underlying systematics are revealed.
We conclude that aliases are present in the crosscorrelation functions used to detect the ionised, atomic and molecular species in exoplanet atmospheres with high S/N. Their impact is that any fitting procedure that treats the region around the signal of the crosscorrelation function solely as noise will bias the detection statistics. K_{p} – V_{sys} maps, if scaled to magnitudes of standard deviations, may in such cases produce detection significance values above the 3σ level away from the planetary signal (Giacobbe et al. 2021; Guilluy et al. 2022). Instead, it is recommended to present K_{p} – V_{sys} maps in terms of their actual numerical output and employ additional methods to infer a signal’s validity. Aliasregression, model injection (Hoeijmakers et al. 2018a, 2019, 2020b), and bootstrapping approaches (Giacobbe et al. 2021; Prinoth et al. 2022) are ways to achieve this. A method that acknowledges and accounts for the systematic effects will likely hold more constraining power.
A profound implication of aliases is that they cannot be mitigated by using larger telescopes or repeated observations. However, improved S/N observations can provide additional insight. With additional observations of the same target, aliases will begin to appear and dominate the noise in the crosscorrelation function. Therefore, this calls for methods such as aliasregression to be applied, in order to discover new species that may be present in the atmospheres of these exoplanets (Kesseli et al. 2022).
We have successfully combined four observations of KELT9 b and used aliasregression to confirm the previous detections of Hoeijmakers et al. (2018a, 2019); Yan et al. (2019); Turner et al. (2020), with the exception of Y II, for which we obtained a significance value of 2.7. The discrepancy in the detection significance of Y II is likely due to the fact that Hoeijmakers et al. (2019) constrained their detection significance by a leastsquares Gaussian fit without considering the alias component of the crosscorrelation function. Thus, they record a detection significance of 6.3 σ. In a leastsquares approach, the omission of an explanatory variable leads to bias in the remaining coefficients. In this specific case, an alias peak for Mg I was found to be in the same position as the signal. Mg I has an alias coefficient value of 2.3 indicating that this species is the likely source of contamination. In this situation, the variance of the estimated coefficients decreases, which can lead to inflated significance values. The reduction in the significance of Y II suggests that the signal was likely constraining some alias contribution.
Furthermore, we have obtained eight discoveries of species unobserved in KELT9 b’s atmosphere, with varying degrees of confidence. Notable detections include Tb II which has not yet been found on any exoplanet, and the tentative detection of Ba II at a significance level of 4.5σ. Ba II has only recently been discovered in the ultrahot gas giants of WASP121 b and WASP76 b (Azevedo Silva et al. 2022). Tb was not predicted to be present in the atmospheres of these planets, nor has detection been proposed. The detection significance of Tb II is 5.0σ, and the extracted orbital velocity is consistent with the positions of the H I, Ca II, and V I detections. Taken together, the significance of the detection and the consistency of the orbital velocities suggest that this detection is sound. However, followup observations should verify the repeatability of this detection. Furthermore, Ca I and Cr I are also detected, increasing the known number of neutralion species pairs from one to three. Additional neutralion pairs will allow us to potentially unlock the ability to empirically constrain the ionisation rate of the same species, which can be used to determine the temperature profile and diagnose nonLTE effects. Moreover, apart from being detected by its own template, Fe II has also been observed via the aliases of Mg I and Ba II with an 8σ level of confidence, this was observed but not statistically verified in Hoeijmakers et al. (2019). It is important to note that while this is the first confirmed detection of a species via the alias structures in the crosscorrelation function, the use of aliasregression does not replace direct signal detections. The implication of the detection of Fe II via aliases is that the alias structures are present in crosscorrelation functions and must be accounted for when detecting fainter species.
The K_{p} – V_{sys} maps in this study vary considerably from clear bright features such as the Fe II and Ca II detections to signals that are less clear and contain many substructures such as the Ni I and Tb II detections. However, there seems to be a discrepancy between what the detection appears as and what the statistics claim. The discrepancy likely lies in the fact that statistical tests perform hypothesis testing under the assumption that there is no underlying systematics at play. This assumption is tenuous with K_{p} – V_{sys} maps, since aliases have been shown to exist. These systematics are likely to compromise the statistical rigour associated with the 3 and 5 σ rulesofthumb used for detection, and do not represent the confidence generally associated with them. However, there is still a clear trend between detection significance and signal clarity, suggesting that we should be more confident with higher statistics. Due to the arbitrary nature of specifying a statistical threshold in the first place, it would be better to use them as a guide rather than the rule when claiming detections, favouring a conservative approach where possible.
The detection of heavier isotopes is important because they can provide a constraint on the age of a planet. This has been shown to be the case for Jupiter in our solar system, which is essential for understanding the presentday architecture of the Solar System (Kruijer et al. 2017). Age determination can constrain the formation and migration scenarios of these planets. For KELT9 b, this could lead to a better understanding of the total mass lost due to the atmospheric accretion exhibited by its host star (SánchezLópez et al. 2022), and a better age constraint for the star itself. Detections are the first footfall in determining the abundances of these heavy isotopes. With better constraints on the systematics and additional followup observations, it may be possible to determine abundances and, by proxy, the age of KELT9 b.
Aliasregression works in both theoretical and observational cases, albeit with some limitations. In this study, the models showed some success in constraining the overall detected signals, but do not account for all of the peaks seen in the 1D crosscorrelation function. However, as shown in Fig. 8, aliasregression is successful in predicting some alias peaks. Deviations from the model could be caused by a variety of sources, including unaccounted systematics, and not yet included alias species. Additionally, templates are functionally dependent on many different input parameters, with line strength depending on global parameters such as temperature, abundances, and planetary bulk parameters. Therefore, alias profiles relate how well the templates match the transmission spectrum of the planet in question, similar to the atmospheric retrievals of Brogi & Line (2019); Gibson et al. (2022).
It is also prudent to consider that the photon noise may still contribute a proportionate amount to the crosscorrelation function, adding to the uncertainty when fitting the aliasregression coefficients. For example, species with few lines absorb less flux and therefore sample more noise in their crosscorrelation, reducing statistical robustness. However, each model applied successfully constrains the 1D signals. The results obtained with the reduced model support this case. Observing Table B.1, we see that the statistical significance of the individual coefficient values of the regression is often small, but in combination, they yield a high significance. Therefore, modelling aliases will constrain the CCF more than if they were omitted when determining the significance of a detection.
We conclude that combining the four nights of observations only begins to resolve the alias floor, making it difficult to fit the aliasregression model accurately. Repeated observations with larger telescopes and additional consideration of the shape of the alias profiles should improve these fits.
Fig. 6 Fitted Gaussian parameters for the systemic velocity centres and widths, obtained from the extracted orbital velocity from the K_{p} – V_{sys} plots. The detections appear to cluster into three groups, suggesting they may be probing different atmospheric layers, which have been colourcoded for easier identification. The median systemic velocity is represented by a grey line on the plot. 
ANOVA Fstatistics for the full aliasregression model and the reduced model.
Fig. 7 Summary of all detection statistics and the largest σ values for the alias coefficients of each aliasregression model. The dashed lines mark the 3, 5, and 10 σ detection significance threshold. The values are extracted from Table B.1. 
Fig. 8 Expected aliases for the 1D crosscorrelation function for Mg I, those with a < 3σ coefficient significance are plotted in grey. The measured 1D crosscorrelation function is shown in teal, and the model fit in dashedblack. The significant alias structure of Fe II is plotted in pink, and clearly contributes to the 1D crosscorrelation function. The remaining predicted nonsignificant aliases is shown in grey. 
6 Conclusion
This work has provided a new approach to the crosscorrelation technique, which has led to the detection of eight new chemical species in the atmosphere of KELT9 b. The method combines the information from multiple spectrographs with different wavelength ranges to produce strong signals in high contrast to the background noise. In addition, by attenuating the noise, alias profiles are resolved, revealing a systematic noise floor in the crosscorrelation function. aliasregression accounts for these aliases and provides a statistical measure of the detection significance, and has led to novel detections of eight undiscovered species in KELT9 b’s atmosphere, as well as reconfirmed detections from previous studies.
Despite its strengths, the approach has some limitations. While the aliasregression models have been shown to have the constraining power of the 1D crosscorrelation signals, they leave much of the crosscorrelation structure unexplained for species with fainter signals. This may be due to the presence of alternative aliasproducing species that are not accounted for, or to other sources of systematic error that have yet to be discovered. As with any statistical inference (model fitting) technique, the effectiveness of the aliasregression technique is limited by the selection of all the relevant alias contributors. In reality, relevant species needed to fully constrain the crosscorrelation function may have been excluded.
Acknowledgements
N.W.B. and B.P. acknowledge partial financial support from The Fund of the Walter Gyllenberg Foundation. B.T. acknowledges the financial support from the Japan Society for the Promotion of Science as a JSPS International Research Fellow. R.F. acknowledges support from the Göran Gustafsson Foundation for Research in Natural Sciences and Medicine and from the Royal Physiographic Society in Lund through The Märta and Eric Holmberg Endowment. D.K. acknowledges financial support from the Center for Space and Habitability (CSH) of the University of Bern. This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. Additionally, this article is based on observations made in the Observatorios de Canarias del IAC with the Telescopio Nazionale Galileo operated on the island of La Palma by the Fundación Galileo Galilei – INAF in the Observatorio del Roque de los Muchachos, and observations collected at the GermanSpanish Astronomical Center, Calar Alto, jointly operated by the MaxPlanckInstitut für Astronomie Heidelberg and the Instituto de Astrofísica de Andalucía (CSIC). We would also like to thank Y. BokhariFriberg for their help in editing and typesetting the manuscript, and the anonymous referee for their helpful contributions.
Appendix A Weight results
An overview of the weight calculation process is presented in Fig. A.1 for two species of Fe I (top panel) and Na I (bottom panel). Atomic iron contains many lines, with a larger flux absorption in the blue wavelengths. From Table 1, the HARPSN nights have larger S/N values, and this explains why the crosscorrelation peaks are much more significant in the two HARPSN nights, and has led to the weights favour these nights. For atomic sodium, however, the CARMENES spectrograph captures a deep triplet in the sodium spectrum at a larger wavelength range. In this case, while the CARMENES data is still poor compared to HARPSN, it contributes a proportional amount to the signal. Table A.1 and Table A.2 decomposes the weights of each night into the pseudoSNR values and the relative flux summation to better clarify how the weights are constructed. Fig. A.2 contains the line lists for all the species detected in this report to give the reader a better idea of which spectrograph performs better with each species.
Fig. A.1 Summary of the weight extraction process. The first two rows illustrate the weight construction for Fe I. Top panel: Line list of neutral iron used for crosscorrelation. The wavelength ranges of the HARPSN and CARMENES visible are plotted in boxes to illustrate the difference in wavelength coverage. Bottom panel: 1D crosscorrelation functions for the four nights of observations. The two HARPSN nights have much higher SNR and contain less background noise, producing larger weights than the CARMENES observations. The second two rows demonstrate the weight constructions for Na I. Top panel: Line list for neutral sodium. The proportion of flux obtained by both spectrographs is comparable to each other, therefore data quality dominates the weight calculation when combined. Bottom panel: Crosscorrelation functions for each night, CARMENES nights obtain a higher proportion in the calculation of weights. 
Fig. A.2 Generated line lists used in the crosscorrelation of this report. The y axis denotes transit depth, with zero equal to the stellar outoftransit flux. The species order ascends first with atomic number, then atoms before ionised species. The included wavelength region span the HARPSN and CARMENES visible arms. 
Computed values deriving the weights for combining the flux of each spectral observation for the species Fe I
Computed values deriving the weights for combining the flux of each spectral observation for the species Na I
Appendix B Detection attributes
Table B.1 presents a summary of all tscores of the statistically significant searches and aliasregression coefficients. The majority of the alias coefficients lie below a significance of 3σ; however, a few obtain values above this threshold. The species which include significant detections above this threshold are Mg I, Sc II, Cr II, V I, and Ba II.
Individual tstatistics for the coefficient values for each multipleleast squares aliasregression function.
Appendix C Model verification for known species
All species were subjected to verification tests outlined in Section 3. Fig. C.1 identifies clear Gaussian profiles which reside well away from the measured line depths, implying the signals come from the planet as opposed to some underlying systematic structure. Fig. C.3 depict the response plots and residual plots for the previously detected species in KELT9 b’s atmosphere. The model’s explanatory power varies, from almost obtaining a onetoone relationship with H I, to Tb II, which only can explain approximately 8.9 % of the variation in the crosscorrelation map. However, the error bars on all the gradient values obtain values above zero, implying that these models do have explanatory power.
Fig. C.1 Bootstrap results for the confirmed detections. Each panel represents the bootstrap results of each new detection in the respective order of: H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. The dark blue histogram represents the random Gaussian sample obtained with the 10 km s^{−1} width, while the light transparent blue histogram is for the distribution width of 20 km s^{−1}. The solid blue vertical line indicates the measured line depth of the signal, with a standard deviation of the value present as a transparent blueshaded region. The legend in the right panel plots labels the inin, inout and outout residual distributions calculated during the bootstrapping process. 
Fig. C.2 Bootstrap results for the new detections. Each panel represents the bootstrap results of each new detection in the respective order of Ca I, V I, Ti I, Cr I, Ni I, Sr II, Ba II, and Tb II. The plots follow the same construction as Fig. C.3 
Fig. C.3 Response and residual plots for each confirmed species in KELT9 b’s atmosphere. Each panel represents a different species. The order of the species is as follows (left to right): H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. Each plot contains a line of best fit plotted through the data of the response plot, spanned by two lines of worst fit. The identity line y = x is plotted to show how closely the model matches the data. The residual plot contains the residuals and the region spanning one standard deviation of the residual data. 
Fig. C.4 Continued: Response and residual plots for each confirmed species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Se II, Ti II, Cr II, and Fe I. 
Fig. C.5 Continued: Response and residual plots for each confirmed species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Fe II and Y II. 
Fig. C.6 Response and residual plots for each newly detected species in KELT9 b’s atmosphere. Each panel represents a different species. The order of the species is as follows (left to right): Ca I, Ti I, V I, and Cr I. The format of the plot follows the same structure as Fig. C.3. 
Fig. C.7 Continued: Response and residual plots for each newly detected species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Ni I, Sr II, Ba II, and Tb II. The format of the plot follows the same structure as Fig. C.3. 
References
 AlonsoFloriano, F. J., SánchezLópez, A., Snellen, I. A. G., et al. 2019a, A&A, 621, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 AlonsoFloriano, F. J., Snellen, I. A. G., Czesla, S., et al. 2019b, A&A, 629, A110 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Azevedo Silva, T., Demangeon, O. D. S., Santos, N. C., et al. 2022, A&A, 666, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 BelloArufe, A., Buchhave, L. A., Mendonça, J. M., et al. 2022, A&A, 662, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Birkby, J. L., de Kok, R. J., Brogi, M., Schwarz, H., & Snellen, I. A. G. 2017, AJ, 153, 138 [NASA ADS] [CrossRef] [Google Scholar]
 Booth, R. A., & Owen, J. E. 2020, MNRAS, 493, 5079 [NASA ADS] [CrossRef] [Google Scholar]
 Borsa, F., Rainer, M., Bonomo, A. S., et al. 2019, A&A, 631, A34 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Brogi, M., & Line, M. R. 2019, AJ, 157, 114 [Google Scholar]
 Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2012, Nature, 486, 502 [Google Scholar]
 Brogi, M., Snellen, I. A. G., de Kok, R. J., et al. 2013, ApJ, 767, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Brogi, M., de Kok, R. J., Birkby, J. L., Schwarz, H., & Snellen, I. A. G. 2014, A&A, 565, A124 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Budde, G., Burkhardt, C., Brennecka, G. A., et al. 2016, Earth Planet. Sci. Lett., 454, 293 [CrossRef] [Google Scholar]
 Budde, G., Burkhardt, C., & Kleine, T. 2019, Nat. Astron., 3, 736 [NASA ADS] [CrossRef] [Google Scholar]
 Burbidge, E. M., Burbidge, G. R., Fowler, W. A., & Hoyle, F. 1957, Rev. Mod. Phys., 29, 547 [NASA ADS] [CrossRef] [Google Scholar]
 Burkhardt, C., Kleine, T., Oberli, F., et al. 2011, Earth Planet. Sci. Lett., 312, 390 [CrossRef] [Google Scholar]
 Cabot, S. H. C., BelloArufe, A., Mendonça, J. M., et al. 2021, AJ, 162, 218 [NASA ADS] [CrossRef] [Google Scholar]
 CasasayasBarris, N., Pallé, E., Yan, F., et al. 2019, A&A, 628, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cauley, P. W., Shkolnik, E. L., Ilyin, I., et al. 2019, AJ, 157, 69 [Google Scholar]
 Chambers, J., & Wetherill, G. 1998, Icarus, 136, 304 [NASA ADS] [CrossRef] [Google Scholar]
 Charbonneau, D., Brown, T. M., Noyes, R.W., & Gilliland, R. L. 2002, ApJ, 568, 377 [NASA ADS] [CrossRef] [Google Scholar]
 Cont, D., Yan, F., Reiners, A., et al. 2021, A&A, 651, A33 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Dauphas, N., Remusat, L., Chen, J. H., et al. 2010, ApJ, 720, 1577 [NASA ADS] [CrossRef] [Google Scholar]
 de Kok, R. J., Brogi, M., Snellen, I. A. G., et al. 2013, A&A, 554, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Lecavelier Des Etangs, A., Hébrard, G., et al. 2008, A&A, 483, 933 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ehrenreich, D., Lovis, C., Allart, R., et al. 2020, Nature, 580, 597 [Google Scholar]
 Ek, M., Hunt, A. C., Lugaro, M., & Schönbächler, M. 2020, Nat. Astron., 4, 273 [NASA ADS] [CrossRef] [Google Scholar]
 Fisher, C., Hoeijmakers, H. J., Kitzmann, D., et al. 2020, AJ, 159, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Gaudi, B. S., Stassun, K. G., Collins, K. A., et al. 2017, Nature, 546, 514 [NASA ADS] [Google Scholar]
 GharibNezhad, E., Iyer, A. R., Line, M. R., et al. 2021, ApJS, 254, 34 [NASA ADS] [CrossRef] [Google Scholar]
 Giacobbe, P., Brogi, M., Gandhi, S., et al. 2021, Nature, 592, 205 [NASA ADS] [CrossRef] [Google Scholar]
 Gibson, N. P., Nugroho, S. K., Lothringer, J., Maguire, C., & Sing, D. K. 2022, MNRAS, 512, 4618 [NASA ADS] [CrossRef] [Google Scholar]
 Grimm, S. L., & Heng, K. 2015, ApJ, 808, 182 [NASA ADS] [CrossRef] [Google Scholar]
 Guilluy, G., Giacobbe, P., Carleo, I., et al. 2022, A&A, 665, A104 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Heng, K. 2016, ApJ, 826, L16 [NASA ADS] [CrossRef] [Google Scholar]
 Heng, K., & Kitzmann, D. 2017, MNRAS, 470, 2972 [Google Scholar]
 Hoeijmakers, H. J., de Kok, R. J., Snellen, I. A. G., et al. 2015, A&A, 575, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hoeijmakers, H. J., Ehrenreich, D., Heng, K., et al. 2018a, Nature, 560, 453 [CrossRef] [Google Scholar]
 Hoeijmakers, H. J., Schwarz, H., Snellen, I. A. G., et al. 2018b, A&A, 617, A144 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hoeijmakers, H. J., Ehrenreich, D., Kitzmann, D., et al. 2019, A&A, 627, A165 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hoeijmakers, H. J., Cabot, S. H. C., Zhao, L., et al. 2020a, A&A, 641, A120 [EDP Sciences] [Google Scholar]
 Hoeijmakers, H. J., Seidel, J. V., Pino, L., et al. 2020b, A&A, 641, A123 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 John, T. L. 1988, A&A, 193, 189 [NASA ADS] [Google Scholar]
 Kesseli, A. Y., Snellen, I. A. G., CasasayasBarris, N., Mollière, P., & SánchezLópez, A. 2022, AJ, 163, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Kitzmann, D., Heng, K., Rimmer, P. B., et al. 2018, ApJ, 863, 183 [Google Scholar]
 Kitzmann, D., Hoeijmakers, H. J., Grimm, S. L., et al. 2023, A&A, 669, A113 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kruijer, T. S., Burkhardt, C., Budde, G., & Kleine, T. 2017, Proc. Natl. Acad. Sci. U.S.A., 114, 6712 [Google Scholar]
 Kurucz, R. L. 2018, ASP Conf. Ser., 515, 47 [NASA ADS] [Google Scholar]
 Lampón, M., LópezPuertas, M., SanzForcada, J., et al. 2021, A&A, 647, A129 [EDP Sciences] [Google Scholar]
 Landman, R., SánchezLópez, A., Mollière, P., et al. 2021, A&A, 656, A119 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Liu, F., Yong, D., Asplund, M., et al. 2020, MNRAS, 495, 3961 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, F., Bitsch, B., Asplund, M., et al. 2021, MNRAS, 508, 1227 [NASA ADS] [CrossRef] [Google Scholar]
 Lowry, R. 2014, Concepts and Applications of Inferential Statistics, 1st edn. (Lowry, Richard) [Google Scholar]
 Lowson, N., Zhou, G., Wright, D. J., et al. 2023, AJ, 165, 101 [NASA ADS] [CrossRef] [Google Scholar]
 Lugaro, M., Pignatari, M., Ott, U., et al. 2016, Proc. Natl. Acad. Sci. U.S.A., 113, 907 [NASA ADS] [CrossRef] [Google Scholar]
 Mayor, M., Pepe, F., Queloz, D., et al. 2003, The Messenger, 114, 20 [NASA ADS] [Google Scholar]
 McLaughlin, D. B. 1924, ApJ, 60, 22 [Google Scholar]
 Meléndez, J., Asplund, M., Gustafsson, B., & Yong, D. 2009, ApJ, 704, L66 [Google Scholar]
 Merritt, S. R., Gibson, N. P., Nugroho, S. K., et al. 2021, MNRAS, 506, 3853 [NASA ADS] [CrossRef] [Google Scholar]
 Nugroho, S. K., Kawahara, H., Masuda, K., et al. 2017, AJ, 154, 221 [Google Scholar]
 Nugroho, S. K., Kawahara, H., Gibson, N. P., et al. 2021, ApJ, 910, L9 [CrossRef] [Google Scholar]
 Olive, D. J. 2017, Linear Regression, 1st edn. (Springer Cham) [CrossRef] [Google Scholar]
 Pai Asnodkar, A., Wang, J., Eastman, J. D., et al. 2022a, AJ, 163, 155 [NASA ADS] [CrossRef] [Google Scholar]
 Pai Asnodkar, A., Wang, J., Gaudi, B. S., et al. 2022b, AJ, 163, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Pelletier, S., Benneke, B., DarveauBernier, A., et al. 2021, AJ, 162, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Pino, L., Désert, J.M., Brogi, M., et al. 2020, ApJ, 894, L27 [Google Scholar]
 Pino, L., Brogi, M., Désert, J. M., et al. 2022, A&A, 668, A176 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icarus, 124, 62 [NASA ADS] [CrossRef] [Google Scholar]
 Pont, F., & Eyer, L. 2004, MNRAS, 351, 487 [Google Scholar]
 Prinoth, B., Hoeijmakers, H. J., Kitzmann, D., et al. 2022, Nat. Astron., 6, 449 [NASA ADS] [CrossRef] [Google Scholar]
 Quenouille, M. H. 1949, Ann. Math. Stat., 20, 355 [Google Scholar]
 Quirrenbach, A., Amado, P. J., Caballero, J. A., et al. 2014, SPIE Conf. Ser., 9147, 91471F [Google Scholar]
 Richard, C., Gordon, I. E., Rothman, L. S., et al. 2012, JQSRT, 113, 1276 [NASA ADS] [CrossRef] [Google Scholar]
 Rossiter, R. A. 1924, ApJ, 60, 15 [Google Scholar]
 Ryabchikova, T., Piskunov, N., Kurucz, R. L., et al. 2015, Phys. Scr, 90, 054005 [Google Scholar]
 Sahlholdt, C. L., Feltzing, S., Lindegren, L., & Church, R. P. 2019, MNRAS, 482, 895 [Google Scholar]
 Saji, N. S., Wielandt, D., Holst, J. C., & Bizzarro, M. 2020, Geochim. Cosmochim. Acta, 281, 135 [NASA ADS] [CrossRef] [Google Scholar]
 Salz, M., Czesla, S., Schneider, P. C., et al. 2018, A&A, 620, A97 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 SánchezLópez, A., Lin, L., Snellen, I. A. G., et al. 2022, A&A, 666, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Sedaghati, E., MacDonald, R. J., CasasayasBarris, N., et al. 2021, MNRAS, 505, 435 [NASA ADS] [CrossRef] [Google Scholar]
 Seidel, J. V., Ehrenreich, D., Bourrier, V., et al. 2020, A&A, 641, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Seidel, J. V., Cegla, H. M., Doyle, L., et al. 2022, MNRAS, 513, L15 [NASA ADS] [CrossRef] [Google Scholar]
 Smette, A., Sana, H., Noll, S., et al. 2015, A&A, 576, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Nature, 465, 1049 [Google Scholar]
 Stangret, M., CasasayasBarris, N., Pallé, E., et al. 2022, A&A, 662, A101 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Stevenson, K. B., Bean, J. L., Seifahrt, A., et al. 2016, ApJ, 817, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Stock, J.W., Kitzmann, D., Patzer, A. B. C., & Sedlmayr, E. 2018, MNRAS, 479, 865 [NASA ADS] [Google Scholar]
 Stock, J. W., Kitzmann, D., & Patzer, A. B. C. 2022, MNRAS, 517, 4070 [NASA ADS] [CrossRef] [Google Scholar]
 Turner, J. D., de Mooij, E. J. W., Jayawardhana, R., et al. 2020, ApJ, 888, L13 [Google Scholar]
 van Sluijs, L., Birkby, J. L., Lothringer, J., et al. 2023, MNRAS, 522, 2145 [NASA ADS] [CrossRef] [Google Scholar]
 Virtanen, P., Gommers, R., Oliphant, T. E., et al. 2020, Nat. Methods, 17, 261 [Google Scholar]
 Wyttenbach, A., Ehrenreich, D., Lovis, C., Udry, S., & Pepe, F. 2015, A&A, 577, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Wyttenbach, A., Lovis, C., Ehrenreich, D., et al. 2017, A&A, 602, A36 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yan, F., & Henning, T. 2018, Nat. Astron., 2, 714 [Google Scholar]
 Yan, F., CasasayasBarris, N., Molaverdikhani, K., et al. 2019, A&A, 632, A69 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Žák, J., Kabáth, P., Boffin, H. M. J., Ivanov, V. D., & Skarka, M. 2019, AJ, 158, 120 [Google Scholar]
All Tables
Computed values deriving the weights for combining the flux of each spectral observation for the species Fe I
Computed values deriving the weights for combining the flux of each spectral observation for the species Na I
Individual tstatistics for the coefficient values for each multipleleast squares aliasregression function.
All Figures
Fig. 1 Overview of the crosscorrelation cleaning process. Panel A: the raw crosscorrelation function timeseries applied after the outoftransit spectra have been divided and the residual broadband features removed. Both the Doppler shadow and the alias structures are present. Panel B: the planet mask covers the expected planet signal, the halfwidth of the planet mask is 20 km s^{−1}. Panel C: the fitted Dopplershadow profile is a twocomponent Gaussian. Panel D: the crosscorrelation function after removal of the Doppler shadow, residual stellar aliases are still present. Panel E: the resulting crosscorrelation map after removing the stellar aliases. The signature of the planet appears as the dark slanted trace. In this plot, the sign is negative; however, it is flipped in preceding figures to denote absorption. 

In the text 
Fig. 2 Comparison plot demonstrating the improved signal from applying the combination method for the species Fe I. Left panel: K_{p} – V_{sys} map of a single transit observation, this specifically is the first HARPSN transit. The red border marks the window where the average signal strength is computed, and the region outside where the standard deviation is computed. Right panel: K_{p} – V_{sys} map when combining all four transit observations. The result is a much higher S/N. 

In the text 
Fig. 3 Crosscorrelation results from correlating species C with Species A, Species B, and the whole toy spectrum. Top left panel: the generated K_{p} – V_{sys} map using the toy spectrum. In the absence of noise, clear signals occur at the same orbital velocity but various systemic velocities, which cannot be related to the search species signal. Top right panel: the extracted 1D CCF, while a peak has formed at the expected location of the signal, the resulting crosscorrelation function is highly structured with two large peaks forming close to 50 km s^{−1} and 100 km s^{−1}. Bottom panel: the decomposition of the 1D crosscorrelation into its signal and alias components, we see that species A is largely responsible for the signal structure away from the planet’s signal, including the two offset peaks. The regression model specified in Eq. (8) has successfully reconstructed the 1D crosscorrelation function. 

In the text 
Fig. 4 Confirmed species detections in KELT9 b. Each panel presents each species’ 2D and 1D crosscorrelation signatures in the following order: H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. In the K_{p} – V_{sys} maps, the white dashed lines locate the maximum peak of the signal. Onedimensional crosscorrelation functions are extracted from the K_{p} – V_{sys} map by selecting the row of pixels at the maximal location of the signal. The dashed dot represents a Gaussian fit to the crosscorrelation function, and the vertical dotted lines mark the location of the peak for each example. The blueshaded region represents the standard deviation of all points away from the peak of each function. 

In the text 
Fig. 5 Newly detected species detections in KELT9 b. Each panel presents each species’ 2D and 1D crosscorrelation signatures in the following order: Ca I, V I, Cr I, Ni I, Sr II, Ba II, and Tb II. The plots are constructed in the same manner as Fig. 4. 

In the text 
Fig. 6 Fitted Gaussian parameters for the systemic velocity centres and widths, obtained from the extracted orbital velocity from the K_{p} – V_{sys} plots. The detections appear to cluster into three groups, suggesting they may be probing different atmospheric layers, which have been colourcoded for easier identification. The median systemic velocity is represented by a grey line on the plot. 

In the text 
Fig. 7 Summary of all detection statistics and the largest σ values for the alias coefficients of each aliasregression model. The dashed lines mark the 3, 5, and 10 σ detection significance threshold. The values are extracted from Table B.1. 

In the text 
Fig. 8 Expected aliases for the 1D crosscorrelation function for Mg I, those with a < 3σ coefficient significance are plotted in grey. The measured 1D crosscorrelation function is shown in teal, and the model fit in dashedblack. The significant alias structure of Fe II is plotted in pink, and clearly contributes to the 1D crosscorrelation function. The remaining predicted nonsignificant aliases is shown in grey. 

In the text 
Fig. A.1 Summary of the weight extraction process. The first two rows illustrate the weight construction for Fe I. Top panel: Line list of neutral iron used for crosscorrelation. The wavelength ranges of the HARPSN and CARMENES visible are plotted in boxes to illustrate the difference in wavelength coverage. Bottom panel: 1D crosscorrelation functions for the four nights of observations. The two HARPSN nights have much higher SNR and contain less background noise, producing larger weights than the CARMENES observations. The second two rows demonstrate the weight constructions for Na I. Top panel: Line list for neutral sodium. The proportion of flux obtained by both spectrographs is comparable to each other, therefore data quality dominates the weight calculation when combined. Bottom panel: Crosscorrelation functions for each night, CARMENES nights obtain a higher proportion in the calculation of weights. 

In the text 
Fig. A.2 Generated line lists used in the crosscorrelation of this report. The y axis denotes transit depth, with zero equal to the stellar outoftransit flux. The species order ascends first with atomic number, then atoms before ionised species. The included wavelength region span the HARPSN and CARMENES visible arms. 

In the text 
Fig. C.1 Bootstrap results for the confirmed detections. Each panel represents the bootstrap results of each new detection in the respective order of: H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. The dark blue histogram represents the random Gaussian sample obtained with the 10 km s^{−1} width, while the light transparent blue histogram is for the distribution width of 20 km s^{−1}. The solid blue vertical line indicates the measured line depth of the signal, with a standard deviation of the value present as a transparent blueshaded region. The legend in the right panel plots labels the inin, inout and outout residual distributions calculated during the bootstrapping process. 

In the text 
Fig. C.2 Bootstrap results for the new detections. Each panel represents the bootstrap results of each new detection in the respective order of Ca I, V I, Ti I, Cr I, Ni I, Sr II, Ba II, and Tb II. The plots follow the same construction as Fig. C.3 

In the text 
Fig. C.3 Response and residual plots for each confirmed species in KELT9 b’s atmosphere. Each panel represents a different species. The order of the species is as follows (left to right): H I, Na I, Mg I, Ca II, Sc II, Ti II, Cr II, Fe I, Fe II, and Y II. Each plot contains a line of best fit plotted through the data of the response plot, spanned by two lines of worst fit. The identity line y = x is plotted to show how closely the model matches the data. The residual plot contains the residuals and the region spanning one standard deviation of the residual data. 

In the text 
Fig. C.4 Continued: Response and residual plots for each confirmed species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Se II, Ti II, Cr II, and Fe I. 

In the text 
Fig. C.5 Continued: Response and residual plots for each confirmed species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Fe II and Y II. 

In the text 
Fig. C.6 Response and residual plots for each newly detected species in KELT9 b’s atmosphere. Each panel represents a different species. The order of the species is as follows (left to right): Ca I, Ti I, V I, and Cr I. The format of the plot follows the same structure as Fig. C.3. 

In the text 
Fig. C.7 Continued: Response and residual plots for each newly detected species in KELT9 b’s atmosphere. The order of the species is as follows (left to right): Ni I, Sr II, Ba II, and Tb II. The format of the plot follows the same structure as Fig. C.3. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.