How optimisation techniques influence the credibility of the high-resolution cross-correlation spectroscopy detection of HD 143105 b

(1)

University of Amsterdam

MSc Physics and Astronomy

Track: Astronomy & Astrophysics

Master Thesis

How optimisation techniques influence the

credibility of the high-resolution

cross-correlation spectroscopy detection of

HD 143105 b

by

Dion Linssen

10814868 (UvA)

August 21, 2020 60 ECTS 09/02/19 - 08/21/20 Supervisor: dr. J.L. Birkby Examiners dr. J.L. Birkby dr. J.M.L.B. D´esert

Anton Pannekoek

Institute

(2)

Abstract

High-resolution cross-correlation spectroscopy is a relatively new technique for character-ising exoplanetary atmospheres. Using the planet’s changing Doppler shift to separate its spectrum from that of the host star, the technique has proven very useful with detections of molecules such as H2O, CO and TiO. Key to a successful detection is the cleaning of

stellar light, telluric absorption and other contaminating effects in order to get down to the planet’s thermal emission. The SYSREM algorithm in particular has been used frequently to this end. Choosing optimal parameter values for SYSREM and other steps of the pipeline are common practices in these analyses. In this thesis, the technique has been employed for two nights of observational data of the Hot Jupiter HD 143105 b, taken with iSHELL at the NASA Infrared Telescope Facility. A 7.3σ match arose with an equilibrium atmo-sphere model that included an inversion layer, due to the presence of TiO and VO, at an orbital inclination of ∼ 37◦. However, the detection strength was highly dependent on the optimisation parameters used and it was found that the method was susceptible to false detections. Most importantly, weighing different spectral orders based on the recovered strength of injected planet signals proved to produce misleading results when the Kp− vsys

parameter space was not sampled well. In addition, the significance values returned by Welch’s t-test could only be interpreted correctly when compared to a distribution of ex-pected values for randomly generated data. Considering these and other pitfalls introduced by the optimisation process, no detection ‘claim’ was made.

(3)

Introduction

In 2019, the Nobel Prize in Physics was awarded for “New perspectives on our place in the universe”, in half to Mayor and Queloz for “the discovery of an exoplanet orbiting a solar-type star”, which was published in Mayor & Queloz (1995). Three years earlier, two exoplanets orbiting a pulsar had already been found by Wolszczan & Frail (1992). In the last few decades following those early discoveries, the field of exoplanetary research has seen an immense growth, with now over 4000 confirmed exoplanets according to the NASA Exoplanet Archive.

These planets exhibit a wide range of sizes, compositions, temperatures and atmospheric structures, and many of them have no analog in our own Solar System (Seager 2013). Fundamental questions in exoplanetary science deal with the formation, evolution and habitability of those planets. The atmospheres of exoplanets are important in all these areas of research.

¨

Oberg et al. (2011) first proposed how the C/O ratio in gas giant atmospheres contains valuable clues about the formation sites of these planets. As Fig. 1.1 shows, H2O, CO2

and CO have different snowlines, which are the distances from the protostar at which these gas molecules freeze out into icy grains. The position of the (proto-)planet in the disk will thus determine the availability of C and O in gaseous and solid form, which should then manifest itself in the C/O ratio of the atmosphere. Naturally, planet migration and chemical evolution of the disk complicate this picture, and recent works such as Cridland et al. (2020) have incorporated these elements into their analyses, showing that planetary atmospheres still contain useful information about the formation and evolution of the planets themselves. In addition, it should be no surprise that atmospheres play a crucial role in determining the habitability of exoplanets. Not only that, but the atmosphere also seems to be the best indicator as to whether a planet is actually inhabited. Under the assumption that life elsewhere is more or less similar to life on Earth, i.e. it produces biosignature gasses, whichever they may be for different life forms, a molecular species in strong dis-equilibrium with the atmosphere can provide an indication of life (Seager 2013). Earth-sized, terrestrial planets in the habitable zone presumably make the best targets, but the ability to detect and characterise them with the current techniques and instrumentation is limited.

Characterising the atmospheres of massive planets on short orbits is more accessible from an observational standpoint and helps the development of observational techniques that could eventually be used for terrestrial planets. These planets are called ‘Hot Jupiters’, after their Jupiter-like masses and extremely hot temperatures as a result of small orbital radii. They make for interesting targets of study since they are predicted to be rare and their formation mechanisms are still a topic of debate. It does seem likely, however, that migration plays an important role and that these planets do not form where they reside currently (Dawson & Johnson 2018).

(5)

ex-Chapter 1. Introduction

Fig. 1.1 – The expected C/O ratios of gas and grains in a protoplanetary disk around a Sun-like star. The positions of the H2O, CO2 and CO snowlines are marked. Since a protoplanet accretes planetesimals

and its atmosphere from the grains and gas in the disk, the C/O ratio of the planet’s atmosphere should reflect the disk composition at the distance where the planet formed. More recent studies such as Cridland et al. (2020) have complicated the picture by including planet migration and complex accretion flows, but the core idea that the C/O ratio is a result of all these processes and thus contains valuable information still stands. Figure taken from ¨Oberg et al. (2011).

plains the observational bias leading to an abundance of Hot Jupiter detections despite their predicted rareness. Concerning the characterisation of their atmospheres, the transit spectroscopy technique has proven very useful. It makes use of the wavelength-dependent opacity of the molecules and atoms present in the atmosphere. At some wavelengths, molecules in the atmosphere will absorb more strongly and the transit depth will be deeper, allowing for the extraction of a transmission spectrum of the planet’s atmosphere, which can also include information about the presence of clouds or hazes (see e.g. Sing et al. 2011, 2016).

On the other end of the orbital distance spectrum, there are far-out planets and their atmospheres that can be observed directly (see e.g. Macintosh et al. 2015). These techniques favour massive planets with high thermal emission to mitigate the contaminating flux of their host stars. Lower resolution spectroscopy observations help gather more photons per wavelength bin, but Snellen et al. (2014) detected spectral signatures of CO in the atmosphere of β Pictoris b at R ∼100,000, showing that high-resolution spectroscopy also has the potential to observe widely separated exoplanets directly.

So where does all this leave the observational potential for detecting lower mass, terrestrial planets that are expected be more suitable for harbouring life? In principle, the transit spectroscopy method has the potential to characterise their atmospheres, but from a sta-tistical point of view, finding transiting Earth-sized planets, preferably close to Earth, is a difficult task. The relatively new high-resolution cross-correlation spectroscopy technique, first employed by Snellen et al. (2010), is less affected by this. It takes away the main limiting factor by not requiring the planet to transit, and thus vastly expands the pool of observable planets.

The technique aims to detect the Doppler shift of the planet’s spectrum as it orbits its host star. Whereas the star will be approximately stationary, the planet will have a radial velocity semi-amplitude which is usually on the order of ∼100 km/s. This radial velocity still corresponds to only a tiny shift in wavelength, but at high resolution (in this context usually R >30,000), it is observable. The planet’s Doppler shift thus allows the disentangling of the planetary and stellar spectra (Fig. 1.2). In addition, observing at high resolution allows otherwise broad molecular absorption or emission bands to be resolved into individual spectral lines. In doing so, it allows the unique spectral pattern of each molecular species to be detected. The individual planetary spectral lines are usually buried even beneath the

(6)

Chapter 1. Introduction

Fig. 1.2 –The core principle of the high-resolution spectroscopy technique. As the planet orbits its host star, its constantly changing radial velocity will Doppler shift its spectrum, while stellar and telluric features are stationary. The derivative of the radial velocity peaks at phases φ = 0 and φ = 0.5, which means that the Doppler shifting of the planet spectrum is best observed close to superior and inferior conjunction, for day side and night side thermal flux, respectively. Figure taken from Birkby (2018), top panel made by Ernst de Mooij.

photon noise of the star, however, and therefore never really reveal themselves, even after extensive removal of stellar and telluric contamination. Cross-correlation with a model spectrum is an excellent way of solving this problem, since it practically combines the simultaneous shifting of many very faint planetary spectral lines together, resulting in an observable signal.

It should also be said that low-resolution spectroscopy and high-resolution cross-correlation spectroscopy are highly complementary as they focus on different properties of the atmo-sphere (Brogi & Line 2019). While high-resolution cross-correlation spectroscopy does well at confirming the presence of molecular species, it is quite insensitive to the spectral line-depth due to the removal of stellar light that normally acts as a continuum calibration. Low-resolution spectroscopy on the other hand retains the continuum and therefore cur-rently does better at constraining molecular abundances, so that a combination of the two can lead to a more efficient characterisation of exoplanet atmospheres.

Although the high-resolution spectroscopy technique in principle has the potential to ob-serve terrestrial planets in our local surroundings, with the current generation of telescopes, it is difficult to reach the sensitivity to do so. Over the last decade, Hot Jupiters have typ-ically been the targets of observations, with detections of molecules such as CO, H2O and

TiO (e.g. Snellen et al. 2010; Birkby et al. 2013; Brogi et al. 2016; Nugroho et al. 2017). Searching for the presence of TiO and VO in particular have drawn the attention of the community in recent years (e.g. Hoeijmakers et al. 2015; Nugroho et al. 2017; Merritt et al. 2020). Hubeny et al. (2003) showed that these molecular species can in some circumstances cause thermal inversions in stellar or planetary atmospheres. Building on their work, Fort-ney et al. (2008) proposed a transition around ∼ 2000K above which the presence of TiO and VO in a planetary atmospheres cause a thermal inversion in its stratosphere, and used this temperature to make a distinction between ‘Hot Jupiters’ and ‘Ultra-Hot Jupiters’ which would display such an inversion. This in turn would lead to a very bright thermal day side flux, and perhaps more important for the current work, molecular spectral features in emission rather than absorption.

(7)

Chapter 1. Introduction

iSHELL data of the Hot Jupiter HD 143105 b, in an attempt to find a match with an equilibrium model spectrum, either including or excluding TiO and VO. A credible detection could therefore both distinguish between these two cases, while the value of the radial velocity semi-amplitude that a detection necessarily constrains could provide the inclination and true mass of the planet, which until now have been unknown. Confirmation of a thermal inversion due to the presence of TiO and VO would place HD 143105 b in the ‘Ultra-Hot Jupiter’ category, while it currently sits at the boundary with a day side equilibrium temperature of ∼2300K.

In addition to analysing the data of HD 143105 b, the implementation of the high-resolution cross-correlation spectroscopy technique was investigated. The data cleaning processes, or more importantly the algorithms to ‘optimise’ them to retrieve the planet signal, are known to influence these analyses (Cabot et al. 2019). Moreover, in literature, detections are typ-ically presented as a probability value, but its exact meaning and implication are often left unstudied. In this work, the results for HD 143105 b were examined in depth by investi-gating the influence of different aspects of optimisation techniques, such as weighing the contributions from different spectral orders. The resulting probability values were inter-preted more accurately by comparing them to the values obtained for randomly generated data.

In Chaps. 2 and 3, the observations, data reduction and implementation of the high-resolution spectroscopy technique are explained. Chap. 4 presents the results, where an apparent detection of the model including TiO and VO is shown. Chap. 5 discusses these results at length by focusing on the side-effects a typical analysis introduces by choosing fundamental parameters in an ‘optimal’ way, while tools that can help make future high-resolution spectroscopy analyses more robust are proposed. Chap. 6 concludes the thesis.

(8)

Chapter 2

Observations and data reduction

The target of observations was HD 143105 (V = 6.75 mag Høg et al. 2000, K = 5.52 mag Cutri et al. 2003). HD 143105 is an F5 star (Cannon & Pickering 1993) with a planetary companion. It was observed on the nights of April 28th and June 14th, 2017, for approxi-mately four hours each, using the iSHELL high-resolution spectrograph at the 3.2m NASA Infrared Telescope Facility, located at the Maunakea Observatories on Hawai‘i. iSHELL’s cross-dispersed echelle spectrograph uses a 2048x2048 H2RG array detector, spanning the 1.06 − 5.3µm range over its different observing formats (Rayner et al. 2016). In these ob-servations, the K3 mode was used, spanning the 2.26 − 2.55µm range over 27 overlapping spectral orders. The slit width of 0.375” corresponded to a resolving power of R ∼ 75, 000. During the first night, 82 observations each lasting 179.8 seconds were taken, with gaps of 20 seconds between exposures, of which 15 seconds readout. During the second night, 80 observations of the same duration were taken. The planet was observed after superior conjunction in the orbital phases 0.55 < φ < 0.63 during the first night, while it was observed before superior conjunction in the orbital phases 0.38 < φ < 0.46 during the second night. The wavelength range and orbital phases were chosen to receive the maximum thermal emission from the planet’s day side. Furthermore, the chosen phases correspond to the maximum change in radial velocity apart from superior conjunction itself, which is unavailable due to the coincidence of stellar and planetary lines. Observing both before and after superior conjunction then manifests itself as a typical ‘cross’ pattern in the final t-test matrix (Sec. 3.5) and helps to constrain the inclination and true planet mass. The planetary companion of the star, HD 143105 b, was first detected using the radial velocity method by H´ebrard et al. (2016). They found a period of P = 2.1974 ± 0.0003 days, no significant evidence of an elliptical orbit (e = 0.021 ± 0.016), a semimajor axis of a = 0.0379 ± 0.0009 AU and a mass of M sini = 1.21 ± 0.06MJ. Its close separation

from the host star and super-Jovian mass make HD 143105 b a Hot Jupiter. This type of exoplanet has a relatively high probability of transiting, which can be approximated as P ≈ R∗/a = 0.12 assuming a solar radius for HD 143105. H´ebrard et al. (2016) made

several attempts to detect a transit signal. One of their six light curves showed a transit-like feature, but the authors concluded that the planet is non-transiting based on follow-up spectroscopy showing no Rossiter-McLaughlin anomaly.

This means that the inclination of the planet’s orbit, its radial velocity semi-amplitude and true mass are unknown. In order to reasonably detect the planet’s spectral features, they are required to shift in wavelength by ∼10 pixels over the course of the observations. For a resolving power of R = 75, 000 = c/∆v, with three pixels per slit width, the radial velocity difference between pixels is ∼1300 m/s. Given the stellar mass, semi-major axis and assuming a circular orbit, the orbital velocity of the planet can be calculated to be ∼188 km/s. Using the orbital phases of the observations presented here then puts a lower limit of i ≥ 10◦ on the inclination to reach the required radial velocity shift of ∼13 km/s.

(9)

Initial data reduction

This inclination corresponds to a radial velocity semi-amplitude of Kp ≈ 33 km/s, which

is therefore the minimum value the observations are sensitive to. An inclination of 90◦, in which case the planet is transiting, means that Kp= 188 km/s, which would lead to a shift

of ∼60 pixels.

2.1 Initial data reduction

Extracting one-dimensional spectra from the raw detector images was done before the start of this project, using the publicly available data-reduction tool Spextool v5 (Cushing et al. 2004; Vacca et al. 2003). Spextool takes as input the raw detector images, flat fields, darks and the ThAr arc lamp wavelength calibration files. For each of the ∼80 observations per night, Spextool was run using the default settings for iSHELL’s K3 observing mode. The tool returned the wavelength solution, flux values, flux uncertainties and error-indicating flags for each spectral order. It was at this stage that I received the data. An example of a raw spectrum is shown in Fig. 2.1.

In order to work with a large number of spectra simultaneously in a clear way, for each night and spectral order, the ∼80 one-dimensional spectra were first combined into a two-dimensional matrix, or “spectrum time series” (Fig. 2.3 top panel). This resulted in two times 27 matrices, containing the flux values for each wavelength channel and frame. The signal-to-noise of the data varied per frame and night, shown for the 21st spectral order in Fig. 2.2. In general, the data of the second night had a higher S/N. The last two frames of the first night were discarded because of very poor quality (due to unknown reasons), and thus the spectrum time series of both nights consisted of exactly 80 frames.

Initially, the data was completely dominated by the stellar flux and telluric absorption. The spectral features of the planet’s atmosphere were buried even beneath the photon noise of the star. In order to eventually cross-correlate with atmospheric models, the data needed to be cleaned of as much of the stellar spectrum and telluric absorption features as possible. The goal was to get down to the Gaussian photon noise, when removing the star any further was not possible anymore, but the planet’s atmospheric spectrum remained relatively intact. The first steps in doing so were correcting the residual bad pixels in the data and aligning the different frames to the same wavelength grid.

(10)

Initial data reduction

Fig. 2.1– A single frame of the 5th spectral order. In further analysis, all 80 frames were combined into a two-dimensional matrix as shown in Fig. 2.3. Telluric absorption features completely obscure parts of the stellar spectrum. The overall continuum shape is not physical but rather due to a combination of the optics, resulting in a shape called the blaze function, and detector efficiency. During later data processing steps, the continuum shape was subtracted out.

Fig. 2.2 –Signal-to-noise values for the 21st spectral order, calculated as the square root of the detector electron counts. The blaze function has propagated into a generally lower S/N at the edges of the spectrum. The average S/N varied from frame to frame, probably as a result of changes in seeing and clouds. The second night had an overall higher and more constant S/N than the first, which is the likely explanation for the consistently better results for the second night in Chap. 4. The last two frames of the first night were discarded because of their low S/N.

(11)

Bad pixel corrections

2.2 Bad pixel corrections

The first step in reducing the data was correcting the residual bad pixels, which are pixels with too high values in the extracted one-dimensional spectra (such as the peak around pixel number 850 in Fig. 2.1). They can be a product of bad extraction from the raw detector images by Spextool, or have a physical origin, which is usually cosmic ray hits. Pixels can also be ‘dead’, in which case their values are too low. These dead pixels were found to be much rarer than the overly bright, ‘hot’ pixels. These processes can be spontaneous, when pixels are bad in one frame but behave properly after readout between frames. In some cases, they stay bad for the entire night, which is visible as a ‘bad column’ in the spectrum time series (see e.g. the bad column visible as a white vertical line around 1100 in Fig. 2.3 top panel).

Bad pixels are a problem when left uncorrected because they become even more pronounced after the detrending algorithms that remove stellar and telluric contamination, described in Secs. 3.2 and 3.3. This would then lead to unwanted effects in the cross-correlation stage (Sec. 3.4), since the bad pixel could match any spectral line in the model and increase the overall noise. The bad pixels were corrected by first identifying their locations and then interpolating between neighbouring pixels. The first step in doing so was normalising the spectra, because this made the bad pixels easier to spot.

The median continuum value for each frame (row) of the spectrum time series was divided out and the mean of each column was subtracted. This step highlighted most of the variation due to telluric absorption features, as well as differences in overall flux per frame. Then, a Laplacian of the Gaussian blob detection algorithm was applied (see e.g. Kong et al. 2013). This algorithm consisted of two steps. The first was to smooth the data on small scales by convolving it with a two-dimensional Gaussian filter with width parameter t (Eq. 2.1). The second step was to find sharp features by applying the Laplacian operator (Eq. 2.2).

g(x, y, t) = 1 2πte −x2+y2_2t _(2.1) ∆f (x, y) = ∂ 2_{f (x, y)} ∂x2 + ∂2f (x, y) ∂y2 (2.2)

Here, x and y are the coordinates of the two dimensions of the normalised, column-mean subtracted data (i.e. pixel number and flux). These two steps resulted in a residual matrix that was smoothed out, but where sharp features on the scale of r =√2t, were emphasized (Fig. 2.3 middle panel). To optimize this algorithm for identifying bad pixels, t = 0.5 was chosen such that it found blobs on the scale of r = 1 pixel.

Any pixels at > 5σ deviation from the median of the surrounding 10 × 10 pixels were taken to be bad pixels. In addition, any pixels next to these bad pixels at > 3σ deviation were also identified as bad. This was done because bad pixels were often found to come in small groups. Not all of these neighbouring pixels turned up at 5σ, however, and it was for this reason that a more lax threshold of 3σ was taken for any pixels next to identified bad ones. After letting the described algorithm find bad pixels, any remaining bad pixels were pin-pointed by eye. All bad pixels were then corrected in the spectrum time series by linear interpolation from horizontally neighbouring pixels. Some bad pixels (on average two small groups per order) were located in an absorption feature, in which case interpolating from neighbouring pixels in their own frame led to outliers when applying the detrending al-gorithms described in Secs. 3.2 and 3.3. Therefore, these pixels were interpolated from the corresponding pixels in the two adjacent frames, while correcting for differences in continuum level between frames.

(12)

Alignment of spectra

Fig. 2.3 –Residual bad pixel correction steps for a part of the 5th spectral order of the first night. Top: The 80 spectra combined into a spectrum time series. The colour map represents flux values. Middle: The data after normalising, subtracting the mean of each column and applying the Laplacian of Gaussian blob detection algorithm. Features on scales larger than a singular pixel are smeared out while small features, such as bad pixels, are enhanced and visible as black dots. Any 5σ outliers were identified as bad pixels. Bottom: The data after correcting the bad pixels and bad columns. Note for example that the bad column around 1100 has been corrected.

Bad columns did not show up in this process, because it involved subtracting the mean of each column during normalisation. Therefore, before normalising, bad columns were iden-tified by eye, and all pixels in those columns were interpolated together with the individual bad pixels as described above. An example of a bad pixel corrected spectrum is shown in the bottom panel of Fig. 2.3.

2.3 Alignment of spectra

The fact that the absorption bands in Fig. 2.3 are slightly tilted means that the spectra drifted in wavelength over the course of the night, typically by roughly 2 pixels between the first and the last frame. This effect can be caused by thermal drift of the instrument (Brown 1990), but it could in theory also be due to the change in barycentric velocity which causes a Doppler shift. The change in barycentric velocity was on the order of ∼160 m/s (∼0.1 pixel) for these observations and was therefore not the main drift source.

Given the precision that is required to successfully match models to the data, it was im-perative that the spectra were properly aligned. In addition, the misalignment of spectra leads to artifacts in the data when applying normalisation steps or detrending algorithms. These become visible at the edges of telluric absorption bands, where due to misalignment, a part of the column is still part of the continuum, while the other part falls in a telluric absorption line. This means that subtracting the mean of those columns results in a typical negative-positive pattern as visible in Fig. 2.4, hindering the ability to adequately clean the data to the photon noise in later steps. These problems were avoided by aligning the spectra. Throughout the project, three different approaches for doing so were tried.

(13)

Fig. 2.4 – The effect of aligning the spectra becomes clear when normalising and subtracting the mean of each column. Top left: Part of the fifth spectral order of the first night, before aligning the spectra. Bad pixels have been corrected. Bottom left: Same as top left, but with each row normalised to the same continuum value and the mean of each column subtracted. Because the spectra are not aligned yet, typical vertical black and white patterns emerge at the edges of telluric absorption features. The horizontal white features around frame number 60 are the result of different continuum shapes of the spectra. Top right: Same as top left, but now the data has been aligned by cross-correlation and a constant shift. Bottom right: Same as top right, but with each row normalised to the same continuum value and the mean of each column subtracted. Because the spectra are aligned, the typical vertical black and white patterns visible in the bottom left image are absent.

(14)

Fig. 2.5 – A zoom-in on the peak of the cross-correlation function between the 70th row and first row of the fifth spectral order. The highest cross-correlation value was found at an offset of -1 pixel, but when allowed for sub-pixel shifts, the Gaussian fit revealed the offset to be -1.23 pixels (indicated by the dotted green line). The 70th row was thus shifted by 1.23 pixels and linearly interpolated back on the same pixel grid as the first row.

2.3.1 Aligning by cross-correlation and a constant shift

This was the simplest way of aligning the data. Every individual spectrum was correlated with the first spectrum of the night. The location of the peak of the cross-correlation function indicated the offset at which the spectra matched best. However, the cross-correlation function had an increment of one pixel, which would therefore only allow finding the peak and shifting the spectra at integer pixel precision. In order to get a more stringent, sub-pixel value for the offset, a Gaussian was fitted to the peak of the cross-correlation function and the maximum of the Gaussian fit was extracted at 0.01 pixel precision (Fig. 2.5). The spectrum was shifted by the offset and linearly interpolated onto the same pixel grid as the first spectrum of the night. An example of part of a spectrum aligned by cross-correlation is shown in the right panels of Fig. 2.4.

The advantage of this alignment approach is that it is fast to calculate. The disadvantage is that it shifts each spectrum in its entirety, while the offset at one end of the spectrum is not necessarily the same as the offset at the other end. More advanced alignment approaches allow for these non-constant shifts, typically by stretching the spectra using a polynomial offset function (such as described in Secs. 2.3.2 and 2.3.3).

2.3.2 Aligning by fitting telluric line positions

Aligning the spectra by the telluric lines tracked the positions of multiple telluric absorption lines throughout the night. For each order, between five and seven features spread as evenly and as widely as possible were selected manually. Only unsaturated lines were used where possible, since the positions of their minima were well-defined. For each spectrum, the five pixels around the minimum of each line were interpolated with cubic splines to a grid of 1/20 pixel resolution. The sub-pixel position of each line was then set to the minimum of this interpolated line, similar to the wavelength calibration approach in Brogi et al. (2018). In orders with insufficient unsaturated lines present, an occasional saturated line was used, with its position set manually, roughly halfway through the saturated region.

Generally, the line positions shifted to redder wavelengths over the course of the night, but there were differences in this shift between lines at different ends of the spectrum (Fig. 2.6). For each individual spectrum, the line positions were compared to the line positions of the

(15)

Fig. 2.6–The positions of five telluric absorption lines spread out over the 19th spectral order of the first night. The colour map is for clarity only and mimics the y-axis in representing the spectrum number. The line positions generally drifted to redder wavelengths over the course of the observations, but the extent and curve of this drift varied between different lines at different ends of the spectrum.

first spectrum of the night (Fig. 2.7). A second order polynomial was fitted through the line positions. Evaluating that parabola at every wavelength bin of the first row gave the corresponding wavelength bins of the row that was to be aligned. In order to work with the data after aligning, however, it was necessary that all spectra were at the same wavelength grid. Therefore, after obtaining the new wavelength bins of the aligned row, the spectrum was interpolated back onto the wavelength grid of the first row.

The advantage of this alignment method is that it allows not only for shifting of a row, but also for stretching, which is achieved by fitting a parabola to the line positions rather than for example averaging out their offsets from the first row. The disadvantage is that it focuses solely on matching the positions of a few lines. The improvement of the alignment on the rest of the row is not assessed via e.g. cross-correlation with other rows, as was done in alignment by a constant shift (Sec. 2.3.1). This means that the fitted lines may be very well aligned, but other parts of the spectrum less so.

2.3.3 Aligning by cross-correlating a grid of wavelength solutions

This alignment method in principle overcame the limitations of the previous two alignment methods: the spectra were allowed to stretch, and were also cross-correlated with the first row to find the best alignment for the spectrum as a whole. The principal steps of the wavelength calibration of Brogi et al. (2016, Sec. 3.1) were followed, although now applied to aligning the data rather than finding a wavelength solution.

An array of ten possible shifts between −2 and +2 pixels in steps of 0.44 pixels was made for the first, middle and last pixel of each row. Together, this resulted in a grid of 1000 trial combinations of shifts per row. A parabola was fitted through the shifted positions of the first, middle and last pixel for each combination. The parabolic fit was used to stretch

(16)

Wavelength solution

Fig. 2.7 – The five selected line positions in the 60th row compared to their positions in the first row, for the 19th spectral order of the first night. The red dotted line indicates a perfectly aligned row. For visual purposes, the offsets from this dotted line have been magnified by a factor 500. A parabola was fitted through the line positions and returned the new wavelengths for the 60th row. After putting the row on the new wavelength grid, it was interpolated back onto the wavelength bins of the first row to get all the rows on the same wavelength grid.

the row by mapping every old integer pixel number to a new pixel value that would no longer necessarily be an integer. All 1000 trials were cross-correlated with the first row, and the best match was used as the resulting aligned row. This process was iterated three additional times, using five instead of ten possible shifts while decreasing the pixel step size every iteration, to eventually align at 0.01 pixel precision. The final shifted and stretched row was eventually interpolated back onto the original integer pixel numbers to place it in the same pixel space as the first row.

The advantage of this alignment approach is that it allows the rows to stretch, while cross-correlation assesses the improvement on the full row (which was not done in Sec. 2.3.2). The disadvantages are that it is slow to calculate, and caution for ‘over-fitting’ is needed. As will be described in more detail in Chap. 5, the three different alignment algorithms were benchmarked against each other by running the full pipeline and comparing the results (Fig. 5.1). Judging which performed best in principle was difficult, given that the inclination of the planet’s orbit and thus where it appeared in the t-test matrix (see Sec. 3.5) were unknown. However, the three plots all showed only one significant feature at the expected system velocity and a radial velocity semi-amplitude of ∼110 km/s. It was found that depending on the detrending method used (see Secs. 3.2 & 3.3), aligning by either a constant shift or a grid of wavelength solutions performed best for this potential planet signal.

2.4 Wavelength solution

Even though the initial data reduction with Spextool returned a wavelength solution for each order (Sec. 2.1), these solutions turned out to be inaccurate. This became apparent at parts of the spectra where neighbouring spectral orders overlapped, but the same spectral features did not coincide (Fig. 2.8). By matching telluric absorption features to a synthetic telluric spectrum, in a fashion similar to the line-fitting alignment as described in Sec. 2.3.2, the wavelengths were recalibrated. The synthetic telluric spectrum was made using

(17)

Wavelength solution

Fig. 2.8–The overlapping region of the 14th and 15th spectral orders, alongside a telluric model spectrum. Not only do the observed spectra show an overall offset from the telluric model, the blue end of the two observed spectra matches up well, while the red end shows significant discrepancies, indicative of a non-constant offset from the true wavelength solution. The new wavelength solutions were produced by fitting the positions of a few of these absorption features and matching them to the telluric model spectrum.

the web-based input form for the software tool ATRAN1. It uses the telescope location and angle of observations, as well as water-vapour and ozone models to calculate atmospheric transmission spectra in the near- and mid-infrared (Lord 1992).

The same absorption features that were used for alignment by telluric line-fitting were used here. For each order, an averaged spectrum of all 80 spectra was made and matched to the synthetic telluric spectrum. Identical to Sec. 2.3.2, the positions of the selected telluric absorption features in both the data and telluric model spectrum were calculated to 1/20 pixel precision and compared to each other (similar to Fig. 2.7). A third order polynomial was fitted through the line positions and evaluated at every ‘old’ wavelength bin to find the newly calibrated wavelength solution. This concluded the first data reduction steps and the residual bad pixel corrected, aligned and wavelength calibrated spectrum time series were used as input to the pipeline that will be described in Chap. 3.

1

(18)

Chapter 3

Methods

HD 143105 completely outshines its planetary companion. This means that at first, all features in the spectra are either from the star and Earth’s atmosphere, or a result of weather conditions and instrument behaviour. The atmospheric spectral features of HD 143105 b are present in the data, but they are so faint that they are buried even beneath the photon noise of the star. In order to ‘probe’ these spectral lines, the data needs to be cleaned and processed rather intensively in a process called ‘detrending’. Even then, the planetary lines only reveal themselves when combined all together using cross-correlation. There are two common ways of detrending high-resolution spectroscopic data. When the technique was only just emerging, analyses usually fitted the variation in the vertical, time direction with a functional form of the airmass. Airmass is defined as the column density of Earth’s atmosphere along the line of sight of observations, normalised to 1 towards the zenith (Green & Center 2006). As the observing direction changes over the course of the observations, the airmass changes and causes variation of flux in the time-direction. The very first high-resolution spectroscopy detection of CO absorption signatures in the atmosphere of HD 209458 b by Snellen et al. (2010) was achieved using airmass detrending. In recent years, the community has started to rely more heavily on SYSREM, which is a Principle Component Analysis algorithm that allows for unequal error bars on each data point (Tamuz et al. 2005). It was first used in high-resolution spectroscopy context by Birkby et al. (2013), who discovered water in the atmosphere of HD 189733 b. Cabot et al. (2019) compared the performance of airmass and SYSREM detrending for a few detections in literature and found that SYSREM generally performed better in retrieving planetary signal on data with high telluric absorption, whereas airmass detrending performed slightly better on data less affected by telluric absorption.

In this research, both airmass detrending and SYSREM have been used. Airmass detrending requires the use of less hyperparameters, and was therefore used first to explore the data. The more optimizable SYSREM was then used to detrend the data better and produce more constraining results.

3.1 Masking telluric absorption bands

Regardless of the detrending method used, telluric absorption features posed problems in analysing the data. Fig. 2.1 shows that absorption by the Earth’s atmosphere caused the flux to approach zero, which means that these parts of the spectrum contained minimal spectral information on the host star and planet. The techniques described in this chapter detrended the spectrum time series such that the columns of telluric absorption features ended up at comparable levels to the continuum. The result was that these columns added nothing but noise to the detrended spectrum time series, hindering the effectiveness of the cross-correlation stage. To prevent this, the telluric absorption features were masked by

(19)

Masking telluric absorption bands

Fig. 3.1–Top: The column-means of the the tenth order of the first night. A few points of the continuum were picked by eye, indicated as red points. A Gaussian was fitted through these points. Bottom: The normalised spectrum after dividing the top panel by the Gaussian fit.

setting their flux to zero, which excluded them from influencing the cross-correlation. Picking the columns to be masked was done with the following algorithm. For every order, the mean was taken along the columns, reducing the spectrum time series to one dimension again. A few points in the continuum were selected manually, and a Gaussian was fitted through them, since it approximated the continuum shape well (Fig. 3.1 top panel). The column-mean spectrum was subsequently divided by this fit, resulting in one continuum normalised spectrum per order (Fig. 3.1 bottom panel).

All 27 continuum normalised spectra were then stitched together, resulting in a single spectrum spanning the full 2.26 − 2.55µm range. The masked columns were chosen based on two parameters, in line with the masking parameters proposed by Cabot et al. (2019). A fraction pm = 0.15 of all columns with the lowest mean were masked. These were the

telluric absorption features, where the flux was low.

In addition to masking columns with low mean, they were also masked based on their variance. The spectrum time series of each order was continuum normalised and column-mean subtracted in the same manner as explained in Sec. 2.2 and shown in the bottom panels of Fig. 2.4. ‘Continuum normalised’ in this sense therefore only refers to setting the highest continuum flux of each row to one, rather than fitting the shape of the continuum and dividing it out, like explained above for low flux masking. The variance was calculated for each column of this normalised spectrum time series. This led to a (one-dimensional) array of variances for each order. Again, all 27 variance arrays were stitched together, and the median and standard deviation of this long array were calculated. Any columns with variance pm = 2 standard deviations higher than the median variance were masked.

These were noisy columns, often located at the edges between the continuum and telluric absorption features.

The standard masking values used in this project were pm = 0.15 and pv = 2, meaning that

the 15% of columns with lowest mean were masked, as well as any columns with variances 2σ above the median of all variances. These parameters were chosen such that the telluric absorption bands were completely covered. Spurious unmasked columns in otherwise large masked regions still remained, and to prevent this, gaps in the mask smaller than 5 pixels were closed. The values of all pixels in the masked columns were set to zero.

(20)

Airmass detrending

Fig. 3.2 – A two-dimensional spectrum with a planet atmosphere model spectrum injected at 3000 times the expected strength. The injected model included the molecular species TiO and VO, causing the spectral lines to appear in emission (Fortney et al. 2008). As time progressed, the planet’s radial velocity changed as it moved along its orbit, and the lines are seen to shift in wavelength. It is this effect that ensures the planet signal remains largely unaffected while trends in the horizontal and vertical directions are removed.

3.2 Airmass detrending

The products of the data reduction steps of Chap. 2 were well-aligned and masked spectra time series for each night and spectral order. The planet’s spectral lines were shifting in wavelength over the course of the observations, which means they did not appear at the same pixel number in all rows, unlike the telluric bands and stellar lines (Fig. 3.2). Detrending the data by fitting the correlation with airmass removed trends in the vertical direction, while the division by the fit normalised each column, thereby also removing variation in the horizontal direction. The planet’s spectrum shifted during the time series and was thus left intact by the detrending steps.

Dark rows in Fig. 3.2 show that the overall flux level varied greatly from frame to frame. This trend was removed first. The continuum of each frame was normalised by fitting a Gaussian through the same continuum points as those used in the masking process (Sec. 3.1) and dividing the row by the fit. The result is shown in the top panel of Fig. 3.3. The division by the Gaussian fit normalised the overall flux level of each row, which already removed part of the variation in flux due to changing airmass.

The next step was to remove the remaining correlation between airmass and flux over time. The airmasses can be considered an array of length 80 for each night, with one value of the airmass for each integration of the observations. Using linear regression, a line was fit through the values of each column as a function of their corresponding airmass value (Fig. 3.4). The column was divided by this line, effectively normalising it and removing any correlations with airmass (middle panel of Fig. 3.3).

At this stage, close to telluric absorption bands, residual patterns still persisted (e.g. those marked by red arrows in Fig. 3.3). These residuals tended to have approximately the same variation through time across a single order, which made it possible to remove them. For a given order, a few columns with prominent residuals were selected and their average was used as independent variable in linear regression for every column, thus repeating the previous step but now for the strongest residual rather than airmass. The result was a detrended spectrum time series as shown in the bottom panel of Fig. 3.3.

(21)

Airmass detrending

Fig. 3.3 – The steps taken in detrending data by fitting the airmass. The spectrum time series shown is the same as that of Fig. 3.2, but without the injected planet model. Masked columns are shown in black. Top: The spectrum time series has been normalised by dividing a Gaussian fit to the continuum for each row. Middle: Every column has been plotted against the airmasses and a linear regression fit has been divided out (see Fig. 3.4). Red arrows indicate columns with strong residuals, used to detrend the data once more. Bottom: Every column of the spectrum time series has been plotted against the average residual column and a linear regression fit has been divided out, producing the final detrended spectrum time series.

Fig. 3.4 – Left: The values of the airmass for all observations. The correlation between every individual column of every normalised spectrum time series (Fig. 3.3 top panel) and the corresponding night’s airmass trend shown here was fitted. Right: Linear regression with airmass as independent variable and the values of a single column as dependent variable. The column values were subsequently divided by the fit, resulting in the middle panel of Fig. 3.3 after iterating over all columns. Repetition of this process with the strongest residual as independent variable resulted in the final detrended spectrum time series (Fig. 3.3 bottom panel).

(22)

SYSREM detrending

3.3 SYSREM detrending

Detrending data by removing correlations with airmass has a clear physical motivation. However, airmass is not the only component of variation in the data. Other effects such as seeing can cause variation as well, and there could even be processes for which there is no functional form available. Ideally, these other components are removed from the data while considering their variation patterns, as was done for airmass in Sec. 3.2. When this is not possible, for example if there is no known functional form, detrending algorithms that do not require knowledge about the sources of variation, such as SYSREM, provide a useful solution.

At the core of this detrending technique lies the SYSREM algorithm, as described in Tamuz et al. (2005). The algorithm is similar to Principal Component Analysis, more specifically Singular Value Decomposition (SVD, see e.g. Golub & Kahan 1965), but allows for unequal errors on the data. Originally designed with the purpose of removing systematic effects in a set of photometric lightcurves, the algorithm is commonly used in high-resolution cross-correlation spectroscopy as well (e.g. in Birkby et al. 2013; Nugroho et al. 2017; Merritt et al. 2020), owing to the fact that the nature of the removed trends does not need to be known.

Before starting, the median continuum level of each row was divided out and the mean of each column was subtracted, which removed some of the variation that would likely be the first trend removed by SYSREM otherwise. SYSREM was then applied up to ten iterations, the principal steps of which will be outlined here. Running more than ten iterations would risk removal of the planet signal, and the major physical components of variation should be removed after ten iterations. For a more detailed walkthrough of the algorithm, see Tamuz et al. 2005.

Consider the input spectrum time series to be a m × n matrix with entries rij and

uncer-tainties σij. The aim is to find the orthogonal vectors {aj; j = 1, ..., m} and {ci; i = 1, ..., n}

that minimise the following expression: S_i2 =X

j

(rij − ciaj)2

σ_ij2 (3.1)

That is, the vectors ci and aj represent the biggest variance components of the matrix in

the x and y-direction, while respecting the errorbars. When some initial values for the components of aj are assumed, differentiation and equating to zero gives a solution for ci:

ci = P j rijaj σ_ij2 P j a2 j σ2 ij (3.2)

These values can now be used to update the assumed components of aj by minimising an

expression similar to Eq. 3.1: S_j2 =X i (rij − ciaj)2 σ2 ij (3.3) which comes down to calculating

a(1)_j = P i rijci σ2 ij P i c2 i σ2 ij (3.4)

and updating the initially assumed aj. This process can be repeated to find c(1)_i , then

(23)

Cross-correlating with a planet atmosphere model

Fig. 3.5–Detrending of a spectrum using SYSREM. For clarity, the values of pixels in masked columns have been set to zero, which in reality were masked by setting their uncertainties sufficiently high. Top: The output of the initial data reduction steps. Second: The spectrum time series after normalisation, used as input to SYSREM. Each row has been divided by its continuum level and the mean of each column has been subtracted. Third: The residual spectrum time series after one iteration of SYSREM. Fourth: The residual spectrum time series after a second iteration of SYSREM. It starts to look like Gaussian noise, which is the desired end product because it means that the stellar light has been cleaned. Bottom: The residual spectrum time series after two iterations of SYSREM, with a very strong injected model. This shows the spectral lines have not been removed due to their shifting in wavelength over the course of the observations.

optimal values have been found. Tamuz et al. (2005) have shown that this process always converges on the same vectors, irrespective of the initially assumed aj.

The implementation of this algorithm was as follows. After normalisation, the spectrum time series and its uncertainties were fed into SYSREM, while the airmass values served as initial aj. They were usually somewhat close to the optimal a∗j since airmass was expected

to be one of the biggest variance components, and therefore helped speed up the calculation time. After the optimal a∗_j and c∗_i had been found, their cross-product was an m × n matrix that approximated the spectrum time series, and was thus subtracted from the data. The residual spectrum time series was saved, and then fed back into SYSREM for the second iteration, using the same uncertainties on the data (which did not change upon subtraction of a constant value) and the same airmass values as initial aj. Higher iterations tended to

take longer to compute because the initial values for aj were now less suited, since

airmass-related trends were usually removed in the first few iterations. The end result was ten detrended spectra time series, one for each number of SYSREM iterations (Fig. 3.5). During this whole process, masked columns did not have their values set to zero but instead have their uncertainties set to 1,000,000, which was so high compared to the normalised flux values around 1 that these columns practically did not influence the detrending process. The advantage of this was that it had to be done only at the start, since the uncertainties were constant throughout, while setting the flux values to zero would require repeating this between each iteration.

3.4 Cross-correlating with a planet atmosphere model

Both airmass and SYSREM detrending aimed to remove as much of the stellar and telluric contamination as possible, while leaving the planet’s signal intact. Only then would it

(24)

Cross-correlating with a planet atmosphere model

Fig. 3.6 – Top: A cross-correlation matrix containing a single cross-correlation value for each radial velocity shift and frame. The matrix resembles noise since most radial velocity shifts did not match the data and resulted in random values. Only one RV shift (and its neighbouring RV bins to a lesser extent) matched the real radial velocity of the planet for each frame, with its value dependent on the orbital phase of the planet and the bulk radial velocity of the host system. Bottom: Same as top panel, but with an artificially injected model spectrum at 30x the expected strength, showcasing the typical white ‘trail’ in the cross-correlation matrix indicative of a detection. In non-injected data, identifying a much fainter trail by eye can be hard, so Welch’s T-test was used to assess the statistical significance of the trails.

become useful to cross-correlate with a model spectrum of the planet’s atmosphere, as the model would match on stellar or telluric features less frequently.

Before cross-correlation with a model spectrum, noisy columns in the detrended spectrum were down-weighed by dividing each column by its variance. This was done to prevent the worst columns (those that had not been detrended as effectively and thus had high variance) from dominating the cross-correlation.

After this initial step, each row of the spectrum time series was cross-correlated with the model spectrum. This cross-correlation was performed for a grid of radial-velocity shifts ranging from -330 km/s to +330 km/s, in steps of 1250 m/s, matching the spectral resolution of the iSHELL instrument. Each frame was thus cross-correlated 528 times, each time with a different RV shift for the model. These shifts were performed by multiplying the wavelength values of the model spectrum by the factor

s

(1 + RV/c)

(1 − RV/c) (3.5)

After applying the shift, the model was binned to the spectral resolution of the data and then cross-correlated. This process resulted in a ‘cross-correlation matrix’ (CCF matrix), which contained the cross-correlation value for each row and radial velocity shift (Fig. 3.6). Such a CCF matrix was calculated for each detrended spectrum time series. This means that each spectral order resulted in ten CCF matrices for SYSREM detrended data. The output of airmass detrending was a single detrended spectrum time series and thus resulted in only one CCF matrix.

(25)

Shifting into the planet’s rest-frame and assessing cross-correlation significances

Fig. 3.7 – The same cross-correlation matrix as the bottom panel of Fig. 3.6, but now shifted into the rest-frame of a planet with Kp= 110 km/s, causing the white trail corresponding to those radial velocities

to appear in the center. Shifting the rows of the CCF matrix into this rest-frame has led to a section at the left where cross-correlation values have not been calculated. The strong model injection also makes visible ‘aliasing’ of the signal, where the model matched itself not only at the right radial velocities, but at constant offsets as well. This can be seen by the light vertical trails at relative systemic velocities away from 0, for example around -90 km/s. Anti-correlation aliasing occured as well, visible as dark vertical trails around +300 km/s.

3.5 Shifting into the planet’s rest-frame and assessing

cross-correlation significances

The white trail in the bottom panel of Fig. 3.6 shows high cross-correlation values at a changing radial velocity through time, equal to that of the planet, and thus indicates a detection of the injected signal. It is useful to quantify the strength of such a trail, especially when analysing non-injected data such as the top panel of Fig. 3.6, where a trail can be quite hard to see by eye. The position of the trail in the CCF matrix is determined by the planet’s radial velocity, which depends on the radial velocity semi-amplitude Kp, the orbital

phase φ(t), the radial velocity of the host system vsys, and the barycentric correction for

Earth vbary(t):

vp(t) = Kpsin(2πφ(t)) + vsys+ vbary(t) (3.6)

in which the radial velocity semi-amplitude for a circular orbit can be approximated by

Kp ≈

r GM∗

a sin(i) (3.7)

with G the gravitational constant, M∗ the host star mass, a the semi-major axis or radius

of the orbit, and i the inclination of the system. The orbital phase is given by φ(t) = t − T0

P (3.8)

where T0 is the time of inferior conjunction and P the orbital period. Given a value for Kp,

the radial velocity of the planet vp(t) can be calculated, assuming that T0 and P are known

from e.g. RV measurements of the host star or transit observations. Tracing out this radial velocity value for each frame in the CCF matrix allows shifting it into the rest-frame of the planet for this given value of Kp (Fig 3.7). When doing so, the cross-correlation value with

the radial velocity equal to vp(t) is positioned exactly at zero for all t. A white trail in the

center of the rest-frame CCF matrix thus indicates a match between model and data. After shifting the CCF matrix into some Kprest-frame, it was useful to quantify the strength

of a trail or its absence. In literature, multiple approaches to doing this have been proposed, and it is a topic of debate and active research. Some convert their cross-correlation values to log-likelihoods and derive statistical significance from there (see e.g. Brogi & Line 2019;

(26)

Buzard et al. 2020). Others simply sum the cross-correlation values along the frame-axis and divide by the standard deviation to get a value for the significance (see e.g. Birkby et al. 2013; de Kok et al. 2013). In this work, Welch’s t-test was used to tackle this problem, previously employed by e.g. Birkby et al. (2013, 2017); Nugroho et al. (2017).

Welch’s t-test assesses the probability that two normally distributed samples have the same mean, where in contrast to the well-known Student’s t-test, it allows the samples to have unequal variances (Welch 1947). The t-value is calculated as

t = qµ1− µ2 σ2 1 N1 + σ2 2 N2 (3.9)

with µ, σ and N the sample mean, standard deviation and size. If both samples have the same mean, the distribution of t-values approximates a chi-squared distribution with the degrees-of-freedom given by the Welch–Satterthwaite equation (Satterthwaite 1946):

ν ≈ _σ2 1 N1 + σ2 2 N2 2 σ4 1 N2 1(N1−1)+ σ4 2 N2 2(N2−1) (3.10)

This allows converting the t-value to the probability that the two samples have the same mean, which can in turn be converted to the commonly used significance values in terms of the standard deviation σ of a normal distribution (e.g. 2σ ⇔ p = 0.05).

After shifting a CCF matrix into the rest-frame of some value of Kp, Welch’s t-test was used

to compare one sample containing the values of a three pixel wide column to another sample containing the values outside of this trail (Fig. 3.8). A three pixel wide trail allowed the in-trail sample to be larger than only the center column, while the radial velocity offset of 1250 m/s to both sides this introduced was low enough to still match the model to the data, as can be seen by the width of the trail in Fig. 3.7. If the in-trail cross-correlation values were a lot higher than out-of-trail, the t-test returned a low probability that they were drawn from the same distribution, which was interpreted as a high significance ‘detection’. Reversely, if the in-trail values were lower than out-of-trail, the t-test also gave a low probability, but the significance values were made negative in this case, indicating anti-correlation between the data and model.

One could imagine calculating the significance in this manner only for the true value of Kp

while centering the three-pixel column around the true vsys, quoting it as a single value.

However, HD 143105 b is non-transiting, which means that its inclination and thus radial velocity semi-amplitude is unknown. The only option therefore was to try a grid of values for Kp to see which one matched the data best. This was done not only for the known vsys,

but also for a grid of vsys-values that were known to be wrong, which is common practice

in high-resolution cross-correlation spectroscopy. This part of the parameter space acted as a sanity check to verify that the best match between data and model appeared at relative vsys = 0, and revealed the structure and distribution of the noise.

In practice, this was realised by shifting the original CCF matrix into a different rest-frame for every value of Kp and sliding the three pixel column along each rest-frame matrix,

checking different values of vsys. The result was a two-dimensional matrix, hereafter referred

to as “t-test matrix”, containing the significance value for each combination of Kp and vsys

(Fig. 3.9). When using SYSREM-detrended data, a CCF matrix and t-test matrix were made for all ten detrended spectra per order. Airmass detrending produced only one detrended spectrum, CCF matrix and t-test matrix per order.

All orders were finally combined into one t-test matrix. This was done by summing the CCF matrices of the 27 orders for each night separately and then concatenating the resulting two matrices along the frame-axis into a single matrix with 160 rows. Even though each order

(27)

Fig. 3.8 – The in- and out-of-trail distributions of cross-correlation values belonging to the peak value of Fig. 4.3. The out-of-trail values are distributed like a Gaussian centered on 0, while the in-trail values show a higher mean, indicative of a match between data and model. Welch’s t-test rejects the two populations having the same mean with 7.3σ confidence.

spanned a different wavelength region, the CCF matrix radial velocity space was the same for each order, allowing addition. The final CCF matrix containing all nights and orders was then processed into a single t-test matrix.

(28)

Injecting models and weighing orders

Fig. 3.9–‘T-test matrix’ showing significance values from Welch’s t-test for the planet HD 189733 b, with a cross-correlated model of H2O and CO. Negative values indicate anti-correlation between data and model.

Birkby et al. (2013) reported a detection for this planet with this exact model, and it was reproduced in this work to verify that the pipeline worked as expected. Comparison with Fig. 3 of that paper shows a very similar structure, but the exact significance values differ, due to Birkby et al. (2013) assessing the significance by different means than Welch’s t-test in their Fig. 3. They afterwards used Welch’s t-test for the known planet parameters (indicated by the black crosshairs) and found a 4.9σ signal. The reason for the difference with the value here is unknown, but could be due to interpolation differences.

3.6 Injecting models and weighing orders

Adding the CCF matrix of every order as explained above naturally introduced the problem of choosing which of the ten SYSREM-detrended spectra to use. At the same time, some orders had higher data quality, less telluric absorption, or more expected planetary spectral lines than others. These factors made it favourable to weigh the contribution of every order to the total. Weighing was done by multiplying a single-order CCF matrix by a constant before adding it to the final, summed CCF matrix.

In order to determine the weight and number of SYSREM detrends to use for every order, model-injected data was used. The model that was to be cross-correlated was injected into the data at the start of the pipeline (after wavelength calibration and before detrending). Since the flux values of the data were in arbitrary units, it was necessary to scale the model to an appropriate level before injecting it in the data. Equation 12 of Brogi & Line (2019) was used to this end, while introducing an additional scaling factor to allow changing the injection strength: Finj(λ) = F (λ) 1 + (a − 1) Fp(λ) B(λ, T∗) R_p R∗ 2 (3.11) Here, Finj are the injected flux values, F are the original flux values of the data, Fp are

the planet model flux values, B(T∗) are flux values for a blackbody with temperature equal

to the surface temperature of the host star, Rp is the planet radius, R∗ is the host star

radius, and a is a scaling factor. The planet model flux values Fp were Doppler shifted with

the appropriate radial velocity for each row of the data. Note that this therefore required choosing values for Kp and vsys to use for injection. The injection effectively assumed that

the original flux was dominated by the stellar flux and that it could be approximated by a blackbody spectrum. The planet model values were scaled to this flux, multiplied by a scaling factor, and added to the data. The definition of the scaling factor reveals the assumption that the planet spectrum is already in the data, since choosing a = 1 results in no injection, while setting a = 2 results in an addition of one time the expected planet flux

(29)

Weighing orders based on model properties

to the data. Running the pipeline with an injected dataset ensured that the cross-correlated model was present in the data. This made it possible to examine how well the planet signal could be recovered in every order, and which number of SYSREM detrends was optimal for doing so.

The general steps were as follows. First, every order was injected at scaling a = 2, for some values of Kp and vsys, referred to as the “optimised location”. The pipeline then

produced ten t-test matrices for every order (one for each number of SYSREM detrends). Since the values of Kp and vsys were known from injecting, the significance value at that

location in the t-test matrix could be read out for all ten t-test matrices and the number of detrends that gave the highest significance was taken as the optimal number. However, if this significance value was below 5σ, it did not constitute a detection of the planet signal and thus the optimal detrends found were in fact not necessarily optimal for detecting the injected planet signal. In this case, the process was repeated for a higher scaling factor a. Only when the significance of the optimal detrends exceded 5σ, the optimal number of detrends was stored.

It was found that often the optimal number of detrends was only one or two, but that the corresponding detrended spectrum time series had remaining large-scale structure due to insufficient cleaning. This remaining structure propagated into patches of high and low values in the CCF matrix, which then propagated into an unreliable t-test matrix with large areas of high and low significances, which could cover the location of the injected planet signal by chance and exceed the significance values of higher, more reliable, numbers of detrends. The algorithm would in this case choose the optimal number of detrends to be only one or two, but the corresponding CCF matrices would be unreliable. To prevent this from happening, a minimum of three SYSREM-detrends was set.

Once the optimal number of detrends for every order had been found, the weights were found by comparing the significance value of every order at the lowest scaling a at which at least one of the orders reached 5σ. The weight for the order with the highest significance was set to 1, the weights of the other orders were scaled proportionally to their significance values, and weights for orders with negative significance were set to 0. This approach ensured that orders that did better at recovering the planet signal weighed more heavily than orders that failed at recovering the planet signal well.

These weights and optimal number of detrends were then used for the non-injected data. Out of the ten cross-correlation matrices per order, the one corresponding to the optimal number of detrends was chosen. The matrix was then multiplied by the weight and added to the total CCF matrix. This summed and weighed CCF matrix was used to make a t-test matrix, as described in Sec. 3.5. An important aspect of this approach is that it ‘optimised’ the t-test matrix to the location in Kp− vsys-space that was chosen for the injected model.

Even though the weights and detrends were in principle used on the non-injected data, they were chosen such that the pipeline performed best at retrieving a planet signal at a particular location. The discussion section will go into detail about the effect this had on the overall data analysis.

3.7 Weighing orders based on model properties

When combining orders into a single t-test matrix for airmass detrended data, the weights were determined differently than for SYSREM detrended data as described in Sec. 3.6. Instead of injecting the data and letting the recovered significances guide the weights, the weights were determined based on the properties of the cross-correlated model. Assuming the cross-correlated model is a good match, the depth and the number of spectral lines in the wavelength range of each order determine the potential for a significant detection. The weights were found by encapsulating exactly these two properties of the model in a single

(30)

Weighing orders based on model properties

Fig. 3.10 – Assessing the number and depth of model spectral lines in order to create weights for the different spectral orders (whose boundaries are not shown here to avoid cluttering). Top: Self-consistent atmosphere model spectrum for HD 143105 b, including TiO and VO, causing the spectral lines to be in emission. The continuum fit is a linear interpolation between the minima of bins of 100 data points. Middle: After subtracting the continuum fit. The black dotted line indicates the manually set threshold for a feature to qualify as an emission line. For every order, the height of its emission lines above this threshold were summed, in this way quantifying their combined depth. Bottom: The weight for each order, calculated as a normalisation of the summed spectral line depth multiplied with a measure for telluric transmittance.

number. This does mean that the weights were no longer dependent on the data (as was the case in Sec. 3.6) but rather on the model solely, and thus changed when switching to a different model.

The model was first Doppler shifted with an approximation of the radial velocity of the planet. For both nights separately, the orbital phase halfway through the observations was used in combination with an average value of Kp = 150 km/s, resulting in the averages

of vrad = −64 km/s for night 1 and vrad = 88 km/s for night 2. Then, the continuum

shape of the model was subtracted. For the model without TiO and VO (resulting in absorption lines), the continuum was taken to be a third order polynomial fit through the maxima of bins of 100 data points. The model including TiO and VO (resulting in emission lines) showed a ‘bump’ in the continuum level around 2.48µm, making fitting a polynomial unpractical. For this model, the continuum was taken to be a linear interpolation between the minima of bins of 100 data points (Fig. 3.10).

After subtracting the continuum, a threshold flux value for spectral lines was chosen manu-ally. A peak-finding algorithm was used to find every spectral line in the model surpassing this threshold value. Then, for every order, the depth of the qualified spectral lines within

How optimisation techniques influence the credibility of the high-resolution cross-correlation spectroscopy detection of HD 143105 b

University of Amsterdam

MSc Physics and Astronomy

Track: Astronomy & Astrophysics

Master Thesis

How optimisation techniques influence the

credibility of the high-resolution

cross-correlation spectroscopy detection of

HD 143105 b

by

Dion Linssen

10814868 (UvA)

Anton Pannekoek

Institute

Abstract

Contents

Chapter 1

Introduction

Chapter 2

Observations and data reduction

2.1

Initial data reduction

2.2

Bad pixel corrections

2.3

Alignment of spectra

2.4

Wavelength solution

Chapter 3

Methods

3.1

Masking telluric absorption bands

3.2

Airmass detrending

3.3

SYSREM detrending

3.4

Cross-correlating with a planet atmosphere model

3.5

Shifting into the planet’s rest-frame and assessing

cross-correlation significances

3.6

Injecting models and weighing orders

3.7

Weighing orders based on model properties