A highly precise shear bias estimator independent of the measured shape noise

We present a new method to estimate shear measurement bias in image simulations that significantly improves the precision with respect to current techniques. Our method is based on measuring the shear response for individual images. We generated sheared versions of the same image to measure how the galaxy shape changes with the small applied shear. This shear response is the multiplicative shear bias for each image. In addition, we also measured the individual additive bias. Using the same noise realizations for each sheared version allows us to compute the shear response at very high precision. The estimated shear bias of a sample of galaxies is then the average of the individual measurements. The precision of this method leads to an improvement with respect to previous methods concerned with the precision of estimates of multiplicative bias since our method is not affected by noise from shape measurements, which until now has been the dominant uncertainty. As a consequence, the method does not require shape-noise suppression for a precise estimation of shear multiplicative bias. Our method can be readily used for numerous applications such as shear measurement validation and calibration, reducing the number of necessary simulated images by a few orders of magnitude to achieve the same precision.


Introduction
Upcoming weak-lensing surveys have the goal to measure cosmology with unprecedentedly high precision. Their very high statistical power requires systematic errors to be very well understood and calibrated. One of the main sources of systematic error for weak gravitational lensing is the bias in the measurement of the galaxy shear, which carries the cosmological information about the large-scale structure and its evolution. For upcoming experiments such as Euclid (Laureijs et al. 2011), LSST (LSST Science Collaboration et al. 2009), or WFIRST (Spergel et al. 2013), we need to calibrate shear biases to sub-percent precision. Traditional calibration methods create large suites of galaxy image simulations and estimate the shear bias for a given shape measurement method, PSF, galaxy population, noise level, etc. Shear bias estimation to date is dominated by the intrinsic ellipticity dispersion of the simulated galaxies. To reach the desired precision requires a simulation volume that exceeds the actual observational data that are to be calibrated, with billions of simulated galaxies. This has a dramatic impact on the computation load of both generating the simulations and shape measurement methods, therefore on limiting the complexity, storage, and re-usability (e.g. Hoekstra et al. 2017).
Existing methods to reduce the number of simulations made use of rotated galaxies with the same shear such that their mean intrinsic ellipticity cancels. The first proposed method was the so-called ring test (Nakajima & Bernstein 2007), using galaxies evenly distributed in their orientation of intrinsic ellipticity with constant modulus, and with constant shear. Massey et al. (2007) reduced the number of objects to a pair of orthogonally oriented galaxies. Both approaches result in a zero net intrinsic ellipticity. This reduces the error of the estimated shear, but does not entirely cancel the contribution from intrinsic ellipticity shape noise. Stochasticity in the measurement, e.g. due to pixel noise and the PSF, will perturb the exact shape-noise cancellation. In addition, systematic biases break the input ellipticity symmetry. Sources of these systematic effects are ellipticity bias (the response of the measurement to intrinsic ellipticity) when it depends on the orientation of the galaxy with respect to the coordinate system, the PSF, or the shear, and also selection effects, including non-equal galaxy weights. The selection-induced shear bias calibration for the HSC survey was performed without shape-noise suppression ).
In case of simulating fields with non-constant shear, noise suppression can be achieved by simulating the intrin-Article number, page 1 of 14 arXiv:1806.10537v1 [astro-ph.CO] 27 Jun 2018 A&A proofs: manuscript no. shear_bias_measurements sic shape noise distribution of galaxies as a pure B-mode field. Using an estimate of the E-mode power spectrum or real-space E-mode correlation is then insensitive to to the intrinsic shape noise (Kitching et al. 2010). However, simulating a realistic intrinsic ellipticity distribution inevitably leads to a power leakage from B to E (Mandelbaum et al. 2014). In addition, measurement stochasticity and biases, causing imperfect noise cancellation as described above, also apply in the case of variable shear.
In this paper we propose a new method to estimate the shear bias from simulations that is insensitive to the shape noise coming from both the intrinsic and the measured ellipticity dispersion. This reduces the number of required simulated galaxies by three orders of magnitude. Even though we estimate a bias for each galaxy individually, for calibration purposes we only use the mean bias averaged over a sample of galaxies. This avoids unstable ratios of two noisy quantities.
Our method is inspired by the metacalibration technique (Sheldon & Huff 2017;Huff & Mandelbaum 2017), where the shear bias is estimated as the shape estimator response to a small shear applied directly to the individual images. However, as our method is simulation-based, we do not need to de-convolve and de-noise observed images with an estimated PSF, which is notoriously difficult. We can apply any shear to the simulated image before the PSF convolution step.
This paper is organised as follows. In Sect. 2 we define the basic required concepts. Sect. 3 introduces our new shear bias estimation method and contrasts it with existing ones. In Sect. 4 we analytically compute the precision of our method, and compare it to existing shear bias estimation techniques. In Sect. 5 we describe the galaxy image simulations, which we use in Sect. 6 to test our analytical descriptions, and to compare our method with existing ones. We discuss potential applications of our method in Sect. 7, and give a summary in Sect. 8.

Shear bias
We define multiplicative and additive shear bias for a population of galaxies. Let g α be the shear of a given galaxy, and e obs α its observed ellipticity, where α = 1, 2 stands for the two components of the complex shear and ellipticity. If the mean intrinsic ellipticity e I α of the galaxy sample is zero, we can estimate the mean reduced shear g obs α as the average of the observed ellipticities 1 , This estimator is biased in general, and c α , m α are the ensemble-average additive and multiplicative shear biases (Huterer et al. 2006;Heymans et al. 2006). Here and in the following, we ignore non-linear contributions to the bias.
We also neglect higher-order terms, e.g. the term g * e I in the denominator of the shear estimator introduced in Seitz & Schneider (1997). Further we set the convergence κ = 0 such that the (observable) reduced shear g = γ/(1 − κ) equals the shear γ. Alternatively, we can describe a shear bias for individual galaxies as the response of the observed ellipticity to a small shear distortion (Huff & Mandelbaum 2017;Sheldon & Huff 2017): R is a 2 × 2 matrix whose diagonal (off-diagonal) terms represent the response of the ellipticity measurements to shear changes of the same (opposite) component. An additive shear bias for individual galaxies is defined as when g = 0.
A perfect shear estimation corresponds to R being the unit matrix and a α = 0. If the shape measurement conserves the spin-2 property of ellipticity and shear, R αβ needs to be a combination of a scalar and a spin-4 tensor. If we neglect the latter, the response collapses to a single non-zero number R 11 = R 22 , with R 12 = R 21 = 0.

Shear calibration
A simulation-based calibration of measured shear estimates typically measures ensemble biases m α , c α from a large number of image simulations with different galaxy properties and shear g α via (1). A calibrated shear estimate is then obtained by correcting the observed ellipticities by the ensemble biases, (e obs α − c α )/(1 + m α ), provided the simulated galaxy population matches the data in all relevant properties.
When measuring the individual responses using (2) and (3), an unbiased shear estimator is given as R −1 e obs − a ≈ R −1 Rγ , see Sheldon & Huff (2017). Below we present our method to compute R for each simulated galaxy without being sensitive to shape noise.

Our method: shape-noise-free shear bias estimation
We measure the shear response (2) for individual simulated galaxy images as follows. For each simulated galaxy with given properties, intrinsic ellipticity e I and given shear g (which can but need not be zero), we create additional, sheared versions of the same galaxy. The galaxy is analytically sheared before PSF convolution and noise addition, so that the differences between the images are only coming from the shear. We then approximate the shear response (2) by finite differences, following Huff & Mandelbaum (2017), galaxy. The galaxy is analytically sheared before PSF convolution, pixellisation, and noise addition, thus the differences between the images are only coming from the shear. We then approximate the shear response (2) by finite differences, following Huff & Mandelbaum (2017), where e obs,± α is the measured ellipticity of the image with additional small shear ±∆g α . We therefore create additional sheared images for each original one. With four sheared images we can estimate all components of R for each galaxy. To determine the shear response averaged over a sample of galaxies, we only require two appropriately chosen shear values, or three images in total, see Sect. 5 and App. A for more details.
To further reduce the stochasticity of our response estimator, we use the same noise realisation for all image copies for each galaxy. This guarantees that intrinsic ellipticity cancels exactly for our bias estimator.
When randomising the noise for each image, we obtain the same mean but noisier response values. Keeping the same noise realisation of our five images is not an artificial noise reduction in the bias estimate, it only helps us to obtain a noise-free numerical derivative. The noise properties will be sufficiently well sampled by the different simulated galaxies.
The additive shear bias for each galaxy is measured via (3), on the original, non-sheared image.
In Fig. 1 we show an example of the estimated component of the response, R 11 , for one galaxy image. The finite-difference estimate is insensitive to the shear value as long as it is small, |∆g α | < ∼ 0.05 for α = 1, 2. More details about the robustness of our new estimator are presented in App. A.
From the measurements of individual galaxy shear biases, we estimate the ensemble multiplicative and additive bias of a galaxy population as the average of the individual estimates, respectively R αα and a α . This can be a weighted average if galaxies have different weights.
We ignore the non-diagonal terms of R, as we have found that their contribution averages out to zero if the shear values are symmetrical around zero, see App. A.
We emphasise again that our new bias estimator is not affected by shape noise coming from the intrinsic galaxy ellipticity. This is true not only for the estimated mean, but for the bias distribution. In our method, the intrinsic ellipticity can be considered as just as another property of the galaxy (such as the flux, radius, etc.) and as such affects the shear bias in a deterministic way, but does not contribute to the statistical uncertainty. Therefore, we can obtain a much more precise bias estimation compared to methods that average over observed galaxy ellipticities. Consequently, our method requires a much smaller number of simulated images. This will be quantified in Sect. 4. Note also that our simulations do not require a vanishing mean intrinsic ellipticity, which can be a hurdle when dealing with selection biases or galaxy weights.
In the following two subsections, we review two commonly used calibration methods to estimate the shear bias.

Linear fit estimation
The most common methods to estimate the shear bias in the literature is to perform a linear fit of (1) to simulated sheared galaxy images (e.g. Heymans et al. 2006;Miller et al. 2013;Zuntz et al. 2013;Mandelbaum et al. 2015;Fenech Conti et al. 2017;Huff & Mandelbaum 2017;Hoekstra et al. 2017;Pujol et al. 2017;Zuntz et al. 2017;Mandelbaum et al. 2017). For each galaxy population (e.g. for each bin of given galaxy properties) we obtain the additive and multiplicative biases c α and m α from a linear fit of the measured ellipticities as function of simulated input shear, as illustrated in the top panel of Fig. 2. The error of the parameter estimation can then be obtained by jackknife resampling, and obtaining the distribution of best-fit parameters for each resample.
Alternatively, the straight line can be fitted to the average measured ellipticities for each input shear, e obs α , as shown in the bottom panel of Fig. 2. Both fitting schemes provide consistent values and error bars for the shear bias parameters.

Linear fit estimation with shape-noise suppression
The precision of the linear fitting technique to measure shear bias is limited by shape noise stemming from the intrinsic ellipticity distribution. To beat down this noise requires the use of a very large number of galaxy images. An alternative method to reduce the shape-noise contribution is to force the mean ellipticity to cancel, by simulating orthogonal pairs of galaxy images (Massey et al. 2007;Mandelbaum et al. 2014), As described in Massey et al. (2007), the estimated shear of a pair of orthogonal objects is where e obs α,A and e obs α,B are the observed ellipticities of respectively two orthogonal galaxies, whose intrinsic ellipticities cancel exactly, e I α,A = −e I α,B for both α = 1, 2. The shear bias is then estimated from a linear fit of g obs α as a function of g α . This estimator is an improvement over the simple linear fit reviewed in the previous section, with reduced contribution from shape noise. However, the observed ellipticities in the absence of shear do not cancel in general, due to various effects. First, the stochasticity of the two (assumed independent) ellipticity measurements means that e obs α,A + e obs α,B is a random variable with non-zero dispersion. We model this dispersion in Sect. 4.3. Second, the response of the measurement to ellipticity, or ellipticity bias, that depends on the galaxies orientation, either with respect to the pixel coordinate system or to the PSF, can cause the estimated shear to be biased with respect to g α (Kacprzak et al. 2012;Pujol et al. 2017). Third, selection effects can break the symmetry if one of the two galaxies is missed. This selection can occur at the detection level or the shape measurement stage, both of which can fail for one of the two objects. This could be due to a dependence on the relative orientation of the galaxy with respect to the PSF, or random noise fluctuations in particular in the low-SNR range. Fourth, when accounting for galaxy weights the ellipticity cancellation is broken.
A generalisation of this method consist on simulating sets of n galaxies on a ring with constant |e I |, rotated uniformly such that their mean intrinsic ellipticity is zero (Nakajima & Bernstein 2007). The case with n = 2 corresponds to the case of orthogonal pairs discussed above. In Sect. 4.3 we show that increasing n beyond n = 2 does not reduce the shape noise contribution to the shear bias estimator.

Error estimation
In this section we study and compare the precision and precision of the different shear bias estimators. In this section, a latin index of shear, ellipticity, bias, etc. serve to indicate a galaxy number from a population.

Our method: shape-noise-free shear bias estimation
Each galaxy i with properties P i has a shear response R i estimated as described in Sect. 3.1, from different sheared versions of the original simulated galaxy image with the same noise realisation. We assume that the statistical uncertainty of this measurement given P i is negligible. This is based on the results shown in App. A. The response R i depends deterministically on P i , given by the input parameters of the simulated image, the PSF, and the stochastically of the random processes of the image realisation. The latter in our case is a simple Gaussian pixel noise realisation, but we can easily include other effects such as Poisson noise, cosmic rays, etc. The effects on R from this stochasticity can be measured by repeatedly estimating R i for fixed P i with different noise realisations. This provides us with samples from the probability density function (PDF) of R i (P i ). This PDF defines the uncertainty σ N,α for both components of the estimated shear response due to stochastic effects.
In Fig. 3 we show two examples of this stochasticity coming from noise. We have measured R 10, 000 times for 10, 000 different noise realisations for the two galaxies shown in the figure (see Sect. 5 for details on the simulated images and shape measurement). As before, for each realisation we do not change the noise for the original and the 4 sheared versions of the image. The mean responses R i depend on the galaxy properties P i . In general, the response is further from 1 for small galaxies (the top panel) and closer for large galaxies (bottom panel), and the two response components can be different as in the top panel. These results are consistent with the bias results from Pujol et al. (2017).
The dispersion for each component σ N,α of the response depends on the noise level and on the properties P of the object. The dispersion is generally larger for smaller objects. For our shear estimation method we only measure R αα once per galaxy, which means that each shear response R ααi (P i ) has a stochasticity σ N,αi .
Quantifying σ N,α allows us to estimate the number of galaxies we need to simulate such that the stochasticity is subdominant in the final bias uncertainty originating in  the distribution of galaxy properties P i . To meet an allowed shear bias uncertainty of σ req,α , assuming that all galaxies have the same stochasticity σ N,α (alternatively one can use the mean, or a worst-case value), we would need at least N min ∼ σ 2 N,α /σ 2 req,α image simulations not to be dominated by pixel noise.
In the following, for the calculation of the precision of our estimator, we do not try to disentangle the contributions from noise and galaxy properties.
Our bias estimator m α for a sample of N galaxies is the average of the individual shear responses, The uncertainty of the estimated response is where σ R,α is the standard deviation of the distribution of Analogously, the additive bias is estimated as with uncertainty where now σ a,α corresponds to the dispersion of the additive bias over the galaxy population. Fig. 4 shows the distributions of the of R 11 and a 1 for our sample of simulated images (see in Sect. 5). Note that only the multiplicative bias is insensitive to intrinsic ellipticity noise. The additive bias estimated via (9) is still affected by shape noise.

Linear fit estimation
The observed ellipticity of a galaxy i with properties P i can be defined as where g αi is the shear, and S αi is the stochasticity around the linear regression of the measurement for galaxy i that will be dominated by the intrinsic ellipticity e I αi . We write the dependence of observed to intrinsic ellipticity as S αi = f (e I αi ) with some generic function f . In general, f is not the identity that would represent a perfect measurement. Nor is this relation S αi = R ααi e I αi , since the ellipticity response has been shown to be different from the shear response. Because ellipticity is typically larger than shear, this relation is likely to be non-linear. When comparing the predictions with results from data, we will only make the weak assumption that S α is dominated by e I α . For the linear fit to (10) we use a set of values of g α and e I α , whose distributions have dispersions σ g,α and σ e,α , respectively. In Fig. 5 we show these distributions measured on our simulated images, which we describe in more detail in Sect. 5.
The best values of (1+m α ) and c α obtained from a linear regression fit from equation (10) are given by (Kenney & Keeping 1962) as, (11) Assuming g α = 0, these relations become: (13) We assume that R αα and g α are not correlated, which is a very good approximation since the shear bias is linear with g α . Then, with we find Note that the estimated m α is consistent with our method if S α g α = 0. A correlation between these two quantities would effectively modify the slope of the distribution of equation (10), resulting in a biased estimate of m α . Note that for our method this condition does not need to be fulfilled.
We can estimate the error σ m,α on m α via simple Gaussian error propagation assuming that the uncertainties in R ααi and S αi are uncorrelated. This assumption would be violated if the shape estimator has a shear bias that depends on ellipticity. We test our assumptions and approximations in Sect. 6, where we compare the numerical predictions with measurements from simulated images. The sensitivity of the bias with respect to these two quantities is Replacing for simplicity the individual galaxies' dispersions σ R,αi and σ S,αi by the mean values, we get Compared to (7) this expressions shows the additional term σ 2 S,α /σ 2 g,α . In most scenarios this is indeed the dominant term for the bias dispersion, which is the main reason why the linear fit achieves a much lower precision in bias estimation compared to our method.
The uncertainty on the additive bias comes directly from the dispersion in the stochasticity, 4.3. Linear fit with shape-noise suppression Here we estimate the uncertainty of the shape-noise suppression estimator (5), which we write similar to (10) as The difference to equation (10) is that the index i now denotes a pair of orthogonal galaxies. The stochasticity depends on the sum of the observed ellipticities of orthogonal pair, Under the scenario of a perfect shape estimator the sum vanishes exactly. However, a shape estimator typically has a non-zero ellipticity bias, for X = A, B, and g α = 0. If the ellipticity bias depends on the galaxy orientation, or the relative orientation between galaxy and PSF or shear, the two bias values b α,A and b α,B are in general not equal, and we find The shear bias uncertainties σ m,α and σ c,α are computed via equations (19) and (20) derived in the previous section, but with σ S,α given by the dispersion of (24). This is a clear improvement, since the pre-factor |b αi,A − b αi,B | can be expected to be smaller than unity. In addition, if the noise realisation is different for each of the objects A and B, this measurement is stochastic even if b αi,A = b αi,B . This stochasticity contributes to σ m,α and σ c,α , which we denote with σ e obs ,α . In the general ring estimator case where we simulate n rotated copies of each galaxy to suppress shape noise, with n j=1 e I j = 0, we can write: Keeping the total number of galaxies used in the linear fit constant, which is now N/n, we get and We can see that forcing shape noise suppression gives a more precise m than the simple linear fit as far as σ e obs ,α σ e,α . We also see that σ c,α does not depend on the number of galaxies used for the shape noise suppression, but σ m,α increases with n. However, in our derivation we neglected the higher-order contributions in the shear estimator (1), which decrease with n. In this paper we do not quantify the optimal n that minimises both contributions, since the second term in (26) dominates over these other two quantities. In conclusion, n does not affect significantly σ m,α .

Simulations
For this analysis we used the public software package Galsim  to generate isolated images of two million galaxies, corresponding to the Control-Space-Constant branch of the GREAT3 challenge (Mandelbaum et al. 2014. The images are organised in 200 fields, each field with a unique PSF and shear (both constant for each field). The galaxy light distribution follows either a single Sérsic profile or a de Vaucouleurs bulge plus exponential disk. Each galaxy is simulated twice, the second one being rotated by 90 degree with respect to the first one to achieve shape noise suppression. For more details about the simulated images, we refer the reader to Pujol et al. (2017) as well as Mandelbaum et al. (2014). This set of simulations are used for the linear fit methods, with (Sect. 3.3) and without (Sect. 3.2) shape-noise suppression. For the latter, we average the observed ellipticity of all galaxies for a given shear g, not specifically accounting for the orthogonal pairs. This results in a mean ellipticity each bin close to zero, but does not reduce the scatter due to the intrinsic shape noise.
For our estimations of R as described in Sect. 3.1, we simulate the 2 million galaxies three times, with two sheared values drawn from the cases g = (±0.02, 0), g = (0, ±0.02).
The two values chosen have to be different in the both components in order to be able to estimate R 11,22 . Since both components of g change for each of the shear versions, the estimation of R αα is affected by the non-diagonal terms as follows: Our estimation of R αα is unbiased as far as ∆g β = 0, for β = α, over the entire sample. This is the case since we choose the sign of the shear changes at random. We can also measure the non-diagonal terms of R by using three images with shear values from g = (0, ±0.02), g = (±0.02, 0) and g = (0, 0) (see Sect. A for more details).
Galaxy shapes are obtained with the KSB method, using the publicly available code shapelens (Viola et al. 2011). This method estimates the ellipticity of the objects from the surface brightness moments defining the ellipticity as: The implementation details of the shape measurement algorithm are not very relevant for this paper, and we refer the reader to Pujol et al. (2017) where we used the same methodology.

Results
In Fig. 6 we compare the shear bias obtained with our method to the linear fit technique. As example of galaxy property we use the input disk flux F d of the simulated bulge+disk galaxies. As shown in the top panel of that figure, both methods give consistent results when using all two million galaxies. However, our method estimates the biases with a significantly better precision. The location of the points on the x-axis corresponds to the centre of the F d bins. In addition to a small shift that we apply for an easier visual comparison, the bin centres for our method in the lower panel are modified, since the galaxies are now a random subsample. It is remarkable that when using all two million galaxies, the curves of m 1 and m 2 for our method are almost identical. We quantify the precision of the different shear bias estimation methods in Fig. 7. as a function of the number of simulated galaxies N sim . We create different random subsets of galaxies with size N sim , and measure for each subset the shear bias for the three methods as described in Sect. 3. Our method Equation (7) Linear fit Equation (19) Shape noise suppression Equation (26) 10 2 10 3 10 4 10 5 10 6 10 7 N sim 10 -4 10 -3 10 -2 10 -1 σ c, 1 Our method Equation (9) Linear fit Equation (20) Shape noise suppression Equation (27)  We compute the RMS for each sub-set by jackknife resampling of the input galaxies for all methods, using 50 subsamples (other numbers of subsamples have given the same results).
We compare these uncertainties as measured from the simulations to the numerical predictions derived in Sect. 4. For the latter, we measure the parameters σ R,α , σ a,α , σ S,α , σ e obs ,α , and σ g,α directly from the simulations, as illustrated in Figs. 4 and 5. The amplitude and N −1/2 sim -dependence of the uncertainty measured from the data shows excellent agreement with the analytical calculations for all three methods. This suggests that the assumptions we made to derive these expressions are valid for the system and regime studied here. For the linear fit predictions, we set σ S,α = σ I e,α , assuming that stochasticity S α is entirely determined by the intrinsic ellipticity. For the linear fit with shape noise suppression, we measure σ e out ,α directly from the distribution of the sum of observed ellipticities of the orthogonal pairs, (e obs A,α + e obs B,α )/2. Our method has a much higher precision on the multiplicative shear bias estimation. Compared to the lin-ear fit, σ m,α for our method is smaller by a factor of 35.9. This means that for this study our method requires 35.9 2 /n ∼ 1300/n times fewer simulated images to obtain the same precision, where n is the number of sheared versions used for each object. In Fig. 7 we show our method with n = 2, where we used shear values symmetrically distributed around zero, but similar results have been found for n = 4.
We demonstrate the high precision of our method in the bottom panel of Fig. 6, where we estimate the shear bias as in the top panel, but now for our method with only a fraction of 1/1300 of the objects, chosen at random. The results are consistent in both mean and error bar, demonstrating that our method reaches precision and precision of existing methods with three orders of magnitude fewer simulations. Note that some of the noise in the data points for our method comes from the more sparely sampled galaxy properties in each bin.
In the case of Euclid, with a global requirement of σ m,α < 2 × 10 −3 one needs at least 2 × 10 7 images for the linear fit method, but only ∼ 10 4 for our method according to this study.
The ratio of the RMS between the two methods is approximately σ e,α /(σ R,α σ g,α ) > 1. The quantities σ e,α and σ g,α in the simulation need to be chosen to match expectations from cosmology and galaxy morphology. Given some basic survey characteristics such as redshift and wavelength coverage, and the survey selection function, these fundamental quantities are fixed. The dispersion σ R,α however strongly depends on instrumental effects such as the PSF size, and on the shape estimator. In this study study we used a KSB method to measure the shapes on GREAT3-CSC-like images, and we expect σ R,α to change when using other simulations and shape estimators. This provides a strong motivation to chose or develop a shape measurement method that minimises this dispersion, and therefore minimises the number of required simulations for calibration.
Applying the shape-noise suppression with orthogonal pairs improves the precision with respect to the simple linear fit by a factor of ∼ 2.8 for the measurements of both the multiplicative and additive bias. This improvement reduces by a factor of ∼ 8 the number of simulated images required for the same level of precision. Note, however, that each galaxy needs to be simulated twice for the shape-noise suppression. This is consistent with the factor of ∼ 9 found in Fenech Conti et al. (2017), where they used n = 4 for the shape noise suppression.
Comparing our method to the linear fit with shape-noise suppression, we obtain an improvement of a factor 12.8 for the multiplicative shear bias. This implies that for the same level of precision we can reduce the number of simulated images required by a factor of 12.8 2 /n ∼ 164/n . When comparing the additive bias precision, our method shows a factor of 2.26 improvement with respect to the linear fit, and a factor 0.56 with respect to the shapenoise suppression. Shape-noise suppression performs better because the additive bias is the average ellipticity over all simulated images, while for our method only 1/n images are used. In principle we could estimate c α with n = 1 (e.g. using only the original image). This would however unfairly not count the n > 1 images we have to simulate to measure m α . We conclude that a similar precision are obtained for both our method and shape-noise suppression when estimating additive shear bias.
For the method with shape noise suppression we have assumed b αi,A = b αi,B , but we can have different ellipticity bias for the orthogonal pairs, since their alignments with respect to the shear or PSF is different in general, and thus their observed ellipticity modulus will be different We have measured the differences in the ellipticity bias as a function of the relative orientation with the shear, finding that orthogonal pairs can have differences in their ellipticity biases up to 2 percent. This would also reduce the precision of this method. Our numerical results suggest that this effect is negligible, although its impact depends on the image specifications and the shape estimator used.
Assuming the ellipticity biases to be the same as the shear bias would also be incorrect. For our image simulations we have found an average shear bias of m 1 = −0.049 and an average ellipticity bias of b 1 = −0.124 (similar results are found for the 2nd component). This large difference between those biases confirms the findings from Pujol et al. (2017).
We also remind the reader that any selection effects or weights affecting differently the orthogonal pairs (which can be expected since their orientation relative to the shear is different and then so their global shapes) would increase the uncertainty of the method significantly.

Discussion and applications
The method presented here is a clear improvement on the precision of the shear bias estimation in simulations with respect to the standard linear fit of equation (1). It is also more precise compared to the linear fit with shape-noise suppression via pairs of orthogonally aligned galaxies. In the following we discuss potentially useful applications to improve shear bias analyses with simulated images.

Shear bias validation and calibration
One of the interests of measuring shear bias in simulations is to validate or calibrate the performance of a shear estimation algorithm. In the case upcoming surveys such as Euclid, LSST or WFIRST the requirements on the knowledge of the additive and multiplicative bias imply the generation of a very large volume of simulations which is computationally very challenging. Our method allows to save significant computational efforts to reach these requirements. In our case of study we require 2 − 3 orders of magnitude fewer images to reach the same precision as common approaches, although the exact factor depends on the shear estimator algorithms and the image and survey specifications.

Selection biases
Shear bias from selection effects have been found to be of the same order of magnitude than those induced by the shape measurement process (Fenech Conti et al. 2017;Mandelbaum et al. 2017). Such biases arise when the galaxy selection function depends on the shear. This is for example the case when detection or shape measurement fails for galaxies that are very elliptical, or aligned with the PSF. Such selection effects also arise by imposed, necessary cuts on galaxy properties such as SNR or size, which can favour certain shear values. The resulting shape catalogue then samples the underlying shear field in a non-representative way, which induces biases on the estimated shear if uncorrected.
Our method does not require shape-noise suppression via tuples of galaxies, and is therefore particularly useful when selection effects and weights are to be simulated and studied. Weights can be applied to the simulated galaxies following an arbitrary distribution, to study the impact on shear bias.
Selection effects that are correlated with the shear can be studied as proposed by Sheldon & Huff (2017): we calculate the mean response of (4) by first averaging the ellipticities of both sheared samples before taking the difference and dividing by the small shear, Now, the two galaxy samples giving rise to the mean observed ellipticities e obs,+ α and e obs,− α , respectively, are not only different by their shear. In addition, a given selection criterium (e.g. a minimum SNR) is applied to the two sheared samples. If the applied shears modify the selection, this results in different mean sample ellipticities, and the selection-induced shear bias translates into the shear response (31).
This shear bias estimator however does not account for selection effects that affected the original galaxy sample. These can only be studied by creating different samples with different selection criteria, and then comparing the shear responses between those samples. In that case, selection effects that act differently on the sheared galaxy images are undesired, since they create additional, spurious selection biases. Such differential detection or shape measurement success or failure is however rare, since the shear is very small and thus the images are very similar. Such occurrences can be further reduced: If an image sheared by a value g can not be measured, e.g. because its increased observed ellipticity pushed it under the SNR threshold, the opposite shear −g (e.g. making the galaxy rounder) should not affect the measurement success. In the case of n = 2 images per original galaxy, we are free to choose the sign of the shear as long as the average shear is zero, largely avoiding such selection-induced measurement failures. A detailed study of selection effects and shear bias is left for future investigation.

Individual shear responses
Studying the shear response as function of galaxy properties for individual galaxies without the need to bin or average can have advantages. For calibration, the shear bias as function of galaxy properties is typically modelled as a smooth function, either parametric, e.g. by fitting an analytical, multi-variate function, or non-parametric, e.g. by interpolating the (smoothed) measured bias values.
For linear fit methods to estimate shear bias we need to compute such a function from data binned in galaxy properties. Then the average shear biases are measured for each bin. However, this average values depend on the galaxy population inside the bin, whose shear responses might not only depend on the binned properties, but also on the properties that have not been used in the binning. As a consequence, our measured shear bias dependencies are sensitive to the property distribution of the galaxy population used. Individual shear responses and biases of simulated galaxies can further serve to further learn shear calibration as a complex non-linear function of galaxy and image properties (e.g. using machine learning techniques), where no binning is needed and we can use a larger set of properties so that the function can be less dependent to the population used.

Variable shear and response on shear statistics
So far we have simulated galaxies either with a random or with no shear, before applying a small constant shear. Switching from constant to variable shear can be carried out by imposing a shear field to our simulation. Similarly to the constant shear case, the shear bias is derived by computing the shear response to a small shear power spectrum perturbation. For example with the shear drawn from a Gaussian random field with a certain power spectrum C , The small shear values applied to each galaxy with an intrinsic shear g can be arbitrary, e.g. it can be proportional to g (see discussion in Sect. 7.5). Going a step further, the same methodology can be applied to derive the influence of shear bias on any shear statistics: the shear two-point correlation function, the shear power spectrum, peak counts, mass maps, higher-order statistics, etc. As an example we hereafter illustrate this with the shear two-point correlation function ξ ± .
First, as described above, we apply a shear field to the simulated galaxies, where each Fourier-space shear coefficientsγ is drawn from a normal distribution with zero mean and variance C . Next, we perturb the shear field by drawing new coefficientsγ ∼ N (0, C + δC ), where we change the power spectrum by a small amount δC . From the original and perturbed shear field we compute the statistics of our choice, e.g. the correlation functions ξ ± and ξ + ± , respectively. The difference between both divided by the perturbation is then the response due to the multiplicative shear bias on the correlation function, which would give us information about the spatially varying shear bias.

Non-isolated images
This analysis has been done using isolated galaxy images. To create more realistic simulations with blended galaxy images leads to the problem that the shape of many galaxies is measured in the presence of one or more closeby galaxy at different redshift and therefore different shear, if the simulation presents a realistic cosmic shear field as described in the previous section. The same issue arises if shapes of blended galaxies are estimated jointly. The presence of nearby isophotes of other objects are known to affect the shear bias (Hoekstra et al. 2015. A common procedure to study these effects is by simulating many combinations of blended objects and close neighbours, and measure the impact on the shear bias statistics over different populations. We claim that we can more efficiently account for these effect, since our shapenoise insensitive method does not require to sample the large space of the distribution of N ellipticies and shears p(e n , . . . , e N , g 1 , . . . g N ). One of the questions to address in this situation is how to produce the different sheared versions of the same images. Here we discuss two possibilities: -We change the shear of only one of the galaxies (the target) from the N -tuple of blended images. The inconvenience is that we need to generate N times more images compared to isolated galaxies. -Alternatively, we can shear every member i of the Ntuple. This shear could be a small additive shear, ∆g = const, as applied to isolated images in this paper. Or it could be a function of g, such as a multiplicative factor, ∆g i = Cg i with C 1 = const. In this case we would preserve the proportions between shears for galaxies at different redshifts. This function of g can be chosen taking into account the statistics or cosmological analysis that we want to do, as discussed in Sect. 7.4.
We leave a comparison of these two approaches, and estimation of shear bias for blended objects in general, to future work, which is beyond the scope of this paper.

Future simulation challenges
Adopting future versions of simulation challenges such as the GREAT series (Bridle et al. 2009;Kitching et al. 2010;Mandelbaum et al. 2014) to our method of shear bias estimation can result in a significant decrease of required image simulations. For GREAT3 the total data volume that had to be downloaded by the participants was 6.5 terabyte (10, 000 simulated galaxies × 200 fields × 20 branches). Reducing this number could result in a more accessible and faster to process challenge.
To use our method, for each galaxy two additional, sheared versions of the same galaxy would need to be simulated with the same noise realisation. The challenge organisers would estimate the shear response for each original galaxy via (2), and apply some metric on the distribution, e.g. the mean, to evaluate the submissions.
To guarantee the blind aspect of the challenge, either all codes have to be run on the organisers' server without direct access to the simulation by the participants. (For testing, smaller training sets of simulations could be provided to the teams for download.). Alternatively, the shear values applied to each galaxy have to be random and kept hidden from the users. Note that it would be trivial for participants to identify the sheared versions of each original galaxy since the noise is the same for the sheared versions, even if the image order was randomised.
To take the example of GREAT3, a similar challenge using our bias estimator could reduce the 10, 000 × 200 simulated galaxies for one branch to a few thousand.
If a variable shear field is to be used in the challenge, with a metric operating on the shear correlation function or power spectrum, a similar method as described in Sect. 7.4 can be employed to measure the response to a small and variable shear.

Summary
In this paper we present a new method to estimate shear bias from image simulations. Our estimator of the multiplicative shear bias is not affected by shape noise, removing the dominant uncertainty in bias estimation. Previous methods constrain the multiplicative and additive bias from a linear fit of the observed average ellipticity as a function of shear. The uncertainty of this parameter estimation is dominated by the intrinsic ellipticity distribution. Shapenoise suppression techniques using matched sets of galaxies with net zero intrinsic ellipticity improve the precision of the measurements, but are affected by selection effects, weights, and ellipticity bias that can break the shape-noise suppression.
Our method consists in measuring the shear response and additive bias of individual galaxy images. To that end, we simulate different sheared versions of the same galaxy, and measure the shear response of the image from the numerical derivative of the measured ellipticity with respect to the shear. We also measure the additive bias for the individual images. For each galaxy the sheared version has the same noise realisation, allowing us to determine the response at the level of machine precision. Then, the multiplicative and additive bias of a sample of galaxy images are obtained from the average shear response and additive bias, respectively. This method improves the precision on the estimation of shear bias significantly because it is not affected by shape noise or by the stochastic uncertainty of the measured ellipticity.
Using numerical simulations as well as analytical predictions, we quantified the uncertainty of the shear bias estimation for our method as well as for linear fits. For the multiplicative shear bias, our methods provides a significant decrease in the shear bias error of a factor of ∼ 36 compared to the linear fit, and a factor of ∼ 12 if the latter is used in combination with shape-noise suppression. The additive bias uncertainty improves by about 2.3 over the linear fit, and under-performs only compared to shape-noise suppression, by a factor ∼ 0.5.
This implies that we can reduce the number of simulated images by a factor of ∼ 1, 300 and ∼ 150, respectively, to measure the shear multiplicative bias with the same precision. Our method has the further advantage that it does not need to impose shape-noise suppression, and hence it can easily be applied for analyses where selection biases or weights play an important role. Appendix A: Robustness of Method 1 To measure R for each individual galaxy, we have generated copies of the same images with the different shears specified in Sect. 3.1. In reality, galaxies can have other values of shear, where both shear components can be different than 0 at the same time. Fixing the other component to 0 can be a simplification of the estimation of R. For this reason, here we measure the impact of different estimations of R using different shear values. In particular, we compare the following estimators: -R A 11 : we obtain R 11 from the fit using the shear values g = (±0.02, 0). -R B 11 : we obtain R 11 from the fit using the shear values g = (±0.02, 0) and g = [0, 0]. -R C 11 : we obtain R 11 from the fit using the shear values g = (±0.02, 0) and a random value of g, with both components random.
-R E 11 : we obtain R 11 from the fit using all the previous values and also g = (0, ±0.02).
If the non-diagonal terms of R are non-zero, R A 11 and R B 11 should behave differently than the rest. In Fig. A.1 we show 10 cases of R estimations, where each case is represented with a different colour. The solid lines correspond to the fits of R A 11 (so they connect the points at g 1 = −0.02 with those at g 1 = 0.02). We can see that all the points, even the random ones, tend to be well adjusted to the fitting line, although not always. These cases are an indication of non-diagonal terms of the shear response R, causing changes in e obs 1 when changing g 2 . In this section we focus on the first component of shear and response, but the same holds for the second.
We have found that R A 11 and R B 11 are identical at the level of machine precision. This means that the method is 0.7 0.8 0.9 1.0 1. and R C 11 . The differences with R D 11 and R E 11 are very similar. We see that in most of the cases the differences are negligible, uncorrelated with R A 11 and they average out because of symmetry. We have found that, when the second component of the shear affects e out , it does it in a symmetric way so that positive g 2 give opposite effects than negative g 2 . For a random distribution of g 2 , the differences cancel out. In the top panel of Fig. A.3 we show that the differences between the different estimations are consistent with 0 and independent of the random g applied. In the bottom panel we show that the mean response as a function of the disk flux of the galaxies is consistent for all the estimators. This is actually the case as a function of all the properties studied, which means that the method is quite deterministic and that the non-diagonal terms of the shear response do not affect the shear response estimation. As a consequence, our method does not depend on the different shear values used for the fit as far as the shear values used are symmetric or homogeneous.
In principle our method also suffers from ellipticity bias, since the observed ellipticities of the different sheared images are not equal. They are however very close, unlike in the case of the orthogonal pairs. We thus expect this ellipticity bias to be very small and negligible compared to the 2 percent maximum ellipticity bias for the orthogonal pairs. 0.05 0.00 0.05 g 1, 2 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0 2.5 and R C 11 as a function of the random g1 applied. Bottom panel : resulting average responses from the different estimations as a function of the disk flux for galaxies with a bulge an a disk. Similar results are found for R22.