Covariance of lucky images: performance analysis

The covariance of ground-based lucky images is a robust and easy-to-use algorithm that allows us to detect faint companions surrounding a host star. In this paper, we analyse the relevance of the number of processed frames, the frames’ quality, the atmosphere conditions and the detection noise on the companion detectability. This analysis has been carried out using both experimental and computer-simulated imaging data. Although the technique allows us the detection of faint companions, the camera detection noise and the use of a limited number of frames reduce the minimum detectable companion intensity to around 1000 times fainter than that of the host star when placed at an angular distance corresponding to the few ﬁrst Airy rings. The reachable contrast could be even larger when detecting companions with the assistance of an adaptive optics system.


I N T RO D U C T I O N
Atmospheric refractive index fluctuations affecting the image quality of a ground-based telescope have been a common topic in astronomy for years.Different techniques, like Speckle Interferometry (Weigelt & Wirnitzer 1983) or Adaptive Optics (AO; Hardy 1998), have been developed to reach the ideal telescope angular resolution.As the atmospheric fluctuations are random, one expects these fluctuations to be occasionally arranged in such a way that they produce a diffraction-limited image.Consequently, David Fried (Fried 1978) suggested to take a series of short exposures images and then select the best ones, i.e. those images with best Strehl ratio.This technique is named lucky imaging (LI).Since then the technique has remarkably evolved (Mackay 2013).
The LI technique seems to be very promising for medium-sized telescopes because of its low complexity and costs in terms of hardware.Furthermore, LI works with reference stars fainter than those required for the natural guide star AO technique.A key parameter when using LI is the image exposure time.The decorrelation timescale of the atmosphere in the case of LI is about 30 ms.Hence, if we want to freeze the atmospheric evolution, the LI exposure time must be around this coherence time.
The image of a point source (the point spread function, PSF) obtained by a perfect optical system with circular aperture can be described by the Airy pattern.However, in ground-based telescopes, where the atmosphere refractive index inhomogeneities distort the incoming wavefront, this image consists of a central peak surrounded by a number of speckles whose temporal average is commonly known as halo.The central peak, corresponding to the coherent part of the incoming wavefront, is added to an incoherent halo due to the incoherent wavefront energy.The shape of the distorted PSF that we obtain will depend on D/r 0 , the ratio between the telescope diameter D and the Fried parameter r 0 , which is the atmospheric coherence length.The number of speckles contained in the PSF halo is roughly given by (D/r 0 ) 2 and they are randomly distributed over an area with angular diameter λ/r 0 .The number of speckles and the area covered by them are strongly dependent on the wavelength since r 0 depends on the detection wavelength (or band).In general, a good balance between the central peak height and incoherent halo height is found for a diameter telescope up to 2.5 m, when observing at I band (700-800 nm wavelength).
Among the different aspects involved in an actual experiment, the most important is, probably, the atmospheric conditions.In a previous paper (Cagigas, Valle & Cagigal 2013), we have shown that the LI technique can only be effective for D/r 0 ratios around 8, since for D/r 0 values over this limit the speckled halo is higher than the coherent peak and the image selection is not possible.
Recently, we have proposed an algorithm based on the temporal evolution of the pixels intensity along a series of lucky images COELI2 681 (Cagigal et al. 2016).The intensity of all those pixels where a faint companion is placed will fluctuate in phase with the main-star intensity along the image series.However, the pixels containing incoherent speckles will fluctuate in counterphase.The finding of pixels in the image series, which are fluctuating in phase, with pixels gathering light from the main star is accomplished by evaluating the normalized covariance between the main star and the rest of the image pixels along the selected LI series.The result is a kind of bi-dimensional covariance map in which the pixel value is the normalized covariance value.This technique can be applied for extracting undetected faint companions from the speckled background in a region that is around the host star and that has an inner radius of 1.22λ/D and outer radius of 1.22λ/r 0 .
Although we have already introduced the principles of the covariance of ground-based lucky images (COELI) algorithm (Cagigal et al. 2016), in this paper we present, along with a brief description of the COELI fundamentals, an analysis of the detectability of a companion as a function of its intensity and distance to the host star, detection noise and atmosphere conditions.The dependence on the number of frames used in the covariance estimate and the effect of the LI frame selection are also included.
We tested the COELI technique using a set of experimental LI images of the object GJ822 taken at the I band by the 2.2 m diameter Calar Alto telescope.Besides, the analysis of the effect of the atmosphere conditions on the companion detectability has been carried out by computer simulation.
Finally, we have compared images obtained by applying the COELI technique with those obtained by the Hubble Space Telescope (HST) for the object GJ822.

C O E L I F U N DA M E N TA L S
Before performing a deeper analysis of the algorithm, it is interesting to state the context in which the algorithm will be applied.We will consider a stack of frames so that the frame exposition time is close to the coherence time of the atmosphere fluctuations.The intensity covariance between the zero mean intensity i(r) of two pixels at positions r 1 and r 2 can then be defined as (1) We will always take τ = 0 since we are working with temporally incoherent frames.
In addition, for the zero delay case, the covariance will also tend to zero when the pixels do not contain an object (halo pixels).When one of the pixels contains an object, covariance will be negative and only when both r 1 and r 2 pixels contain an object, the covariance will be positive and its value will depend on its relative intensity.

COELI algorithm
The algorithm we named COELI is composed of the following steps.
(1) We first perform the image stack centreing.As a criteria for an efficient centreing, we have chosen the superimposition of the most intense pixel of every frame.
(2) The second step consists of the image pre-processing.Detection noise can be softened by convolving a Gaussian filter with every frame of the stack (Gonzalez & Woods 2002).This step can be omitted or even performed after the step 3.
(3) Then, we estimate the normalized covariance (correlation) between the most intense peak of the reference star and the rest of the frame pixels along the frame series [i(r 1 ) = i sp , i(r 2 ) = i(r) and τ = 0, in equation ( 1)].
All the preceding steps have been included in an ImageJ plugin named COELI (ImageJ Plugins).
The procedure we followed was to estimate the expression: where Cov[,] stands for covariance, C[,] stands for normalized covariance or correlation, σ is the standard deviation and is the ensemble average (frame series average).The intensity i sp is the star peak intensity defined as the addition of the coherent peak intensity (i cp ) plus the intensity star halo at the centre (i h ): where detection noise has been considered negligible compared to the star peak intensity.The intensity detected at a position r = (i, j) is the addition of the star halo background plus detection noise, i(r) = i h (r) + i n .However, in those pixels where there is an object, it would be necessary to add the object intensity, i o :

Covariance contrast
To evaluate the capability of this technique for detecting objects, we estimate the contrast of the object against the background in the covariance map.The contrast is defined as the ratio between the covariance value of the pixel containing an object and the covariance value of the background pixels.For this purpose, covariances are estimated using equation ( 2), with and without the term i o in equation ( 4).We first define the object peak intensity i o , which is proportional to i sp : (5) On the other hand, the intensity star halo is a function of the distance to the main star and it is related to i cp through the expression: where k(r) states the halo intensity spatial dependence.In this simplified model, we will consider k as a constant.Finally, the readout noise intensity is given by i n .
Introducing equations ( 3)-( 6) into equation ( 2), the covariance contrast between pixels containing an object and those without an object for pixels inside the halo can be expressed by where σ 2 x is the variance corresponding to the intensity i x .To obtain equation ( 7), we have considered that the value of σ 2 ir only depends on the speckle and detection noises and it does not change significantly when a companion is found in that position.Covariance terms Cov[i cp , i n ] and Cov[i h , i n ] will tend to zero and Cov[i h , i h (r)] will also tend to zero for a distance r greater than the speckle size.Under these considerations, the contrast can be approximated by Hence, from equations ( 7) and ( 8), we can see that the contrast will increase as the companion intensity (k o ) increases and will decrease when the covariance between the different noises takes high values.

Covariance map estimate
In an actual experiment, the theoretical covariance expression given by equation ( 2) has to be evaluated using the following estimator where N is the number of frames used for the covariance estimate and Hence, to achieve a good contrast estimate, the number of samples has to be high enough to allow the covariance terms containing speckle or detection noise to tend to zero.A detailed analysis of the errors affecting the covariance estimate can be found in (Stuart & Ord 1994).

E X P E R I M E N TA L C H E C K I N G
To check the COELI technique, a series of experimental measurements were completed.The observations were carried out in 2013 September with the LI camera Astralux at the 2.2 m telescope of the Calar Alto Observatory in Almería/Spain.This instrument incorporates a fast readout electron multiplying CCD chip (EMCCD) that is able to acquire images with a very low readout noise, thanks to internal charge amplification before conversion to voltage by an output amplifier.Astralux allows the acquisition of a large number of images, typically several thousand for each target, with exposure times of about a few tens of milliseconds, which allows us to freeze atmospheric speckles.Conventional long integration time averages these speckles producing a seeing-limited PSF.
The observations were done in Sloan Digital Sky Survey I band with a pixel scale of 47 mas pixel −1 and 7000 images were acquired, each with exposure time of 30 ms.The internal electron multiplying gain was adopted to work in the EMCCD linear regime and therefore determined by the luminosity of the target.To carry out a precise calibration of the pixel scale and camera rotation, we observed the core of the globular cluster M15, and correlated the astrometric data of detected point sources with their precise astrometry, obtained with the HST.This provided an accurate astrometric calibration for each observing night with plate-scale precision as good as 0.1 mas pixel −1 in general.The main source of astrometric error for a given star results from the uncertainty in the measurement of the barycenter, which is in turn mainly determined by the signalto-noise ratio.In our data, this uncertainty is typically 0.1 pixel and reaches 0.2 in the faintest stars.This leads to typical errors in separation down to about 10 mas.

C O E L I P R E P RO C E S S I N G
The performance analysis of the technique is quite difficult because many different parameters are involved in an actual experiment.Hence, in this section we will analyse separately the aspects we consider more relevant.
A key point is the way in which the COELI technique is applied.In our case, there are different blurring or sharpening filters that can be applied before carrying out the covariance estimate.The use of a particular filter or of a combination of them will depend on the case to be analysed.In particular, before applying the COELI algorithm, we may convolve the stack of lucky images with a Laplacian filter to improve the point objects detectability or even with a Gaussian filter to smooth the image noise before applying the Laplacian filter.
Fig. 1 shows the comparison between the results obtained using different processing algorithms.Fig. 1(a) shows the results of applying the shift-and-add (SAA) algorithm to a stack of 100 images of the star GJ822, obtained as described in Section 3. Fig. 1(b) shows the result of processing the same stack as in Fig. 1(a) but with the COELI algorithm.The result of convolving the stack image series with a one-pixel-radius Laplacian mask before applying COELI is shown in Fig. 1(c).For obtaining Fig. 1(d), the stack was convolved with a one-pixel-radius Gaussian mask before using the COELI algorithm.
It can be seen that the speckled halo surrounding the host star (Fig. 1a) disappears, whereas the companion is clearly seen when the COELI algorithm is applied (Fig. 1b).When a Laplacian filter is convolved with the stack, the halo area is noisier (Fig. 1c).Finally, the use of a Gaussian mask has a double effect.On one hand, the object signal is stronger than before but, on the other hand, there is a loss of resolution.

C O N T R A S T D E P E N D E N C E O N T H E N U M B E R O F F R A M E S
Since the COELI algorithm is based on the experimental estimate of the intensity pixels covariance, as described by equation ( 9), the number of frames (N) used for the estimate plays an important role.In general, one can expect that larger the number of frames, better the results are obtained.However, in an actual experiment, the image quality may fluctuate according to the instantaneous atmosphere seeing.Hence, to estimate the covariance map, we have to choose between using a short number of excellent-quality frames or a large number of frames including the poor-quality ones.In general, it is necessary to find a balance depending on the number of frames and the seeing value during the frames acquisition.

Image quality
Fig. 2 shows the result of applying LI frame selection according to the intensity of the brightest spot to the stack containing 1000 images of the star GJ822.Fig. 2(a) shows the result of applying COELI to the best 100 images of the stack.In Fig. 2(b) appears the result of applying COELI to the worst 100 images of the stack.It can be seen that the effect of using the 100 worst images is that the size of the object clearly increases when compared to the object size obtained from the 100 best images stack.Hence, we can conclude that not only the number but also the image quality have to be taken into account for obtaining the best results.

Number of frames
To check the effect of the number of frames on accuracy of the covariance map, we started from a 1000 frame stack of the star GJ822.From this stack, we have created a new stack containing the 100 best images and another one containing the 300 best images.
As it can be seen in Fig. 3, the noise decreases as the number of frames, included in the stack, rises.The result was expected given the statistical nature of the estimator.We have already commented that covariance terms in equation ( 7) containing the noise intensity will tend to zero as the number of frames increases.However, bearing in mind the results of Fig. 2, it may not always be useful to include a large number of frames since, although the result improves with the number of frames, the use of poor-quality images may ruin the whole processing.It means that it is necessary to find a balance between number of frames and their quality.

Effects of atmospheric seeing
Among the different parameters affecting the companion detectability, like its distance or magnitude difference to the host star, the instantaneous conditions of the atmosphere are of crucial importance.Given the difficulty of finding experimental measurements covering the whole range of cases, we have carried out a computer simulation analysis.
The telescope outcome is generated from the Fourier Transform of the entrance pupil field described by random phase screens.We have simulated a stack of 1000 phase screens following the Kolmogorov statistics (Roddier 1990).The phase-screen size was 1024 × 1024 pixels, whilst the phase information was contained in an area of 128 × 128 pixels in order to have a good sampling balance.By rescaling the original phase-screen stack, we have obtained three different stacks of images corresponding to the cases of D/r 0 = 3, 6 and 9, since the limit of applicability of the LI technique is around D/r 0 = 8.In this particular analysis, the set of simulated images has not been affected by detection noise, although the noise due to the Fast Fourier Transform algorithm is present.
We have taken a stack containing 1000 frames of a point source, duplicated the stack, translated the number of pixels equivalent to the corresponding distance (in units of Airy rings), multiplied it by a reducing factor and added it to the original stack.The resulting stack simulates the presence of a companion fainter than the host star placed at a particular distance.This simulation technique has been widely used for simulating binary stars, see Bagnuolo (1982), Lee & Yee (2003) and Gladysz & Christou (2009).
Fig. 4 shows that the minimum detectable companion relative intensity decreases as the distance to the host star increases.This behaviour is repeated for all the D/r 0 values analysed.It is worth to note that the kind of image obtained after the COELI processing  is just a covariance map.Hence, the Airy rings, whose intensity fluctuates in phase with the main peak intensity, may clearly appear affecting the companion detectability.This fact is particularly noticeable for low D/r 0 values, as shown in Fig. 4.
Fig. 5 shows the covariance map obtained by applying COELI and a Laplacian mask to a 1000 frames stack containing a simulated star plus a companion 250 times fainter.In Fig. 5(a), the companion is placed at an angular distance corresponding to the sixth Airy ring of the host star, for the case D/r 0 = 3.The companion can clearly be detected and the host star Airy rings can also be seen.
However, for higher D/r 0 values (D/r 0 = 9 in Fig. 5b), the Airy rings disappear from the covariance map, and the distance at which an object 250 times fainter than the host star can be detected, increases up to the ninth Airy ring.Therefore, as D/r 0 increases, it is more difficult to detect faint objects.This result agrees with equation (8) so that an increase of D/r 0 , which implies a σ 2 h increase, reduces the object contrast.The contrast improvement would require to increase the number of samples or, when possible, the assistance of an AO system (Basden 2014).

Effects of the detection noise
To perform the analysis of the effect of the detection noise, we have used peak intensity and noise standard deviation values similar to those experimentally obtained for the object GJ822.Hence, we will use in our analysis noise standard deviation values around the experimental one.
The detection noise clearly affects the covariance map obtained by applying COELI, since it affects cross-terms in equation ( 7).However, the resulting map not only depends on the detection noise standard deviation but also on the D/r 0 value.As an example, Fig. 6(a) shows the covariance map obtained from a simulated stack for D/r 0 = 3, so that the highest peak takes an averaged value of 1 and the detection noise has a standard deviation σ = 0.2.When the standard deviation decreases to σ = 0.04 (Fig. 6b), the halo area clearly appears.
The same simulation has been carried out for D/r 0 = 9.  result for σ = 0.04.In both cases, D/r 0 = 3 and D/r 0 = 9, when the background noise is high (Figs 6a and c), the improvement provided by COELI is really limited since the dominant noise comes from the detection process, not from the speckle statistics.For low-noise background (Figs 6b and d), a dark region appears surrounding the central star when applying COELI, since it is able to cancel out the light following the speckle statistics.In Fig. 6(d), the dark area is larger than in Fig. 6(b) since an increase of the D/r 0 value leads to a larger speckled area.

F I E L D O F V I E W
Usually we have restricted ourselves to the analysis of companions around a host star.However, the technique can be implemented to images containing a series of objects distributed across the whole field of view.All the speckled haloes surrounding every star located in the isoplanatic angle will disappear by applying COELI.As an example, we show the result of using SAA (Fig. 7a) and COELI (Fig. 7b) to a stack of 100 frames of globular cluster M15 taken with Astralux at the 2.2 m telescope of the Calar Alto Observatory.It can be seen that the covariance map contains not only the same objects as the SAA image, but many other clearly appear.
The image contrast was set in both cases to allow dim objects to be clearly visible.Distances and relative position are very similar in both images.The main difference is that in the COELI image the photometry is lost, given that it is a covariance map.However, the resolution is much better in the COELI image, which could be useful for detecting possible faint companions located at short distances of the host star.
Hence, the technique can be successfully applied not only to a small area surrounding the host star, as in the case of searching for faint companions, but also to a wide field to perform other kind of analyses.Besides, the peak intensity of a bright star can be used for estimating the covariance map around any faint star within the same isoplanatic angle.

C O M PA R I S O N TO S TA N DA R D T E C H N I Q U E S
In a recent paper, Ginski et al. (2016) show images obtained with Astralux at the Calar Alto 2.2 m telescope of a series of newly detected comoving companions of exoplanet host stars.We paid attention to Kepler-21 that exhibits a comoving companion, Kepler-21B, which is located about 0.8 arcsec south-east of the exoplanet host star and whose intensity is about 100 times fainter than the host star.With standard LI techniques and high-pass filtering, the companion was recovered with a signal-to-noise ratio of about SNR = 1.5 using the best 10 per cent of all images (50 000) of the star, taken with AstraLux (see Fig. 8a).With COELI processing, we clearly recover the companion at much higher signal-to-noise ratio (SNR = 8) using fewer images (500), i.e. a selection rate of the best 1 per cent of the total stack.The COELI algorithm allows us now distinguishing it clearly from the background.Fig. 8(b) shows the COELI map obtained from a stack containing the 1 per cent of the 2013 data.Experimental contrast curves (Fig. 9; which have been estimated according to the procedure stated in Ginski et al. 2016) confirm that the contrast inside the halo area is much higher when using COELI.This improvement makes COELI algorithm extremely useful to search for companions close to the host star.
It can also be seen that the COELI performance outside the halo region is clearly outperformed by the standard data reduction.The reason is that the COELI algorithm has been particularly developed for the removal of light, whose temporal statistics is that of the star halo.

C O N C L U S I O N S
After the analysis shown in the preceding paragraphs, we can conclude that the evaluation of the covariance of lucky images is a robust  algorithm that allows us to detect faint companions surrounding a host star from ground-based images.Since the technique requires a previous lucky image selection, it can only be applied for D/r 0 < 9, which implies a maximum telescope aperture of roughly 2.5 m.
The object detectability provided by the technique depends on different aspects like the number of frames, the frame quality, the seeing during the frame series detection or the camera detection noise.
We have shown that objects 1000 times fainter than the host star can be detected when located at an angular distance equivalent to the fifth Airy ring.However, reducing the camera detection noise and using AO for improving image quality, this contrast can even be increased.

Figure 1 .
Figure 1.Result of processing a stack of 100 images of the object GJ822 using: (a) the SAA algorithm, (b) COELI, (c) Laplacian mask plus COELI (d) Gaussian mask plus COELI.

Figure 2 .
Figure 2. (a) Result of applying COELI to a stack containing the best 100 images obtained from a stack of 1000 images of the star GJ822.(b) As in (a) but selecting the worst 100 images of the stack.

Figure 3 .
Figure 3. Result of processing a stack containing the object GJ822 using COELI.The number of frames in each stack is (a) 100, (b) 300 and (c) 1000.

Figure 4 .
Figure 4. Minimum detectable relative intensity (in log scale) obtained from a 1000 frames stack as a function of the companion position (in Airy rings) for D/r 0 = 3 (blue), 6 (red) and 9 (grey).

Figure 5 .
Figure 5. Result of processing a stack of 1000 frames containing a simulated star plus a companion 330 times fainter using COELI.(a) The companion is placed at the sixth Airy ring for D/r 0 = 3.(b) The companion is placed at the 11th Airy ring for D/r 0 = 9.

Figure 6 .
Figure 6.Result of processing a stack of 1000 frames containing a simulated star using COELI once the stack has been convolved with a Laplacian mask.(a) The detection noise has a standard deviation of 0.2 for D/r 0 = 3.(b) The detection noise has a standard deviation of 0.04 for D/r 0 = 3. (c) The detection noise has a standard deviation of 0.2 for D/r 0 = 9.(d) The detection noise has a standard deviation of 0.04 for D/r 0 = 9.

Figure 7 .
Figure 7. Result of processing a stack of 100 frames of globular cluster M15.(a) Applying SAA.(b) Using COELI.

Fig. 8
Fig.8(c) shows standard LI using the best 1 per cent of all images (50 000).In this case, the companion and background intensities are similar.Experimental contrast curves (Fig.9; which have been estimated according to the procedure stated inGinski et al. 2016) confirm that the contrast inside the halo area is much higher when using COELI.This improvement makes COELI algorithm extremely useful to search for companions close to the host star.It can also be seen that the COELI performance outside the halo region is clearly outperformed by the standard data reduction.The reason is that the COELI algorithm has been particularly developed for the removal of light, whose temporal statistics is that of the star halo.

Figure 8 .
Figure 8. Result of processing a stack of 500 frames of Kepler-21: (a) standard reduction and high-pass filter (10 per cent image selection), (b) COELI result (1 per cent image selection), (c) standard reduction (1 per cent image selection).

Figure 9 .
Figure 9. Contrast curves experimentally obtained of processing a stack of 500 frames of Kepler-21, (red) standard reduction and (blue) COELI processing.