Propagating photo-z uncertainties: a functional derivative approach

Photometric redshifts are a key ingredient in the analysis and interpretation of large-scale structure (LSS) surveys. The accuracy and precision of these redshift estimates are directly linked to the constraining power of photometric surveys. It is hence necessary to define precision and accuracy requirements for the redshift calibration to not infer biased results in the final analysis. For weak gravitational lensing of the LSS, the photometry culminates in the estimation of the source redshift distribution (SRD) in each of the tomographic bins used in the analysis. The focus has been on shifts of the mean of the SRDs and how well the calibration must be able to recover those. Since the estimated SRDs are usually given as a normalized histogram with corresponding errors, it would be advantageous to propagate these uncertainties accordingly to see whether the requirements of the given survey are indeed fulfilled. Here we propose the use of functional derivatives to calculate the sensitivity of the final observables, e.g. the lensing angular power spectrum, with respect to the SRD at a specific redshift. This allows the propagation of arbitrarily shaped small perturbations to the SRD, without having to run the whole analysis pipeline for each realization again. We apply our method to an EUCLID survey and demonstrate it with SRDs of the KV450 data set, recovering previous results. Lastly, we note that the moments of the SRD of order larger than two will probably not be relevant when propagating redshift uncertainties in cosmic shear analysis


INTRODUCTION
Cosmic shear, the weak gravitational lensing effect imprinted on distant galaxies by the large-scale structure (LSS), is one of the primal science goals for EUCLID and Rubin-LSST.The blueprint for these missions has been set by current stage-3 surveys, including the Kilo-Degree Survey (Hildebrandt et al. 2017b;Asgari et al. 2021, KiDS) 1 , the Dark Energy Survey (Abbott et al. 2018;Amon et al. 2022, DES) 2 and the Subaru Hyper Suprime-Cam (Takada 2010;Hamana et al. 2020, HSC) 3 , yielding tight constraints on the matter distribution in the late Universe.
The cosmic shear signal is estimated by measuring the coherent distortion of background galaxies.Since the intrinsic ellipticity of galaxies is much larger than the lensing effect, millions (or even billions) of galaxies are required to measure a significant signal.This makes a complete spectroscopic survey unfeasible.There are two main techniques to obtain an estimate of the true redshift for the background galaxy sample: The first one is called photometric redshifts (see e.g.Lima et al. 2008;Bonnett et al. 2016; Hildebrandt ⋆ E-mail: reischke@astro.ruhr-uni-bochum.de 1 https://kids.strw.leidenuniv.nl/ 2 https://www.darkenergysurvey.org/ 3 https://hsc.mtk.nao.ac.jp/ssp/ et al. 2021), that is using broadband photometry which is then calibrated using a significantly smaller spectroscopic reference sample.Recently self-organising maps have become a standard technique for photometric redshifts (Masters et al. 2015;Wright et al. 2020;Hildebrandt et al. 2021;Myles et al. 2021).The second approach is clustering redshifts (see e.g.Newman 2008;Matthews & Newman 2010 or van den Busch et al. 2020;Gatti et al. 2022 for more recent works), where the redshift distribution of the galaxy sample is estimated by an angular cross-correlation measurement.There also exist hybrid methods, combining photometry and clustering measurements into a Bayesian hierarchical model (Sánchez & Bernstein 2019;Alarcon et al. 2020;Rau et al. 2020;Rau et al. 2022).
All of the above methods yield an as unbiased as possible estimate of the distribution of background galaxies in redshift (sourceredshift-distribution, SRD).The statistical precision of the experiments sets limits on the required accuracy of the SRD estimation technique.Likewise, the question arises of how the residual uncertainties in the SRD itself affect the inference process.Most works use simple shift parameters in the mean of the SRD which are then marginalised over in a given (informed) prior range.There are, however, other approaches using different parametrisation for these uncertainties which itself can be informative priors or even be self-calibrated using different two-point function measurements.This includes non-parametric approaches (Rau et al. 2020), using higher moments beyond the mean (McLeod et al. 2017), explicit parametrisation of outliers in photometric redshifts (Schaan et al. 2020) or Gaussian mixture models (Stölzner et al. 2021) Recently full shape techniques, i.e. methods beyond the mean shifts of the SRD, have been applied to real data.In Amon et al. (2022) used Hyperrank (Cordero et al. 2022) to propagate the uncertainties of the SRD, finding that simple shift parameters are sufficient for DES.For the next generation of surveys, stage-4, this might no longer be true, however.Including the full-shape SRD might also require more efficient techniques for the marginalisation over the nuisance parameters, something which is for example analytically done in Stölzner et al. (2021) thanks to the Gaussian mixture model used.Another possibility was presented in Zhang et al. (2023) where the SRD was sampled from the full shape and the resulting Markov chains of the cosmological parameters were combined using Baysian model evidence.
In summary, it is not entirely clear to date how accurate the error propagation of the residual SRD uncertainties have to be in order to obtain unbiased cosmological results, i.e. uncertainties on the cosmological parameters together with their maximum posterior values.While unbiased results are obtained with simple shifts at the moment, EUCLID and Rubin-LSST will change this and the full shape of the SRD posterior distribution becomes important for cosmological inference.It is hence vital to investigate the most general sensitivity of cosmic shear observables to the underlying SRD, without assuming a specific functional form of the SRD itself.Therefore, in this paper, we calculate the functional derivative of the cosmic shear angular power spectrum with respect to the SRD.That is, we investigate arbitrary (but small) perturbations to the SRD at a certain redshift and how they propagate into the angular power spectrum.This functional derivative can then be used to calculate the total error in the measurement process by simple error propagation.We take the constraint of the normalisation of the SRD into account when calculating the functional derivative.Therefore we can propagate arbitrary perturbations to the SRDs (subject to some underlying covariance) and propagate them into the C ℓ of cosmic shear.This allows us to estimate the difference in χ 2 induced by the uncertainty in the SRD, without having to run thousands of realizations of the analysis pipeline used.By using a Fisher matrix for the cosmological parameters, this ∆χ 2 can then be mapped to potential biases in cosmological parameters.Here we studied a rather idealised scenario by working in Fourier space, assuming a Gaussian likelihood and ignoring intrinsic alignments.The method, however, easily generalises and including these effects is straightforward.The approach therefore tries to fill the gap between cheap marginalisation over shift parameters and very expensive marginalisation over different Monte Carlo Markov Chains (MCMCs).It furthermore allows us to investigate where uncertainty in the SRD is wreaking the most havoc on cosmological inference.
We structure the paper as follows: In Section 2 we briefly review cosmic shear basics and introduce the methodology used by calculating the functional derivative of the weak lensing angular power spectrum.The results are presented in Section 3, where we apply the procedure to a survey with EUCLID's specifications and to KiDS-VIKING-450 (KV450).We conclude in Section 4. In the appendices, we also investigate the possibility of an Edgeworth expansion of the SRD (Appendix A), discuss photometric galaxy clustering (Appendix B), the distribution of the mean and standard deviation of the SRD in Appendix C, the general relationship to observables (Appendix D), the functional derivative of the non-Limber projection in Appendix E and intrinsic alignments (Appendix F).

METHODOLOGY
In this section, we present the basic methodology of our analysis.In particular, we describe the basics of cosmic shear and derive the functional derivative of the lensing angular power spectrum with respect to the SRDs.

Cosmic shear basics
The equation for the cosmic shear power spectrum in tomographic bins i and j in the Limber projection is (Limber 1954;Loverde & Afshordi 2008) where P δ is the matter power spectrum, for which we use the emulated spectrum from Mead et al. (2015).W (i) κ (χ) is the lensing weight of the i-th tomographic bin as given by: Here χ is the co-moving distance, a the scale factor, Ω m0 the matter density parameter today, χ H the Hubble radius and n (i) s is the SRD in the i-th tomographic bin which builds on photo-z measurements and its calibration.It is normalized in each tomographic bin such that dz n (i)  s (z Since photo-z is just an estimate of the true redshift, the estimated source-redshift distribution, n (i) s , is not exactly known.Here we investigate two approaches: i) Use functional derivatives to evaluate the change of the lensing power spectrum when perturbing the n (i)  s at different redshifts.Given specific survey settings and precision goals, limits on the allowed change of the n (i)  s can be determined, which in turn can be mapped to changes in the cumulants or moments of the underlying distribution (see Section 2.2).
ii) We expand the underlying source-redshift distribution in an asymptotic Edgeworth series and investigate the requirements on the cumulants directly in a Fisher analysis.The second approach is not feasible for realistic SRDs (see Appendix A).

Functional derivative of the lensing power spectrum
Here we wish to investigate the sensitivity of the weak lensing power spectrum to the full shape of the source-redshift distribution using functional derivatives.In particular we start by perturbing n (i)  s (χ(z)) at a certain redshift z 0 , such that χ 0 = χ(z 0 ).The corresponding perturbed lensing weight is thus s (χ 0 ) This expression evaluates, how the lensing weight changes if the source-redshift distribution is perturbed by an amount ∆n (i) s at the co-moving distance χ 0 corresponding to the redshift z 0 .
Ultimately, we are interested in the change of the lensing power spectrum, Equation (1).First, by applying the Leibniz rule The missing ingredient is the functional derivative of the lensing kernel, for which we find Θ(x) is the Heaviside function to ensure that the functional derivative vanishes if the SRD is perturbed outside the integration bounds.
Using Equation (4) and Equation ( 5) we can write the change in angular power spectrum ∆C (i j) ℓ (χ ′ ) due to a change in the sourceredshift distribution at co-moving distance χ 0 as Integrating the perturbed lensing spectrum then gives the total perturbation: So far we have treated the function n (i) (z) as being completely free.However, the functional derivative needs to respect the constraint given in Equation (3), thus limiting the possible variations of n (i) (z).
The normalization condition itself is again a function and we write this constraint can be implemented by first defining which will be normalized by construction.n (i) s (z) is a functional of f and we can now evaluate the functional derivative of C[n[ f ]] as an unconstrained derivative but evaluated at f = n.To avoid clutter we ignore the sub-and superscripts in this part one finds where we denote that we want to keep the normalization fixed by the variation δ 1 .This is a very intuitive expression: the first term evaluates the standard functional derivative, while the second term corrects this variation to respect the normalization.

Fisher forecast
The next step is to set some requirements for the lensing power spectra.Here we will look at the difference in the χ 2 , assuming a Gaussian likelihood and thus setting a lower limit on the required accuracy of n (i)  s (z).For modes a ℓm with zero mean and covariance C ℓ , the ∆χ 2 between multipoles ℓ min and ℓ max can be written as note that C ℓ is the matrix with the components C (i j) .The factor f sky takes into account the observed sky fraction.Using Equation ( 8) we rewrite the previous equation as a Riemann sum with the measure Dχ r .If we define the Fisher matrix in this case as: where we labelled n (i) (χ r ) → n α , we recover for a difference in χ 2 using a scalar product on the finite-dimensional Hilbert space of shifts in the redshift distribution where the Fisher matrix acts as a norm-inducing metric where ∆n is the vector containing shifts of the components n α .The Fisher matrix, Equation ( 16), describes, how well the shifts n α can be determined by a measurement of the angular power spectra C α given certain survey settings.If one would try to measure all possible perturbations, neighbouring δn(χ) are strongly correlated.This is, however not the question we would like to ask in this work.Instead, we want to look at the situation which we allow any perturbation ∆n, irrespective of the correlation.Therefore, by turning this argument around, we only use the diagonal part of the Fisher matrix.Lastly one should note that the functional derivative is strictly defined as a limiting process for infinitesimally small perturbation to the function at hand.The relation in general can be non-linear, but as long as relative perturbations to the function are small with respect to unity, these non-linear contributions are sub-dominant.Especially for surveys with tight requirements on the SRDs this is essentially always fulfilled.

Allowed Perturbations to the Source Redshift Distribution
First, we will look at the allowed perturbations to the SRD by allowing for a total ∆χ 2 of unity, corresponding to a one σ shift of a linear model parameter.Clearly, there are many different solutions ∆n that satisfy ∆χ 2 = 1 subject to Equation (17).To show the structure of the Fisher matrix we therefore distribute the allowed ∆χ 2 per ∆n α equally.We will assume EUCLID specifications for the survey as given in Blanchard et al. (2020) and assume n tomo = 10 tomographic bins, a sky fraction of 0.3.Furthermore, we will collect multipoles between ℓ min = 10 and ℓ max = 3000.We then calculate the diagonal Fisher matrix from Equation ( 16) and distribute the errors equally as described above.This results in a possible realisation of ∆n yielding ∆χ 2 = 1 subject to the constraint Equation (3).
Figure 1 shows the resulting perturbed SRDs.The solid lines show the fiducial SRD, while the shaded areas show the allowed perturbations to not cause a bias of more than 1 σ for a linear model parameter.Lastly, the tomographic bin index is shown as a colour bar.The general trend is very clear, the allowed perturbations become very large around a small interval ∆χ around the mean of the distributions.For most tomographic bins this coincides with the peak of the distribution as they are very close to Gaussian.Only for the first and the last bin, these spikes are a bit offset since the distributions are a bit more asymmetric.This already confirms that the most important part of the SRDs in cosmic shear measurements is to calibrate the mean redshift of each tomographic bin very well.Furthermore, we observe that the spikes tend to be narrower at higher redshifts, indicating that the uncertainty on the mean of the SRD is more important at higher redshifts (see also fig.2).We want to stress again, that this is just one realization of ∆n that produces a ∆χ 2 = 1, but by distributing the errors equally, it is possible to see, which perturbations the final measurement is most sensitive to.However, the uncertainties should not be taken at face value and are extreme values, they just give a general trend.Furthermore, it is also important to notice that we treat the variation at each co-moving distance (or redshift) as independent.While the exact spacing does not affect the results as long as the Riemann sum, Equation (15), converges, an issue can arise, however, when inverting the Fisher matrix, Equation ( 16), for finely sampled data points as it can become degenerate.This issue was discussed as well in Kuntz et al. (2023) and requires the definition of an equivalence class by restricting possible variations of ∆n to the quotient space of those functions (SRDs) which the data can distinguish.Having said this, Figure 1, should be understood as a sensitivity scan of the angular power spectra to the SRD only.
Next, we use the perturbed SRDs to calculate their central moments µ n : for a probability distribution function p(x) with mean µ.The perturbed SRDs are used to calculate the change in the central moments relative to the fiducial SRD. Figure 2 shows the resulting relative change for all tomographic bins as a function of the order of the central moment.Clearly, the first moment is most important and while the second one still needs to be known at a 10% level, all higher-order moments are essentially unimportant.This is of course reminiscent of the behaviour observed in Figure 1, where the perturbations are such that they essentially fix the mean.It is of course entirely possible, that we alter the shape of the distribution differently, but still achieve the desired accuracy.As mentioned before, there are perturbations to the SRD which can source larger changes in the moments, what we would like to show here is the relative importance of the moments.How an ensemble of perturbations to the SRD will affect an actual cosmological analysis is studied in the next section, where the correlations at different redshifts are taken into account.
Nonetheless, the results show that for the SRD for cosmic shear, only the mean redshift and the width are important with the former influencing the result way more (by over an order of mag- nitude).In Appendix C we sample from the allowed changes in the SRD and show the relative difference of the first two moments to illustrate their scatter.

Propagating Redshift Errors
In this section, we will revisit the KV450 data for the SRD (Hildebrandt et al. 2020).This data set is used since it includes a covariance matrix from the direct calibration (DIR), for the clustering redshifts (van den Busch et al. 2020) or the self-organising maps (Wright et al. 2020) no bootstrap covariance was estimated so far.
For completeness, the allowed perturbations are shown in Figure 3. Due to the lower signal-to-noise ratio of the measurement, the allowed perturbations are much larger than in the previous case.The features, however, are very similar.
Since we are expressing everything in co-moving distance, the covariance matrix needs to be transformed accordingly.Let C KV450 n(z) be the covariance matrix in n(z) space, the transformed covariance is then where J is the Jacobian with components J i j = δ i j dz/dχ.Alternatively, the Fisher matrix of the SRD perturbations can be expressed in redshift space by the inverse transform.
Perturbations ∆n are now sampled from C KV450 n(χ) and propagated to obtain ∆χ 2 according to Equation ( 14).If the redshift errors as given in C KV450 n(χ) are sufficiently small to not produce a significant bias in the cosmological parameters such as S 8 we expect most realisations (i.e.68% Hildebrandt et al. 2020) to yield ∆χ 2 < 1. Figure 4 shows the resulting distribution in ∆χ 2 for the 10 6 realizations of ∆n for KV450.The vertical dashed lines show the 50th, 68th and 95th percentile.It is clear from this plot that the precision of the SRD used in KV450 is high enough to not yield any spurious detection in the final parameter constraints since the 68th percentile is still well below unity.This is also in agreement with Figure 6 in Hildebrandt et al. (2020), showing almost identical S 8 results when ignoring the uncertainties in the SRD (the shift parameters in this case).
Before moving on, two comments about the limitations of the method are in order: (i) The produced realisations from the KV-450 covariance matrix produce perturbations ∆n(z) which are not necessarily small.Therefore the functional derivative is just an first-order approximation to the non-linear dependency on the SRD fluctuation.However, as surveys become more constraining the requirements on the SRD precision become tighter as well, making the method more applicable in the future.Furthermore, using just a linear model is a fair comparison to what has been done in KV-450, where also Gaussian error propagation of the mean shifts was used to estimate the induced bias on cosmological parameters.
(ii) As mentioned earlier we collect multipoles from 10 to 3000 for the signal of EUCLID.KV-450 on the other hand only measures real space correlation functions over a finite angular range.In principle correlation functions in configuration space get contributions from all multipoles at every angular scale (see appendix D).Here we choose ℓ min and ℓ max such that they match the inverse angular scales used for the KV-450, [0.5, 300] arcmin, effectively amounting to the multipole range used for the KiDS bandpower analysis ℓ ∈ [100, 1500] (Joachimi et al. 2021).In general, the choice of the scales involved in the analysis will change the requirement on SRD uncertainties slightly since different scales obtain most of their signal from different redshifts.
One could now further propagate these uncertainties into cosmological parameters using the corresponding Fisher matrix.For a given shift in the SRD ∆n, the corresponding shifts in the cosmological parameters, ∆θ can be calculated: where Greek indices run over the perturbations in the SRD, while Latin indices label cosmological parameters.Here we assumed the sum convention.F i α hence is the mixed pseudo Fisher matrix: and its inverse is a pseudo inverse.Since the inversion of this matrix is not necessarily stable we choose to go another route here.
Since the distribution of ∆χ 2 is known, we are interested in samples of cosmological parameters with the same ∆χ 2 with respect to the best-fit value.For a Gaussian posterior in one dimension this would amount to a distribution such that the absolute value of each sample is fixed to √ ∆θ 2 .We sample from a standard Gaussian distribution and modify its width by √ ∆θ 2 .This Gaussian is then mapped into the frame of the cosmological parameters under consideration via the Cholesky decomposition of the Fisher matrix of the cosmological parameters.In Figure 5 we apply this procedure to the ∆χ 2 distribution of KV450 (Figure 4).Each dot represents one sample of the ∆χ 2 distribution with its value shown as a colour bar.It can be seen as the geodesic distance to the fiducial value for the cosmological parameters in the parameter manifold (Giesel et al. 2021).The red contours depict the expected 1, 2, 3σ confidence regions from the Fisher forecast for KV450.Since in the original analysis more than the two parameters here were used, we re-scale the ∆χ 2 accordingly, in particular by the χ 2 quantile function χ 2 k (p), where k = 10 is the number of parameters in the actual analysis (Hildebrandt et al. 2017a) and p = 0.68.This is done to obtain a fair comparison.It is clear from the plot, that all samples for the photometric redshift distribution lie well within the 1σ contour.Furthermore, it should be noted that we are considering a very idealised forecast with two free parameters and no systematics here.The procedure, however, can be generalized to any number of parameters.Furthermore, one can apply the same analysis to a full Monte-Carlo-Markov-Chain (MCMC) by matching those samples which are ∆χ 2 away from the maximum likelihood of the MCMC.Lastly, the samples from Figure 5 can be mapped to S 8 = σ 8 √ Ω m0 /0.3.Figure 6 shows the resulting histogram of the scatter due to the photo-z uncertainties.Comparing this to ∆S 8 = 0.076 at 68% confidence (Hildebrandt et al. 2020) shows that the scatter induced by the redshift uncertainties as sampled from the KV450 SRD covariance have a small effect on the overall error budget.In Hildebrandt et al. (2017a) a Fisher matrix method for the shifts of the mean of the SRDs was investigated as a source of systematics, which found similar results to the ones presented here.The main difference between the two methods is that we allow for general perturbations to the redshift distribution (provided their correlation is given).Generalizing the procedure in Hildebrandt et al. (2017a) to moments higher than the variance is bound to fail (see Appendix A).However, we would also conclude that even for EUCLID, the analysis of the first two moments is probably sufficient.
In appendix C the mean and standard deviation of each SRD in the five tomographic bins are shown for the realisations used in this section as sampled from the DIR covariance matrix.Figure C1 shows a very similar behaviour to what we found in Figure 2. In particular we find that the mean scatters less at higher redshifts, while the standard deviation scatters roughly equally for most of the bins.
We close the section with a general discussion about the usage of ∆χ 2 or directly assessing uncertainties in cosmological parameters.It is in general advantageous to make accuracy assessments for the SRD using the ∆χ 2 and not by inverting the Fisher matrix for the parameters of interest to obtain the shift values for those.The reason for this is that ∆χ 2 is an invariant quantity, while shifts in parameter space are dependent on the specific model choice.The only caveat in the ∆χ 2 is that the number of parameters must be taken into account, this is, however, much easier than calculating the Fisher matrix.

CONCLUSIONS
In this paper, we have analysed the dependence of the cosmic shear angular power spectrum on the SRD.This has been done by employing functional derivatives of the cosmic shear C ℓ with respect to the SRD at a fixed co-moving distance χ 0 .By integrating over the introduced error we estimated the ∆χ 2 introduced by arbitrary uncertainties in the SRD.We applied our method to a cosmic shear √ Ω m0 /0.3 parameter.This is directly derived from the samples of Figure 5.The scatter is roughly 15 per-cent of the statistical error budget reported in Hildebrandt et al. (2017aHildebrandt et al. ( , 2020)).
survey with EUCLID specifications and KV450 since a covariance of the SRD estimate was given.Our main findings can be summarised as follows: (i) Allowed perturbations of the SRD are such that they preserve the mean of the underlying distribution.If they do, they can be rather larger, even for a survey like EUCLID.This is in line with the common practice of using only shifted means of the underlying redshift distribution.
(ii) To achieve the accuracy required for EUCLID, the mean of the redshift distribution needs to be determined between 1 and 0.01 per-cent, depending on the tomographic bin under consideration.The variance of the SRD is still important at the 10 per-cent level.There is still some sensitivity left in the skewness, but all other moments are not relevant.
(iii) We performed a simplistic analysis of the KV450 SRDs to check whether they fulfil the requirements and found that the uncertainties, in this very idealised scenario, only yield biases up to 1σ in the final constraints.In a full analysis, this bias would be even smaller.Thus confirming the redshift calibration used in KV450.
(iv) Even for EUCLID it is most likely not necessary to investigate moments of the redshift distribution n > 2. This conclusion could change for different settings and self-calibration methods.
(v) The procedure outlined here has the advantage of being very cheap computationally since the functional derivatives only need to be computed once.It is then only a matter of sampling from the underlying SRD and propagating these perturbations with the previously calculated functional derivative.It is hence not necessary to push thousands of realisations of the SRD through the analysis pipeline.
The method outlined here can thus be used to analyse whether a perturbation in the SRD still fulfils the requirements of a given experiment so that no biases of model parameters are introduced.It allows for arbitrary perturbations to the SRD without requiring a fit to the actual distribution.We intend to apply the presented method to the updated SRDs of KiDS in the future.
For the interested reader the appendices Appendix A -Appendix E discuss various aspects of the analysis which could be refined in future work.In particular, we look at the Edgeworth expansion of the SRD in Appendix A, i.e. an expansion in the cu-mulants of the underlying SRDs.However, we find that, even for a realistic setting, the Edgeworth expansion cannot reproduce the original SRDs if cumulants n > 2 are considered.If a distribution is then expanded as an asymptotic series around a normal distribution one finds where λ n κ n /κ n/2 2 .We are now interested in the sensitivity of the distribution with respect to its cumulants.Here the cases n = 1, 2 are a bit special: and for κ 2 where where we also defined the product: For all cumulants with n ≥ 3 one finds: It should be noted, however, that the Edgeworth expansion is not a convergent series but rather an asymptotic expansion.One therefore needs to check whether the expansion is a good approximation of the underlying distribution.
In this case, one can define the ordinary Fisher matrix using partial derivatives: where κ (i) m is the m-th cumulant of the source-redshift distribution in the i-th tomographic bin.
Figure A1 shows the fiducial redshift distributions for EU-CLID and their Edgeworth expanded approximations as solid and dashed lines respectively.The top plot uses the expansion up to κ 3 , while the bottom plot sums contributions up to κ 6 .For all but the first and last tomographic bin, the Edgeworth series is a good approximation.This is expected as they are essentially Gaussian and therefore κ n ≈ 0 for n > 2. The first tomographic bin experiences boundary effects at z = 0 and is therefore slightly skewed.This effect is even larger for the last tomographic bin, which has a very long tail to high redshifts.While the first bin can still be described by the Edgeworth expansion and the series converges, the 10th bin shows negative probability in the Edgeworth series already at third order.The situation becomes worse if higher-order cumulants are included.This goes to show that even for such an idealized case as the EUCLID forecast, the use of the Edgeworth expansion can be very dangerous.
For the case n = 2 we show the Pearson correlation coefficient of the joint covariance matrix between the first three cumulants in each tomographic bin and four cosmological parameters in Figure A2.We observe some correlations between the first and second moments of each tomographic bin.There is a very strong correlation between the first and second moment of two different redshift bins.Furthermore, one can see that parameters controlling the amplitude of the lensing spectrum are anti-correlated with the mean.We want to stress again, however, that the expansion, even in this case, is not convergent and results obtained with n > 2 have thus to be taken with care.

APPENDIX B: PHOTOMETRIC GALAXY CLUSTERING
For photometric galaxy clustering, the procedure can be simply adopted by changing the weight function (up to galaxy bias, which we absorb in the power spectrum).Again by using the Limber projection: with the galaxy power spectrum P gg and corresponding weights given by: therefore the functional derivative takes the very simple form We show the corresponding allowed perturbations (cf. Figure 1) in Figure B1 where the different shape compared to the cosmic shear case is clearly visible.While the latter had a distinctive peak around the mean of the distributions, the allowed perturbations of galaxy clustering have two maxima, indicating that also the width of the distribution is important.Overall the allowed perturbations are much smaller, this is, however, because the clustering signal is much larger in general.

APPENDIX C: DISTRIBUTION OF THE MEAN AND VARIANCE
We show the relative difference between the mean redshift and the standard deviation of the SRD for each tomographic bin.As before we distinguish between the EUCLID's survey settings and KV450.In particular, we sample from the diagonal covariance obtained from the functional Fisher matrix as described in Section 3 for the former, while we use the DIR covariance for the latter.
The top plots of Figure C1 show the distribution of the mean and the standard deviation and show generally good agreement with Figure 2, that is that the mean must be known below the per-cent level for most bins, while the standard deviation needs to be determined by roughly 10 per-cent.It should be noted that Figure 2 considers the extreme case where we exactly look at the envelope shown in Figure 1.
Finally, the bottom two plots show the same for KV450, where we find much wider errors on mean and standard deviation, a few per-cent and a few ten per-cent respectively.The general trend, however, is the same -high redshift bins are more important than lower redshift bins.

APPENDIX D: RELATIONSHIPS TO OBSERVABLES
Real surveys usually do not use the angular power spectra as a final statistic.This is for example due to incomplete sky coverage, masking effects, variable depth or simply the dimensionality of the data vector.All these factors require a sufficient summary statistic.Very commonly used ones are the correlation function or band powers (or similarly pseudo-C ℓ ).All of these are essentially linear transformations of the pure angular power spectrum C ℓ and assume the following general form: where O is some observable of interest and W O (ℓ) is the associated kernel defining the transformation.Again by the chain rule, the functional derivative of this new observable with respect to the SRD is readily available: where we dropped all the indices for less clutter.For band powers, C l , this would for example assume the following form: where S ℓ is the band power response function and N l is the normalisation.For the two-point correlation function ξ ± one finds: δξ ± (θ) δn(χ 0 ) = 1 2π dℓℓJ 0,4 (ℓθ) δC ℓ δn(χ 0 ) . (D4)

APPENDIX E: NON-LIMBER C ℓ
The Limber projection used for the C ℓ is not valid on large angular scales, where it must be replaced by the full expression.In full generality, for any tracers i and j of the matter density where the functional I k,i [n i ] is given by Here P ii is the auto power spectrum of the tracer i and W i is its associated weight.Thus we find: δC i j ℓ δn a = dk k 2 I ℓ,k,i δI ℓ,k, j δn a δ K ja + I ℓ,k, j δI ℓ,k,i δn a δ K ia , (E3) where δI ℓ,k,i δn i = dχ δW i δn i P ii (k, χ) j ℓ (χk) . (E4) The derivative of the weight function is calculated as before.

APPENDIX F: INTRINSIC ALIGNMENTS
In this work, we have ignored intrinsic alignments (IA).Its inclusion is, however, straightforward by noting that the IA angular power spectrum is simply given by where P II is the IA power spectrum, which summarises the reaction of galaxy shapes to the ambient LSS on the two-point level.
The functional derivative there proceeds in the same way as in Appendix B. For the GI term of intrinsic alignments, one proceeds as before for cosmic shear (compare Section 2).

Figure 1 .
Figure 1.Allowed perturbation for EUCLID to the SRD of the ten tomographic source bins.Solid lines show the fiducial SRD, while the bands show the allowed perturbation to it, corresponding to a ∆χ 2 = 1.

Figure 2 .Figure 3 .
Figure 2. Allowed relative change in per-cent of the central moment of the SRD in each tomographic bin.The changes are calculated from the perturbed SRD distributions as shown in Figure 1.

Figure 5 .
Figure 5.The scatter points show the induced shifts by the photo-z uncertainty in the Ω m0 -σ 8 -plane, derived from the ∆χ 2 (colour-coded) of Figure 4.In red we show the contour from the Fisher matrix for KV450 enclosing the 1σ confidence interval.

Figure A1 .
Figure A1.SRD for EUCLID in all 10 tomographic bins.Solid lines represent the fiducial SRD, while dashed lines represent their respective Edgeworth expansion.Cumulants up to order n = 3 , 6 are used respectively.

Figure B1 .
Figure B1.Allowed perturbations for the redshift distribution for photometric galaxy clustering.

Figure C1 .
Figure C1.Distribution of the relative deviation of the mean and the variance of the SRD, n (i)s (z).Top: For the EUCLID survey settings with realisations from the inverse of the diagonal Fisher matrix used in Figure1.Bottom: For KV450 using the samples generated from the DIR bootstrap covariance.