Issue 
A&A
Volume 586, February 2016



Article Number  A74  
Number of page(s)  10  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/201526933  
Published online  27 January 2016 
Evaluating the effect of stellar multiplicity on the point spread function of spacebased weak lensing surveys
Laboratoire d’astrophysique, École Polytechnique Fédérale de Lausanne (EPFL), Observatoire de Sauverny, 1290 Versoix, Switzerland
email: thibault.kuntzer@epfl.ch
Received: 10 July 2015
Accepted: 25 November 2015
The next generation of spacebased telescopes used for weak lensing surveys will require exquisite point spread function (PSF) determination. Previously negligible effects may become important in the reconstruction of the PSF, in part because of the improved spatial resolution. In this paper, we show that unresolved multiple star systems can affect the ellipticity and size of the PSF and that this effect is not cancelled even when using many stars in the reconstruction process. We estimate the error in the reconstruction of the PSF due to the binaries in the star sample both analytically and with image simulations for different PSFs and stellar populations. The simulations support our analytical finding that the error on the size of the PSF is a function of the multiple stars distribution and of the intrinsic value of the size of the PSF, i.e. if all stars were single. Similarly, the modification of each of the complex ellipticity components (e_{1},e_{2}) depends on the distribution of multiple stars and on the intrinsic complex ellipticity. Using image simulations, we also show that the predicted error in the PSF shape is a theoretical limit that can be reached only if large number of stars (up to thousands) are used together to build the PSF at any desired spatial position. For a lower number of stars, the PSF reconstruction is worse. Finally, we compute the effect of binarity for different stellar magnitudes and show that bright stars alter the PSF size and ellipticity more than faint stars. This may affect the design of PSF calibration strategies and the choice of the related calibration fields.
Key words: gravitational lensing: weak / methods: data analysis / binaries: close
© ESO, 2016
1. Introduction
Weak gravitational lensing is one of the main cosmological probes to study dark energy and dark matter (e.g. Weinberg et al. 2013; Frieman et al. 2008). In practice, this statistical method requires measuring the shape of billions of faint and small galaxies with the least possible contamination by the instrumental point spread function (PSF). Spacebased wide field surveys provide the highest quality check to perform these challenging measurements. Two missions are either under construction or in project to conduct these wide field optical and nearIR surveys: the ESA Euclid satellite^{1} (Laureijs et al. 2011) to be launched in 2020 and the NASA Wide Field InfraRed Survey Telescope (WFIRST; e.g. Spergel et al. 2015).
Even with the high spatial resolution of a space telescope, any weak lensing experiment must account for the small but nonnegligible instrumental PSF, which needs to be removed accurately. The weak lensing distortions are an order of magnitude smaller than the PSF effects on the shape of the observed galaxies (e.g. Kilbinger 2015). To achieve the required accuracy, the PSF must be characterised using stars in the field of view. However, in practice, the high spatial resolution needed to resolve faint galaxies also means that double or multiple stars may significantly affect the PSF measurement, i.e. its size and ellipticity. The goal of the present work is to assess the bias and error introduced by stellar multiplicity on the PSF determination, which is crucial for weak lensing experiments (e.g. Cropper et al. 2013; Massey et al. 2013).
Binary and multiple stellar systems are very common (see review by Duchêne & Kraus 2013, hereafter DK13); they formed in and their characteristics depend on the initial environment and mass (Reipurth et al. 2014). Several volume and magnitudelimited surveys were dedicated to the search of binary stars, allowing a rigorous construction of catalogues with reliable distribution inferences beginning with the work of Duquennoy & Mayor (1991) and compiled by DK13. Nevertheless, large uncertainties remain in the distribution of the multiple system parameters, even in the solar neighbourhood, reflecting the difficulty of such surveys in accounting for the properties of multiple stellar systems in detail. The frequency of multiple system varies drastically depending on the mass of the primary object, with almost 100% of the most massive stars being multiple systems and on average 35–40% of mainsequence stars being at least binary (DK13).
Fig. 1
Reconstruction errors for a single binary star with the main component centred on (0,0). Each point on the map shows the measurement of the quantity labelled in each panel when a single companion star is located at that position. The labels are in arcsecond. The first column on the left shows the effect for a PSF aligned along e_{1} and a contrast of Δm = 0 between the main star and its companion. The next column is for the same PSF shape and orientation, but for a contrast of Δm = 1. The block of six graphs on the right is for the same two contrasts, but for a PSF aligned along e_{2}. The red line represents the orientation of the PSF. From top to bottom, the total ellipticity (star + companion) of the PSF, the relative error on the ellipticity due to the companion star, and the relative error on the size are shown. 
In this paper, we simulate the error in the reconstruction of the PSF arising from stellar multiplicity. The information on the distribution of multiple systems and their properties (separation, contrast) is fairly scarce but remains sufficient to simulate a spacebased survey. For this purpose, we use the Besançon model of the Galaxy (BMG; Robin et al. 2003) in combination with the binary star fraction characteristics drawn from observed stellar distributions summarised in DK13. These data allow us, first, to build mock catalogues with realistic stellar properties and, second, to simulate spacelike images of stars. With these simulations, mocks, and images, we derive the bias and reconstruction errors on the PSF width and ellipticity. We also derive analytical predictions and show that the analytical predictions agree very well with the results of the simulations, even for a realistic spacelike PSF with complex structures and spikes.
We first study the effect of one single companion as a toy model of the problem in Sect. 2. In Sect. 3, a formalism is developed to analytically compute the expected PSF reconstruction error and bias. In Sect. 4, we discuss the current knowledge of the parameters of multiple systems in the Milky Way, describe how we construct the stellar population, and detail the simulations that support our analytical findings. We present the results of the simulations and their parameters in Sect. 5 and conclusions are drawn in Sect. 6.
2. Reconstruction error of PSF for a single binary star
The typical physical separations in binary stars are of the order a = 50 AU and typical light contrast is between 0 and 1 mag. With the resolving power of a Euclidlike telescope, 0.1′′, binary stars can be resolved up to distances of d = 2.5 kpc. These resolved binaries can be identified and removed from samples of stars to be used for the PSF reconstruction. However, most binaries cannot be explicitly identified and still affect the PSF determination, which is an effect that we illustrate in this section.
We generate a simple simulated image that uses a spacelike PSF with the GalSim software (Rowe et al. 2015). This simulation is composed of a main star located in the centre of the image to which we add a companion star at a given position within the central pixel. We then measure the error in the reconstruction of the PSF of the double star and compare its ellipticity and size to the values found for a single star.
The resulting errors are presented in Fig. 1 for two different ellipticities and for two different contrasts in magnitude, Δm = 0 and Δm = 1. Each PSF has its major axis along e_{1} and e_{2}, respectively. The colour maps in Fig. 1 show the ellipticity of the reconstructed PSF, the relative errors on the ellipticity, Δe/e_{0}, and on the size, , where e_{0} and are the intrinsic ellipticity and size of the PSF in the absence of a companion.
Two effects are striking. First, the relative errors can easily reach 1% in ellipticity when the companion star is as close as 0.01′′ from its host. This figure is reached when the contrast in magnitude is Δm = 0. Even for a larger contrast of Δm = 1, the error remains of the same order of magnitude. Similarly large effects are seen on the size measurement. Second, the shape of the error maps for the ellipticity depend on the PSF orientation and ellipticity, while the error maps for the size show no obvious dependence on the intrinsic PSF ellipticity.
The above test suggests that the presence of binary stars can have an important impact on the PSF reconstruction. They either must be avoided or their effect must be accounted for in the error budget of the mission. While in principle the Gaia catalogue is expected to be complete up to V_{lim} ~ 20−25 mag (Robin et al. 2012), binary and multiple stars catalogues will not be complete for the small separations in which we are interested in this work (de Bruijne et al. 2015). In the following, we investigate, both analytically and with image simulations, the effect of unidentified binary stars on the PSF reconstruction for a Euclidlike telescope.
3. Expected deviations and reconstruction error
3.1. General case
Fig. 2
Sketch of the system with the main object at , with total intensity I_{m} and the companion located (δx,δy) away with a total intensity of I_{m}f = I_{c}. 
We investigate in an analytical manner how the PSF is altered by the presence of multiple stellar signals. We show that (i) the bias on the complex ellipticities is a function of both the characteristics of the multiple stars distribution and of the intrinsic complex ellipticities, and that (ii) similarly, the bias on the size is a function of the multiple stars distribution characteristics and of the intrinsic size.
Let q_{xx} be the quadratic moments along the x axis of the image of a star, q_{yy} along the y axis and q_{xy} along the x−y diagonal, (1)where j,k are either x or y. The parameter denotes the position of the centroid for the axis a. The image is centred on the main star, which is at . In a noisefree image, the complex ellipticity components are obtained via (Kaiser et al. 1995): (2)We compute the expectation values for the quadratic moments to mimic the result of the stacking of many binary systems. We restrict this analysis to the case of the same intensity for all main stars, while each companion may have a different intensity ratio f_{k}. All companion stars are randomly distributed around the main star. To get the expectation value of the quadratic moments, we use the limit case when the number of systems in the stack N_{⋆}− → ∞. We assume that the light profile of the companion star is given by I_{c}(x,y) = fI_{m}(x−δx,y−δy) and total intensity of I_{c} = fI_{m}.
For a binary system, the centroid in the x axis is located at (3)The computations are analogous for the y axis. We compute the total quadratic moment of the system k, q_{xx,k}, in the xx direction, via the parallel axis theorem. The quadratic moment of the binary system k along xx is denoted q_{xx,m}, i.e. (4)The derivation for q_{yy,k} is similar.
3.1.1. Biases for complex ellipticities
The position of the companion stars are random, which implies that ⟨ δx^{2} ⟩ = ⟨ δy^{2} ⟩, where the notation ⟨·⟩ denotes the expectation value. The expectation value for the complex ellipticity component e_{1} is (5)As the term ⟨ fr^{2}/ (1 + f)^{2} ⟩ is close to zero in Eq. (5), its Taylor expansion is taken and we get the value of the bias for the system k,(6)for i = 1,2 and e_{i,0} the ith component of the complex ellipticity, such that the expectation value for the bias is (7)where is the intrinsic size of the PSF defined as the sum of the quadratic moments in the xx direction plus the moments in the yy direction. This quantity is referred to as T in Kaiser et al. (1995). Similarly, for e_{2}, (8)The above set of equations shows that the bias, δe_{i}, is a function of the intrinsic ellipticity of the PSF, e_{i,0}, and of the companion distribution through the quantity ⟨ r^{2}f/ (1 + f)^{2} ⟩. It also shows that the magnitude of the bias is the same both in e_{1} or e_{2} components. The intrinsic size of the PSF can mitigate the effect of the binaries. If we use the same binary population for small and large PSFs (i.e. ensuring that the same binaries are excluded in both cases), the bias is smaller than for small PSFs. Another conclusion drawn from the presence of a minus sign is that the presence of binaries in the stack tends to make the reconstructed PSF rounder than the intrinsic PSF shape.
3.1.2. Bias for size
The size of the PSF is R^{2} = q_{xx} + q_{yy}. The bias on the size is written . Computing the bias gives (9)Thus changes in size depend on the binary star fraction distribution and not on the intrinsic PSF ellipticity.
3.2. Ideal reconstruction error
The reconstruction error can be bounded by a requirement for future weak lensing surveys. We write the reconstruction error as (10)where s = 1...N_{s} runs over the N_{s} stacked images of stars and e_{i,s},e_{i,0} are the ellipticity components of the resulting stacked image and of the ideally reconstructed PSF, respectively. If N_{s} is large enough, we can write (11)Similarly, for the size (12)
3.3. Suboptimal reconstruction error
The above equations are derived for an image that contains a very large number of coadded stellar systems N_{⋆}. If this number gets too small, shotnoise error dominates, introducing an extra source dispersion in the reconstruction error σ. When N_{⋆} is too low, the expectation value ⟨ r^{2}f/ (1 + f)^{2} ⟩ will no longer be representative of the stellar binary distribution in the image. We discuss the minimal value of N_{⋆} to escape this suboptimal regime in Sect. 5.3.
The separation r and ratio of the intensity f distributions are typically lognormal as described in DK13. The bias δe is the multiplication of those two lognormal distributions, thus the bias is also an exponentiallike distribution. The standard deviation of such exponentiallike distributions has a dependence (e.g. Hoogenboom et al. 2006). At low N_{⋆}, this σ_{N} is more important than the reconstruction error σ(x), where x is either the ellipticity component e_{i} or the size R^{2}. The parameter σ(x) is the binary PSF reconstruction error signal, while σ_{N} is a statistical artefact. The total measured error σ_{meas} is a combination of σ_{N} and σ(x). For low N_{⋆} or low binary PSF reconstruction errors, the statistical error σ_{N} dominates the reconstruction error. On the other hand, when enough binaries are stacked or when the PSF properties exhibit a large value of the size or the ellipticity, the binary PSF reconstruction error σ(x) dominates.
4. Binary stars in the Milky Way
Fig. 3
Mainsequence stellar population breakdown for different Galactic latitudes. The Galactic longitude on the left panel is ℓ = 0° (towards the centre of the disk) and the right panel ℓ = 180° (Galactic anticentre). The bin of massive stars (>1.5 M_{⊙}) is not visible as it represent less than 0.2%. We stress the different behaviour of the binary fraction for massive stars (blue vs. green bins) for the opposite directions in Galactic longitude. 
We now test our predictions for the PSF reconstruction errors. To accomplish this, we draw a realistic population of stars from the Besançon Model of the Galaxy (BMG; Robin et al. 2003). To this population of single stars, we add companions according to the properties of binary stars described in Duchêne & Kraus (2013, hereafter DK13), including the orbital and mass ratio parameters. We simulate the PSF reconstruction process by stacking many observations of a mix of stellar systems with and without companions. We then measure the ellipticity and size of the stacked light profile and compare them with the intrinsic values in the absence of any stellar companion.
In the following section, we present the properties of the multiple systems in the Milky Way and discuss the demographics of binary stars with respect to their position on the sky before we describe the simulations themselves.
4.1. Characteristics of multiple stellar systems
Duchêne & Kraus (2013) found that multiple systems are characterised by two main quantities: (i) the frequency of multiple system, MF, i.e. the number of occurrences of a binary system in a given population; and (ii) the companion frequency, CF, the average number of companions per star, which can exceed one. We restrict ourself to binary stars and only use the frequency of multiple systems, MF, to draw populations of binary stars. The physical separation distributions follow a lognormal curve. The detailed representation of the lognormal distribution function suggests different formation processes for multiple systems (e.g. Heacox 1996). The mass ratio of the companion mass to the primary mass, q ≤ 1, can be determined from the observation of the flux ratio assuming massluminosity relations. The distributions of the mass ratio follow an empirical power law. Table 1 summarises important trends in the population of binary systems. Some of the main trends are:

1.
The more massive the primary star, the larger the probability of acompanion,

2.
The more massive the primary star, the lower the mass ratio powerlaw index γ. This implies that massive stars tend to host companions of lower relative masses than lowmass stars,

3.
The more massive the primary star, the broader the physical separation with the companion.
4.2. Spatial distribution of binaries
The demographics of stellar populations depend on the Galactic coordinates and most notably on the Galactic latitude, and since different regions of the Milky Way have different stellar populations, they have different fractions of multiple systems. This fraction is mainly controlled by the fraction of massive stars. An informed choice of the stellar population used for the PSF reconstruction can improve the purity of the sample by minimising the fraction of multiple systems.
We illustrate this effect in the following with several realisations of the BMG mainsequence stars at different Galactic coordinates. For each simulation, the number of binary stars is evaluated. Figure 3 shows two different examples for Galactic longitude ℓ = 0°, near the Galactic centre and ℓ = 180°, in the direction of the Galactic anticentre.
Towards the Galactic centre, the stellar population demographics remain approximately constant for latitudes between b = 15° and b = 75°. Thus, there is no significant variation in the fraction of binary stars. On the other hand, pointing close to the anticentre results in an increasing fraction of massive stars with increasing Galactic latitude. Stars of masses between 0.5 M_{⊙} and 1.5 M_{⊙} represent less than 30% of the stellar population at b = 15°, while they account for almost half of the population at b = 75°. When this observation is combined with the binary star information of Table 1, this leads to a statistically significant increase of the fraction of binary stars with increasing Galactic latitude, as shown in Fig. 3. This means that the impact of the multiplicity of stars on the PSF determination is not the same at different locations on the sky, not only because of the projected number density of stars but also because of the varying stellar populations. It also means that PSF calibration fields must be carefully chosen and the fraction of multiple stars in the selected fields must be accounted for. In the following, we choose the pointing ℓ = 180°,b = 45°, as the population demographics and fraction of binary stars seems to be representative of all available data and because this pointing is located at a portion of the sky surveyed by Euclid (Laureijs et al. 2011).
4.3. Synthetic images of single and binary stars
Our goal is to measure the PSF reconstruction errors for a given number of stars that can be coadded but that are sometimes binaries. To this end, several simulations with different fractions of binary stars are generated. The stars potentially useful for the PSF reconstruction of future spacebased surveys have magnitudes in the range i(AB) ~ 16−24, depending on the exposure duration and detector saturation levels. Thus, we only keep stars falling into this range in the sample. Mainsequence stars more massive than 1.5 solar masses are discarded as well owing to the lack of data about their multiplicity. However, these stars represent less than 0.2% of the whole dataset. We also discard stars that have more than one companion. All our simulations are therefore on the optimistic side in this respect.
The stellar populations generated by the BMG do not contain stellar companion information. We therefore add companions to the BMG stars using the information in DK13 and summarised in Table 1. We replace missing values with the value in the lowest bin of mass, which represents a conservative estimate.
A massluminosity relation for lowmass mainsequence stars is used to turn the mass of the companion into a magnitude (Habets & Heintze 1981). For simplicity, we assume that all the stars have the same flat spectral energy distribution. In other words, our simulated images are for flat SEDs but the allocation of companions to the stars follows the prescription of DK13 applied to the stellar populations of the BMG.
The position of a companion star along the assumed circular orbit is drawn from a uniform distribution. The angle between the plane of the orbit and the line of sight, as well as the longitude of the ascending nodes, are also uniformly distributed. The light profile of the companion is normalised according to the massluminosity relation, and interpolated and shifted by the required vector. Magnitudes of both objects are modified such that the total apparent magnitude of the system is the original apparent magnitude given in the BMG catalogue. We account for the large uncertainties in the stellar multiplicity parameters and explore the behaviour of the PSF reconstruction at different binary star fraction, by randomly modifying the total binary star fraction of the stellar population between each realisation of the catalogue to yield a fraction of binary stars between 0 and 1.
When building images of stars and their companions, we adopt a spatial sampling 12 times better than the expected survey resolution. We simulate Euclidlike images that have the same resolution as the VIS camera, using Gaussian PSFs with a FWHM of 0.137′′. In doing so, we adopt a pixel size of 0.0083′′ and a filter bandpass that matches the VIS instrument of Euclid. Downsampling of the images to the same pixel size as the actual data is not necessary and has little impact on the results in the absence of noise. The results we present in the rest of the paper do not depend on the details of the shape of the PSF, i.e. Gaussians give the same results as more realistic diffractionlimited PSFs generated with the GalSim software (Rowe et al. 2015). This is mainly because the effect of unresolved stellar companions affects the very centre of the PSF.
When dealing with real (noisy) data, the PSF is best recovered by simultaneously considering many different point sources to improve the signal to noise. To mimic this process, we stack many images of stars (some containing binaries) generated as described above. In the present case, all images are (i) noisefree and (ii) centred on the measured centroid of the binary system or the centroid of the single star. Stars are only stacked within a narrow magnitude range. Binaries separated by 0.05′′are removed from the sample, i.e. we assume that in a real survey these binaries can be identified as such and removed from the stars used to build the PSF. Each image of a stack is composed of a mix of single and binary stars with realistic separations and contrasts. Measurements of the ellipticity and the size of the stacks are computed with adaptive moments by Hirata & Seljak (2003), as implemented in GalSim (Rowe et al. 2015).
A large number of stacks are measured and divided into bins of similar binary fractions. The biases and random errors on the reconstructed size and complex ellipticities are then estimated by comparing with the input values in the absence of binaries. For each group of stars with a similar binary fraction we measure the standard deviation for the sizes and ellipticities of all objects in the group. This gives an estimate of the precision achieved on a given parameter. The bias is the standard deviation of the difference between the ellipticity of the stack and the input values in the absence of binaries. In the rest of the paper we use the following terminology:

A stack of stars is an ensemble of stars that all have the same PSFshape and are drawn from the BMG using the binary properties ofDK13.

N_{⋆} is the number of stars in each stack; a binary system is counted only once.

N_{s} is the number of realisations of the stacks generated for each simulation. Each stack contains N_{⋆} stars drawn N_{s} times. N_{s} is always large enough to make the numerical errors on the simulation negligible in front the parameters that are reconstructed, i.e. the size and ellipticity of the PSF. Typical values for N_{s} are hundreds of thousands of realisations.

PSF reconstruction refers to the measurement of the PSF profile of the stack (size and ellipticity).

The error on the reconstructed parameters are defined in Eqs. (10) and (12).

The bias on the reconstructed parameters are defined in Eqs. (7)−(9).
5. Results
Fig. 4
PSF reconstruction errors (left column) and biases (right column). From top to bottom, plots for the size and two ellipticity components, e_{1} and e_{2}, are shown. In each panel, the values are given as a function of the binary fraction in the stellar population considered and for the three different PSF shapes described in Table 2. The solid lines represent the predicted values using the formalism in Sect. 3. The latter is valid for N_{⋆} → ∞ while the simulations used in this figure are for N_{⋆} = 40 stars with apparent magnitude in the range 18 ≤ i(AB) ≤ 19. 
In the following, we analyse the simulations described in Sect. 4.3. We focus on how the error in the PSF reconstruction due to multiple star blending varies with the intrinsic PSF ellipticity, i.e. in the absence of blends, and we show our results as a function of the total fraction of binary stars. We carry out the experiment both for Gaussian PSFs and for diffraction limited PSFs. We then explore the impact of the number of stars, N_{⋆}, available to reconstruct the PSF. Finally, we show the effect of changing the stellar populations, i.e. the magnitude distribution of the stars contained in each stack.
5.1. Varying the fraction of binaries for different PSF
We design a simple numerical experiment using three different Gaussian PSFs with the properties shown in Table 2. For each PSF, we generate ensembles of stars with different fractions of binaries and we produce noisefree stamp images on which we measure the ellipticity and size using a simple moment calculation.
Our results are summarised in Fig. 4, where we choose N_{⋆} = 40 as an example. In practice, N_{⋆} can change with different survey characteristics, in particular, the temporal and spatial stability of the PSF. In other words, N_{⋆} reflects the number of stars that can be used simultaneously to reconstruct the PSF.
Figure 4 indicates the reconstruction error for the size, and for the two ellipticity components, σ(e_{1}) and σ(e_{2}), as a function of the fraction of binary stars. We also show the biases on the size and on the ellipticity. In each panel, the figure compares the errors and biases, as derived from our simulations (symbols), to the predictions (solid lines) derived in Sect. 3.
Shape of the three PSFs used in the simulations.
The two upper panels of Fig. 4 show that the error and bias on the size as a function of the binary star fraction is independent of the PSF ellipticity, as suggested by Eq. (9). As one would expect, increasing the fraction of binary stars in the stacks increases the size of the PSF. The increase with the binary fraction is linear as the characteristics of the binary stars are the same in each bin. The only difference between the bins is the probability of observing a binary system. The figure also shows excellent agreement between the theoretical prediction and the simulation, even for the small number of stars in the stacks, here N_{⋆} = 40.
The reconstruction errors σ(e_{1}) and σ(e_{2}) on the complex ellipticity, as obtained from the simulations, appear very similar. They seemingly do not depend on the value of the intrinsic ellipticity as Eq. (11) suggests. However, with N_{⋆} = 40 stars, we fall in the suboptimal reconstruction scheme described in Sect. 3.3. Those cases are shotnoiselimited because of the small number of stars in the stack, i.e. the PSF reconstruction error of the ellipticity is dominated by noise due to the small number of random positions for the companions. We show in Sect. 5.3 that, as is intuitively expected, increasing N_{⋆} reconciles the errors found from the simulations and theoretical prediction: σ(e_{1}) and σ(e_{2}) do depend on the ellipticity of the PSF following Eq. (11) in the limit of large values for N_{⋆}.
We now turn to the biases due to stellar binarity. The bias on the size follows the same curve as a function of stellar binarity as the error on the size. The equations describing the bias and the error, as reported in Sect. 3, are (13)which is also well illustrated with our Gaussian simulations in the upper right panel of Fig. 4. The bias, both on the size and ellipticity, is easily measured from simulations with N_{⋆} = 40 in contrast to the measurement of the errors, which require much larger values for N_{⋆}.
The behaviour of the ellipticity is not the same as for the size. For a round PSF with e_{1} = e_{2} = 0, both the bias and reconstruction error equals zero in the limit where N_{⋆} → ∞, as the net effect of adding round PSFs at random angular positions around a star is null. The remaining bias for a round PSF, both on e_{1} and e_{2}, is of the order of 2 × 10^{6}. This can be attributed to the effect of a large but finite number of stacks N_{s} and to numerical errors.
For noncircular PSFs, the bias on e_{1} and e_{2} shows a linear correlation to the intrinsic complex ellipticity of the PSF, as does the error. The bias and reconstructions on the ellipticity have the same analytical form with opposite signs, and both are proportional to the intrinsic ellipticity of the PSF, i.e. The simulated measurements fall very close to the predictions in all panels of the right column of Fig. 4, i.e. the predictions are validated by simple Gaussian simulations.
The slight departure from linearity, for very large binary star fractions, is due to a changing value of ⟨ r^{2}f/ (1 + f)^{2} ⟩, i.e. for large binary fractions the stellar population is modified with respect to populations with a lower fraction of binaries.
The simple linear relation between the bias and intrinsic ellipticity of the PSF hints at a possible calibration to reduce the effect of binary stars on the PSF shape or, at the very least, on a reliable control of the bias. However, that the lack of knowledge of the distribution of the binary star fraction (⟨ r^{2}f/ (1 + f)^{2} ⟩) in each specific field used to reconstruct the PSF may complicate the task in practice.
5.2. Synthetic spacebased PSF
Up to now, we have only considered Gaussian PSFs. In order to test the effect of a more complex PSF, we use the GalSim software to produce a synthetic PSF similar to future spacebased telescopes. More specifically, our PSF is the same size as the Euclid VIS PSF and has six struts. No PSF aberration such as coma or trefoil is applied. The strut angle, with an arbitrarily chosen orientation, is kept constant.
Similar to the simulations in Sect. 5.1, we build three PSFs with the characteristics shown in Table 2. We choose the strutrotated validation case to have the same e_{1} and e_{2} input as the PSF along e_{1}.
The results show the same behaviour as for the Gaussian toy model. No significant departure from the predictions are detected. This is likely because we consider here very narrow angular separation binaries, with a companion star well within the central pixel of the PSF. It is therefore not surprising that the main change in shape, due to the binarity, occurs in the very centre of the PSF. The details of the PSF shape on larger scales do not matter much. This suggests that (i) simple Gaussian models are sufficient to assess the impact of binarity of the PSF shape; and (ii) our analytical predictions based on Gaussian PSFs are a good approximation of the change in the PSF shape. These predictions can therefore be safely used to predict the PSF reconstruction error as a function of the binary fraction and mean angular separation between the stars and their companions.
Fig. 5
PSF reconstruction errors on the two ellipticity components, e_{1} (left) and e_{2} (right). These plots correspond to the two panels in the lower left corner of Fig. 4, but we now show the effect of changing the number of stars, N_{⋆}, used in the stack. Depending on the intrinsic ellipticity of the PSF (see Table 2), thousands of stars are needed to avoid being dominated by the shot noise. 
5.3. Effect of the number of stars N_{⋆} in the stack
As seen in Sect. 5.1 and illustrated in Fig. 4, varying the number of stars in the stack N_{⋆} used to reconstruct the PSF changes the reconstruction error. Increasing the number of stars increases the sampling of the coordinates for the companion stars. Eventually, for very large values of N_{⋆}, the predicted errors as computed in Sect. 3 are realised by the simulations.
For N_{⋆} = 40 stars, the reconstruction error on the size, as well as the values for the biases on the size and ellipticity, are compatible with the predictions, as in Fig. 4. However, the errors on the ellipticity components require larger values for N_{⋆}, which also depends on the intrinsic ellipticity of the PSF.
This is tested in Fig. 5 for N_{⋆} = 40,750,2100,10 000 and for two PSF ellipticities. For low values of N_{⋆}, the error is large and independent of the ellipticity of the PSF and the error is dominated by the shot noise due to the number of companions. As can be expected, it takes a much larger N_{⋆} to converge to the theoretical predictions, and the smaller the ellipticity, the larger N_{⋆} needed. For example, N_{⋆} = 10 000 stars in the stack are barely sufficient to converge to the theoretical estimate of the error for PSF2, which has a small e_{1} component of −0.021.
Figure 5 shows that when N_{⋆} is very large, the theoretical values for the error estimates on the ellipticity can be modelled by simulating binary populations and the equations in Sect. 3. However, the figure also shows that if for any reason the prediction we make in Sect. 3 with Gaussian PSF does not hold true, then the simulations required to replace the theoretical estimates involve tens of thousands of companions. Given that we also need large numbers of realisations N_{s}, estimating the effect of binary stars on the PSF reconstruction may quickly become computationally expensive.
5.4. Dependence on the apparent magnitude of the stars
Fig. 6
Theoretical reconstruction error according to the formalism developed in Sect. 3. The stellar distribution is drawn from realisations of the BMG for pointing coordinates to ℓ = 180°,b = 45°. Top: value of the ⟨ r^{2}f/ (1 + f)^{2} ⟩ parameter as a function of apparent i(AB) magnitude. Middle: reconstruction error as a function of the apparent magnitude i(AB) of the stars in the stack and as a function of the complex ellipticity components e_{1} or e_{2} (colour code). Bottom: reconstruction error on the PSF size. The red line indicates typical requirements for stage IV experiments. 
Our ability to resolve a binary star system depends mostly on its physical distance to the Earth, i.e. on its apparent magnitude. The PSF reconstruction errors due to the blend of binaries are correlated with the distance to the binary system through the ⟨ r^{2}f/ (1 + f)^{2} ⟩ parameter. We now turn to the impact of stars in different magnitude bins on the reconstruction of the PSF. In doing so, we adopt a binary star fraction of 35%, which is representative of the mean binary star fraction on the sky (DK13). The resulting reconstruction errors as a function of the i(AB) magnitude are shown on Fig. 6. In the figure, we also indicate, as a reference, the total reconstruction error budget for a stage IV weak lensing experiment: σ(e_{i}) ≤ 2 × 10^{4} for the ellipticity components and σ(R^{2})/ ⟨ R^{2} ⟩ ≤ 1 × 10^{3} for the size (PaulinHenriksson et al. 2008). These values are the total error budgets allowed for these parameters.
For the brightest stars with 16 <i(AB) < 18, the error in reconstruction takes the full error budget both on the size and ellipticity, even for small intrinsic PSF ellipticities. Stars in the magnitude range 18 <i(AB) < 19 still suffer from a relatively large blending bias, taking up to 60% of the total error budget on the size and ellipticity. Requiring stars fainter than magnitude 20 lowers the blending bias to ≲ 10% of the error budget, which is still large enough to be accounted for.
Stars brighter than roughly i(AB) ≃ 19 may not be the best choice for PSF calibration strategies, as brighter stars are much more likely to alter the PSF both in size and ellipticity. Instead, longer observations of faint and most distant stars are preferred and minimise the effect of binary stars. Finally, we emphasise that the impact of binary stars on the PSF shape becomes more important with higher spatial resolution if the stellar population remains the same. Surveys like WFIRST may therefore be more impacted by binary stars than the Euclidlike telescope used in this work. On the other hand, because WFIRST has deeper limiting magnitudes, more faint stars can potentially be used to reconstruct the PSF.
6. Conclusions
Multiple stellar systems are ubiquitous in the Milky Way, with a total fraction of multiple stars of about 35% (DK13). These multiple stars, especially when they are not resolved by the instrument, may significantly impact the PSF measurement given the strong requirements on the modelling of the PSF.
In this paper, we present analytical predictions for the PSF reconstruction error introduced by binary stars based on stellar multiplicity parameters and on the intrinsic PSF size and ellipticity. We then verify our predictions by means of realistic numerical simulations. To do so, we use the BMG (Robin et al. 2003) to generate mock stellar distributions. We then add companion stars to produce binary systems following the most recent observational data in DK13. We produce noisefree images used to measure the error and bias on the PSF ellipticity and size.
In our simulations, the parameters of the detector and of the PSF are taken to be Euclidlike, but the PSF profile is Gaussian. We also check that our results are valid for more realistic diffraction limited PSFs. Binary stars with a separation larger than r = 0.05′′ are assumed to be easily identified as multiple and eliminated from the star catalogue to build the PSF. In other words, we only consider double stars that cannot be identified as such but that statistically affect the PSF reconstruction. Our main findings can be summarised as follows.

Binary stars can significantly alter the reconstruction of the PSF.They introduce both a bias and an error on the size and ellipticity ofthe reconstructed PSF.

The PSF binary reconstruction errors can be predicted analytically as a function of characteristics of the stellar population considered, i.e. the separations and intensity ratio of the stars, and as a function of the intrinsic PSF size and ellipticity.

The effect of the binaries cannot be suppressed by averaging a large number of stars.

The analytical predictions are supported by numerical experiments. The PSF used in the simulations are Gaussian and diffractionlimited.

Bright stars of magnitudes i(AB) ≲ 19 are most affected by binarity. These stars are either massive stars with very likely companions or they are nearby stars, with potentially large angular separations to their companion(s).

Fainter stars in the range 18 < i(AB) < 19) may cause significant biases on the PSF shape. If no mitigation scheme is found, these stars can take on their own up to ~ 60% of the total budgeted error on the size and ellipticity described in PaulinHenriksson et al. (2008) and Cropper et al. (2013).

Blending errors in systems fainter 20 < i(AB) < 21 can still contribute up to 10% to the typical error budget on the size and ellipticity.

If any PSF calibration fields are chosen for space missions, targeting fields of faint stars can be preferable to bright stars to avoid nearby binary systems that have larger effects on the PSF shape.

The total fraction of binaries varies with Galactic longitude. There are hints that the binary star fraction increases with Galactic latitude. Thus, the variations in the underlying stellar populations can have a varying effect on the PSF shape as a function of Galactic coordinates. This may also imply a bias on the shear power spectrum that must be controlled.

Different stellar types differ in their fraction of multiple systems and also in the distribution of the parameters of the orbits. Regions of the sky with lower binaries fractions could be exploited to quantify the effect of binary stars on the PSF determination differentially across the sky.
Given the above conclusions, it will be important that future spacebased weak lensing surveys account for the effect of binary stars on the PSF determination. In particular, their contribution to the error budget on the PSF shape (both the error and the bias) must be estimated and any PSF reconstruction method must account for the extra error caused by binarity. The above conclusions depend on the ability to clean star samples from multiple systems. This can be done in many ways and with different efficiency depending on signaltonoise ratio and stellar type. Gaia may provide the spectral type of the stars and, in some cases, identify binary stars from their apparent motion on the plane of the sky. Alternatively, specific image processing techniques using the Euclid images themselves may be used, but these remain to be devised.
Acknowledgments
We are indebted to the anonymous referee whose valuable comments increased the quality of this work. The authors would like to thank Laurent Eyer, David Harvey, Henk Hoekstra, Matthew Nichols, Pierre North, and Thomas Kitching for useful discussions. This work is supported by the Swiss National Science Foundation (SNSF).
References
 Cropper, M., Hoekstra, H., Kitching, T., et al. 2013, MNRAS, 431, 3103 [NASA ADS] [CrossRef] [Google Scholar]
 de Bruijne, J. H. J., Allen, M., Azaz, S., et al. 2015, A&A, 576, A74 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Duchêne, G., & Kraus, A. 2013, ARA&A, 51, 269 [Google Scholar]
 Duquennoy, A., & Mayor, M. 1991, A&A, 248, 485 [NASA ADS] [Google Scholar]
 Frieman, J. A., Turner, M. S., & Huterer, D. 2008, ARA&A, 46, 385 [NASA ADS] [CrossRef] [Google Scholar]
 Habets, G. M. H. J., & Heintze, J. R. W. 1981, A&AS, 46, 193 [NASA ADS] [Google Scholar]
 Heacox, W. D. 1996, PASP, 108, 591 [NASA ADS] [CrossRef] [Google Scholar]
 Hirata, C., & Seljak, U. 2003, MNRAS, 343, 459 [NASA ADS] [CrossRef] [Google Scholar]
 Hoogenboom, J. P., den Otter, W. K., & Offerhaus, H. L. 2006, J. Chem. Phys., 125, 204713 [NASA ADS] [CrossRef] [Google Scholar]
 Kaiser, N., Squires, G., & Broadhurst, T. 1995, ApJ, 449, 460 [NASA ADS] [CrossRef] [Google Scholar]
 Kilbinger, M. 2015, Rep. Prog. Phys., 78, 086901 [Google Scholar]
 Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, Euclid Study Report, ArXiv eprints [arXiv:1110.3193] [Google Scholar]
 Massey, R., Hoekstra, H., Kitching, T., et al. 2013, MNRAS, 429, 661 [Google Scholar]
 PaulinHenriksson, S., Amara, A., Voigt, L., Refregier, A., & Bridle, S. L. 2008, A&A, 484, 67 [Google Scholar]
 Reipurth, B., Clarke, C. J., Boss, A. P., et al. 2014, Protostars and Planets VI, 267 [Google Scholar]
 Robin, A. C., Reylé, C., Derrière, S., & Picaud, S. 2003, A&A, 409, 523 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Robin, A. C., Luri, X., Reylé, C., et al. 2012, A&A, 543, A100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Rowe, B., Jarvis, M., Mandelbaum, R., et al. 2015, Astron. Comput., 10, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Spergel, D., Gehrels, N., Baltay, C., et al. 2015, ArXiv eprints [arXiv:1503.03757] [Google Scholar]
 Weinberg, D. H., Mortonson, M. J., Eisenstein, D. J., et al. 2013, Phys. Rep., 530, 87 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
All Tables
All Figures
Fig. 1
Reconstruction errors for a single binary star with the main component centred on (0,0). Each point on the map shows the measurement of the quantity labelled in each panel when a single companion star is located at that position. The labels are in arcsecond. The first column on the left shows the effect for a PSF aligned along e_{1} and a contrast of Δm = 0 between the main star and its companion. The next column is for the same PSF shape and orientation, but for a contrast of Δm = 1. The block of six graphs on the right is for the same two contrasts, but for a PSF aligned along e_{2}. The red line represents the orientation of the PSF. From top to bottom, the total ellipticity (star + companion) of the PSF, the relative error on the ellipticity due to the companion star, and the relative error on the size are shown. 

In the text 
Fig. 2
Sketch of the system with the main object at , with total intensity I_{m} and the companion located (δx,δy) away with a total intensity of I_{m}f = I_{c}. 

In the text 
Fig. 3
Mainsequence stellar population breakdown for different Galactic latitudes. The Galactic longitude on the left panel is ℓ = 0° (towards the centre of the disk) and the right panel ℓ = 180° (Galactic anticentre). The bin of massive stars (>1.5 M_{⊙}) is not visible as it represent less than 0.2%. We stress the different behaviour of the binary fraction for massive stars (blue vs. green bins) for the opposite directions in Galactic longitude. 

In the text 
Fig. 4
PSF reconstruction errors (left column) and biases (right column). From top to bottom, plots for the size and two ellipticity components, e_{1} and e_{2}, are shown. In each panel, the values are given as a function of the binary fraction in the stellar population considered and for the three different PSF shapes described in Table 2. The solid lines represent the predicted values using the formalism in Sect. 3. The latter is valid for N_{⋆} → ∞ while the simulations used in this figure are for N_{⋆} = 40 stars with apparent magnitude in the range 18 ≤ i(AB) ≤ 19. 

In the text 
Fig. 5
PSF reconstruction errors on the two ellipticity components, e_{1} (left) and e_{2} (right). These plots correspond to the two panels in the lower left corner of Fig. 4, but we now show the effect of changing the number of stars, N_{⋆}, used in the stack. Depending on the intrinsic ellipticity of the PSF (see Table 2), thousands of stars are needed to avoid being dominated by the shot noise. 

In the text 
Fig. 6
Theoretical reconstruction error according to the formalism developed in Sect. 3. The stellar distribution is drawn from realisations of the BMG for pointing coordinates to ℓ = 180°,b = 45°. Top: value of the ⟨ r^{2}f/ (1 + f)^{2} ⟩ parameter as a function of apparent i(AB) magnitude. Middle: reconstruction error as a function of the apparent magnitude i(AB) of the stars in the stack and as a function of the complex ellipticity components e_{1} or e_{2} (colour code). Bottom: reconstruction error on the PSF size. The red line indicates typical requirements for stage IV experiments. 

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