Keplerian-Stacker: A new technique for coronagraphic images to increase the signal to noise ratio and detect exoplanets with direct imagning instruments

KTH Royal Institute of Technology in Stockholm

Department of Engineering Physics

Degree Project in Engineering Physics track Subatomic and Astrophysics

Keplerian-Stacker : A new technique for coronagraphic images to increase the signal

to noise ratio and detect exoplanets with direct imaging instruments

Dimitri Estevez

Under the direction of Herv´ e Le Coroller

TRITA-FYS 2017:05

ISSN 0280-316X

ISRN KTH/FYS/–17:05-SE

Contents

Abstract 1

Introduction 2

1 Framework of direct imaging detection 4

1.1 Coronagraphy . . . . 4

1.1.1 Lyot coronagraph . . . . 4

1.1.2 Four Quadrant Phase Mask (FQPM) . . . . 5

1.1.3 Vortex Phase Mask . . . . 6

1.1.4 Adaptive Optics . . . . 7

1.1.5 Spectral Differential Imaging (SDI) . . . . 8

1.1.6 Angular Differential Imaging (ADI) . . . . 8

1.2 Telescopes . . . . 9

1.2.1 SPHERE at VLT . . . . 9

1.2.2 NaCo at VLT . . . . 10

2 Keplerian-Stacker (K-Stacker) 11 2.1 The K-Stacker principle . . . . 11

2.1.1 The orbital model . . . . 11

2.1.2 Recombining the images : an optimization problem of the Signal to Noise Ratio . . . . 13

2.1.2.1 First part: The brute force algorithm . . . . 13

2.1.2.2 Second part: The gradient descent optimization . . . . 14

2.1.2.3 Third part: Determining the detection of a companion . . 14

3 Adaptation of K-Stacker to work on real observations 16 3.1 Determining the grid for the free parameters . . . . 17

3.2 Number of reoptimized orbits . . . . 18

3.3 Modification of the calculation of background and noise profiles . . . . 19

3.4 Detection near the edge of the inner mask . . . . 20

3.5 True North default . . . . 21

3.6 Tolerance on the mass of the star . . . . 24

3.7 Accuracy on the distance of the stars . . . . 25

3.8 Influence of several companions in the images . . . . 26

4 Results on real astronomical observations 28 4.1 Images of GJ3331A with NaCo . . . . 28

4.1.1 Search for a companion with K-Stacker . . . . 29

4.2 Images of HD95086 with SPHERE/IRDIS . . . . 30

4.2.1 Search for a companion with K-Stacker . . . . 30

4.2.2 Injection of false planets in the images . . . . 31

Discussion 34 Conclusion 35 Appendix A Reduction of the inclination range 36 Appendix B K-Stacker user guide 37 B.1 Shell files (.sh) and Python files (.py) . . . . 37

B.1.1 Path of the files . . . . 37

B.1.2 Shell files . . . . 37

B.1.3 Python files . . . . 37

B.1.3.1 K-Stacker program . . . . 37

B.1.3.2 Useful files for K-Stacker . . . . 38

Appendix C Readme text file 39

Bibliography 42

Abstract

Direct imaging detection of exoplanets has been one of the main techniques to be devel- oped and improved for the last decade using high contrast imaging facilities such as Sub- aru/HiCIAO, VLT/SPHERE or Gemini/GPI. Combined with a data reduction method like Angular Differential Imaging (ADI) which uses the field rotation induced in 1 or 2 hours of observation to reduce static speckle noise, the contrast ratios of an instrument like SPHERE can reach 10 ⁻⁶ . However, young Jupiter-like planets remain difficult to detect because of their very low contrast ratios of 10 ⁻⁵ to 10 ⁻⁷ .

We report here a new method of direct detection called Keplerian-Stacker (K-Stacker) based on stacking and recentering images with respect to the orbital motions of exo- planets taken during different nights of observation to make a potential faint companion visible although the latter is not detected in a single frame. This technique combined with an ADI reduction on SPHERE images could help to reach contrast ratios of 10 ⁻⁷ and it gave optimistic results in a previous work with sequences of 25 SPHERE/IRDIS simulated images.

The main goal of my work was to adapt K-Stacker to real data. We make an adaptation of the program to real observations by first simulating sequences of 5 SPHERE/IRDIS images which is closer to the number of real observations the SPHERE consortium has so far on several stars. We then inject planets at the limit of detection in each frame corresponding to a Signal to Noise Ratio SNR ∼ 3 which grows as a factor close to √

n where n is the number of images and reach in our case (n = 5) a total recentered SNR

∼ 7. We study the theoretical limits of our optimization algorithm to detect companions in those sequences of images by discussing tolerances on different parameters which have an effect on the positions of the planets projected onto the detector. Eventually we try K-Stacker on a few real observations.

Our simulations show that K-Stacker can detect planets at low SNR (SNR ∼ 3 in each

image) on sequences of a few images (typically 5 images) with different tolerances on the

accuracy of the True-North, the mass of the star, the distance observer/star, and the

number of planets in the images. The orbital parameters found to stack a planet are not

accurate if the planet did not move enough to reduce consequently the number of orbits

passing through the orbital positions but more observations spread over months or years

could be done if K-Stacker detects a new exoplanet in order to constrain its orbital param-

eters. Eventually, the work we did on real observations (NaCo, SPHERE/IRDIS) with

the adaptation of the program is a preliminary work aiming at proving that K-Stacker

can work on real images. The main issue we encountered was that there were very few

observations on a given star (< 5) in the database with the telescopes we used and those

observations were not always taken with a coronagraph. However at the end we had the

chance to work on 5 images of the star HD95086 with a good quality of reduction and

we report the results of the HD95086 b detection and also the detection of a false planet

injected in those real observations.

Introduction

The first exoplanet, 51 Pegasi b, has been found by Mayor et al. 1995 [12]. Since this epoch, detecting exoplanets has been one of the most challenging activity in astronomy.

Studying exoplanets allows the observation of exoplanets at different stages of evolution and radically challenges planet formation models which had been elaborated only based on the observations in the Solar System. Also, with the improvement of the detection methods, astronomers are trying to bring out Earth-like planets orbiting in the habitable zone of their host star where it is possible for liquid water to exist on the surface which is a prerequisite for life on Earth.

In order to detect exoplanets astronomers can use indirect detection methods or direct detection methods which will give different information on the observed exoplanets. Up to now, the main part of the ∼ 3500 exoplanets found results from indirect detection including radial velocity measurements, which use the movement of a star due to gravita- tional influence of a planet; photometric transit observations, which detect the variation of the stellar flux due to a companion passing through the line of sight to its host star;

astrometry, which also uses stellar motion; and gravitational microlensing, which involves unrepeatable observations biased toward planets with short orbital periods. However, these methods provide limited information about the planets and do not allow spectral characterization of the companions which could reveal atmospheric physics and chemistry as well as internal structure. These parameters are needed in the formation and evolution of planetary systems and only direct detection methods are able to reach those missing information. Consequently, direct detection is really complementary to the indirect de- tection methods and it allows to explore a comparative approach in exoplanetary science.

Thus, new facilities have been recently developed to implement direct detection methods using adaptive optics and coronagraphy. As planets are extremely faint light sources com- pared to stars with contrasts from 10 ⁻⁶ to 10 ⁻¹⁰ , it is very difficult to detect and resolve them directly from their host star. High contrast and high angular resolution are thus needed and these requirements urge the scientific community to develop instruments with very high performance such as VLT/SPHERE [1], Gemini/GPI [6] or Subaru/HiCIAO [2].

The improvement in correcting atmospheric turbulence in telescopes have reached such a high level with the eXtreme Adaptive Optics (XAO) that the detection imaging limita- tion is mainly due to intrinsic optical defaults in the instruments leading to static speckle patterns which prevent the detection of faint planets.

Techniques called Angular Differential Imaging (ADI) [10] and Spectral Differential Imag-

ing (SDI) [9] enabled to substantially attenuate the speckles and to obtain a high gain in

contrast. With these methods and their improvements (LOCI) [3], (TLOCI) [8] used on

SPHERE data, a contrast of 10 ⁻⁶ has been reached [20] [22] and a few more exoplanets

have been able to be detected. However the number of planets expected to be seen has

not been reached yet and astronomers keep trying to develop techniques to increase the

contrast and remove speckles from atmospheric turbulence and instrumental defaults.

A new way of direct imaging detection of exoplanets called K-Stacker was proposed in a previous paper [4]. The idea is to combine different observations of a star, spread over months or years, in order to increase the contrast limit of direct imaging instruments. A previous work done on 25 simulated SPHERE/IRDIS observations showed that K-Stacker could indeed be used to combine several high-contrast images taken in different nights to detect exoplanets (Nowak et al. 2017 submitted, A&A). Obviously, the next step is to deal with real observations and try to find an exoplanet.

My work at LAM focused on the adaptation of the program to real observations and is split in two big parts. The first one is a preliminary work on 5 simulated SPHERE/IRDIS images to be closer to the number of observations the consortium has so far, and to set the limits and tolerances on different parameters impacting the convergence of K-Stacker.

The second part of my work focuses on adapting the program and searching exoplanets

in real data. Eventually, I also simplified the use of K-Stacker for any user that would

work on it.

Chapter 1

Framework of direct imaging detection

1.1 Coronagraphy

In the first years of coronagraphy, a coronagraph was an optical instrument used in as- tronomy designed to reproduce a solar eclipse and study the solar corona. Nowadays new types of coronagraphs have been made for exoplanets detection with different techniques and masks to block the light coming from the host star of faint companions without block- ing the light of the companions.

1.1.1 Lyot coronagraph

The coronagraph was invented in 1931 by the french astronomer Bernard Lyot to study the solar corona [5]. This optical instrument is principally based on an amplitude mask and a diaphragm also called Lyot stop (Figure 1.1).

Figure 1.1: Basic scheme of the Lyot coronagraph

The mask has the size of the image of the Sun formed by the objective of the telescope and is placed at the focal plane of the latter to mask the Sun. But due to diffraction, there remains diffracted light which is not stopped by the mask and has to be blocked by a diaphragm. This diaphragm has a slightly smaller size than the image of the telescope’s objective (image of the pupil) formed by the intermediate lens and is placed in the plane of the telescope’s objective image. Finally a last lens forms an image of the solar corona on a detector.

However, when using this coronagraph for stellar coronagraphy, diffraction is a bigger

issue since the image of a star formed by a telescope is theoretically an Airy pattern

(Figure 1.2). The Airy pattern depends on the diameter of the entrance pupil D, the focal length of the telescope f and the wavelength λ. The radius of the central spot can be computed as follow R _Airy = 1.22 ^λ _D ^·f or expressed in radian for an aperture far from the image θ ∼ 1.22 _D ^λ .

Figure 1.2: Image of a star through a Lyot coronagraph

Moreover, the diameters of the mask and the diaphragm limit the field of view and may hide exoplanets one would like to detect. As a result, the mask should block the light of the star without hiding a companion but enough not to lost the planet in the flux of the star (Figure 1.3). Nevertheless, a simple Lyot coronagraph is able to reach a maximum contrast of ∼ 10 ⁻⁴ (in perfect conditions without atmospheric turbulence). The light of the planets (10 ⁻⁵ to 10 ⁻¹⁰ fainter than their host star) is still embedded in the diffraction light of the star even with a classical Lyot coronagraph.

Figure 1.3: Image of a planet with a Lyot coronagraph

An improvement of this set-up called Apodized Pupil Lyot Coronagraph (APLC) has been developed in a few papers [18] [16] [19] using an apodized entrance pupil. The apodizer profile can be defined as a prolate spheroidal function which exists for any aperture geometry and focal mask diameter to reach quasi-achromatic solutions.

1.1.2 Four Quadrant Phase Mask (FQPM)

Another way to improve the coronagraph is to use a phase-mask [14] instead of an am-

plitude mask to block the light. A phase-mask coronagraph uses a transparent mask to

shift the phase of the stellar light in order to create a self-destructive interference.

We describe briefly here the principle of the Four Quadrant Phase Mask (FQPM) [15]

and how it improves the set-up aiming to detect exoplanets.

The scheme is the same as Figure 1.3 but the amplitude mask is replaced by a FQPM (Figure 1.4) which gives a resulting image as in Figure 1.5.

Figure 1.4: Four-Quadrant Phase-Mask.

Two quadrants have a 0 phase shift and the two others have a π phase shift.

Figure 1.5: On the left, the Airy pattern of a star and a companion before the coro- nagraph. On the right, the image of the star and the companion after the coron- agraph.

In order for this phase-mask to work correctly, the star has to be centered very well on the mask so that the star light gets the different phase shifts and the resulting four beams combine into destructive interferences at infinity with an optimal attenuation of 20 mag for a perfect system. Note that the Lyot stop still remains to perform a spatial filtering of the diffracted light. As a result, an exoplanet will be probably localized in only one quadrant and its flux will not be attenuated. However, if a planet lies along an axis of the FQPM, it lights will experience a π phase shift and will be attenuated by a fac- tor of about 4mag which is non negligible. In this case, the planet may remain undetected.

Advantages of the FQPM:

-Strong attenuation of the star.

-No attenuation of the planet out of the four transitions phase in the mask.

Drawbacks of the FQPM:

-Not easy to make it achromatic because the phase shift depends on the wavelength.

The phase shift between two adjacent quadrants is given by : ∆ϕ = ^{2π(n −1)∆d} _λ where n is the refraction index of the mask, ∆d the difference of thicknesses between two quadrants and λ the wavelength.

-The farther the wavelength of observation is from the requested wavelength the mask has been made, the worse the attenuation will be.

-Very sensitive to centering (tip-tilt residuals at the focus of the telescope) in order to have the best destructive interferences.

1.1.3 Vortex Phase Mask

The Vortex Phase Mask stands in the continuity of the FQPM as an improvement of

the latter. It is a phase mask made of a concentric circular subwavelength grating which

induces local phase shifts in any points of the disk and result in a phase singularity in the center (Figure 1.6). This so called ”dark core” forces the intensity to vanish by a total destructive interference. Moreover this mask is achromatic for specific local characteristics [11].

The optical vortex is generated by a phase shift induced by e ^iθ , θ being the azimuthal angle. In the particular case of coronagraphy it is convenient to take ϕ = 2θ where θ is the azimuthal angle to have 2 × 2π phase shift for a full 2π rotation. Thus, when the PSF is centered on this type of vortex, the starlight is theoretically totally rejected outside the original pupil area and is then blocked by a Lyot stop (Figure 1.7).

Figure 1.6: Vortex phase mask made of annular grooves with period Λ and depth h.

Figure 1.7: Rejection of the starlight af- ter the vortex phase-mask. This repre- sents : F T [e ^i2θ · P SF ] where F T stands for the Fourier Transform.

1.1.4 Adaptive Optics

When observing a star with a based-ground telescope, turbulence in Earth’s atmosphere

limits the performance and this effect results in an image of the star broken in speckles

due to wavefront errors. A solution to this problem is to launch telescopes into space in

order to have a better resolution and get rid of the atmosphere’s turbulence. However

another way of skirting this issue keeping based-ground telescopes is to correct in real

time the defaults of the wavefront coming from the star as it passes through the Earth’s

atmosphere. An adaptive mirror reflect the light from the telescope and send it to a

beamsplitter which separates the beam in two; one which goes on a wavefront sensor and

the other one going to the camera. A control system assesses the distorted wavefront

in a smaller time than the wavefront defaults coherence time and adapts the shape of

the mirror to correct the wavefront (Figure 1.8). The difference between an image of a

target in the sky without adaptive optics and an image of the same target after adaptive

optics corrections is shown in Figure 1.9. This set-up must be placed upstream of the

coronagraph to be able to observe a corrected image but atmosphere’s turbulence and

instrument defaults lead to moving speckles and quasi-static speckles respectively. Most

of these speckles (star light not in phase) are not stopped by the coronagraph and limit the reachable contrast. The decorrelation time of the moving speckles is short enough to increase the detection limit by increasing the exposure time. However, there still remain quasi-static speckles arising upstream from the adaptive optics system not well-corrected by the sensor, and others intrinsically linked to the defaults of the instrument arising after the beamsplitter. Those artefacts are caused by the misalignment or the imperfections of the optics, and are also created by the differences of temperature surrounding the instru- ment. They are called ”quasi-static” since they evolve slowly in time with a correlation time τ _corr ∼ a few minutes and limit the detection of faint objects.

Figure 1.8: Adaptive optics system Figure 1.9: On the left, an observation without adaptive optics. On the right, the same observation with adaptive optics.

The target is a binary star. ESO credit.

1.1.5 Spectral Differential Imaging (SDI)

The idea of Spectral Differential Imaging is to remove quasi-static speckles arising from instrumental defaults. It consists in taking several images simultaneously of the same star in different wavelengths and substract a reference PSF after a rescaling of the images. If there is a planet orbiting around the star, it will appear in the images at the same place for every wavelengths whereas the speckles will have the same pattern but shifted for each wavelength. As the position of the speckles changes with the wavelength of observation, one can create a reference PSF without the planet. By rescaling the images, the speckles superimpose and the planet is shifted in each image. Then substracting the reference PSF, most of the speckles cancel out and the signal of the planet remains almost unchanged.

1.1.6 Angular Differential Imaging (ADI)

Another method to reduce consequently quasi-static speckles is called Angular Differential

Imaging. This technique consists in taking advantage of the ”instrumental” rotating field

of view during an observation of a star close to the transit. The instrument field derotator

is switched off so that the field of view of an altitude/azimuth telescope rotates while the

optics stay aligned with the instrument. During a sequence of images taken in this mode,

the planet moves with the field of view whereas the quasistatic speckles stay fixed in the

focal plane. Then, for instance, by taking the median of the images, it is possible to create a reference PSF and substract it from all the other images to gain in contrast ratio planet/star. Since 2006, a lot of more efficient methods have been proposed to create a reference PSF (see for example : LOCI, TLOCI, PCA methods [17]).

1.2 Telescopes

We present here in a nutshell two modules of the Very Large Telescope (VLT) located in Chile at Cerro Paranal : the Spectro-Polarimetric High-contrast Exoplanet-REsearch (SPHERE) and the NAOS-CONICA (NaCo).

SPHERE is directly the instrument managed by the consortium I joined at LAM and we had a collaboration with people of the NaCo team also, who kindly gave us a few images from this instrument to search exoplanets.

1.2.1 SPHERE at VLT

SPHERE is an extreme adaptive optics and coronagraphic facility divided into four sub- systems : CPI, IRDIS, IFS and ZIMPOL (see Figure 1.10).

This instrument is designed to give the highest image quality and contrast performance in a narrow field of view around a star in the visible or near infrared. SPHERE is installed at the UT3 Nasmyth focus on the VLT.

Figure 1.10: SPHERE common path infrastructure and the three science instruments.

The sub-systems included are the following ones :

CPI: The common path infrastructure which receives the light directly from the tele- scope and gives AO-corrected and coronagraphic beams to the three science instruments.

IRDIS: The infrared dual-band imager and spectrograph can provide classical images,

dual-band imaging, dual-polarization imaging, and long slit spectroscopy either between

0.95 - 2.32µm, with resolving power of R ∼ 50 or between 0.95 - 1.65µm with R ∼ 350.

IFS: The integral field spectrograph provides a data cube of 30 monochromatic images either at spectral resolution of R ∼ 50 between 0.95 - 1.35µm or at R ∼ 30 between 0.95 - 1.65µm.

ZIMPOL: The Zurich imaging polarimeter provides diffraction limited classical imag- ing and differential polarimetric imaging at < 30 mas resolution in the visible.

1.2.2 NaCo at VLT

NaCo is a first generation facility (Figure 1.11) which provides adaptive optics, coronagra- phy, imaging polarimetry and spectroscopy working in the 1 − 5µm range. It is composed of the NAOS system and the CONICA camera.

Figure 1.11: NAOS-CONICA facility.

NAOS: This is an adaptive optics system containing visual (0.45 − 1.0µm) and in- frared (0.8 − 2.5µm) wavefront sensors.

CONICA: The infrared camera and spectrometer attached to NAOS. It is equipped

with filters, masks, polarizing elements and grisms.

Chapter 2

Keplerian-Stacker (K-Stacker)

2.1 The K-Stacker principle

The main idea of K-Stacker (Le Coroller et al. 2015 [4], Nowak et al. 2016) is to combine several observations of a star taken on different nights separated by months or years in order to increase the Signal to Noise Ratio (SNR) and the contrast limit. Between the dif- ferent observations, the orbital position of a companion may have move substantially and we need to take its Keplerian motion into account when stacking the images. Obviously if the planet is detected in each image, the recombination is very easy and it only increases the SNR. But it is more interesting to focus on observations where a companion remains undetected in each single frame. In this case, the recombination should increase the SNR and help us to reach higher contrast levels. The images are rotated and translated ac- cording to the orbital motion of the companion so that the flux of a potential exoplanet in each image adds up while the background speckle noise is averaged and bring a detection in the final recombined image.

2.1.1 The orbital model

We use here the two-body problem in Newtonian dynamics to describe the orbital motion of the planets. Two hypotheses are made: the mass of the star is far more important than the mass of the planets orbiting around it, and the planets are not interacting with each other. Moreover, even if there are some interaction between the planets, these effects should be negligible for our problem because the orbits are not very well resolved (typically on orbit is described by 100 − 1000 PSF with SPHERE), and the deviations of positions due to these interactions will be smaller than the size of the PSF (in most of the cases).

The motion of a companion in the orbital frame can be described by three parameters: the eccentricity e, the semi major-axis a and the true anomaly θ which is the angle between the position of the planet and the perihelion seen from the star. However, it is convenient for reasons of calculations to introduce the eccentric anomaly E (see Figure 2.1) which is the angle between the perihelion and the position of the planet projected on the circle of radius a perpendicularly to the semi major-axis of the ellipse seen from the center of the latter.

To solve the two-body problem we use Kepler’s equation which describes the motion of the planet along its orbit:

E(t) − e · sin[E(t)] = n · (t − t 0 ) (2.1)

where t ₀ is the epoch at the perihelion, e the eccentricity, E the eccentric anomaly which depends on the epoch t and n = ^2π _T the orbital pulsation with T the orbital period.

Figure 2.1: Motion of a planet according to Kepler’s second law.

Figure 2.2: Position of the orbital frame (in blue) with respect to the frame of the observer (in red).

By solving equation (2.1) (using an iterative Krylov method [21] since the equation is non-linear) for E at all the epochs t, the positions of the planet can be determined and using geometrical properties of the ellipse one obtains the coordinates in the orbital frame:

 



x = a · (cos(E) − e) y = a · √

1 − e ² sin(E) z = 0

(2.2)

As we want to project the position of the planet onto a CCD plane, we need to compute the position of the planet in the observer’s frame. To do so, we use three Euler’s angles linking the two frames we are interested in (Figure 2.2). The longitude of the ascending node Ω is the angle between the line of nodes and a reference direction in the plane of reference, the inclination i is the angle between the orbital plane and the plane of reference, and the argument of the periapsis θ ₀ is the angle between the perihelion and the line of nodes in the orbital plane. The transformations from a system of coordinates to another can be done using rotations matrices associated to those three angles as follow:



 X Y Z



 = R ₃ (Ω) × R 1 (i) × R 3 (θ ₀ ) ×



 x y z



 (2.3)

x,y,z denoting the coordinates in the orbital frame and X,Y ,Z denoting the coordinates in the reference frame. R _k (α) where k ∈ {1, 2, 3}({x, y, z}) and α ∈ {Ω, i, θ 0 } is the rotation matrix of angle α along the axis k.

Finally, to get the final position (X ^′ ,Y ^′ ) on the camera, we set Z = 0 and divide the X and Y coordinates by the distance of the star:

[ X ^′ Y ^′ ]

= 1 d

[ X Y

]

(2.4)

2.1.2 Recombining the images : an optimization problem of the Signal to Noise Ratio

We simulate images that are close to the one produced by SPHERE/IRDIS, the diameter of the primary mirror is 8.2m, the plate scale is 12.25mas/pixel and the wavelength of observation is 1.6µm. The field of view of the A.O. corrected area (area where the maximum contrast is reached) is about 1.5” which corresponds to 128 ×128 pixels images.

Atmospheric speckles are simulated on the images and are not correlated from an image to an other and we do not simulate an ADI/SDI reduction as we could have with real data. The goal here is just to see if we can detect a faint companion at low SNR in each image (i.e. the companion is not detected in a single frame).

The algorithm we use to find exoplanets is mainly composed of three parts. The first part is mainly a brute force algorithm which will test all the possible orbits and save the q best SNR values on each point of a predefined grid after having loaded all the images to use and computed the background average and noise profiles at each radius in the images.

The second part will perform a gradient descent optimization on the q best values of SNR (the optimization is done on −SNR since it’s a gradient descent) and save the q best of them.

The third part is to assess by eye in the q stacked images if we detect a planet.

By using a brute force algorithm first and then a gradient descent algorithm to converge eventually towards the good final solutions, we expect not to miss any planet once the total SNR is around 6 or above.

2.1.2.1 First part: The brute force algorithm

First we define the SNR =

∑

s

√ ∑

σ

where p is the number of images, s _n and σ _n ² the flux and the noise (variance of the background intensity) respectively, integrated on an orbital position on the image n.

This part of the program loads n images and computes the background and noise profiles assuming a radial symmetry in the distribution of speckles. The background intensity is computed at a given radius r where we draw 1000 values of an angle θ using a random distribution in the range [ −π, π]. For each position the flux is integrated in a circular photometric box with a radius equals to the FWHM of the PSF corresponding to ∼ 3 pixels for SPHERE/IRDIS parameters. Then we can compute the mean background intensity over the 1000 values and take the standard deviation of the sequence to get the noise level. As a result, we obtain radial intensity and noise profiles for each image.

Then we define a grid to perform the brute force algorithm in order to sweep the field of view with a maximum of possible orbits. We cannot take an infinitesimal step for reasons of computation time so we first took the chart done by the previous developer of K-Stacker and then changed a few values for reasons explained in Section 3.1 (Table 2.1, the changed values are in bold characters).

The brute force algorithm is then executed 8 ×10×100×35×7×35 = 6.86×10 ⁷ times on

the n images. It computes the orbital motion of a potential planet in the images according

Param. Unit Min Max Step Points

a a.u. 2.5 7.5 0.7 8

e - 0 0.8 0.08 10

t ₀ year 0 20 0.2 100

Ω rad −π π 0.18 35

i rad 0 π 0.5 7

θ ₀ rad −π π 0.18 35

Table 2.1: Grid used by the brute force algorithm for 5 simulated SPHERE/IRDIS images with a star located at 10 pc and a mass of 1M _⊙ .

to the orbital model explained in Section 2.1.1 and it integrates the flux at a given orbital position determined by the epoch the observation was done in each image. Then it sums the fluxes and save the value. It also gets the noises computed previously on the radial positions and takes the square root of the sum of the variance in each image. As the algorithm tests all the possible orbits, this stage of the program takes the longest time of calculation and is thus split into 100 nodes of calculation running in parallel corresponding to the 100 points of t ₀ . To do so, we use the cluster of calculation at LAM where one node corresponds to a computation time of ∼ 10h for this set of parameters, so it is very convenient to have the nodes running in parallel. Once the brute force algorithm is done, we launch a program that computes the SNR on the 100 nodes running in parallel and save the q best values corresponding to the q best orbits found on each node. As a result, 100 files are saved with q SNR values each and also their corresponding orbital parameters.

The q best values are determined regarding the SNR threshold we choose to have a potential detection and the distribution of SNR values for randomly stacked speckles with p images (see Section 2.2.2).

2.1.2.2 Second part: The gradient descent optimization

When launching the gradient descent algorithm, a first part gathers in one file the 100 (number of nodes) previous files containing the best SNR values and orbital parameters.

Among those 100 × q values the program chooses the q best values and re-optimizes them with the Broyden-Fletcher-Goldfarb-Shanno minimization method (BFGS). It plots the orbits on which a potential planet has been found, saves the stacked images so that we can see if we detect a planet and collects the values of SNR and orbital parameters in a file.

2.1.2.3 Third part: Determining the detection of a companion

In this part we look at the q stacked images in order to determine if our K-Stacker

algorithm has found an exoplanet or has recentered bright speckles. In a paper by Males

et al. 2013 [7], we are told that recombining images with respect to an orbital motion

to increase the SNR will result in a false alarm detection probability of almost 1 because

of the huge number of orbits to test. However, and this is where assessing by eye the

presence of a potential planet is crucial, on every simulations we have done so far we

have never found a single false positive. The result is striking but there might be an

explanation to this. The false positives they talk about in the paper are defined above a

particular SNR threshold of 5 which could be underestimated for our simulations since, for instance, we have a lot of recentered points in the range 4 ≤ SNR ≤ 5.5 in our simulations corresponding to stacked speckles but when looking at the recombined images, we do not see any brighter point than the surrounding noise located where the algorithm recentered the images. When the recentered SNR is 5 ≤ SNR ≤ 5.5 we sometimes see a spread spot but the shape does not look like a PSF and we conclude that this spot is not a detection of a planet. However, when we see a PSF shape with a SNR higher than 5.5 the conclusion is that K-Stacker has found a planet and by looking retrospectively at the injection of the companion in the images we always get a true detection.

In any case, if we see a PSF shape recentered on real observations for a SNR for which

the bright spot is a bit drowned in the noise, it would be an opportunity to have the star

observed once more to confirm or not the detection of a new exoplanet.

Chapter 3

Adaptation of K-Stacker to work on real observations

The main goal of my project is to adapt the K-Stacker algorithm to real data.

All the previous work done on the K-Stacker algorithm was achieved with 25 images with 1 seconde exposure each, simulated with SPHERE/IRDIS characteristics and dominated by speckle noise. However images from real observations are often a sequence of fewer images with a longer time exposure.

In this section, we thus choose to focus on simulations with only 5 images (for now, it is rare to be able to find more than 5 observations on one target with the main SPHERE survey) of 20 secondes time exposure each and we draw conclusions about the theoretical limits of the algorithm on images that looks close to real data in order to then work on real images. Figure 3.1 is a typical simulated SPHERE/IRDIS image used in my work.

Figure 3.1: 128 × 128 simulated SPHERE/IRDIS image of 20 secondes time exposure

with speckle noise. The inner mask has a radius of 20 pixels.

3.1 Determining the grid for the free parameters

As mentioned in the 2.1.2.1 part, I changed a few values in the parameters the brute force algorithm uses compared to the ones used in the previous work.

In order to give relevant steps to the grid we want to use with the brute force algorithm, we plot the distribution of SNR in the images versus each orbital parameter. We explore

Figure 3.2: Distribution of SNR in 5 simulated SPHERE/IRDIS images versus the ec- centricity e. A planet has been injected at e = 0.5 and the other orbital parameters are randomly fixed. The FWHM of the global maxima is 0.08 and defines ∼ 10 steps for the eccentricity e ∈ [0, 0.8].

two possibilities to define our steps. The first one is to inject a planet in the images with a total SNR = 7 on a random orbit and plot the SNR in the recentered image as a function of each parameter while fixing the others. As the SNR for speckles should not be higher than SNR = 5.5 for the highest ones, we should see a clear peak standing out where the planet is located (Figure 3.2), around SNR = 7, and we use the FWHM of this peak as a step for the studied parameter. Taking the FWHM as a step should enable the brute force algorithm to fall close enough to the global maximum for one point of the grid and then be chosen as one of the best local maxima to be reoptimized by the gradient descent algorithm. On the other hand, we can also process a run of K-Stacker on images without any planet and plot the distribution in the same way as before but this time we count the local maxima due to stacked speckles and define the step of the studied parameter as the number of those local maxima.

These two methods seem to give the same results at ±1 step each time we plot the distribution. For the chart presented in Table 2.1 we chose to inject a planet on a random orbit and then plot the distribution of SNR.

The first parameter I changed in the grid is the semi major-axis a. I set the minimum value to 2.5 a.u. instead of 0.9 a.u. because we use a numerical circular mask of 20 pixels radius which corresponds to ∼ 2.45 a.u. for a star located at 10 pc. As a result we chose to divide the range of a into 8 points instead of 10 points.

The maximum value of the epoch at the perihelion t ₀ is 20 years which corresponds to

the period of a circular orbit at 7.5 a.u. corresponding to 61 pixels of radius for a star located at 10 pc, so on the edge of the A.O. corrected area of the field of view. This maximum value resulted in 100 steps of 0.2 years according to the distribution of SNR in the images versus the epoch at the perihelion. To compute this maximum value we used the third Kepler’s law with the distance in Astronomical Units, the time in year and the stellar mass in Solar Mass Unit where the gravitational constant G is then equal to 4π ² :

T ²

a ³ = 4π ²

GM _star ⇔ T =

√ a ³

M _star (3.1)

The inclination range in the previous work was first i ∈ [−π, π] with 13 steps. However, I noticed that the projected orbits on the CCD camera were exactly the same when i swept the range [ −π, 0] or [0, π] and we derive the reduction of domain in Appendix A. Thus, we chose to set i ∈ [0, π] with 7 steps and as a result, the changes we did for the definition of the grid enabled to go from more than 20 hours of computation for each node to ∼ 10 hours.

3.2 Number of reoptimized orbits

When the brute force part is done, we want to reoptimized q orbits with a gradient descent algorithm. However, we do not reoptimized the 6.86 × 10 ⁷ orbits (section 2.1.2.1) because it would take a huge amount of computation time and lots of the reoptimized orbits would be irrelevant. Thus, we need to set a threshold SNR value at which we agree that a planet is detectable and see the minimum number of values we should reoptimize in order not to miss a planet. To do this we run K-Stacker on 5 simulated SPHERE/IRDIS images without planet and plot the distribution of SNR values after the brute force algorithm Figure 3.3.

Figure 3.3: Normalized histogram representing the distribution of SNR values on 5 sim-

ulated SPHERE/IRDIS images without planet after the brute force algorithm. The star

is located at 10 pc, the stellar mass is 1.0M _⊙ and the total time of observation is 3 years.

We look at the number of SNR values found at several percentiles of the distribution to see the number of false positives we would have for different SNR thresholds. The 99.99999-th percentile corresponds to ∼ 7 values located beyond SNR = 4.63 which means that if we had a planet at this SNR in the stacked images we should reoptimized at least 8 orbits to have a chance to detect the planet among the reoptimized speckles. The 99.999999-th percentile is already at ∼ 0.7 value located beyond SNR = 5.01 meaning that we can choose a threshold SNR > 5.01 to avoid false positives in our simulations.

However, with our criterion of seeing a planet in the final stacked image, we fix the threshold at SNR = 5.5 and we are clearly not limited by a minimum number of values to choose regarding the histogram. As the reoptimization process is quite fast, we choose an arbitrary number of orbits which combines a reasonable computation time with a reasonable number of images to check. We thus keep q = 100 orbits corresponding to a computation time around 30min and it enables us to see the behaviour of the gradient descent optimization algorithm when it finds a planet.

3.3 Modification of the calculation of background and noise profiles

The first issue I dealt with was the calculation of noise in the simulated images. The previous work which had been done on K-Stacker computed the background and noise values before injection of the planet in the image and then used them to retrieve the companion.

However, in real data, the planets we would like to find are already in the images and the background and noise profiles can only be computed taking into account the flux of the unknown planets.

Figure 3.4: Mean intensity computed on a ring in a simulated coronagraphic SPHERE/IRDIS image without planet versus the distance in pixels to the center of the image.

Figure 3.5: Mean intensity computed on a ring in a simulated coronagraphic SPHERE/IRDIS image with a planet at 4.5a.u. (36.7 pixels) at SNR ∼ 3 versus the distance in pixels to the center of the image.

Thus we focused on the impact the flux of a planet could have in the resulting background

and noise profiles. We did several profiles for injections of a planet at different positions

Figure 3.6: Noise level computed on a ring in a simulated coronagraphic SPHERE/IRDIS image without planet versus the distance in pixels to the center of the image.

Figure 3.7: Noise level computed on a ring in a simulated coronagraphic SPHERE/IRDIS image with a planet at 4.5a.u. (36.7 pixels) at SNR ∼ 3 versus the distance in pixels to the center of the image.

in simulated SPHERE/IRDIS images with a SNR such that the planet is at the limit of detection in one single frame, typically SNR ∼ 3 in each image.

Note that above SNR ∼ 3, the planet is seeable and K-Stacker has no real interest;

moreover, at this level (SNR ∼ 3) the error on the estimation of the noise is maximum (when K-Stacker is testing an orbit with the planet in the noise). Thus, if the estimation of the noise is not to much perturbated by a planet at SNR ∼ 3 it will also be good for planets at lower SNR.

The curves after injection of a planet do not show a real difference with the profiles computed without any planet (see Figure 3.4, 3.5, 3.6, 3.7). Those results were expected since the companion is injected at a very weak signal, i.e. ∼ 5% of the mean intensity at a given position. The two main bumps of the profiles correspond to the diffraction effect due to the edges of the mask and the Lyot stop.

To convince ourselves more, we also did different runs of K-Stacker with background and noise profiles computed with and without the planet. The results we were given were the same, meaning that when we had (or not) a detection of a planet it occured for both types of profiles.

3.4 Detection near the edge of the inner mask

Close to the mask a bright annulus of noise appears due to diffraction. We thus wonder whether this noise peak would be a problem to detect a planet drowned in this area. A planet orbiting on a circular orbit inside this annulus corresponding to a radius of 3 a.u is simulated. We change the epoch at perihelion for different simulations in order to sweep the entire annulus with the orbital positions and we vary the SNR (see Table 3.1).

The results show that the detection appears when SNR _tot = 7 (Figure 3.8). As the

noise is stronger in this region, the flux of the planet needs to be stronger also to keep

a SNR _tot = 7. We could expect false positives rising when the detection of the planet is

unsuccessful but we did not get any speckles recentered with a total SNR sufficiently high

Epochs at perihelion t 0 and SNR 4 5 6 7 8

t ₀ = 0 year X X X V V

t ₀ = 1.0 year X X X V V

t ₀ = 2.0 years X X X V V

t ₀ = 3.0 years X X X V V

t ₀ = 4.0 years X X X V V

Table 3.1: Detection (V symbol) and non-detection (X symbol) for 5 simulated SPHERE/IRDIS images spread over 2 years of observations of circular orbits at radius 3 a.u. with different epochs at perihelion and different SNR. One period corresponds to 5.2 years for a stellar mass M = 1M _⊙ .

Figure 3.8: Detection for a planet at SNR tot = 7 close to the mask with 5 stacked images.

The planet follows a circular orbit of radius 3 a.u.

to be detected. This can be explained by the fact that speckle statistics do not follow keplerian motions and the algorithm stacks speckles with very different fluxes which always sum up to SNR _speck ≤ 5.5 for the best re-optimized values.

3.5 True North default

A parameter to take into account during astronomical observations is the accuracy of the

alignment with the True North. For K-Stacker, we need the images all oriented in the

same direction. An error on the True North is responsible for misalignment between the

frames. Some frame may undergo a default of a few degrees with respect to the True

North. The True North accuracy for SPHERE should be much better than ±0.5 ^o per

observation. However, the consortium realized that observations they did during a few

months period could have been done with a True North default around ±1 ^o because of a

bad calibration. Thus there was a need to simulate several observations with a random

default angle of a few degrees of amplitude between the images in order to see the impact

it has on the position of a planet injected at different distances from the star and up to which point the algorithm can find the companion.

We first made a calculation to evaluate the displacement in pixels a certain default leads to, with respect to the radius of the field of view.

For our simulated images in 128 × 128 pixels, the maximum radius is 64 pixels which represents about 7.8 a.u. for a star at 10pc since the plate scale of SPHERE/IRDIS is 12.25 mas/pixel. A default amplitude of dθ = ±0.5 ^o on the True North at this distance corresponds to δ = R.dθ = 0.56 pixel = 0.17 FWHM of a PSF at λ = 1.6µm (the FWHM is 3.37 pixels at this wavelength) on the position of a planet at the epoch t.

We simulate a planet orbiting circularly at several radii from its star with random defaults in the 5 stacked images for different ranges of True North default amplitudes (see Table 3.2).

True North default vs Orbits a = 3.4a.u. a = 4.5a.u. a = 6.0a.u. a = 7.2a.u. Flat

±0.5 ^o V V V V V

±1 ^o V V V V V

±2 ^o V V X X X

±3 ^o V X X X X

±4 ^o X X X X X

Table 3.2: True North default in degrees versus different orbits : the four first orbits are circular and the ”flat” column corresponds to a particular situation where the orbital positions are aligned along a radius of the field (orbit parallel to the line of sight, viewed from the observer position). V corresponds to a good detection whereas X corresponds to a non-detection. The total SNR is 6.

Figure 3.9: ”Flat” orbit with parameters [a = 7.0 a.u., e = 0.1, t0 = 10.0 years, Ω = 0.0, i = ^π ₂ , θ ₀ = 0.0], the orbital plane is orthogonal to the observer’s plane. The red crosses are the five positions of the planet evenly spread during three years of observations.

The ”Flat” column corresponds to an orbit which appears to be flat after projection on

the CCD (Figure 3.9) which is the worst case scenario since the points are on a radius of the FoV and will be shifted on both sides of the radius with different pixel displacements for a given True North default.

Figure 3.10: The big dots in blue are centered on orbital positions with True North default of ±2 ^o ran- domly spread (here the sequence is [ −2 ^o , +2 ^o , −2 ^o , +2 ^o , 0]) for a circular orbit of radius a = 3.4 a.u. The red disks are centered on the best orbital position K-Stacker found.

Figure 3.11: Stacked image for a True North default of 2 ^o on a circular orbit.

This stack corresponds to the red points on Figure 2.13.

By converting the degrees into pixels displacement for each configuration we clearly see that in each case the companion is not detected when the positions are shifted more than 1.5pixels = 0.44 FWHM of a PSF. Below 1.5 pixels displacement the planet is found by the algorithm. In Figure 3.11 we show an example of true detection for a circular orbit of radius a = 3.4 a.u. and a True North default of 2 ^o with the best orbit found by the algorithm (Figure 3.10, 3.11).

Thus, the conclusion is that the tolerance accepted on the True North should not be higher than 1.5 pixels displacement on an orbital position which corresponds for our images to

±1 ^o on the edge of the FoV. However, the ”flat” case is more sensitive to displacements

and for 1 ^o of True North default the final image looks like Figure 3.13. Here the stacked

points result in a bit blurry signal which is a mixture of a PSF shape and a random

speckle shape stuck on it. By looking at the points the algorithm stacked on Figure 3.12,

the integration of the fluxes are quite accurate despite the non-orbital positions of the

injected companion but the stacking leaves a sort of tail on the left side of the true PSF.

Figure 3.12: The big dots in blue are centered on orbital positions with True North default of ±1 ^o ran- domly spread (here the sequence is [ −1 ^o , 0, −1 ^o , +1 ^o , +1 ^o ]). The red disks are centered on the best orbital position K-Stacker found.

Figure 3.13: Stacked image for a True North default of 1 ^o on the flat orbit. This stack corresponds to the red points on Figure 2.14. The companion is located in the red circle.

3.6 Tolerance on the mass of the star

With today’s devices and techniques, the mass of a star (observed in the surveys of coronagraphy) can be known with a good precision (i.e. less than ±0.1M ⊙ ) but in some case it still remains an uncertainty of a few solar masses. It is thus interesting to simulate the search for a planet injected in the images with a star at a given mass and let the algorithm try to find it with another stellar mass. Those simulations will provide the tolerance on the stellar mass the algorithm could bear to find a companion and to choose relevant stars to observe. However, if we do not want to be limited in the choice of the stars to observe, we could still add the mass of the star as a free parameter in the grid but it would increase the computation time by the number of steps included in the range of search for the stellar mass.

We introduce a planet in the images with a total recentered SNR = 6 at different radii in the field of view and for different stellar masses. Then we launch K-Stacker setting the mass of the star at 1.0M _⊙ (Table 3.3) and at 2.0M _⊙ (Table 3.4) and see if it detects the companion.

Around 1M _⊙ the accepted tolerance on the mass of the star is ±0.1M ⊙ which means that if we are given an observed star with an accuracy on its mass worse than ±0.1M _⊙ we could launch the search for a companion at different stellar masses with a step of 0.1M _⊙ around the given mass.

For a stellar mass of 2.0M _⊙ the results are quite similar as the previous ones except for

a = 3.4 a.u. where the tolerance is less than ±0.1M ⊙ but more simulations would be

required to better constrain statistically this result. For higher radii the tolerance seems

better than around 1M _⊙ showing that K-Stacker manages to compensate for the stellar

mass default with the free orbital parameters. As a result, Table 3.4 informs us that we could search also a planet for different stellar masses with a step of 0.1M _⊙ around heavier masses.

Mass default vs Orbits a = 3.4 a.u. a = 4.5 a.u. a = 6.0 a.u. a = 7.2 a.u. Flat

∆M = −0.2M ⊙ V V X V X

∆M = −0.1M _⊙ V V V V V

M = 1.0M _⊙ V V V V V

∆M = +0.1M _⊙ V V V V V

∆M = +0.2M _⊙ X X X X X

Table 3.3: Mass default around 1.0M _⊙ versus different orbits : the four first orbits are circular and the ”flat” column corresponds to a particular situation where the orbital positions are aligned along a radius of the field. V corresponds to a good detection whereas X corresponds to a non-detection. The total SNR is 6.

Mass default vs Orbits a = 3.4 a.u. a = 4.5 a.u. a = 6.0 a.u. a = 7.2 a.u. Flat

∆M = −0.2M ⊙ V V V V X

∆M = −0.1M _⊙ V V V V V

M = 2.0M _⊙ V V V V V

∆M = +0.1M _⊙ X V V V V

∆M = +0.2M _⊙ X V V V X

Table 3.4: Mass default around 2.0M _⊙ versus different orbits : the four first orbits are circular and the ”flat” column corresponds to a particular situation where the orbital positions are aligned along a radius of the field. V corresponds to a good detection whereas X corresponds to a non-detection. The total SNR is 6.

3.7 Accuracy on the distance of the stars

Another concern we had to deal with was the accuracy on the distance of an observed star and how it may affect the orbital positions in the images. We looked up in the SPHERE database to see the uncertainties on the distances of the stars and we gather in Table 3.5 the mean percentage of error on the distances for different ranges for a 12.25mas/pixel plate scale (SPHERE/IRDIS) to see the impact it may have on the convergence of K- Stacker.

Range of distances in pc [1, 10] [10, 30] [30, 70] [70, 150] [150, 270]

Average percentage of error ±0.63% ±2.75% ±0.84% ±11.68% ±23.34%

Average error in pixel ±0.10 ±0.11 ±0.03 ±0.10 ±0.10

Table 3.5: Range of distances of the stars observed with SPHERE and how the errors on the distances affect the orbital positions of a planet.

On average, the shift of the image due to errors on distances of the stars does not exceed

±0.1pixel which means that K-Stacker should not be sensitive to this default regarding the tolerance we found for the True-North (1.5 pixel). However if a distance is given with an accuracy resulting in a shift higher than ±1.5 pixel, it should be relevant to launch K-Stacker on different values of distances with a step chosen to have less than 1.5 pixel shift for each run.

3.8 Influence of several companions in the images

Here we consider 5 simulated observations with 3 planets hidden in each image in order to see if K-Stacker succeed in finding all the planets. The companions are injected in the images with a SNR ∼ 3 in each frame over a total time of observation of 3 years. Table 3.6 gives the orbital parameters of the 3 planets.

Orbits a (a.u.) e t ₀ (years) Ω i θ ₀

Orbit 1 32 0.2 50 0 0.3 0

Orbit 2 26 0 10 0 0 0

Orbit 3 20 0.3 50 0 2.5 0

Table 3.6: Orbital parameters of the 3 orbits on which the planets were injected at SNR

∼ 3 in each image. The mass of the star is 1.5M _⊙ and the distance is d = 50pc.

We first reoptimized with the gradient algorithm the 2000 best orbits the brute force algorithm found and only the planet moving on Orbit 2 was detected. Then reoptimizing by packs of 2000 orbits for the next best ones ( ∼ 10h of computation for each pack), the planet moving on Orbit 3 was detected a few times in the pack [8000, 10000] best orbits among the detections of the Orbit 2. We did not go further than the 10000 first best orbits for reasons of time and practicability when dealing with future real observations.

Nevertheless what one can do when a planet has been found is to mask the orbital positions of the planet (i.e. Orbit 2) in each image regarding the best orbital parameters found by the algorithm with a radius of the mask ∼ 1 FWHM. As a result, the best orbits found in the next run of K-Stacker would hopefully be those from another planet (i.e. Orbit 3) that we could mask at the end of this run and launch K-Stacker a last time to find Orbit 1.

Indeed, we found Orbit 3 after having masked the orbital positions for Orbit 2 (Figure 3.14). Eventually we mask the orbital positions of Orbit 2 and 3 to see if K-Stacker finds Orbit 1. Once again the planet is found using this method of masking the previous orbital positions found which is really convenient.

Obviously, on real observations, the number of planets in the images will not be known

before running K-Stacker. However if a companion is just above or clearly above the limit

of detection by eye, we could mask it and then run K-Stacker to find putative planets

that we cannot see in each single frame.

Figure 3.14: Stacked image by K-Stacker for the detected planet on Orbit 3. On the

bottom left part of the image, one can see the masks on the 5 images for the planet on

Orbit 2 shifted following the recentering with respect to Orbit 3.

Chapter 4

Results on real astronomical observations

4.1 Images of GJ3331A with NaCo

We were given 5 observations of the GJ3331A star done on NaCo spread over 6 years without coronagraph. The star is located at 18.328pc and has a mass of 0.5 ± 0.1M ⊙ . Unfortunately the ADI could not be used on these observations and a less effective reduc- tion was done on the images. However, we wanted to try K-Stacker on those data for a first work on real observations.

Figure 4.1: One typical observation of the star GJ3331A located at 18.328 pc after reshaping and recentering the im- age. The scale is in pixel and the total field of view corresponds to 3.48” mean- ing a maximum circular orbit radius of 31.8 a.u.

Figure 4.2: The same observation of GJ3331A after a numerical masking of 10pixels radius ( ∼ 5 a.u.). The quasi- static speckles are clearly visible despite the data reduction and the diffraction rings are not removed since there is no coronagraph.

First of all there was a need to recenter and reshape the images to only keep the useful

field of view we wanted to use (Figure 4.1). The plate scale on NaCo is 27.19mas/pix and

we chose to keep an image of 128 × 128 pixels which corresponds to a total field of view of 3.48”. In astronomical units the outer part of the field is located at 31.8a.u.

We masked numerically the central disk of stellar light in order not to have K-Stacker searching for a companion in this area where the flux of the star is too high to have a chance to detect a planet (Figure 4.2).

As the observations were done in the L ^′ band ( ∼ 4µm) we computed the size of the PSF and the FWHM we had to use to define the radius of the circular boxes to integrate the fluxes on orbital positions. At this wavelength, the FWHM ∼ 4 pixels for NaCo plate scale which is thus the radius of the circular boxes.

4.1.1 Search for a companion with K-Stacker

Param. Unit Min Max Step Points

a a.u. 6.0 31.0 2.7 10

e - 0 0.8 0.08 11

t ₀ year 0 250.0 1.7 147

Ω rad −π π 0.24 27

i rad 0 π 0.6 6

θ ₀ rad −π π 0.24 27

Table 4.1: Grid used by the brute force algorithm for NaCo with the star GJ3331A located at d = 18.33 pc and a mass of 0.5M _⊙ .

Figure 4.3: One of the best stacked image of GJ3331A reoptimized by K-Stacker.

Figure 4.4: The same stacked image with

different signals encircled. The red circle

is the sum of the signals K-Stacker recen-

tered and stacked whereas the green cir-

cles stand for unfortunate similar spots

due to remaining quasi-static speckles

looking like the one K-Stacker found.

We first gather the results of the steps to define the grid on those images in Table 4.1 following the methods we explained in Section 3.1. The expectations for these images were not very high since they were taken without coronagraph, meaning a very weak contrast, but it was the first opportunity to have K-Stacker running on real observations and see if the code was well adapted to real images. As one can see on Figure 4.3, we cannot conclude that K-Stacker has found a planet since several spots are similar to the one stacked by the algorithm and Figure 4.4 is needed to know which spot was really found. Moreover among the observations we got, some of them had strong quasi-static speckles mixed with diffraction rings which obviously lead to this type of images. The total recentered SNR for the spot encircled in red is SNR _tot = 5.3 which is right in the range of bright speckles SNR _speck ≤ 5.5.

4.2 Images of HD95086 with SPHERE/IRDIS

We worked on 5 observations of HD95086 taken with SPHERE/IRDIS spread over 1 year.

The star is at a distance of 90.4 ± 3.4pc and has a mass of 1.7 ± 0.1M ⊙ . We know that the planet HD95086 b was discovered in 2013 by direct detection with VLT/NaCo [13]

so we would like to see if K-Stacker is able to find it. Then we use the same data with a false injected planet before PCA reduction and try to find it after reduction for several magnitude differences between the star and the planet.

4.2.1 Search for a companion with K-Stacker

Once again we first determined the grid on which we would launch K-Stacker but now for the IRDIS detector. The steps are gathered in Table 4.2.

Param. Unit Min Max Step Points

a a.u. 11.5 87.5 3.8 20

e - 0 0.8 0.08 10

t ₀ year 0 627 10.5 60

Ω rad −π π 0.3 20

i rad 0 π 0.3 10

θ ₀ rad −π π 0.3 20

Table 4.2: Grid used by the brute force algorithm for SPHERE/IRDIS with the star HD95086 located at d = 90.4 pc and a stellar mass of 1.7M _⊙ .

The planet HD95086 b is seeable in 3 frames over 5 but we want to see if K-Stacker increases the total recentered SNR by a factor close to √

5. Thus, we launch the program on the 5 images and the 100 best orbits found correspond to HD95086 b with a maximum total recentered SNR _tot = 7.66 (Figure 4.5). By looking at the best orbital parameters stacking the orbital positions we compute the SNR in each frame for the corresponding orbital positions. The mean SNR in each frame is 3.48 which would lead theoretically to a total SNR _tot = 3.48 × √

5 = 7.78 (this law is true for a normal random distribution,

however speckles statistics do not follow a normal law but we show here that the result

does not differ too much from what we could have for a normal distribution). This result

is consistent with the value found by K-Stacker and strengthens our hypothesis about