© ESO 2020

1. Introduction

Over the past century advances in observational techniques in cosmology have led to a number of important discoveries of which the accelerating expansion of the Universe is perhaps the most surprising. Moreover, a wide range of detailed observations can be described with a model that requires a remarkably small number of parameters, which have been constrained with a precision that was unimaginable only 30 years ago. This concordance model, however, relies on two dominant ingredients of the mass-energy content of the Universe: dark matter and dark energy. Neither of these ingredients can be described satisfactorily by our current theories of particle physics and gravity. Although a cosmological constant/vacuum energy is an excellent fit to the current data, the measured value appears to be unnaturally small. Many alternative explanations have been explored, including modifications of the theory of general relativity on large scales (see e.g. Amendola et al. 2013, for a review), but a more definitive solution may require observational constraints that are at least an order of magnitude more precise.

The concordance model can be tested by studying the expansion history of the Universe and by determining the rate at which structures grow during this expansion. This is the main objective of the Euclid mission (Laureijs et al. 2011), which will carry out a survey of 15 000 deg² of the extragalactic sky. Although Euclid will enable a wide range of science topics, it is designed with two main probes in mind: first, the measurement of the clustering of galaxies at z >  0.9 using near-infrared, slitless spectroscopy; and second, the direct measurement of the distribution of matter as a function of redshift using weak gravitational lensing, the effect whereby coherent shear distortions in the images of distant galaxies are caused by the differential deflection of light by intervening large-scale structures. The two-point statistics of the weak gravitational lensing caused by large-scale structure is known as cosmic shear (see Kilbinger 2015, for a recent review). In this paper we explore the impact of instrumental effects and scanning strategy on the accuracy and precision with which dark energy parameters w₀ and w_a (Chevallier & Polarski 2001; Linder 2003) can be measured using the cosmic shear from Euclid.

The challenge of measuring the cosmic shear signal is that the typical change in polarisation, i.e. third flattening or eccentricity, caused by gravitational lensing is approximately one percent. This value is much smaller than the intrinsic, unlensed ellipticities of galaxies. To overcome this source of statistical uncertainty, cosmic shear is measured by averaging over large numbers of galaxies pairs. For the result to be meaningful, sources of bias caused by systematic effects need to be sub-dominant. Systematic effects can be mitigated through instrument design, but some need to be modelled and removed from the data. In order to determine how such systematic effects can bias cosmic shear measurements – and cosmological parameter inference – a series of papers derived analytic expressions that represented the measurement and modelling processes involved. Following an initial study by Paulin-Henriksson et al. (2008) that focussed on point spread function (PSF) requirements, Massey et al. (2013, M13 hereafter) presented a more general analytic framework that captures how various systematic effects affect the measurements of galaxy shapes. This study provided the basis for a detailed breakdown of systematic effects for Euclid by Cropper et al. (2013, C13 hereafter), which has been used in turn to derive requirements on the performance of algorithms and supporting data. Another approach, based on Monte Carlo Control Loops (MCCL), has also been presented (Bruderer et al. 2018; Refregier & Amara 2014), which uses a forward modelling approach to calibrate the shear measurement.

Although these previous studies provide a convenient way to compare the impact of various sources of bias, their analytic nature means that particular assumptions are made, and they cannot capture the full realism of a cosmic shear survey. Therefore we revisit the issue in this paper for a number of reasons:

(1) To avoid an implicit preference for implementation, the derivations in M13 are scale-independent, i.e. they do not explicitly depend on angle θ or multipole ℓ. In more realistic scenarios, such as those we consider in this work, spurious signals are introduced on specific spatial and angular scales on the celestial sphere. For example, the PSF model is determined from the full instrument field of view, whereas detector effects, such as charge transfer inefficiency (CTI) occur on the scale of the region served by a single readout register on a CCD. In addition, the biases may depend on the observing strategy or time since launch. This is particularly true for CTI, which is exacerbated by radiation damage, and thus increases with time (Massey et al. 2014; Israel et al. 2015). An initial study of the implications of scale-dependent scenarios was presented in Kitching et al. (2016), who find that survey strategy can play a critical role in the case of time-dependent effects. Their results suggest the expected biases in cosmological parameters may be reduced if the correct scale dependences are considered.

(2) The residual systematic effects may depend on the region of the sky that is observed. For example, the model of the PSF can be constrained to a higher precision when the density of stars is higher. On the other, hand these effects may also have an adverse effect on the galaxy shape measurement of the shear (Hoekstra et al. 2017). The impact of CTI depends on the sky background level, and thus is a function of ecliptic latitude, whereas Galactic extinction may introduce biases in the determination of photometric redshift that depend on Galactic latitude (and longitude). These subtle variations across the survey should be properly accounted for, and their impact on the main science objectives of Euclid evaluated.

(3) In the analytic results of works such as C13, a distinction was made between convolutive (i.e. caused by PSF) and non-convolutive contributions. The impact of the former, such as the PSF, is relatively easy to propagate because it is typically clear how these contributions depend on galaxy properties. The latter, however, which include biases introduced by CTI, are more complicated to capture because their dependence on galaxy properties such as size and flux can be non-linear. Moreover, the allocations implicitly assume that residual errors are independent because correlations between effects could not be easily included. Hence, the impact of a more realistic error propagation needs to be examined.

(4) The interpretation of the requirements presented in C13 is unclear, in particular whether they should be considered as values that are never to be exceeded, the mean of a distribution of possible biases, or upper limits corresponding to a certain confidence limit. As shown below, we expect our limited knowledge of the system to result in probability density distributions of biases that should be consistently combined to evaluate the overall performance.

(5) Finally, in Kitching et al. (2019) we show that these previous studies made simplifying assumptions with regard to the analytic relationship between position-dependent biases and the cosmic shear statistics, where the correct expression involves second and third order terms. This motivates our study in two ways. Firstly the correct expression involves previously unstudied terms. Secondly, the correct expression is computationally demanding, meaning its calculation is intractable for realistic cosmic shear measurements.

In this paper we present a general framework for investigating systematic effects that addresses all these issues, but does not require full image-level, end-to-end simulations, which would require fully realistic mock data and data processing stages. Instead our approach starts at the object catalogue level and systematic effects are propagated through a chain of processes on an object-by-object basis. This does not mean that systematic effects are not in common between galaxies, but it assumes that the measurement process is. This is a reasonable assumption for weak lensing studies in which the shape measurement itself is confined to a narrow angular region about the vicinity of the galaxy on the sky. This allows us to create scenarios where systematic effects are calculated in a more realistic fashion, starting from a catalogue of sources with appropriate parameters, and propagated all the way to the evaluation of cosmological parameters. This approach may not capture all correlations between systematic effects because this can only be achieved through a full end-to-end simulation of the pipeline; however, it does present a major advance over the initial studies presented in M13 and C13. The pipeline as presented also stops a number of steps short of realistic parameter inference. This was an explicit design choice such that we start from a very well defined and minimal-complexity modular baseline that then allows us to expand this on a modular-by-modular basis.

We describe the general framework in more detail in Sect. 2, where we also discuss the properties of the input catalogue, sky parameters, and observational characteristics. Results are presented in Sect. 3. A more complete exploration of the many possible sources of bias for Euclid is deferred to future work, but in Sect. 4 we consider a few case studies: in Sect. 4.1 the residuals in the PSF correction and in Sect. 4.2 the impact of imperfections in the correction for CTI. Although the performance analysis in this paper takes Euclid as a reference mission, the framework is sufficiently general that it can be applied to any future Stage IV weak gravitational lensing survey (e.g. Large Synoptic Survey Telescope Science Collaboration 2009; Spergel et al. 2015).

2. General framework

The general framework we present is a causally connected pipeline or transfer function-like methodology. This pipeline modifies the values of quantities associated with each individual galaxy according to the effects of the instrument and measurement processes. These are in turn used to compute cosmic shear power spectra to evaluate the impact on cosmological parameter inference. The general framework is captured in Fig. 1, which we summarise in Sect. 2.6.

Fig. 1. Overall structure of the concept as described in the main text. The quadrupole moments Q are initiated with intrinsic moments and then modified by incorporating the shear, PSF, and detector effects. Survey characteristics such as dither pattern, slew pattern, and observation time are entered in the initial catalogue. A measurement process M subsequently converts the observed moments to polarisations. The estimation of the galaxy polarisation is then made (as described in Eqs. (10) and (11)). This is done per object. Next a power spectrum for the reference and perturbed scenarios is computed. For the perturbed line the PSF and detector moments are drawn from distributions that represent the measurement uncertainty as described in the text. This process is repeated for 150 random realisations for the set of galaxies that are in the input catalogue. Finally the residual power spectrum is computed per realisation, and the statistics of each of the realisations is passed onto the Fisher matrix, from which uncertainties and biases of dark energy parameters are calculated. White circles indicate moment space, where modifications are performed on an object-by-object basis. Grey circles indicate ensemble average in the harmonic space. Diamonds show cosmological parameter space.

2.1. Causally connected pipeline

As light propagates from a galaxy, several processes occur that act to transform a galaxy image. We represent this as a series of sequential processes, or a pipeline, which are causally ordered, i.e.

(1)

where I is a surface brightness and the subscripts i refer to an object (a galaxy in our case); the other subscripts refer to the addition of an effect (labelled as a +), where in this example “shr” denotes shear, “PSF” denotes the PSF, det denotes the detector, and so on. The pipeline is initiated by a projected initial (intrinsic) surface brightness distribution for object i that is modified/transformed via a series of processes, that is, shearing by large-scale structure and convolution by the PSF, which depend on the preceding step. The last step M represents a measurement process that converts the observed surface brightness distribution into quantities that can be used for science analyses. Equation (1) is an example, which includes shear and PSF effects, of a more general framework that we define here in this equation

(2)

where α is some general process that modifies the surface brightness distribution of object i that precedes process α + 1, and so forth. In this paper we only focus on the impact of PSF and detector effects on cosmic shear analyses, but emphasise that the approach is much more general. This approach can be readily extended to include more effects, such as photometric errors, spectral energy distribution (SED) dependent effects, or the impact of masking. These will be explored in future work.

The objects in question for weak lensing measurements are stars, which are used for PSF determination, and galaxies. The primary quantities of interest for these galaxies are the quadrupole moments of their images, which can be combined to estimate polarisations and sizes. The unweighted quadrupole moments Q_i, mn of a projected surface brightness distribution (or image) I_i(x) are defined as

(3)

where F is the total observed flux, m and n are (1, 2) corresponding to orthogonal directions in the image plane, and we assumed that the image is centred on the location where the unweighted dipole moments vanish. We can combine the quadrupole moments to obtain an estimate of the size , and shape of a galaxy through the complex polarisation, or third eccentricity¹

(4)

Therefore, the pipeline process for the cosmic shear case is similar to that given by Eq. (1), but for the quadrupole moments of the surface brightness distribution. In this case each process acting on the surface brightness distribution is replaced by its equivalent process acting on the quadrupole distribution, and the final measurement process is the conversion of quadrupole moments into polarisation,

(5)

where we suppress the mn subscripts for clarity. In this expression χ_obs, i is the observed polarisation for object i that is a function of , where these quantities are related by Eq. (4) in the general case. The result is then used for cosmic shear analysis. Importantly, at each stage in the pipeline, the relevant quantities that encode the intrinsic effects of the ellipticity, shear, PSF, and detector, instead of being fixed for all objects, can be drawn from distributions or functions that capture the potential variation owing to noise in the system and the natural variation of object and instrumental properties.

2.2. Reference and perturbed scenarios

Next we introduce the concept of a reference scenario, representing the ideal case, and a perturbed scenario, which results in biased estimates caused by mis-estimation and uncertainty in the inferred values of the quantities that are included in the set of causally linked processes as described in Eq. (1). We define these below.

Reference. In this scenario the systematic effects that have been included in the pipeline are perfectly known, so that in the final measurement process their impact can be fully accounted for and reversed. In this case the distribution of parameter values that are used to undo the biases are all delta functions centred on the reference values, that is there is no uncertainty in the system.

Perturbed. In this scenario systematic effects that have been included in the pipeline are not known perfectly. As a consequence the corrections result in biased measurements. In this case relevant quantities that are used to undo the systematic effects are drawn from probability distributions that represent the expected level of uncertainty.

We can then define the elements in a pipeline for each scenario. The difference between the observed reference polarisation for a given object, and the observed perturbed polarisation is a realisation of expected polarisation uncertainty caused by a semi-realistic treatment of systematic effects in a data reduction scenario. We explain this further using the specific example with which we are concerned in this paper: the assessment of cosmic shear performance.

In our case, the output of the pipeline process, Eq. (5), leads to a set of measured polarisations and sizes that represent the true response of the system, that is an ellipticity catalogue that includes the cumulative effects of the individual processes as they would have occurred in the real instrument and survey. As detailed in M13, we can compute how PSF and detector effects change the polarisation and size of a galaxy² as follows:

(6)

where χ_obs, i is the observed polarisation, χ_ini, i is the intrinsic/unlensed polarisation, χ_shr, i is the induced polarisation caused by the applied shear γ, χ_PSF, i is the polarisation of the PSF, and χ_det, i is the detector-induced polarisation; the same subscripts apply to the R² terms (, see Eq. (3)). The relation between the applied shear, γ, and the corresponding change in polarisation, χ_shr, is quantified by the shear polarisability P^γ so that

(7)

(Kaiser et al. 1995). The shear polarisability depends on the galaxy morphology, but it can be approximated by the identity tensor times a real scalar (where I is the identity matrix) in the case of unweighted moments (Rhodes et al. 2000). We simplify this equation, in terms of notation, to

(8)

where χ_gal, i = χ_ini, i + χ_shr, i (the polarisation that would be observed given no PSF or detector effects), and

(9)

These quantities are constructed from the corresponding quadrupole moments in Eq. (5).

Given a set of observed galaxy polarisations and sizes and perfect knowledge of the systematic effects Eq. (8) can be inverted, yielding an estimate for the galaxy shape in the reference case given by

(10)

where the superscript R denotes the reference case. In this case the quantities , , and are known exactly and constructed from the quadrupole moments in Eq. (5), and we obtain (trivially) the underlying true . Even though this is a trivial inversion we nevertheless perform this step since in general the measurement process may not be exactly invertable.

In the perturbed case, the uncertainties in the measurement and modelling process result in a set of estimated values that include residual effects of the PSF and detector

(11)

where the superscript P denotes the perturbed case. In this equation, , , and are constructed from the quadrupole moments drawn from relevant probability distributions that represent uncertainties in the system. The resulting polarisation estimates correspond to a realisation of the system that encodes the expected uncertainty in our understanding of PSF and detector effects. Each of these steps is then repeated for realisations of the probability distributions present in the perturbed quantities. The implementation of these probability distributions for the PSF and CTI cases are detailed in Appendices A and B.

To convert the estimated reference and perturbed polarisations to their corresponding shear estimates we use

(12)

which provides a noisy, but unbiased estimate of the shear γ (M13). We note that P^γ does not change between the reference and perturbed cases. In practice, shape measurement algorithms use weighted moments to suppress the noise in the images, which changes the shear polarisation compared to the unweighted case. The correction for the change in shape caused by the weight function depends on the higher-order moments of the surface brightness (Melchior et al. 2011) and is a source of shape measurement bias that can be quantified using image simulations (e.g. Hoekstra et al. 2017). This also leads to a sensitivity to spatial variations in the colours of galaxies if the PSF is chromatic (Semboloni et al. 2013; Er et al. 2018). However, for this study this complication can be ignored as we implicitly assume that the biases in the shape measurement algorithm have been accounted for to the required level of accuracy (C13).

In future work, we will include more effects in the perturbed scenario. Observable quantities can be generalised to a function of redshift and wavelength, i.e. . We will then explore the effects of masking, shape measurement errors, photometric errors, and SED variations within a galaxy.

2.3. Shear power spectrum estimation

The estimated polarisations contain

(13)

where is the change in polarisation. We note that χ_gal, i is in general redshift-dependent, leading to tomographic power spectra, but that in this initial proof of concept we do not include redshift-dependent effects in ; the consequence of this is shown in Eq. (18). We assume higher-order terms are subdominant, i.e. terms involving for n >  1. The quantity is caused by the uncertainty in systematic effects, that is defined by expanding the denominator in Eq. (11) to linear order, and substituting Eq. (8) as follows:

(14)

where the denominator in Eq. (11) is expanded by assuming .

The polarisations in Eq. (13) can be converted to estimates of the corresponding shears using Eq. (7) and Eq. (12), and . These can be subsequently used to calculate shear power spectra, and the residual between the reference and perturbed spectra,

(15)

where

(16)

where are the spherical harmonic coefficients of the perturbed shear field, i.e.

(17)

In the above expressions (θ_i, ϕ_i) is the angular coordinate of galaxy i, the ₂Y_ℓm(θ_i, ϕ_i) are the spin-weighted spherical harmonic functions, and a ^* refers to a complex conjugate. Similarly for the reference case . is a realisation n of one that may be observed given the limited knowledge of uncertainties in systematic effects. We can split the residual power spectrum into three terms: δ − δ quantifies the auto-correlation of the systematic uncertainties and gal-δ and δ-gal are the cross-correlation power spectra between the systematic uncertainties and the true cosmological signal (i.e. the signal that would have been observed if all systematic effects were perfectly accounted for).

Although selection effects can result in a correlation between the shear and systematic effects, we stress that we are interested in residual effects, and thus implicitly assume that such selection effects have been adequately accounted for. Hence, when taking an ensemble average over many realisations, we are left with ⟨δC_n(ℓ)⟩ = C^{δ − δ}(ℓ) as the mean of these additional terms should reduce to zero and any variation is captured in the error distribution of the ensemble of δC(ℓ). Hence we can determine the power spectrum caused by uncertainties in systematic effects.

We sample from all parameter probability distributions in the perturbed case, and compute the mean and variance over the resulting ensemble of {δC_n(ℓ)}. In the cases in which random numbers are required for the reference case, care must be taken to ensure that the seed is the same in the reference and perturbed cases.

In this initial proof of concept we do not investigate redshift-dependent systematic effects, that is the change in polarisations is applied to all galaxies regardless of their redshifts. This then means that we can assume that PSF and CTI effects have an equal impact on all tomographically binned cosmic shear power spectra, and we can generalise the discussion above such that, for example

(18)

for all redshift bin labels α and β.

2.4. Comparison to previous work

To compare to previous work, in M13 generic non-parametric realisations of δC(ℓ) are generated and used to place conservative limits on a multiplicative and additive fit to such realisations δC(ℓ) = ℳ C^R(ℓ)+𝒜, where ℳ and 𝒜 are constant so that biases in the dark energy parameters, using Fisher matrix predictions, were below an acceptable value. This represents a worst case because the residual power spectra are assumed to be proportional to the cosmological signal (apart from the additive offset). In Kitching et al. (2016), simple models for systematic effects are used to create simplified but realistic δC(ℓ) values. In Taylor & Kitching (2018) the constant multiplicative and additive formulation is generalised to include the propagation of real-space multiplicative effects into power spectra as a convolution. The full expression for the analytic propagation of constant and scale-dependent multiplicative and additive biases is derived in Kitching et al. (2019). This reveals that the analytic propagation of biases into cosmic shear power spectra involves second- and third-order terms that result in an intractable calculation for high-ℓ modes. Our approach, therefore, differs from the earlier works in that it captures any general scale and redshift dependence on an object-by-object level, and, very importantly, creates δC(ℓ) values that correctly incorporate the uncertainty in the system. This procedure enables a complete evaluation of the performance, which differs from a true end-to-end evaluation only in that we do not use the images and image-analysis algorithms that are used to analyse the real data.

These catalogue-level simulations have the major advantage that they are much faster than full end-to-end image simulations, allowing for realisations of systematic effects to be computed so that a full probability distribution of the effect on the cosmological performance of the experiment can be determined. This allows us to explore various survey strategies and other trade-off considerations, whilst capturing most of the complexities of the full image-based analysis. The catalogue-level simulations include survey-specific features, such as the detector layout, survey tiling, and PSF pattern (see Sect. 3.2). It also allows for foreground sky models to be included to account for variations in Galactic extinction, star density, and Zodiacal background. Calibration uncertainties can be incorporated by adjusting the probability density distributions of the relevant parameters accordingly.

2.5. Propagation to cosmological parameter estimation

To assess the impact of the power spectrum residuals on cosmological parameter inference, we used the Fisher matrix Euclid Collaboration (2019a), and bias (Kitching et al. 2008; Amara & Réfrégier 2008; Taylor & Kitching 2018) formalism. We used the w0waCDM Fisher matrix from Euclid Collaboration (2019a); all code and files can be found on the associated repository for that paper³.

We very briefly summarise the Fisher matrix and bias formalism, which is based on the notation in Euclid Collaboration (2019a). In general, a change in the power spectrum caused by a residual systematic effect can influence the size of the confidence region about any parameter as well as the maximum likelihood location. In this paper we only consider the change in the maximum likelihood position.

The expected confidence regions for the cosmological parameters can be expressed using the Fisher matrix, which is given by

(19)

where (α, β) and (ω, ρ) are redshift bin pairs and (μ, ν) denote cosmological parameter pairs and (ℓ_min, ℓ_max) are the minimum and maximum angular wavenumbers used. The are the true cosmic shear power spectra; we note that due to sample variance. The covariance is given by

(20)

where f_sky is the fraction of the sky observed. We note that we assume a Gaussian covariance in this case and do not include non-Gaussian terms. The Euclid Collaboration (2019a) find that the signal-to-noise ratio of the power spectra decreases by ∼30% at ℓ_max = 5000 when the non-Gaussian contributions are included and that this loss of information content corresponds to an effective cut at ℓ_max = 1420 in a forecast that only uses the Gaussian covariance given by Eq. (20).

The noise power spectrum is defined as , where N_g, β is the total number of galaxies in bin β for full sky observation and δ_αβ is a Kronecker delta. The intrinsic shape noise is quantified by σ_χini = 0.3, the dispersion per ellipticity component. This can be used to compute the expected marginalised, cosmological parameter uncertainties σ_μ = [(F⁻¹)_μν]^1/2.

The changes in the maximum likelihood locations of the cosmology parameters (i.e. biases) caused by a change in the power spectrum can also be computed for parameter α as

(21)

where the vector B for each parameter β is given by

(22)

We note that the biases computed in this section are the one-parameter, marginalised biases and that this may result in optimistic assessments for multi-dimensional parameter constraints. The quantity n refers to the realisation number as previously discussed and we note that δC_n(ℓ) is not redshift dependent in this case. A multi-dimensional constraint may be biased by more than 1σ along a particular degenerate direction, and yet the marginalised biases may both be less than 1σ.

The fiducial cosmology we use in the Fisher and bias calculations is a flat w₀w_aCDM cosmology with a redshift-dependent dark energy equation of state, defined by the set of parameters Ω_m, Ω_b, σ₈, w₀, w_a, h, n_s; these are the matter density parameter; baryon density parameter; the amplitude of matter fluctuations on 8 h⁻¹ Mpc scales, which is a normalisation of the power spectrum of matter perturbations; the dark energy equation of state parameterised by w(z) = w₀ + w_az/(1 + z); the Hubble parameter H₀ = 100 h km s⁻¹ Mpc⁻¹; and the scalar spectral index of initial matter perturbations, respectively. The fiducial values are defined in Euclid Collaboration (2019a). The uncertainties and biases we quote on individual dark energy parameters are marginalised over all other parameters in this set. The survey characteristics we used are defined in Euclid Collaboration (2019a) with area of 15 000 deg² a galaxy number density of 30 arcmin⁻², and ten equi-populated tomographic bins. We used the weak lensing only “optimistic” Fisher matrix from Euclid Collaboration (2019a), where further details can be found; for a flat w₀w_aCDM cosmology the marginalised 1σ errors from that paper (Table 11) are: σ(Ω_m) = 0.034, σ(Ω_b) = 0.42, σ(w₀) = 0.14, σ(w_a) = 0.48, σ(h) = 0.20, σ(n_s) = 0.030, σ(σ₈) = 0.013 for an optimistic setting (defined in that paper). In this paper we only quote biases on dark energy parameters, relative to the expected parameter uncertainty. We chose the optimistic Fisher matrices from Euclid Collaboration (2019a) since these yield smaller expected errors and hence biases are more sensitive to systematic effects. The ℓ range used to compute these optimistic Fisher matrices is 2 ≤ ℓ ≤ 5000; for further discussion of the range we used in the computation of the δC_n(ℓ), see Sect. 3.4.

2.6. Summary of the pipeline

In Fig. 1 we summarise the overall architecture of the current concept. This propagates the changes in the quadrupole moments, converts these to observed polarisation, determines the estimated galaxy polarisation, and then determines the power spectra and the residuals. The steps are listed below.

• Survey: specifies input positional data for each galaxy, for example the position, dither pattern, slew pattern, and observation time.

• Q_ini, i: initial, intrinsic quadrupole moments are assigned to a galaxy.

• Q_shr, i: shear effects are included for each galaxy in the form of additional quadrupole moments.

• Q_PSF, i: PSF effects are included for each galaxy; these can be drawn from a distribution representing the variation in the system.

• Q_det, i: detector effects are included for each galaxy; these can be drawn from a distribution representing the variation in the system.

• Q_obs, i: observational effects are included for each galaxy such as the impact of shape measurement processes. In this paper these are not included, but we include them in the pipeline for completeness.

• M: moment measurements are converted into polarisations χ_obs, i. At this step, where the systematic effects are removed, the reference and perturbed lines separate.

• : a reference polarisation is computed, from Eq. (10), which includes , , and , which are the same values used in the construction of χ_obs, i.

• : a perturbed polarisation is computed, from Eq. (11), which includes , , and constructed from quadrupole moments drawn from relevant probability distributions that represent uncertainties in the system.

• C^R(ℓ): computes the power spectrum of .

• C^P(ℓ): computes the power spectrum of .

• δC(ℓ): computes the residual power spectrum for realisation n.

• F_αβ: computes the Fisher matrix and biases given the perturbed power spectrum, which can be used to derive uncertainties σ_α and biases b_α.

3. End-to-end pipeline

Having introduced the general formalism, we now describe the details of the current pipeline. As we work at the catalogue level, we have full flexibility over the steps that are included in or excluded from the pipeline. Furthermore, the approach (and code) is modular, giving us full flexibility in terms of developing the pipeline further. As certain steps in the pipeline mature, the relevant modules can be updated with increasingly realistic performance estimate.

3.1. Input catalogue

To evaluate the performance we need an input catalogue that contains galaxies with a range of sizes, magnitudes, and redshifts⁴. It is also important that the catalogue captures spatial correlations in galaxy properties, such as clustering, because the morphology and SED of a galaxy correlate with its local environment.

3.1.1. Mock catalogue: MICE

We used the Marenostrum Institut de Ciències de l’Espai (MICE) Simulations catalogue to assign galaxy properties, such as magnitude, right ascension (RA), declination (Dec), and shear. It is based on the DES-MICE catalogue and designed for Euclid (Fosalba et al. 2015a,b; Crocce et al. 2015). This catalogue has approximately 19.5 million galaxies over a total area of 500 deg² (11 arcmin⁻²), with a maximum redshift of z ≃ 1.4. The catalogue is generated using a halo occupation distribution (HOD) to populate friends of friends (FOF) dark matter haloes from the MICE simulations (Carretero et al. 2015). The catalogue has the following observational constraints: the luminosity function is taken from Blanton et al. (2003); the galaxy clustering as a function of the luminosity and colour follows Zehavi et al. (2011); and colour-colour distributions are taken from COSMOS (Scoville et al. 2007).

A model for galaxy evolution is included in MICE to mimic correctly the luminosity function at high redshift. The photometric redshift for each galaxy is computed using a photo-z template-based code, using only Dark Energy Survey (DES) photometry; see Fosalba et al. (2015a,b), Crocce et al. (2015) for details of the code. Our magnitude cut is placed at 20.0 ≤ m_VIS ≤ 25.0 in the Euclid visible (VIS) band. We used a 10  ×  10 deg² area of the catalogue, containing approximately 4 million galaxies.

3.1.2. Intrinsic polarisations

The MICE catalogues contain the information about the position, redshift, and (apparent) magnitudes of the galaxies and we wished to assign each galaxy an initial triplet (Q₁₁, Q₂₂, Q₁₂) of unweighted quadrupole moments. The Cauchy-Schwartz inequality for quadrupole moments implies that |Q₁₂| is bounded by . Thus, the distributions of the moments are not independent of each other and cannot be sampled independently from a marginal distribution as was done in Israel et al. (2017a). Moreover, the shapes and sizes of the galaxies depend on parameters such as redshift, magnitude, and morphology. Faint galaxies are more likely to be found at higher redshifts and thus may have smaller angular sizes; see for example M13. The polarisation distribution can have a mild dependence on the local environment as well (Kannawadi et al. 2015).

To learn the joint distribution of the quadrupole moments from real data, we used the galaxy population in the COSMOS field as our reference and assigned shapes and sizes that are consistent with the observed distribution in the COSMOS sample. Since the unweighted moments are not directly available from the data, we have to rely on parametric models fitted to the galaxies. We used the publicly available catalogue of best-fit Sérsic model parameters for COSMOS galaxies as our training sample (Griffith et al. 2012). The catalogue consists of structural parameters such as Sérsic indices, half-light radii, and polarisation prior to the PSF convolution, in addition to magnitudes and photometric redshifts for about 470 000 galaxies. In their paper, Griffith et al. (2012) model the PSF at each galaxy position.

We modelled the six-dimensional multivariate distribution of magnitude, redshift, polarisation, half-light radius, and Sérsic index using a mixture of 6D Gaussians. A generative model such as this has the advantage that we can arbitrarily generate large mock catalogues that are statistically similar to the catalogue we begin with without having to repeat the values in the original catalogue. We find that with 100 Gaussian components, we are able to recover the one-dimensional and two-dimensional marginal distributions very well. We obtain a mock catalogue, sampled from the Gaussian mixture model, with three times as many entries as the MICE catalogues have. We remove from the mock catalogue any unrealistic values (such as polarisation above 1 or redshift less than 0), caused by over-extension of the model into unrealistic regimes. We then find the closest neighbour for each galaxy in the MICE catalogues in magnitude-redshift space using a kd-tree and assign the corresponding polarisations. The orientations of the galaxies are random and uncorrelated with any other parameter, thus any coherent, intrinsic alignment among the galaxies is ignored. The model is hence too simplistic to capture the environmental dependences on shapes and sizes.

Using the knowledge of circularised half-light radii along with their Sérsic indices, the R² = Q₁₁ + Q₂₂ values assigned to the galaxies are second radial moments computed analytically for their corresponding Sérsic model. Additionally, with the knowledge of polarisation and position angle, which are in turn obtained from the best-fit Sérsic model, we obtained all three unweighted quadrupole moments (Q₁₁, Q₂₂, Q₁₂).

3.2. Survey

A key feature of our approach is that survey characteristics are readily incorporated. Having assigned the galaxy properties, we simulated a 10 × 10 deg² survey with a simple scanning strategy. We tiled the VIS focal plane following the current design; see Sect. 3.3.2.

To fill the gaps between its CCDs, Euclid will observe in a sequence of four overlapping exposures that are offset (or “dithered”) with respect to each other; a re-pointing between the sets of overlapping exposures, i.e. dither, is called a slew. The nominal pattern of offsets for exposures i = 1, …, 4 creates an “S”-shaped pattern (see Markovič et al. 2017, for more details), where the angular shifts with respect to the previous field positions are (Δx₁, Δy₁) = (50, 100); (Δx₂, Δy₂) = (0, 100); (Δx₃, Δy₃) = (50, 100) in arcsec. The code uses Mangle (Swanson et al. 2008) to create the corresponding weight map and tiles this map across the survey patch; the code is flexible enough to incorporate any dither pattern. The weight map for a pointing with four dithers is shown in Fig. 2.

Fig. 2. Coverage of a single slew by VIS. The default dither pattern in Euclid is “S”-shaped (shown as the black lines in the bottom left corner) with displacements (Δx,Δy) = (0,0; 50,100; 0,100; 50,100)″. The weights show the number of times an area has been observed. In each field of view there are 6 × 6 non-square CCDs that have asymmetric spacing between them in the vertical and horizontal directions, which results in a non-square field of view.

The propagation of the PSF and CTI stages of the pipeline, and the inverse relations described in Eqs. (10) and (11), are performed on a per exposure basis. The resulting polarisations are then averaged over all of the exposures that each galaxy receives, subject to the dither pattern; some areas of sky have fewer than four exposures, and this is captured by the dither pattern described in this work.

We also simulated a simple scanning strategy by ordering the tiling of the survey area in row (right ascension) order followed by column (declination) order, i.e. a rectilinear scanning strategy (see Kitching et al. 2016). In future implementations this will be generalised to match the full Euclid reference survey scanning strategy (Scaramella et al., in prep.).

In this first implementation and presentation of the code we did not include uncertainties in the spatial variation of foreground sources of emission or extinction. However, given the pipeline infrastructure these can be readily included and will be investigated further in future studies.

3.3. Instrumental effects

We limited our analysis to the two main sources of instrumental bias, namely uncertainties in the PSF caused by focus variations and the impact of an imperfect correction for CTI. There are other systematic effects that impact the inference of cosmological parameter using cosmic shear, which could in principle result in larger effects than these (such as photometric redshift uncertainty) but as a proof of concept we limit this study to these instrumental effects.

3.3.1. Point spread function

Correcting the observed shapes to account for their convolution by the PSF is an important step in any weak lensing measurement pipeline, and much effort has been spent on the development of algorithms to achieve this. A critical ingredient for the correction is an accurate model of the PSF itself (Hoekstra 2004). Current cosmic shear studies take a purely empirical approach where the spatial variation of the PSF is captured by simple interpolation functions that are fitted to the observations. In the case of Euclid with its diffraction-limited PSF this is no longer possible: the PSF depends on the SED of the galaxy of interest (Cypriano et al. 2010; Eriksen & Hoekstra 2018). Moreover, compared to current work, the residual biases that can be allowed are much smaller given the much smaller statistical uncertainties afforded by the data. Therefore, a physical model of the telescope and its aberrations is being developed (Duncan et al., in prep.). The PSF model parameters are then inferred using measurements of stars in the survey data, supported by additional calibration observations.

The model parameters however will be uncertain because they are determined from observations of a limited number of noisy stars. Constraints may be improved by combining measurements from multiple exposures thanks to the small thermal variations with time. The PSF will nevertheless vary with time, and thus can only be known with finite accuracy. Moreover, the model may not capture all sources of aberrations, resulting in systematic differences between the model and the actual PSF. Fitting such an incorrect model to the measurements of stars results in residual bias patterns (e.g. Hoekstra 2004), which may be complicated by undetected galaxies below the detection threshold of the algorithms used for object identification (Euclid Collaboration 2019b).

The PSF uncertainties in the pipeline are based on the current Euclid PSF wavefront model and capture one of the main sources of uncertainty, which is the nominal focus position, as detailed in Appendix A. We note that our results are expected to be somewhat conservative for this particular example because we ignore the correlations in focus positions between subsequent exposures. On the other hand, a more realistic scenario is expected to introduce coherent patterns on smaller scales caused by errors in the model itself. This will be studied in more detail in future work.

3.3.2. Detector

The VIS focal plane is comprised of 6 × 6 CCDs that each have dimensions of (2 × 2048)×(2 × 2066) pixels, where we explicitly indicate that each CCD consists of four separate readout circuits (quadrants).

Thanks to their high quantum efficiency and near linear response, CCDs are the most practical devices to record astronomical images. They are, however, not perfect and various detector effects can degrade the images. Examples include the brighter-fatter effect (BFE; e.g. Antilogus et al. 2014; Plazas et al. 2018), which affects bright objects such as stars, detection chain non-linearity, offset drifts, and photo-response non-uniformity. In this work we focus on CTI, caused by radiation damage that accumulates over time in the detectors. The resulting trailing of charge changes the measured shape and has a larger impact on fainter objects and is therefore most damaging for weak lensing studies.

There is an extensive, ongoing, characterisation programme that focusses on CTI for the Euclid detectors, the CCD273 from e2v, (see e.g. Gow et al. 2012; Hall et al. 2012; Prod’homme et al. 2014; Niemi et al. 2015). The results from this on-ground characterisation work, together with calibration measurements acquired in flight with the actual Euclid detectors, will allow the data processing to mitigate the biases caused by CTI, using correction algorithms such as those described in Massey et al. (2014). There is a fundamental floor to the accuracy of CTI correction, even if the model exactly matches the sold-state effect, owing to read noise in the CCD. The model will also have associated systematic errors and uncertainties that will translate into increased noise and residual biases for the shape measurements, with preferred spatial scales corresponding to those of the quadrants (which are approximately in right ascension and 4′ in declination) and the CCDs (which are approximately 7′ × 8′).

As there are more electrons from brighter sources, the relative loss of charge due to CTI is lower. As a result, CTI affects fainter and extended sources more (e.g. see Figs. 10 and 11 in Hoekstra et al. 2011). In our current implementation, which is detailed in Appendix B, we ignore these dependences. Instead we consider a worst case scenario, adopting the bias for a galaxy with S/N = 11 and FWHM of and a trap density that is expected to occur at mid-survey. These parameters are based on the results from Israel et al. (2015) (with updated parameters as presented in Israel et al. 2017b), who adopted the same approach.

As discussed in Appendix B, CTI is expected to increase with time as radiation damage accumulates. To account for this increase, we assume that trap densities grow linearly with time. This gradual trend is further deteriorated by intermittent steps, which are caused by solar coronal mass ejections (CMEs), which largely increase the flux of charged particles through the detectors over the baseline level. This means the estimate of the trap density parameter has to be updated periodically using images acquired in orbit. To investigate this effect in the model we define “reset on” or “reset off” cases. The two cases affect the estimated trap densities, ρ, and the associated errors in the model. In the first case the relative error in the density of species i, δρ_i, is the same throughout the whole patch of sky under study⁵, sampled from a normal distribution with zero mean and standard deviation σ_p. Hence, for each realisation all measurements in the observed patch are affected by the same relative error in trap density; we refer to this case as reset off.

The second case is reset on, in which we model the potential effect of resetting the CCD after a CME event, a so-called CME jump on scales smaller than those of the considered patches. In this case the relative error in trap densities are re-estimated midway through the patch, meaning it has one value in one half of the patch and another in the other half, both drawn from the same distribution as that used in the reset off case. And again these biases are updated (sampled from the same normal distribution) in every realisation. This scenario would correspond to a more frequent, but equally accurate, update of the trap densities than the reset off case and the coherence of the biases across the angular scales is decreased by the jumps, or resets, across the patch halves. The point is that the error is never exactly zero. But we have to re-do the model in the case of a CME jump that causes a different model uncertainty.

3.4. Power spectrum computation

For each realisation we took a spherical HEALPix map of the galaxies to make an estimate of the shear map for both the reference and perturbed catalogues. The unobserved areas are masked, and we apodised this mask with a Gaussian with a standard deviation σ = 1.5π/2048 (0.13 deg) to minimise the effect of the result of leakage due to the boundaries. We then used ANAFAST from HEALPix to calculate the E-mode power spectrum of the masked map.

Throughout we used an ℓ range 13 ≤ ℓ ≤ 4096. The minimum ℓ approximately corresponds to the maximum angular seperation in a 10 × 10 deg² patch (along the diagonal). The maximum corresponds is less than that used in the optimistic Fisher matrices used (see Sect. 2.5), where ℓ_max = 5000; however we note that the apodisation required for the power spectrum computation corresponds to a smoothing above ℓ ≃ 1500, and in we find in practice that δC_n(ℓ)→0 above ℓ ≃ 1000. Therefore we expect this assumption to have a minimal impact on results.

3.5. Pipeline set-up

A key feature of our approach is that we create realisations of the systematic effects, for each galaxy and each pointing, which enables us to determine the expected probability distributions for the changes in the cosmological parameter inferences caused by these systematic effects. This is done by creating 150 random realisations that are propagated through the Fisher matrix and bias calculations as discussed in Sect. 2.5; we chose 150 since this then means the total area is 150 × 100 square degrees which is equal to the total Euclid wide survey. The run in which we combined PSF and CTI residuals took 20 h to compute on a machine with 25 1.8 GHz CPUs and 6 GB RAM. The PSF-only scenario took 14 h, and the CTI-only run took seven hours on the same architecture. As each realisation can be run in parallel, the calculations can be sped up accordingly on a machine with more processors.

4. Results

As a demonstration of the usefulness of our approach, we assessed the impact of two prime sources of bias for the Euclid cosmic shear analysis: PSF and CTI modelling. We computed the expected residual systematic power spectra caused by imperfect removal of systematic effects from realistic uncertainties in the modelling. We then propagated the power spectrum residuals through a Fisher matrix to compute the biases in dark energy parameters.

4.1. PSF

The top left panel of Fig. 3 shows the residual systematic power spectrum caused by uncertainties in the PSF model caused by focus variations. The thick line indicates the mean of the 150realisations, whereas the thin lines delineate the 68% interval. As discussed in Appendix A, we consider only the uncertainty in the PSF model given the assumed nominal focus position, which is the dominant contribution and introduces residuals in the power spectrum on large scales. Other imperfections in the optical system typically introduces residuals on smaller scales.

Fig. 3. Residual power spectra caused by imperfect removal of systematic effects. Thin lines show 68% intervals. Top left panel: residual power spectrum due to PSF, caused by the limited precision with which the nominal focus position can be determined from the stars in the data; it can be seen that residuals have on average been removed. Top right panel: residual power spectrum caused by CTI when the CTI-removal model parameters are updated throughout the survey (reset on case, see text for details). There are residuals on the scales corresponding to half the distance between the CCDs, as shown in the insets. Bottom right panel: results when the CTI-removal model parameters are kept constant during the survey (reset off). As can be seen the residuals have a slightly wider distribution compared to the reset on case. Bottom left panel: residual systematic effects from uncertainties in the modelling of both PSF and detector effects; as shown in the inset the two effects seem to work in opposite directions, where the positive offset present in the PSF-only case has reduced in the combined case. We note that owing to the sensitivity of dark energy parameters to relatively large angular scales ℓ ≃ 50−1000, the deviations on these scales are of more importance.

To understand the relevant scales in the PSF case, it is helpful to look at Fig. 4, where some of the relative correlated scales are indicated. A point in one field of view is correlated with the same point in all the other fields of view; i.e. the angular distances between the field of view are also relevant, not only the scales of field of view itself. Also the field of view is not square, and hence the distances to the same point in the fields of view are not the same in both directions. In our 10 × 10 deg² area, this gives us correlated scales in the range 13 ≤ ℓ ≤ 300. The minimum distance between adjacent fields of view corresponds to ℓ = 300, and the diagonal in our square survey area (the maximum angular separation) corresponds to ℓ = 13. Incidentally this is also the range in which cosmic variance dominates.

Fig. 4. Part of the observed area with 3 slews in each direction and 4 dithers for each slew. The slews are plotted at 1.2× their nominal value for presentation purposes, causing apparent gaps, which are not present in the actual simulated survey. The lines show some of the correlated scales relating to the same point in each field of view. We also note that there are correlations at 2×,3×,n× of these harmonic scales. It should be noted that relevant scales are determined by the distances between the fields of view, not the size of the field of view itself.

The average residual power spectrum in the top left panel of Fig. 3 is close to zero and does not show sharp features, but the residual PSF biases contribute over a range of scales. This is because the averaging over the four dithers for each slew reduces the average induced biases in the polarisations, which in turn reduces the correlations between slews; and the polarisations in the perturbed line for each field of view (i.e. each dither and each slew) are drawn from a distribution, so that the average impact is typically less extreme.

4.2. CTI

The thick line in the top right panel in Fig. 3 shows the average residual power spectrum when we consider the imperfect correction for time-dependent CTI for the reset off case (see Sect. 3.3.2). The amplitude of the residuals are slightly larger than that of the PSF case. Compared to the PSF case, there are additional angular scales on which correlations can occur, namely the distances between the CCDs in the detector. The inset shows a zoom in around ℓ ≃ 3080, which corresponds to half the distance between CCDs. This is because in our setting, CTI systematic effects are induced only in the serial readout direction (see Appendix B), inducing biased polarisation estimates at half the CCD scale (quadrant scale).

In the second case, reset on (see Sect. 3.3.2), the results presented in the bottom right panel of Fig. 3 show that this procedure does not improve the residuals around ℓ ≃ 3080. It does, however, reduce the variance on the largest scales, even though the average residual power spectrum is largely unchanged, except for increased variation for ℓ in the range 150−300.

4.3. PSF and CTI

Rather than considering individual sources of bias separately, we can simultaneously propagate different types of systematic effects and capture their correlated effects. This is demonstrated in the bottom left panel of Fig. 3, which shows the residual systematic power spectrum resulting from both CTI (reset on) and PSF systematic uncertainties⁶. Both features of CTI and PSF systematic effects can be seen in the residual power spectrum. The inset shows the residual power spectrum in the range corresponding to the CCD scales, where CTI contributes most. The residuals on these scales are now dominated by both the CTI and PSF systematic effects.

4.4. Impact on cosmology

For each residual power spectrum we compute the change in the expected maximum likelihood locations for the parameters w₀ and w_a. The tolerable range for biases on dark energy parameters is generically |(b/σ)_w0| ≤ 0.25 (where b is the bias, and σ is the 1σ marginalised uncertainty) as derived in M13 and Taylor et al. (2018), which ensure that the biased likelihood has a greater than 90% overlap integral with the unbiased likelihood. This tolerance is applicable for all systematic effects in an experiment, not per each systematic effect.

The results are presented in Fig. 5 and reported in Table 1. We show results for the PSF-only case, the CTI-only case with resetting on, and the combined case. The panels, respectively, show the biases in w₀ and w_a relative to the statistical uncertainty. In Table 1 we list the mean and its uncertainty for the quantities as well as the standard deviation of the distributions themselves. We also quote the 90% confidence limits of the bias distributions.

Fig. 5. Left panel: ratio of the bias in w0 and the 1σ uncertainty in this parameter for PSF-only (cyan), CTI-only with resetting on (blue), and both PSF and CTI with resetting on (red) scenarios. We note that the two darker purple shades are overlaps between the red and blue and red, blue and green histograms. Right panel: ratio of the bias in wa and the 1σ uncertainty in this parameter. Although the distributions are wide in some scenarios, we find that they are well within limits set in C13; see also Table 1.

We find that the PSF residuals have a minimal impact, which is expected as the amplitudes of the residual power spectra were small. The induced biases b, relative to the uncertainty σ on the dark energy parameters are expected to be (b/σ)_w0 = [−0.024, 0.033] and (b/σ)_wa = [−0.042, 0.015] at 90% confidence interval. These are well within the tolerable range.

For the case in which the CTI model parameters are kept fixed during the simulated observations of a 100 deg² patch (reset off), the impact on the induced biases are (b/σ)_w0 = [−0.328, 0.077] and (b/σ)_wa = [−0.054, 0.281], which are just outside the tolerable range. However for the case in which we resample the CTI model parameters (reset on), the results are improved with (b/σ)_w0 = [−0.078, 0.152] and (b/σ)_wa = [−0.121, 0.067]. The effects seen are very similar to effects seen using the simplified models of CTI in Kitching et al. (2016).

Perhaps most interesting are the results for the case in which we include both CTI and PSF residuals, since this joint case was not captured in the C13 “flow down”. We find that the biases are expected to be (b/σ)_w0 = [−0.046, 0.144] and (b/σ)_wa = [−0.124, 0.032], which is again within the tolerable range.

4.5. Discussion

It is useful to compare our findings to the requirements derived in C13. In the latter study, requirements on systematic effects were set through a formalism that flowed down (i.e. subdivided requirements in progressively finer details via a series of inter-related subsystems) changes in the power spectrum parameterised by

(23)

Requirements on ℳ and 𝒜 were determined for various effects such as PSF and CTI. To compare to this formalism we could naively fit the residual power spectra that we find using such a linear model. However, this would neglect the correct formulation of how to propagate biases into cosmic shear power spectra (Kitching et al., in prep.).

Therefore to assess the difference between the C13 approach and our approach we need to flow up the requirements on the uncertainties set in C13 (referred to as σ in that paper) for individual effects, and compare the outcome of the two approaches at the level of biases in cosmological parameters rather than comparing ℳ and 𝒜 values. We do this by determining the multiplicative and additive biases, ℳ and 𝒜, associated with each systematic effect in C13, constructing Eq. (23) for these values, and then adding this to Eq. (21); a process we refer to as flow up.

Whilst uncertainties are included in this flow-down approach, these are taken to be constant across the survey (both spatially and temporally). They are also assumed to be independent of each other. Our approach does not suffer from these limitations. By modelling biases simultaneously, they also have a chance of acting at different scales, or even cancelling each other out. Hence any comparison with prior work should not be interpreted as there being a margin in previously derived requirements. Nevertheless, such a comparison is useful to show how different the approaches are, and if previous requirements were exceeded this would be of concern.

Assuming PSF modelling errors in the shear power spectrum at the maximum values permitted by the C13 requirements of 𝒜 = 5 × 10⁻⁸ and ℳ = 4.8 × 10⁻⁴, we find biases on cosmological parameters (b/σ)_w0 = 0.25 and (b/σ)_wa = 0.31. Assuming CTI correction biases at the maximum values permitted by C13 of 𝒜 = 1.21 × 10⁻⁸ and ℳ = 0 (CTI contributions to multiplicative bias are subdominant) yields (b/σ)_w0 = 0.14 and (b/σ)_wa = −0.2. In contrast, our flow-up analysis predicts biases on cosmological parameters that are lower by a factor between 2 and 5. None exceed previously derived requirements and all are within acceptable tolerances to meet top-level scientific goals.

We note that the scatter in δC_n(ℓ) and the corresponding scatter in the cosmological parameter biases are smaller in the PSF+CTI case than for the PSF alone. What is happening in this case is a partial cancellation of terms where the CTI is predominately a positive Q₁₁ component at the chip edges, whereas the perturbations in the PSF can have negative Q₁₁ components. Therefore in combination the overall amplitude cancels out in some regions of the field of view. This is the first time that these systematic effects have been combined, and it is not unexpected that systematic effects that act in opposite senses may cancel each other out. However we leave a full investigation of these effects to future work.

Finally we emphasise several assumptions in this analysis that should be relaxed in the future, which may mean the results are either optimistic or pessimistic:

• We do not model intrinsic alignments, the environmental dependence on the intrinsic size and shapes of galaxies.

• The smooth increase in CTI over adjacent pointings may be considered optimistic, if CTI has sudden jumps in reality. Furthermore the choice of 45%−55% end-of-mission radiation dose is average. In a tomographic analysis CTI residuals may also mimic redshift-dependence of cosmic shear, which may mean the results are optimistic.

• The uncorrelated PSF residuals between consecutive exposures may be conservative or optimistic, depending on the final state of the telescope at launch.

5. Conclusions

We have presented an end-to-end approach that propagates sources of bias in a cosmic shear survey at a catalogue level. This allowed the capture of spatial variations, temporal changes, dependences on galaxy properties, and correlations between different sources of systematic and stochastic effects in the pipeline. We use our methodology to revisit the performance of a Euclid-like weak lensing survey. We limit the analysis to quantify the impact of imperfect modelling of the PSF and CTI, as these are two major sources of bias caused by uncertainty in the modelling of instrument(s). Other effects can be readily included, which will be done in future work.

The PSF systematic effects are introduced through the expected uncertainty in fitting the PSF model to noisy data given the assumed nominal focus of the telescope. Additional imperfections introduce residuals on smaller scales, but these should not affect our main conclusions because the dark energy measurements are most sensitive to variation on large scales. We also consider a time-dependent CTI, where the CTI increases with the survey time due to accumulation of radiation damage on the detectors. We consider a conservative scenario because the parameters we adopted apply to the faintest galaxies in the analysis, whereas the biases are smaller for brighter objects. We also do not include intrinsic alignment effects or source blending effects, both of which will be included in future studies.

These effects were propagated through to residual cosmic shear power spectra and cosmological parameters to estimate the expected biases in the parameters w₀ and w_a. Compared to requirements based on a more restricted flow-down approach by C13, we find that the biases on the dark energy parameters from our more realistic performance estimates are well within the requirements. Even for the combined scenario of CTI and PSF, we find the biases on dark energy parameters are well within the required tolerances.

This paper presents the first step towards a more comprehensive study of the performance of a Euclid cosmic shear survey. The same approach however can also be readily applied to other cosmic shear surveys. In future work we will introduce more complexity in the PSF and detector systematic effects, so that the resulting redshift dependences of these effects can be assessed. As alluded to earlier, CTI is dependent on flux and morphology, which implies it will change with redshift. Other systematic effects, such as shape measurement uncertainties, will also be implemented in the pipeline. These improvements will enable us to examine the impact of systematic effects on an increasingly realistic tomographic analysis.

Acknowledgments

PP is supported by an STFC consolidated grant. CW is supported by an STFC urgency grant. TDK is supported by a Royal Society University Research Fellowship. HH acknowledges support from Vici grant 639.043.512 and an NWO-G grant financed by the Netherlands Organization for Scientific Research. LM and CD are supported by UK Space Agency grant ST/N001796/1. VFC is funded by Italian Space Agency (ASI) through contract Euclid – IC (I/031/10/0) and acknowledges financial contribution from the agreement ASI/INAF/I/023/12/0. We would like to thank Jérome Amiaux, Koryo Okumura, Samuel Ronayette for running the ZEMAX simulations. AP is a UK Research and Innovation Future Leaders Fellow, grant MR/S016066/1, and also acknowledges support from the UK Science & Technology Facilities Council through grant ST/S000437/1. MK and FL acknowledge financial support from the Swiss National Science Foundation. SI acknowledges financial support from the European Research Council under the European Union’s Seventh Framework Programme (FP7/2007–2013)/ERC Grant Agreement No. 617656 “Theories and Models of the Dark Sector: Dark Matter, Dark Energy and Gravity”.

References

Appendix A: Details of PSF modelling

In this section, we describe the propagation of uncertainties which result from inaccuracies in the PSF model, using the broadband parametric phase retrieval method from Duncan et al. (in prep). In this method, the PSF variation is modelled in the wavefront domain. The corresponding real-space optical PSF is obtained as the modulus squared of the Fourier transform of the wavefront at the exit pupil, including the effect of telescope distortion, integrated over the bandpass. In the present application, a simple Gaussian model for the distribution of telescope guiding accuracy is used, linear detector effects (including pixelisation and linear charge diffusion) are included, and it is assumed that the detectors lie precisely in the focal plane. More realistic guiding and detector effects will be included in the future.

The wavefront at the exit pupil of the telescope describes the coherent perturbations in the optical path differences of infalling photons, caused by the design and alignment of the telescope optical elements. The wavefront can be split into two parts: an amplitude component, which describes where the light is vignetted by structures in the telescope, and a phase component, which describes the variation of the optical path differences. Both change with position in the focal plane. To capture the amplitude variation, we use a geometric model that describes the projection of intervening structures in the telescope (i.e. the secondary mirror M2 and its struts) at the focal plane. To model the phase variation, we used a suite of simulated wavefronts obtained with the optical design program ZEMAX⁷, configured to the specifications of Euclid. Each phase map was fitted by a sum of Zernike polynomials, and the variation of the corresponding Zernike coefficients with focal plane position was captured by a set of polynomials. Several optical elements in the telescope design were displaced or deformed by turn, and the corresponding effects on the phase maps were captured by so-called telescope modes. As a result, the wavefront can be predicted for any telescope set-up, with a realistic focal plane variation. Given a model wavefront, the real-space PSF is then computed for a range of densely sampled wavelengths. The final PSF is obtained by integrating over the spectral telescope response, weighted by the SED of the source with additional convolution effects of guiding and CCD pixel response included. In this application of the model, detector offset and high frequency contributions such as those arising from surface errors are not included.

As the Euclid VIS PSF model is jointly fitted to stars in the entire field of view, PSF errors are correlated across that field. In order to capture this, we investigated the effect of varying one of the principal model parameters, the wavefront error associated with the defocus of the telescope. Higher-order wavefront errors are also expected to contribute to the PSF uncertainty. But the PSF variations of this mode should be a realistic representation of the actual correlated PSF errors in the absence of possible effects at cryogenic temperatures, which could cause deformation in the y-axis displacement of M2 rather than the z-axis. We modelled the effect of a shift in the focus position resulting in an optical defocus for a given source SED. We chose the source SED to be the template spiral Sbc galaxy of Coleman et al. (1980) with redshift of 1.

We assume a nominal offset focus position that is drawn at random from a normal distribution whose variance matches the expected σ_z ≃ 0.5 μm uncertainty in this model parameter, which can be obtained from fitting the telescope model to the stars that appear in each survey field. Even consecutive exposures are assumed to have independent nominal focus values.

A realisation of errors in the quadrupole moments of the PSF was then obtained in a two-step process. In the first step, the focus position was shifted to give a minimum in PSF size at the centre of the field of view, which was taken as the nominal in-focus position. The mirror was then perturbed in both directions (positive and negative z offsets) until the second-moment measure of the model PSF’s size, R², varied from the nominal value by a tolerance . This value is both the requirement on the knowledge of the mean PSF size across the survey, set by C13, and also is approximately the measurement uncertainty that we expect to obtain from measurements on individual survey exposures.

In a second step, the FoV-dependent behaviour was determined by fitting a sixth order polynomial⁸ to ΔQ_ij(x, y, Δz) = Q_ij(x, y, z)−Q_ij(x, y, z_nom) across x and y, where x and y describe the FoV position, z the focus mirror position, z_nom the nominal focus position and Δz = z − z_nom. The variation with z was modelled by training the coefficients of the FoV position fit across Δz = 0 ± 1.35 μm, using a quadratic form. Finally, it was verified that the fit recovered the expected perturbed quadrupole moments to ≲5% accuracy across the field of view, except where ΔQ was close to zero so that relative differences – defined as – became large owing to numerical inaccuracy. This verification was conducted both on the training data, as well as a coarse grid of field-of-view values and M2 shifts within the range σ_z = 0.5 μm, which was not used to train the fit. The residual quadrupole moment variations that these changes induce is shown in Fig. A.1. Figure A.2 shows the variation in polarisation across the field of view at the nominal focus position (z = z_nom).

Fig. A.1. Change in the quadrupole moments in units of Euclid pixels squared, for focus position shifts of Δz = −1.35 μm (left columns) and Δz = 1.35 μm (right columns), for Q11 (top panels), Q12 (middle panels), and Q22 (bottom panels) as a function of field of view position (in degrees).

Fig. A.2. Stick plot describing the variation in polarisation across the field of view (in degrees) for the nominal focus (z = 0). All values shown are taken from the training data, not the fit.

Appendix B: Details of CTI modelling

Charge transfer inefficiency is caused by the capture and delayed release of photoelectrons by traps, i.e. localised, unintended quantum levels in the silicon lattice of a detector. These defects are created when high energy particles displace silicon atoms and, above the protection of the Earth’s atmosphere, will accumulate throughout the mission. The timescales for capture and release depend upon the type of damage and the proximity to any lattice impurities. During readout, if electrons are captured from a charge packet that is moving through the lattice and released after a sufficient delay, they become part of a later charge packet (or pixel in the resulting image). This creates faint luminous “trails” behind the images of galaxies and stars, which bias measurements of their polarisation and size. Objects farther from the readout nodes gain brighter trails because their electrons must travel farther and are subject to more traps.

We adopted the Israel et al. (2015) model of CTI in the Euclid VIS serial readout direction⁹. This treats charge capture as instantaneous and charge release as a stochastic process governed by exponential decay (Massey et al. 2014). Euclid CCDs will contain three trap species in the serial register. Each species i has a different characteristic release time τ_i, and a time-evolving surface density (abundance) ρ_i. Table B.1 shows the trap properties expected after the radiation dose accumulated by the end of the mission (for 90% of realisations of solar weather).

To account for the accumulation of radiation damage over time, we assume that trap densities grow linearly over the patch of the sky observed, from 45% of the values in Table B.1 at one end to 55% at the other. This is very conservative for large angular scales as it introduces an increase of approximately 10% over roughly a ten-day cycle (i.e. the typical amount of time it would take to observe our 100 deg² patch), which is much larger than the milder increase over the actual mission. On the other hand, we assumed consecutive observations, and hence consecutive increases in CTI at each exposure in this patch. In reality, a 100 deg² patch of the observed area will not have this smooth increase as the exposures that cover it will not, in general, be consecutive. Therefore, this is also a somewhat optimistic approach in that sense.

Euclid’s baseline strategy for CTI mitigation is a pixel-by-pixel movement of flux from trails, back to the pixels it came from. This “back-clocking” approach is limited by read-out noise (RON), model parameter uncertainties, and model inaccuracy. Because RON is added at the amplifier, it is not trailed during readout – but it is spuriously corrected by pixel-level methods as if it had been trailed (see Sect. 5.3 of Israel et al. 2015).

In the reference case (i.e. assuming a perfect CTI model), the model and model parameters are known perfectly and the only source of bias is the RON. Images of galaxies contain a residual shape measurement error

(B.1)

where Δη may refer to either Δχ₁ or ΔR²/R²; N_Tr is the “serial”/“horizontal” distance (in pixels) of the object to the readout amplifier, and = 2099 pixels is the maximum number of serial transfers given the detector design; the function f^res has the form shown in Eq. (14) of Israel et al. (2015) and coefficients listed in rows 7 and 8 of Table 1 in Israel et al. (2015). The linear prefactor reflects the increasing number of transfers (hence increasing number of encountered traps) traps for galaxy images farther from the readout amplifier. Finally, to model the non-deterministic filling history of traps along the readout direction, which can perturb CTI trailing, we add Gaussian noise to the scaled Δη in each exposure with zero mean and 5% standard deviation. This was not considered in the analysis of Israel et al. (2015) and we find its effect to be negligible.

In the perturbed case, we also propagate uncertain knowledge of the model parameters (e.g. accumulated radiation dose) at each point in time. Since all our galaxies have the same flux, we introduce model parameter errors in the assumed density Δρ_i and release time Δτ_i of the traps. In addition to errors from back-clocking the RON, as in the reference case, Israel et al. (2015) find that model parameter errors introduce a bias,

(B.2)

(B.3)

where the function f^deg provides the change (“degradation”) in shape parameters because of CTI, without mitigation, as a function of the model parameters. Its functional form is shown in Eq. (14) of Israel et al. (2015) and it uses coefficients listed in rows 3 and 4 of Table 1 in Israel et al. (2015). The difference of this function evaluated at (ρ + Δρ, τ + Δτ) from the same function at (ρ,τ) in Eq. (B.3) is a reflection of how the iterative mitigation of CTI is, at its root, equivalent to an additional degradation of the images, similar to that caused by CTI, but applied in the opposite sense (hence the name back-clocking). We note that Eq. (B.3) is equal to Eq. (17) in Israel et al. (2015). Halfway through the mission, both terms account for roughly equal levels of residual.

To assign values to the biases in model parameters, we adopt a constant bias Δτ_i = 1% in the release time parameters; and Δρ_i drawn from a Gaussian distribution with zero mean, where the average bias is zero; and standard deviation of 1%, over the true value of ρ at each time. These are both conservative in the sense that they could be derived from Euclid calibration each day (Nightingale et al., in prep.), but τ_i are likely to be constant for the entire mission and ρ_i smoothly increasing, so errors could be reduced by iterative calibration. These values therefore do not necessarily reflect the ultimately achievable uncertainty in the model parameters, but are useful as reference values. Figure B.1 shows the pattern of induced biases due to imperfect CTI mitigation for a random selection of galaxies in one FoV.

Fig. B.1. Random selection of galaxies are shown representing the pattern of the induced polarisations owing to imperfect CTI mitigation in one field of view. The biases are larger with distance from the readout nodes on either side of the CCDs. Biases are only considered in the serial direction.

All Tables

All Figures

Fig. 1. Overall structure of the concept as described in the main text. The quadrupole moments Q are initiated with intrinsic moments and then modified by incorporating the shear, PSF, and detector effects. Survey characteristics such as dither pattern, slew pattern, and observation time are entered in the initial catalogue. A measurement process M subsequently converts the observed moments to polarisations. The estimation of the galaxy polarisation is then made (as described in Eqs. (10) and (11)). This is done per object. Next a power spectrum for the reference and perturbed scenarios is computed. For the perturbed line the PSF and detector moments are drawn from distributions that represent the measurement uncertainty as described in the text. This process is repeated for 150 random realisations for the set of galaxies that are in the input catalogue. Finally the residual power spectrum is computed per realisation, and the statistics of each of the realisations is passed onto the Fisher matrix, from which uncertainties and biases of dark energy parameters are calculated. White circles indicate moment space, where modifications are performed on an object-by-object basis. Grey circles indicate ensemble average in the harmonic space. Diamonds show cosmological parameter space.

In the text

Fig. 2. Coverage of a single slew by VIS. The default dither pattern in Euclid is “S”-shaped (shown as the black lines in the bottom left corner) with displacements (Δx,Δy) = (0,0; 50,100; 0,100; 50,100)″. The weights show the number of times an area has been observed. In each field of view there are 6 × 6 non-square CCDs that have asymmetric spacing between them in the vertical and horizontal directions, which results in a non-square field of view.

In the text

Fig. 3. Residual power spectra caused by imperfect removal of systematic effects. Thin lines show 68% intervals. Top left panel: residual power spectrum due to PSF, caused by the limited precision with which the nominal focus position can be determined from the stars in the data; it can be seen that residuals have on average been removed. Top right panel: residual power spectrum caused by CTI when the CTI-removal model parameters are updated throughout the survey (reset on case, see text for details). There are residuals on the scales corresponding to half the distance between the CCDs, as shown in the insets. Bottom right panel: results when the CTI-removal model parameters are kept constant during the survey (reset off). As can be seen the residuals have a slightly wider distribution compared to the reset on case. Bottom left panel: residual systematic effects from uncertainties in the modelling of both PSF and detector effects; as shown in the inset the two effects seem to work in opposite directions, where the positive offset present in the PSF-only case has reduced in the combined case. We note that owing to the sensitivity of dark energy parameters to relatively large angular scales ℓ ≃ 50−1000, the deviations on these scales are of more importance.

In the text

Fig. 4. Part of the observed area with 3 slews in each direction and 4 dithers for each slew. The slews are plotted at 1.2× their nominal value for presentation purposes, causing apparent gaps, which are not present in the actual simulated survey. The lines show some of the correlated scales relating to the same point in each field of view. We also note that there are correlations at 2×,3×,n× of these harmonic scales. It should be noted that relevant scales are determined by the distances between the fields of view, not the size of the field of view itself.

In the text

Fig. 5. Left panel: ratio of the bias in w0 and the 1σ uncertainty in this parameter for PSF-only (cyan), CTI-only with resetting on (blue), and both PSF and CTI with resetting on (red) scenarios. We note that the two darker purple shades are overlaps between the red and blue and red, blue and green histograms. Right panel: ratio of the bias in wa and the 1σ uncertainty in this parameter. Although the distributions are wide in some scenarios, we find that they are well within limits set in C13; see also Table 1.

In the text

	Fig. A.1. Change in the quadrupole moments in units of Euclid pixels squared, for focus position shifts of Δz = −1.35 μm (left columns) and Δz = 1.35 μm (right columns), for Q11 (top panels), Q12 (middle panels), and Q22 (bottom panels) as a function of field of view position (in degrees).
In the text

	Fig. A.2. Stick plot describing the variation in polarisation across the field of view (in degrees) for the nominal focus (z = 0). All values shown are taken from the training data, not the fit.
In the text

	Fig. B.1. Random selection of galaxies are shown representing the pattern of the induced polarisations owing to imperfect CTI mitigation in one field of view. The biases are larger with distance from the readout nodes on either side of the CCDs. Biases are only considered in the serial direction.
In the text