# Radio weak lensing shear measurement in the visibility domain – II. Source extraction

Radio weak lensing shear measurement in the visibility domain – II. Source extraction Abstract This paper extends the method introduced in Rivi et al. (2016b) to measure galaxy ellipticities in the visibility domain for radio weak lensing surveys. In that paper, we focused on the development and testing of the method for the simple case of individual galaxies located at the phase centre, and proposed to extend it to the realistic case of many sources in the field of view by isolating visibilities of each source with a faceting technique. In this second paper, we present a detailed algorithm for source extraction in the visibility domain and show its effectiveness as a function of the source number density by running simulations of SKA1-MID observations in the band 950–1150 MHz and comparing original and measured values of galaxies’ ellipticities. Shear measurements from a realistic population of 104 galaxies randomly located in a field of view of 1 $$\deg ^2$$ (i.e. the source density expected for the current radio weak lensing survey proposal with SKA1) are also performed. At SNR ≥ 10, the multiplicative bias is only a factor 1.5 worse than what found when analysing individual sources, and is still comparable to the bias values reported for similar measurement methods at optical wavelengths. The additive bias is unchanged from the case of individual sources, but it is significantly larger than typically found in optical surveys. This bias depends on the shape of the uv coverage and we suggest that a uv-plane weighting scheme to produce a more isotropic shape could reduce and control additive bias. gravitational lensing: weak, methods: statistical, techniques: interferometric, cosmology: observations 1 INTRODUCTION Cosmological or targeted surveys of weak gravitational lensing at radio wavelengths may have a relevant role in the next years, when the Square Kilometre Array (SKA)1 radio telescope will start to operate, providing a density of detected galaxies sufficient for shear measurement and a resolution to reliably measure their shapes. They will also be able to probe to higher redshifts given the different galaxy redshift distributions compared to the optical band (Brown et al. 2015). Although the galaxy number density in the radio band will be lower than in the optical, the possibility to observe deeper can make radio weak lensing surveys for cosmology measurements competitive with the corresponding optical surveys, as shown in recent forecasts from SKA simulations (Harrison et al. 2016). Moreover, radio observations have the advantage of a deterministic knowledge of the image-plane point spread function (PSF), being the Fourier Transform of the uv coverage, and will provide unique approaches for mitigating intrinsic alignments, such as concurrent measurements of polarization (Brown & Battye 2011) and galaxy rotation velocities (Blain 2002; Morales 2006). Being subject to different observational systematics, cross-correlation with optical observations of the same field will allow suppression of systematic errors on shear measurement from future large surveys (Patel et al. 2010; Demetroullas & Brown 2016; Camera et al. 2017). This is quite relevant for precision cosmology as these errors may become comparable to, and larger than, the statistical noise. The precursor radio weak lensing survey SuperCLASS2 is already underway and will soon provide data that may be used to test new methods required for accurate galaxy shape measurement in the radio band. A natural approach for such methods is working in the visibility domain where the data originates and the noise is Gaussian, avoiding non-linear data manipulation of the imaging process. SKA simulations have already shown that current imaging methods produce images with structures in the residuals which may dominate the cosmological signal (Patel et al. 2015). Also cross-correlation analysis using real data images shows that no evidence of correlation is found between the optical and radio intrinsic shape of the matched objects (Patel et al. 2010; Tunbridge, Harrison & Brown 2016). This result suggests the presence of systematics in the procedure adopted for turning the visibility data into images, although a significant percentage of AGN sources in the observed population may be another possible explanation, as well as the astrophysical scatter between optical and radio position angle due to the different emission mechanisms in the two bands. Currently, cosmic shear in the radio band has been successfully detected only working in the visibility domain but obtaining sources position from the image (Chang, Refregier & Helfand 2004). Galaxies’ ellipticities from the VLA FIRST survey (Becker, White & Helfand 1995) were measured by decomposing them into shapelets, an orthonormal basis of functions corresponding to perturbations around a circular Gaussian invariant under Fourier transform (Chang & Refregier 2002). Since the data size and the number of resolved sources ($${\sim } 20\hbox{--}30 \deg ^{-2}$$) of each pointing is quite small, a joint fitting of the shapelet coefficients was possible by solving normal equations. Such an approach, computationally convenient, becomes very challenging when dealing with the order of 104 sources per square degree and a very large data set per pointing (order of PetaBytes), as expected from SKA Phase 1 continuum surveys (Brown et al. 2015). Moreover, shapelets introduce a shear bias as they cannot accurately model steep brightness profiles and highly elliptical galaxy shapes (Melchior et al. 2010). In the companion paper (Rivi et al. 2016b, hereafter Paper I), we presented RadioLensfit, an alternative method working in the visibility domain where model fitting is performed on a single source at a time using an exponential profile as model for the galaxy. It is an adaptation of the optical Bayesian lensfit method (Miller et al. 2013) to radio data, where model visibilities are defined analytically and the likelihood is marginalized over uninteresting parameters. The method was tested in the simple case of individual galaxy visibilities simulated adopting the SKA1 uv coverage described in Section 3, and the shear noise bias (Refregier et al. 2012; Melchior & Viola 2012) estimated as a function of the signal-to-noise ratio (SNR). Results compared with requirements (Brown et al. 2015) for the proposed SKA1 radio weak lensing survey (Braun 2014) showed that the multiplicative shear bias may need calibration corrections similar to those for optical surveys, while the additive bias have to be controlled by an isotropic sampling of the visibility plane. In this paper, we extend this work implementing the method for isolating source visibilities from realistic data, i.e. when many galaxies are in the field of view. We estimate its effectiveness in terms of ellipticity fitting and shear measurement by running SKA1-MID simulations as we did in the previous paper. We finally investigate the effect of the shape of the uv coverage on the additive shear bias. This paper is organized as follows. In Section 2, we summarize RadioLensfit and present the extraction algorithm. In Section 3, details of SKA1 simulations are provided, while in Sections 4 and 5 results for galaxy ellipticity and shear measurements are presented respectively. Finally, we discuss the shear additive noise bias in Section 6. 2 OVERVIEW OF RadioLensfit RadioLensfit is a method for measuring radio galaxy ellipticities in the visibility domain. The idea is to adapt the approach used in the optical case to radio data, i.e. extracting from visibilities and model fitting a single source at a time. Source extraction is difficult in the Fourier domain because signals from all sources in the primary beam are mixed altogether in the visibilities and sources are no longer localized. For this reason, a joint analysis with the image domain may be needed: it allows us to identify sources and measure their position and flux with sufficient accuracy. With such information, we can also compute a model of the observed sky and use it to approximate the signal from the other galaxies that must be removed when extracting each source. The extraction is completed using a faceting technique that phase shifts the phase centre to the source position and further reduces its signal contamination by averaging visibilities in a coarse grid. Finally, the model fitting can be performed as for the simple case of a single galaxy in the primary beam located at the phase centre as summarized in Section 2.1. This way we can largely reduce the computational time when a huge number of sources are in the field of view (as for SKA) instead of trying a challenging joint fitting of all sources. A detailed algorithm for the fitting of many sources in the primary beam is presented in Section 2.2. 2.1 Galaxy ellipticity fitting In Paper I, we introduced this method as an adaptation of lensfit (Miller et al. 2013) by performing the chi-square fitting of single source visibilities. They are evaluated at the uv points, that are defined by the baselines formed between two antennas projected on the plane orthogonal to the pointing direction. Model visibilities of a star-forming (SF) galaxy are computed analytically as the Fourier transform of the exponential brightness profile (Sérsic index n = 1):   \begin{eqnarray} V(u,v) = \Big ( \frac{\lambda _\mathrm{ref}}{\lambda }\Big )^\beta \frac{S_{\lambda _\mathrm{ref}}{\boldsymbol {e}}^{2\pi \mathrm{i} (u l + v m)}}{(1+4\pi ^2 \alpha ^2 |\boldsymbol{A}^{-T}\boldsymbol {k}|^2)^{3/2}}, \end{eqnarray} (1)where $$\boldsymbol {k}=(u,v)^T$$ is measured in wavenumber units, β = −0.7 is the assumed spectral index for the synchrotron radiation emitted by the galaxy disc, (l, m) and α are the source position and scalelength, respectively, $$S_{\lambda _\mathrm{ref}}$$ is the source flux at reference wavelength λref. The ellipticity parameters (e1, e2) are contained in the matrix $$\boldsymbol{A}$$ that linearly transforms the circular exponential profile to an ellipse:   \begin{eqnarray} \boldsymbol{A} = \left(\begin{array}{cc}1-e_{1} & -e_{2} \\ -e_{2} & 1+e_{1} \end{array} \right). \end{eqnarray} (2)We assume the following ellipticity definition:   \begin{eqnarray} {\boldsymbol{\it e}} = e_1 +\mathrm{i}e_2 = \frac{a-b}{a+b}\mathrm{e}^{2\mathrm{i}\theta }, \end{eqnarray} (3)where a and b are the galaxy major and minor axes, respectively, and θ is the galaxy orientation. The likelihood is marginalized over non-interesting parameters such as flux, scalelength and position, adopting uniform priors for the flux and position, and a lognormal prior dependent on the flux for the scalelength (see Section 3). This way we obtain a likelihood function of only the ellipticity parameters. The galaxy ellipticity measurement is given by the likelihood mean point and 1D standard deviation (defined as the square root of the covariance matrix determinant), obtained after sampling the likelihood with an adaptive grid covering a neighbourhood around the maximum point. In real observations, the finite channel bandwidth and sampling time introduce smearing effects that are dependent on the source position. These effects may be approximated analytically and included in the visibilities model (Bridle & Schwab 1999; Smirnov 2011). For example, for frequency smearing, assuming a square bandpass filter in the expression of the smeared visibility presented in Bridle & Schwab (1999), we obtain:   \begin{eqnarray} \tilde{V}(u,v)= V(u,v) \mathrm{sinc}[\pi (ul+vm)\Delta \nu /(\nu _0 c)], \end{eqnarray} (4)where uv coordinates are taken at the mid-channel frequency ν0, Δν is the channel bandwidth and sinc(x) = sin (x)/x. Another option is to make the observation with very tiny frequency channels and sampling time intervals. Harrison & Brown (2015) proposed to use ∼30 kHz channel bandwidth and 0.5 s sampling time to make smearing tolerable, but meaning a huge number of uv points. In this case, raw visibilities may be averaged into a single uv grid without jeopardizing ellipticity measurements. In fact, observations from the same pair of antennas at different frequencies (resp. times) correspond to visibilities evaluated at different uv points along a radial (resp. tangential) direction, therefore these visibilities can be treated as the ones evaluated at uv points related to different baselines. Depending on the grid size, data volume and then computational time may be considerably reduced. 2.2 Source extraction We assume flux and source positions are provided. For example, they may be estimated from a cleaned image of the same data that are analysed, or applying MC methods to the visibilities (e.g. using a multimodal nested sampling with a single source model as in Feroz, Marshall & Hobson 2008). From this information, we define an initial sky model where the visibilities of each source s in the field of view are computed according to equation (1) with ellipticity $$\boldsymbol {e}= \mathbf {0}$$, i.e. circular source, and scalelength provided by the linear relation between the log of the median scalelength αmed and flux density S (Rivi et al. 2016a):   \begin{eqnarray} \ln {[\alpha _\mathrm{med}/\mathrm{arcsec}]} = -0.93 +0.33\ln {[S/\mu \mathrm{Jy}]}. \end{eqnarray} (5)The sky visibilities are obtained adding the model visibilities of each source in the beam:   \begin{eqnarray} V_\mathrm{sky}(u,v) = \sum _{s=1}^N V_s(u,v). \end{eqnarray} (6)Starting from this sky model, the source extraction and fitting procedure is performed according to decreasing flux order, i.e. decreasing SNR, as follows: Given the position of the source (l, m), remove the corresponding circular source model visibilities from the sky model and then take the difference between the data and the sky model, so that the visibilities of the current source (with a reduced contamination from the others) are isolated. Apply faceting (Cornwell & Perley 1992): phase shift these visibilities in order to move the phase centre to the location of the source, by multiplying each visibility by the factor exp ( − 2πi(ul + vm)), and average them in a coarse grid (facet). This way, we reduce the field of view to a small patch around the source, with the advantage of reducing the number of visibilities used for the model fitting and therefore accelerating the computation. On the other hand, this procedure limits the maximum wavelength of the Fourier mode that can be measured because of the finite spacing of the facet uv points. Use the source visibilities for measuring the corresponding source ellipticity as in Section 2.1. Use the estimated ellipticity to improve model visibilities of the current source and remove them from the data. Repeat from step (i) until all sources are fitted. Note that in this algorithm the sky model is improved after each source fitting by replacing circular sources with the elliptical source that has been fitted. Moreover, by ordering the source extraction by decreasing flux, the source fitting is performed with a better approximation of the sky model for sources at low SNR. In the case of ‘bad measurements’, the corresponding sources are not removed in the first instance from the data and sky model visibilities, but they are re-fitted at the end of the procedure, when the ellipticities of all the other sources are measured and a better sky model is obtained. These unreliable fits are recognized by a too small standard deviation of the ellipticity likelihood to be realistic. This seems related to errors in the likelihood computation, when the cross-correlation function is not sufficiently smooth to be marginalized over the source position, possibly due to PSF sidelobes or too much noise in the data. Bad measurements are given weight zero in the shear computation. 3 SKA1 SIMULATIONS As in Paper I, we simulate SKA-MID 8-h observations of 1 deg2 at declination δ = −30 deg by using the SKA-MID3 Phase 1 antennas configuration provided in Heystek (2015). This integration time provides a complete circular coverage, i.e. without large gaps (because of the three telescope arms), and allow to reach a sensitivity of 10 μJy at 10σ. It would allow a targeted area of $$800 \deg ^2$$ to be observed with such sensitivity in 10 000 h in a forthcoming SKA1 radio weak lensing survey, sufficient for measuring cluster lensing. We choose the following conservative approximation of the frequency bandwidth: 950–1190 MHz, as proposed in Bonaldi et al. (2016). This seems to be the optimum observation frequency for a weak lensing survey with SKA1-MID, in case only 30 per cent of the full bandwidth of SKA Band 2 is usable (because of RFI problems, other surveys commensality, etc.). Visibilities are sampled every τacc = 60 s and we consider one large channel because smearing effects are not included on shorter time and bandwidth scales. The observations are simulated by using equations (1) and (6) and adding an uncorrelated Gaussian noise whose variance is given by Wrobel & Walker (1999):   \begin{eqnarray} \sigma ^2 = \frac{\mathrm{SEFD}^2}{2\eta ^2\Delta \nu \tau _\mathrm{acc}}, \end{eqnarray} (7)where Δν is the frequency channel bandwidth, SEFD = 400 Jy is the system equivalent flux density for SKA-MID dishes and η = 0.9 is the system efficiency (Braun 2014). For simplicity, we assume that the SKA1-MID core is composed of only SKA dish antennas even if part of it actually contains 64 MeerKAT dishes. SF galaxy populations are generated according to distributions estimated from archival data of VLA radio surveys at 1.4 GHz. As described in Rivi et al. (2016a), we estimated distributions of the flux S and the scalelength α of sources modelled by an exponential profile by the analysis of faint sources (order of tens μJy) catalogue of the SWIRE field survey (Owen & Morrison 2008). This is the radio catalogue with the largest number of SF galaxies, being related to the deepest survey so far in terms of the radio source density (source flux cut ∼ 15 μJy), and it contains source size measurements from the imaged data. The flux distribution is fitted by a power law: p(S) ∝ S−1.34. The scalelength is obtained by the relation with the measured full width at half-maximum:4 FWHM = 2αln (2), and its distribution is fitted dependently on source flux by a lognormal function with mean μ = ln (αmed) and variance σ ∼ 0.3, where αmed is given by equation (5). The variance value is suitably chosen in the middle of a range well-representing scalelength distributions for different flux values. The modulus e of the intrinsic ellipticities is generated according to a function proposed by Miller et al. (2013):   \begin{eqnarray} P(e) = \frac{Ne\left(1-\exp \left[ \frac{e-e_\mathrm{max}}{c}\right]\right)}{(1+e)(e^2+e_0^2)^{1/2}}. \end{eqnarray} (8)The parameter values, c = 0.2298 and e0 = 0.0732 are obtained from fitting the VLA-COSMOS field data. Although this survey is less deep than the SWIRE, but still detecting μJy sources (flux cut ∼75 μJy), we rely on a recent re-analysis of the L-band radio visibility data where the level of systematics in the measurement of the galaxy position angle is significantly reduced (Tunbridge et al. 2016). In fact, the previous analysis (Schinnerer et al. 2007) as well as the one of the VLA-SWIRE were mainly focused on faint source counts. The normalization factor is N = 2.595. We generate galaxy populations with flux densities ranging between 10 and 200 μJy. According to our flux distribution we obtain a source number density of 2.7 gal arcmin−2. To be consistent, we adopt this source number density in our simulation, although more accurate modelling from recent radio continuum surveys suggest that a higher source number density should be detected at such a flux density cut in real observations (Condon et al. 2012; Mancuso et al. 2017). This is the expected source number density for the proposed 2-yr SKA1 radio weak lensing survey covering 5000 $$\deg ^2$$ (Brown et al. 2015). 4 GALAXY SHAPE MEASUREMENT In this section, first we select the facet size to be used in the source extraction by simulating visibilities of individual sources. Then, we simulate populations of galaxies located simultaneously in the field of view in order to show the efficacy of our source extraction algorithm as a function of the source number density. 4.1 Facet size The facet size is affected by the weighting scheme (Briggs, Schwab & Sramek 1999) adopted in the gridding phase. For example, natural weighting optimises the sensitivity for detecting weak sources by emphasizing the data from short baselines. In this case, a relatively large facet size is expected even for covering a single galaxy because a large contribution to the signal is from long wavelength modes which must be adequately sampled by small facet cells. In effect, the source in the image domain turns to be convolved with a large natural-weighted PSF with a broad low-level plateau (see left-hand panel of Fig. 1). Uniform weighting will require instead a much smaller facet size because it emphasizes data from long baselines, where most of the source shape signal is contained. This is reflected by the small uniform-weighted PSF (see right-hand panel of Fig. 1). On the other side, the weighting scheme used in the faceting procedure shouldn't affect the model fitting, provided that the measurement uncertainties are propagated correctly in the likelihood computation (see equation 21 of Paper I for the natural case) and model visibilities are consistent with the observed data. This may not be true for measurement methods in the image domain, as shown in Tunbridge et al. (2016). Figure 1. View largeDownload slide Image of the PSF corresponding to our SKA1 uv coverage Left: natural weighting. The PSF has a broad, low-level plateau because uv points have short spacings close to zero, as they tend to spend more time per unit area near the uv origin. Right: uniform weighting. Figure 1. View largeDownload slide Image of the PSF corresponding to our SKA1 uv coverage Left: natural weighting. The PSF has a broad, low-level plateau because uv points have short spacings close to zero, as they tend to spend more time per unit area near the uv origin. Right: uniform weighting. Since we are interested in the detection of faint sources for radio weak lensing, we adopt a natural weighting scheme. To minimize the number of sources falling in the same facet, we define a facet size dependent on source flux, as it is related to the size of the source. We split the flux total range of the simulated population, i.e. 10–200 μJy, in seven bins as shown in Table 1. Facet uv point coordinates have to be re-computed only once per bin as the model fitting is performed according to source flux order. We chose larger bins at large fluxes because the sizes of such sources increase more gradually with the flux compared to the ones with low flux (see equation 5). To estimate the facet size for each bin, we simulate raw visibilities of a single galaxy in the primary beam in order to avoid any source contamination effects, and vary the noise added to the visibilities (to have SNR ≥ 100) in order to see the effect of the source size only. We measure the galaxy ellipticity after averaging visibilities in the facet. The best-fitting slope5 for the ellipticity measurements of 1000 galaxies is computed for different facet sizes and source flux ranging between the flux bin bounds. We select the facet size when a fixed best-fitting slope threshold of about 0.97 is reached, as listed in Table 1. Table 1. Facet sizes dependent on source flux range and corresponding best-fitting slopes for 1000 sources. S (μJy)  No. cells  Best-fitting slope      e1  e2  200–150  600  0.9774 ± 0.0025  0.9675 ± 0.0025  150–100  550  0.9795 ± 0.0030  0.9664 ± 0.0030  100–80  500  0.9797 ± 0.0032  0.9660 ± 0.0032  80–60  460  0.9774 ± 0.0029  0.9614 ± 0.0029  60–40  420  0.9760 ± 0.0031  0.9631 ± 0.0031  40–20  350  0.9756 ± 0.0028  0.9557 ± 0.0029  20–10  280  0.9765 ± 0.0030  0.9510 ± 0.0030  S (μJy)  No. cells  Best-fitting slope      e1  e2  200–150  600  0.9774 ± 0.0025  0.9675 ± 0.0025  150–100  550  0.9795 ± 0.0030  0.9664 ± 0.0030  100–80  500  0.9797 ± 0.0032  0.9660 ± 0.0032  80–60  460  0.9774 ± 0.0029  0.9614 ± 0.0029  60–40  420  0.9760 ± 0.0031  0.9631 ± 0.0031  40–20  350  0.9756 ± 0.0028  0.9557 ± 0.0029  20–10  280  0.9765 ± 0.0030  0.9510 ± 0.0030  View Large Note that the fitting for the first ellipticity component is better than the second one because of a slightly anisotropy of the PSF as discussed in Section 6. The selected facet sizes are consistent with the relation between the uv grid cell size Δu (in units of wavelengths) and the related field of view (in radians): ψ = 1/Δu (Briggs et al. 1999). We also note that for small sources (low flux), it is actually better to use smaller facets even in the case of a single source in the field of view. This is shown in Fig. 2 where the difference between binned measurements and true values of the first ellipticity component of 1000 galaxies, with realistic flux distribution in the range 10–200 μJy, are plotted both for the case where the facet size is constant and equal to 600 (left-hand panel) and when the facet has a variable size dependent on the source flux (right-hand panel). Similar plots are obtained for the second component. In the latter case we obtained 25 bad measurements (see Section 2.2) and the best-fitting slopes of the two ellipticity components are 0.9552 ± 0.0057 and 0.9426 ± 0.0061, respectively, whereas for the case of 600 × 600 facet the best-fitting slopes for the same galaxy population and noise are 0.9306 ± 0.0054 and 0.9135 ± 0.0056 and the number of bad measurements is three times larger. Figure 2. View largeDownload slide Binned measurements (minus true values) of the ellipticity first component of 1000 individual galaxies randomly located in the field of view and with flux ranging between 10 and 200 μJy. Left: facet size 600 × 600 for all sources. Right: facet size dependent on source flux. Figure 2. View largeDownload slide Binned measurements (minus true values) of the ellipticity first component of 1000 individual galaxies randomly located in the field of view and with flux ranging between 10 and 200 μJy. Left: facet size 600 × 600 for all sources. Right: facet size dependent on source flux. These results are due to the fact that we do not model exactly the primary beam because the model visibilities are directly sampled on the uv facet points. This means that in the image domain the sidelobes of the source model are not suppressed by any apodisation, whereas the gridding of the original uv coverage causes the full image to be apodised by a broad 2D sinc function which has the effect in the data of suppressing background sources that are a long way from the phase centre and the distant sidelobes from the primary source. The grid sampling causes the resulting image domain facet to become a small, but aliased version of the apodised image. The aliased model is an incorrect description of the apodised and aliased data, and the discrepancy will get worse for smaller facets and at large distances from the source. Fig. 2 shows that a suitable facet size dependent on the source flux/size may be a trade-off between these two effects. We could improve the model by applying the same gridding operations as in the data (sampling on the original uv coverage and then averaging in the facet), but this will add a large amount of computational time. Our results show that the adopted model approximation is acceptable, provided that the facet sizes are sufficiently large to not affect the significant sidelobes in the image domain. Otherwise we expect the discrepancy between data and model to become severe and the biases may become less robust and hence less calibratable. 4.2 Dependence on source number density We estimate the efficacy of the source extraction method by measuring the slope of the best-fitting line of 104 galaxy ellipticity measurements as a function of the source number density. For each measurement, we simulate sources located simultaneously in the field of view according to a uniform distribution. Results are plotted in Fig. 3 and show reliable fits, independent of the source density up to 2.8 gal arcmin−2. In this case, the best-fitting slopes of the correlation for each ellipticity component are 0.9365 ± 0.0017 and 0.9262 ± 0.0017 respectively and the number of bad measurements is about 1 per cent. At higher densities galaxy ellipticity measurement starts to deteriorate, as residuals of nearby galaxies affect the model fitting, but may still be good enough for shear measurement because of the improved statistics (as shown in Section 5). Shape measurements of galaxies may be improved by a joint fitting within facets by applying the Hamiltonian Monte Carlo technique (Neal 2011). RadioLensfit results used as starting points should reduce the burn-in phase and accelerate convergence. Since the number of sources in the facet will be relatively small this approach becomes more feasible and preliminary results with this method show a better accuracy in the galaxy ellipticity fitting, although requiring a large computational time (Rivi et al., in preparation). Figure 3. View largeDownload slide Best-fitting slope for both galaxy ellipticity components as a function of the source number density. Figure 3. View largeDownload slide Best-fitting slope for both galaxy ellipticity components as a function of the source number density. Using a single channel, the serial version of RadioLensfit running on an Intel Xeon E5-2650 takes on average about 10 s gal−1 computing time for the model fitting. As discussed in Paper I, the shared memory parallelization with OpenMP allows us to exploit all the computational resources when the amount of memory for the source model fitting requires the usage of the full CPU. Its implementation has been optimized by distributing to each thread the likelihood computation and marginalization over the position parameters for different scalelength values of the model. It doesn't scale linearly with the number of threads because the likelihood mariginalization over the scalelength parameter is not parallel and there is an overhead for the creation and destruction of OpenMP threads at each iteration of the likelihood maximization and sampling. This version running on all the eight cores of the CPU takes on average about 2.4 s gal−1. 5 SHEAR Following Paper I, for shear measurement we simulate galaxy populations as described in Section 3 in a field of view of 1 $$\deg ^2$$. We generate populations free of shape noise (Massey et al. 2007; Nakajima & Bernstein 2007): for each ellipticity modulus, 10 equally spaced galaxy orientations are generated so that the corresponding ellipticity values are distributed uniformly on a circle, and galaxies whose ellipticity is on the same ring are given the same size and flux. We generate sources randomly located according to a uniform distribution. All measurements are performed with facet size dependent on the source flux as defined in Table 1. We measure the reduced shear $$\boldsymbol {g}=g_1 + \mathrm{i}g_2$$ as the weighted average of the galaxies’ ellipticities using weights that approximate the inverse-variance of each ellipticity measurement. Error bars are given by the standard deviation of the shear values estimated from 1000 bootstrap resamples. Fig. 4 shows shear measurements as a function of the source number density up to 8 gal armin−2, when no input shear is applied, i.e. g1 = g2 = 0. At high densities, the larger number of sources compensates the less accurate galaxy shape fitting (Fig. 3), still producing shear values consistent with the results obtained at the SKA1 source density corresponding to a population of about 104 galaxies. Figure 4. View largeDownload slide Shear components estimated from a galaxy population in 1 $$\deg ^2$$ as a function of the source number density for input g = 0. Figure 4. View largeDownload slide Shear components estimated from a galaxy population in 1 $$\deg ^2$$ as a function of the source number density for input g = 0. For this population, we measure the shear also for input reduced shear values with amplitude g = 0.04 and eight different orientations. The input shear $$\boldsymbol {g}$$ action on the intrinsic galaxy ellipticity $$\boldsymbol {e}^s$$ is simulated following Seitz & Schneider (1997):   \begin{eqnarray} \boldsymbol {e} = \frac{\boldsymbol {e}^s + \boldsymbol {g}}{1+\boldsymbol {g}^*\boldsymbol {e}^s}, \end{eqnarray} (9)where $$\boldsymbol {g}^*$$ is the shear complex conjugate. We compare results with the optimal case where each galaxy is at the phase centre and the only one contained in the field of view, considering the same galaxy population. Results are plotted in Fig. 5, both for SNR ≥ 10 and SNR ≥ 25, where measurements from individual source visibilities are green crosses and the ones from the same population but with all sources simultaneously in the primary beam are black crosses. Clearly error bars (cross arms) are larger at SNR ≥ 10 as the galaxy population is dominated by lower flux sources. Figure 5. View largeDownload slide Comparison of shear measurements: input values are blue points, measured values from single sources at the phase centre are green crosses, measured values from sources simultaneously in the f.o.v are black crosses. Figure 5. View largeDownload slide Comparison of shear measurements: input values are blue points, measured values from single sources at the phase centre are green crosses, measured values from sources simultaneously in the f.o.v are black crosses. The measured shear bias, defined as   \begin{eqnarray} g_i^m - g_i = m_i g_i + c_i, \qquad i=1,2, \end{eqnarray} (10)is shown in Table 2. At SNR ≥ 10, the multiplicative biases mi for the two shear components are, respectively, 1.6 and 1.4 times the ideal case of a single source in the field of view, while additive bias ci are almost consistent. Selecting galaxies with SNR ≥ 25 the population reduces to 5810 sources (i.e. 1.6 gal arcmin−2). This should affect the bias uncertainty only, as the bias on the shear measurement should not depend on the number of sources. At this regime, the m values reduce by a factor two instead of three as in the single source case. This is due to the source signal contamination by residuals of nearby galaxies, which is a new contribution to the shear bias that seems to have an effect on the multiplicative terms only. It may be mitigated by refining ellipticity measurements by joint fitting within larger facets, as explained in Section 4.2. Note that ‘neighbour bias’ affects optical lensing measurements also (see Miller et al. 2013; Jarvis et al. 2016; Mandelbaum et al. 2017; Zuntz et al. 2017; Samuroff et al. 2018). Table 2. Shear bias components estimated from a realistic population of ∼104 galaxies randomly distributed in 1 $$\deg ^2$$, corresponding to a source number density of 2.8 gal arcmin−2.     m1  c1  m2  c2  SNR ≥ 10  Single source  − 0.0904 ± 0.0186  0.00655 ± 0.00050  − 0.1297 ± 0.0171  0.00632 ± 0.00046    Multiple sources  − 0.1428 ± 0.0274  0.00872 ± 0.00073  − 0.1864 ± 0.0256  0.00578 ± 0.00067  SNR ≥ 25  Single source  − 0.0352 ± 0.0127  − 0.00006 ± 0.00034  − 0.0506 ± 0.0132  − 0.00070 ± 0.00035    Multiple sources  − 0.0677 ± 0.0242  0.00089 ± 0.00064  − 0.0992 ± 0.0230  − 0.00112 ± 0.00060      m1  c1  m2  c2  SNR ≥ 10  Single source  − 0.0904 ± 0.0186  0.00655 ± 0.00050  − 0.1297 ± 0.0171  0.00632 ± 0.00046    Multiple sources  − 0.1428 ± 0.0274  0.00872 ± 0.00073  − 0.1864 ± 0.0256  0.00578 ± 0.00067  SNR ≥ 25  Single source  − 0.0352 ± 0.0127  − 0.00006 ± 0.00034  − 0.0506 ± 0.0132  − 0.00070 ± 0.00035    Multiple sources  − 0.0677 ± 0.0242  0.00089 ± 0.00064  − 0.0992 ± 0.0230  − 0.00112 ± 0.00060  View Large As discussed in Paper I, the noise bias values exceed SKA1 survey requirements6 except for the additive component at SNR ≥ 25. However they are comparable to the ones obtained from optical surveys using lensfit (Fenech Conti et al. 2017) and im3shape (Zuntz et al. 2017), where a shear calibration correction reduced the multiplicative bias to well below the percent level. Standard approaches in the optical domain derive such calibration by inferring the bias from simulated data matching the observations (Bruderer et al. 2016) or parametrizing the bias as a function of the observed galaxy properties (Kuijken et al. 2015; Jarvis et al. 2016). Recently, a self-calibration approach, implemented in the Metacalibration method (Huff & Mandelbaum 2017; Sheldon & Huff 2017) and used in the analysis of the Dark Energy Survey7 (DES) (Zuntz et al. 2017), proved to be the most efficient, being able to recover the input shear in realistic simulations to better than a part in a thousand. It also isotropises the PSF to remove any additive bias. The key idea of the method is to compute the shear estimator response for a shape measurement directly from observed data perturbed with a small known shear. This way all the features present in real data are already included, which are instead extremely difficult to model accurately in external simulations, and it can be applied to any shear measurement method based on averages of galaxy shapes. A similar approach may then be implemented quite easily in the interferometer data analysis, with the advantage that for the additive bias at radio wavelengths we know the PSF much better than at optical wavelengths and we can make the PSF isotropic directly by weighting the uv-plane (as discussed in Section 6). 6 ADDITIVE BIAS DEPENDENCE ON THE IMAGE-PLANE PSF SHAPE It is well known from weak lensing optical surveys that shear additive bias is dependent on the PSF shape (Miller et al. 2013). For radio interferometers, the PSF is deterministically defined by the uv coverage of the telescope (i.e. antennas locations and pointing direction). For example, if we increase the antenna pointing declination to a larger zenith distance the PSF shape becomes compressed along the y-axis. We measure the additive noise bias for different pointing declinations at various SNR values. Starting from our uv coverage (corresponding to the same declination as the observatory latitude), we simulate the effect on the uv points when the phase centre declination is increased by an angle ϕ. A plot of the image-plane PSF ellipticity components as functions of the zenith distance is given in Fig. 6. The R-squared size of the PSF slightly increases from 14.83 arcsec2 to 15.54 arcsec2. These values are computed from the quadrupole moments of the image domain as follows (Schneider 2006):   \begin{eqnarray} \boldsymbol {e} = \frac{Q_{xx}-Q_{yy}+2\mathrm{i}Q_{xy}}{Q_{xx}+Q_{yy}+2(Q_{xx}Q_{yy}-Q^2_{xy})^{\frac{1}{2}}}, \end{eqnarray} (11)  \begin{eqnarray} R = Q_{xx}+Q_{yy}. \end{eqnarray} (12) Figure 6. View largeDownload slide PSF ellipticity components versus the zenith distance at the starting of the observation. Figure 6. View largeDownload slide PSF ellipticity components versus the zenith distance at the starting of the observation. We simulate individual source visibilities, to avoid nearby source contamination effects, and assume a constant maximum facet size 1000 × 1000 to ensure that the galaxy convolved with the PSF is contained in the facet even when the PSF becomes highly anisotropic. We observe that the additive bias is dependent on source size. In fact measurements at the same SNR obtained by lowering the noise instead of increasing the source flux cut, produce larger bias values, meaning that the additive bias worsens when source sizes decreases. This is consistent with the analysis presented in Massey et al. (2013). Noise bias causes a correlation between measured shear and PSF ellipticity even when we correct for the PSF in the model fitting. This becomes noticeable at low SNR, where the first ellipticity component increases significantly towards larger negative values, as the PSF becomes more compressed along the y direction (see left-hand panel of Fig. 7). At large SNR the additive bias almost disappears independently of the PSF shape (see right-hand panel of Fig. 7). This is in good agreement with what we should expect. Because of our long integration time the PSF is not isotropic even when the declination equals the observatory latitude, as in the baseline simulation, besides the fact that the distribution of the SKA-MID baselines on the ground are not circularly symmetric. Therefore, a large additive bias is still measured at small zenith distances. The PSF anisotropy may be reduced by combining snapshots obtained over a range of hour angles, however this may not be sufficient to reach an additive bias acceptable for weak lensing surveys. To further reduce the noise bias at low SNR, we also need to weight the uv plane to ensure that the PSF is more isotropic. A standard technique in radio imaging to improve PSF shape is to use tapering functions (Briggs et al. 1999) to define uv points weights, although a more specific weighting scheme may be required. Moreover, as for the multiplicative bias, we can calibrate the additive shear bias with simulations. This is more feasible than in optical surveys because the PSF is deterministic at radio wavelengths. However, any such calibration would be strongly dependent on distributions of source properties, so isotropising the PSF is a much better option. Figure 7. View largeDownload slide Shear additive bias as a function of the zenith distance at different lower signal to noise. Facet size fixed at 1000 × 1000. Figure 7. View largeDownload slide Shear additive bias as a function of the zenith distance at different lower signal to noise. Facet size fixed at 1000 × 1000. Note that when using variable size faceting, the facet size must be dependent not only on the source size but also on the PSF shape. In fact, as the PSF becomes anisotropic the facet size may become too small relative to the size of the source convolved with the PSF, modifying the effective shape of the source. For example, the left-hand panel of Fig. 8 shows what happens at SNR ≥ 10 when we maintain the same flux dependent facet size (Table 1) for all pointing declinations: at large zenith distances the shortest baselines occupy smaller uv frequencies and therefore can measure wavelengths longer than the limit imposed by the small facet size used to extract the majority of the galaxy population. If we increase the facet size according to the source declination, e.g. 50 cells per side every 10 deg declination increment, we obtain consistent results with the case of one single large facet (see right-hand panel of Fig. 8). Figure 8. View largeDownload slide Shear additive bias at SNR ≥ 10 as a function of the zenith distance. Left: facet size dependent on source flux according to Table 1 for all declinations. Right: facet size increased with source declination produce similar results as for a constant large facet size. Figure 8. View largeDownload slide Shear additive bias at SNR ≥ 10 as a function of the zenith distance. Left: facet size dependent on source flux according to Table 1 for all declinations. Right: facet size increased with source declination produce similar results as for a constant large facet size. 7 CONCLUSIONS We have extended the presentation of the RadioLensfit method, introduced in Paper I for the simple case of individual galaxies located at the phase centre, to the real case where many galaxies are randomly located in the field of view. This has been done by isolating the visibilities of each source and shifting the phase centre to the source position so that a coarse grid can still be used to reduce nearby galaxy residuals contamination and accelerate ellipticity measurement computation. Source extraction has been performed by removing apart from the data the simulated visibilities of the sky model, but the source of interest, given the positions and fluxes of all sources in the field of view from the image, and down-weighting what remains of nearby source-contamination by averaging visibilities in a coarse grid (facet). For gridding, we adopted a natural weighting to maximize sensitivity and estimated the smallest facet size dependent on source flux thresholds in order to minimize the number of nearby galaxies included in the facet. We tested the source extraction procedure, simulating visibilities of SF galaxy populations observed by SKA1-MID in the first 30 per cent of frequency Band 2. We adopted flux and scalelength parameters distributions estimated from the VLA SWIRE catalogue and used the lensfit ellipticity prior with coefficients fitted from a new version of the VLA COSMOS catalogue optimized on shape measurements. We showed the efficacy of our source extraction algorithm as a function of the source number density, obtaining a reliable galaxy ellipticity fitting for the density expected from the current proposal of the SKA1 radio weak lensing survey. Shear measurements from 8-h observation of 1 deg2 show that the bias due to the extraction procedure mainly affects the multiplicative bias as no significant change has been observed for the additive bias when comparing with the bias obtained for the ideal case of a single source at the phase centre at a time. This bias may be mitigated by a second step in the galaxy ellipticity measurement, where a joint fitting within the facets is performed with HMC, starting from the values obtained with RadioLensfit. However, multiplicative noise bias calibration is also required as for the optical domain. We finally observed that because of our uv coverage the PSF is slightly anisotropic even if pointing close to the zenith, therefore we obtain a large additive bias (on average about 0.0068 ± 0.0006 at SNR ≥ 10). Although a suitable choice of the integration time (split and distributed along a longer period of time) may reduce the PSF anisotropy, a uv weighting scheme may still be required to satisfy weak lensing requirements. It should be optimized to avoid any significant reduction of the signal to noise. ACKNOWLEDGEMENTS We thank Tom Kitching for useful discussions and Ben Tunbridge for providing the distribution parameters of our ellipticity prior, obtained by fitting VLA-COSMOS data as for the prior function presented in his paper. We also thank Sphesihle Makhathini for the support in the simulation of the SKA1-MID uv coverages. MR acknowledges the support of the Science and Technology Facilities Council (STFC) via an SKA grant. LM acknowledges STFC grant ST/N000919/1. Footnotes 1 https://www.skatelescope.org 2 http://www.e-merlin.ac.uk/legacy/projects/superclass.html 3 SKA-MID latitude is −30°49΄48.00΄S. 4 The FWHM is derived from the Gaussian model fit of the source after PSF deconvolution. 5 Consistently with Paper I, we refer to the ellipticity best-fitting slope instead of the multiplicative bias when measuring galaxy shapes. This terminology is used to clearly distinguish it from the shear multiplicative bias, which is obtained from the best fit of shear measurements (each being the weighted average of galaxy ellipticities). 6 For a 2-yr SKA1-MID weak lensing survey over 5000 $$\deg ^2$$ and zmed = 1.0, the requirements for cosmological parameters measurements to be dominated by statistical rather than systematic errors are: multiplicative bias m < 0.0067, additive bias c < 0.00082 (Brown et al. 2015). They are derived using the rules provided in Amara & Refregier (2008). 7 https://www.darkenergysurvey.org REFERENCES Amara A., Refregier A., 2008, MNRAS , 391, 228 https://doi.org/10.1111/j.1365-2966.2008.13880.x CrossRef Search ADS   Becker R., White R., Helfand D. J., 1995, ApJ , 450, 559 https://doi.org/10.1086/176166 CrossRef Search ADS   Blain A., 2002, ApjL , 570, L51 https://doi.org/10.1086/341103 CrossRef Search ADS   Bonaldi A., Harrison I., Camera S., Brown M., 2016, MNRAS , 463, 3686 https://doi.org/10.1093/mnras/stw2104 CrossRef Search ADS   Braun R., 2014, SKA Memo SKA-TEL.SCI-SKO-SRQ-001  Bridle A., Schwab F., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser. Vol. 180, Synthesis Imaging in Radio Astronomy II . Astron. Soc. Pac., San Francisco, p. 371 Briggs D., Schwab F., Sramek R., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser. Vol. 180, Synthesis Imaging in Radio Astronomy II . Astron. Soc. Pac., San Francisco, p. 129 Brown M., Battye R., 2011, MNRAS , 410, 2057 Brown M. et al.  , 2015, Advancing Astrophysics with the Square Kilometre Array, Vol. 2, Proc. Sci. , SISSA, Trieste, p. 1365 Bruderer C., Chang C., Refregier A., Amara A., Bergé J., Gamper L., 2016, ApJ , 817, 25 https://doi.org/10.3847/0004-637X/817/1/25 CrossRef Search ADS   Camera S., Harrison I., Bonaldi A., Brown M., 2017, MNRAS , 464, 4747 https://doi.org/10.1093/mnras/stw2688 CrossRef Search ADS   Chang T., Refregier A., 2002, ApJ , 570, 447 https://doi.org/10.1086/339496 CrossRef Search ADS   Chang T., Refregier A., Helfand D., 2004, ApJ , 617, 794 https://doi.org/10.1086/425491 CrossRef Search ADS   Condon J. et al.  , 2012, AJ , 758, 23 https://doi.org/10.1088/0004-637X/758/1/23 CrossRef Search ADS   Cornwell T. J., Perley R. A., 1992, A&A , 261 Demetroullas C., Brown M., 2016, MNRAS , 456, 3100 https://doi.org/10.1093/mnras/stv2876 CrossRef Search ADS   Fenech Conti I., Herbonnet R., Hoekstra H., Merten J., Miller L., Viola M., 2017, MNRAS , 467, 1627 Feroz F., Marshall P., Hobson M., 2008, preprint (arXiv:0810.0781) Harrison I., Brown M., 2015, SKA ECP150007, preprint (arXiv:1507.06639) Harrison I., Camera S., Zuntz J., Brown M., 2016, MNRAS , 463, 3674 https://doi.org/10.1093/mnras/stw2082 CrossRef Search ADS   Heystek L., 2015, SKA Memo SKA-TEL-INSA-0000537  Huff E., Mandelbaum R., 2017, ApJ , preprint (arXiv:1702.02600) Jarvis M. et al.  , 2016, MNRAS , 460, 2245 https://doi.org/10.1093/mnras/stw990 CrossRef Search ADS   Kuijken K. et al.  , 2015, MNRAS , 454, 3500 https://doi.org/10.1093/mnras/stv2140 CrossRef Search ADS   Mancuso C. et al.  , 2017, AJ , 842, 95 https://doi.org/10.3847/1538-4357/aa745d CrossRef Search ADS   Mandelbaum R. et al.  , 2017, MNRAS , preprint (arXiv:1710.00885) Massey R. et al.  , 2007, MNRAS , 376, 13 https://doi.org/10.1111/j.1365-2966.2006.11315.x CrossRef Search ADS   Massey R. et al.  , 2013, MNRAS , 429, 661 https://doi.org/10.1093/mnras/sts371 CrossRef Search ADS   Melchior P., Viola M., 2012, MNRAS , 424, 2757 https://doi.org/10.1111/j.1365-2966.2012.21381.x CrossRef Search ADS   Melchior P., Böhnert A., Lombardi M., Bartelmann M., 2010, A&A , 510, A75 https://doi.org/10.1051/0004-6361/200912785 CrossRef Search ADS   Miller L. et al.  , 2013, MNRAS , 429, 2858 https://doi.org/10.1093/mnras/sts454 CrossRef Search ADS   Morales M., 2006, ApJ , 650, L21 https://doi.org/10.1086/508614 CrossRef Search ADS   Nakajima R., Bernstein G., 2007, ApJ , 133, 1763 https://doi.org/10.1086/511957 CrossRef Search ADS   Neal R. M., 2011, in Brooks S., Gelman A., Jones G., Meng X.-L., eds, MCMC using Hamiltonian dynamics. Astronomical Journal of the American Astronomical Society , Chapman & Hall/CRC Press Owen F., Morrison G., 2008, AJ , 136, 1889 https://doi.org/10.1088/0004-6256/136/5/1889 CrossRef Search ADS   Patel P., Bacon D. J., Beswick R. J., Muxlow T. W. B., Hoyle B., 2010, MNRAS , 401, 2572 https://doi.org/10.1111/j.1365-2966.2009.15836.x CrossRef Search ADS   Patel P. et al.  , 2015, Advancing Astrophysics with the Square Kilometre Array, Vol. 2, Proc. Sci. , SISSA, Trieste, p. 1427 Refregier A., Kacprzak T., Amara A., Bridle S., Rowe B., 2012, MNRAS , 425, 1951 https://doi.org/10.1111/j.1365-2966.2012.21483.x CrossRef Search ADS   Rivi M., Miller L., Makhathini S., Abdalla F., 2016a, EXTRA-RADSUR2015, Proc. Sci. , SISSA, Trieste Rivi M., Miller L., Makhathini S., Abdalla F. B., 2016b, MNRAS , 463, 1881 (Paper I) https://doi.org/10.1093/mnras/stw2041 CrossRef Search ADS   Samuroff S. et al.  , 2018, MNRAS , 475, 4524 https://doi.org/10.1086/516587 preprint (arXiv:1708.01534) Schinnerer E. et al.  , 2007, ApJS , 172, 46 https://doi.org/10.1086/516587 CrossRef Search ADS   Schneider P., 2006, in Meylan G., Jetzer P., North P., Schneider P., Kochanek C. S., Wambsganss J., eds, Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro. Springer-Verlag, Berlin, p. 269 Seitz C., Schneider P., 1997, A&A , 318, 687 Sheldon E. S., Huff E. M., 2017, ApJ , 841, 24 https://doi.org/10.3847/1538-4357/aa704b CrossRef Search ADS   Smirnov O. M., 2011, A&A , 527, A106 https://doi.org/10.1051/0004-6361/201016082 CrossRef Search ADS   Tunbridge B., Harrison I., Brown M., 2016, MNRAS , 463, 3339 https://doi.org/10.1093/mnras/stw2224 CrossRef Search ADS   Wrobel J., Walker R., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser., Vol. 180, Synthesis Imaging in Radio Astronomy II , Astron. Soc. Pac., San Francisco, p. 171 Zuntz J. et al.  , 2017, MNRAS , preprint (arXiv:1708.01533) © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Monthly Notices of the Royal Astronomical Society Oxford University Press

# Radio weak lensing shear measurement in the visibility domain – II. Source extraction

, Volume 476 (2) – May 1, 2018
10 pages

Publisher
The Royal Astronomical Society
ISSN
0035-8711
eISSN
1365-2966
D.O.I.
10.1093/mnras/sty371
Publisher site
See Article on Publisher Site

### Abstract

Abstract This paper extends the method introduced in Rivi et al. (2016b) to measure galaxy ellipticities in the visibility domain for radio weak lensing surveys. In that paper, we focused on the development and testing of the method for the simple case of individual galaxies located at the phase centre, and proposed to extend it to the realistic case of many sources in the field of view by isolating visibilities of each source with a faceting technique. In this second paper, we present a detailed algorithm for source extraction in the visibility domain and show its effectiveness as a function of the source number density by running simulations of SKA1-MID observations in the band 950–1150 MHz and comparing original and measured values of galaxies’ ellipticities. Shear measurements from a realistic population of 104 galaxies randomly located in a field of view of 1 $$\deg ^2$$ (i.e. the source density expected for the current radio weak lensing survey proposal with SKA1) are also performed. At SNR ≥ 10, the multiplicative bias is only a factor 1.5 worse than what found when analysing individual sources, and is still comparable to the bias values reported for similar measurement methods at optical wavelengths. The additive bias is unchanged from the case of individual sources, but it is significantly larger than typically found in optical surveys. This bias depends on the shape of the uv coverage and we suggest that a uv-plane weighting scheme to produce a more isotropic shape could reduce and control additive bias. gravitational lensing: weak, methods: statistical, techniques: interferometric, cosmology: observations 1 INTRODUCTION Cosmological or targeted surveys of weak gravitational lensing at radio wavelengths may have a relevant role in the next years, when the Square Kilometre Array (SKA)1 radio telescope will start to operate, providing a density of detected galaxies sufficient for shear measurement and a resolution to reliably measure their shapes. They will also be able to probe to higher redshifts given the different galaxy redshift distributions compared to the optical band (Brown et al. 2015). Although the galaxy number density in the radio band will be lower than in the optical, the possibility to observe deeper can make radio weak lensing surveys for cosmology measurements competitive with the corresponding optical surveys, as shown in recent forecasts from SKA simulations (Harrison et al. 2016). Moreover, radio observations have the advantage of a deterministic knowledge of the image-plane point spread function (PSF), being the Fourier Transform of the uv coverage, and will provide unique approaches for mitigating intrinsic alignments, such as concurrent measurements of polarization (Brown & Battye 2011) and galaxy rotation velocities (Blain 2002; Morales 2006). Being subject to different observational systematics, cross-correlation with optical observations of the same field will allow suppression of systematic errors on shear measurement from future large surveys (Patel et al. 2010; Demetroullas & Brown 2016; Camera et al. 2017). This is quite relevant for precision cosmology as these errors may become comparable to, and larger than, the statistical noise. The precursor radio weak lensing survey SuperCLASS2 is already underway and will soon provide data that may be used to test new methods required for accurate galaxy shape measurement in the radio band. A natural approach for such methods is working in the visibility domain where the data originates and the noise is Gaussian, avoiding non-linear data manipulation of the imaging process. SKA simulations have already shown that current imaging methods produce images with structures in the residuals which may dominate the cosmological signal (Patel et al. 2015). Also cross-correlation analysis using real data images shows that no evidence of correlation is found between the optical and radio intrinsic shape of the matched objects (Patel et al. 2010; Tunbridge, Harrison & Brown 2016). This result suggests the presence of systematics in the procedure adopted for turning the visibility data into images, although a significant percentage of AGN sources in the observed population may be another possible explanation, as well as the astrophysical scatter between optical and radio position angle due to the different emission mechanisms in the two bands. Currently, cosmic shear in the radio band has been successfully detected only working in the visibility domain but obtaining sources position from the image (Chang, Refregier & Helfand 2004). Galaxies’ ellipticities from the VLA FIRST survey (Becker, White & Helfand 1995) were measured by decomposing them into shapelets, an orthonormal basis of functions corresponding to perturbations around a circular Gaussian invariant under Fourier transform (Chang & Refregier 2002). Since the data size and the number of resolved sources ($${\sim } 20\hbox{--}30 \deg ^{-2}$$) of each pointing is quite small, a joint fitting of the shapelet coefficients was possible by solving normal equations. Such an approach, computationally convenient, becomes very challenging when dealing with the order of 104 sources per square degree and a very large data set per pointing (order of PetaBytes), as expected from SKA Phase 1 continuum surveys (Brown et al. 2015). Moreover, shapelets introduce a shear bias as they cannot accurately model steep brightness profiles and highly elliptical galaxy shapes (Melchior et al. 2010). In the companion paper (Rivi et al. 2016b, hereafter Paper I), we presented RadioLensfit, an alternative method working in the visibility domain where model fitting is performed on a single source at a time using an exponential profile as model for the galaxy. It is an adaptation of the optical Bayesian lensfit method (Miller et al. 2013) to radio data, where model visibilities are defined analytically and the likelihood is marginalized over uninteresting parameters. The method was tested in the simple case of individual galaxy visibilities simulated adopting the SKA1 uv coverage described in Section 3, and the shear noise bias (Refregier et al. 2012; Melchior & Viola 2012) estimated as a function of the signal-to-noise ratio (SNR). Results compared with requirements (Brown et al. 2015) for the proposed SKA1 radio weak lensing survey (Braun 2014) showed that the multiplicative shear bias may need calibration corrections similar to those for optical surveys, while the additive bias have to be controlled by an isotropic sampling of the visibility plane. In this paper, we extend this work implementing the method for isolating source visibilities from realistic data, i.e. when many galaxies are in the field of view. We estimate its effectiveness in terms of ellipticity fitting and shear measurement by running SKA1-MID simulations as we did in the previous paper. We finally investigate the effect of the shape of the uv coverage on the additive shear bias. This paper is organized as follows. In Section 2, we summarize RadioLensfit and present the extraction algorithm. In Section 3, details of SKA1 simulations are provided, while in Sections 4 and 5 results for galaxy ellipticity and shear measurements are presented respectively. Finally, we discuss the shear additive noise bias in Section 6. 2 OVERVIEW OF RadioLensfit RadioLensfit is a method for measuring radio galaxy ellipticities in the visibility domain. The idea is to adapt the approach used in the optical case to radio data, i.e. extracting from visibilities and model fitting a single source at a time. Source extraction is difficult in the Fourier domain because signals from all sources in the primary beam are mixed altogether in the visibilities and sources are no longer localized. For this reason, a joint analysis with the image domain may be needed: it allows us to identify sources and measure their position and flux with sufficient accuracy. With such information, we can also compute a model of the observed sky and use it to approximate the signal from the other galaxies that must be removed when extracting each source. The extraction is completed using a faceting technique that phase shifts the phase centre to the source position and further reduces its signal contamination by averaging visibilities in a coarse grid. Finally, the model fitting can be performed as for the simple case of a single galaxy in the primary beam located at the phase centre as summarized in Section 2.1. This way we can largely reduce the computational time when a huge number of sources are in the field of view (as for SKA) instead of trying a challenging joint fitting of all sources. A detailed algorithm for the fitting of many sources in the primary beam is presented in Section 2.2. 2.1 Galaxy ellipticity fitting In Paper I, we introduced this method as an adaptation of lensfit (Miller et al. 2013) by performing the chi-square fitting of single source visibilities. They are evaluated at the uv points, that are defined by the baselines formed between two antennas projected on the plane orthogonal to the pointing direction. Model visibilities of a star-forming (SF) galaxy are computed analytically as the Fourier transform of the exponential brightness profile (Sérsic index n = 1):   \begin{eqnarray} V(u,v) = \Big ( \frac{\lambda _\mathrm{ref}}{\lambda }\Big )^\beta \frac{S_{\lambda _\mathrm{ref}}{\boldsymbol {e}}^{2\pi \mathrm{i} (u l + v m)}}{(1+4\pi ^2 \alpha ^2 |\boldsymbol{A}^{-T}\boldsymbol {k}|^2)^{3/2}}, \end{eqnarray} (1)where $$\boldsymbol {k}=(u,v)^T$$ is measured in wavenumber units, β = −0.7 is the assumed spectral index for the synchrotron radiation emitted by the galaxy disc, (l, m) and α are the source position and scalelength, respectively, $$S_{\lambda _\mathrm{ref}}$$ is the source flux at reference wavelength λref. The ellipticity parameters (e1, e2) are contained in the matrix $$\boldsymbol{A}$$ that linearly transforms the circular exponential profile to an ellipse:   \begin{eqnarray} \boldsymbol{A} = \left(\begin{array}{cc}1-e_{1} & -e_{2} \\ -e_{2} & 1+e_{1} \end{array} \right). \end{eqnarray} (2)We assume the following ellipticity definition:   \begin{eqnarray} {\boldsymbol{\it e}} = e_1 +\mathrm{i}e_2 = \frac{a-b}{a+b}\mathrm{e}^{2\mathrm{i}\theta }, \end{eqnarray} (3)where a and b are the galaxy major and minor axes, respectively, and θ is the galaxy orientation. The likelihood is marginalized over non-interesting parameters such as flux, scalelength and position, adopting uniform priors for the flux and position, and a lognormal prior dependent on the flux for the scalelength (see Section 3). This way we obtain a likelihood function of only the ellipticity parameters. The galaxy ellipticity measurement is given by the likelihood mean point and 1D standard deviation (defined as the square root of the covariance matrix determinant), obtained after sampling the likelihood with an adaptive grid covering a neighbourhood around the maximum point. In real observations, the finite channel bandwidth and sampling time introduce smearing effects that are dependent on the source position. These effects may be approximated analytically and included in the visibilities model (Bridle & Schwab 1999; Smirnov 2011). For example, for frequency smearing, assuming a square bandpass filter in the expression of the smeared visibility presented in Bridle & Schwab (1999), we obtain:   \begin{eqnarray} \tilde{V}(u,v)= V(u,v) \mathrm{sinc}[\pi (ul+vm)\Delta \nu /(\nu _0 c)], \end{eqnarray} (4)where uv coordinates are taken at the mid-channel frequency ν0, Δν is the channel bandwidth and sinc(x) = sin (x)/x. Another option is to make the observation with very tiny frequency channels and sampling time intervals. Harrison & Brown (2015) proposed to use ∼30 kHz channel bandwidth and 0.5 s sampling time to make smearing tolerable, but meaning a huge number of uv points. In this case, raw visibilities may be averaged into a single uv grid without jeopardizing ellipticity measurements. In fact, observations from the same pair of antennas at different frequencies (resp. times) correspond to visibilities evaluated at different uv points along a radial (resp. tangential) direction, therefore these visibilities can be treated as the ones evaluated at uv points related to different baselines. Depending on the grid size, data volume and then computational time may be considerably reduced. 2.2 Source extraction We assume flux and source positions are provided. For example, they may be estimated from a cleaned image of the same data that are analysed, or applying MC methods to the visibilities (e.g. using a multimodal nested sampling with a single source model as in Feroz, Marshall & Hobson 2008). From this information, we define an initial sky model where the visibilities of each source s in the field of view are computed according to equation (1) with ellipticity $$\boldsymbol {e}= \mathbf {0}$$, i.e. circular source, and scalelength provided by the linear relation between the log of the median scalelength αmed and flux density S (Rivi et al. 2016a):   \begin{eqnarray} \ln {[\alpha _\mathrm{med}/\mathrm{arcsec}]} = -0.93 +0.33\ln {[S/\mu \mathrm{Jy}]}. \end{eqnarray} (5)The sky visibilities are obtained adding the model visibilities of each source in the beam:   \begin{eqnarray} V_\mathrm{sky}(u,v) = \sum _{s=1}^N V_s(u,v). \end{eqnarray} (6)Starting from this sky model, the source extraction and fitting procedure is performed according to decreasing flux order, i.e. decreasing SNR, as follows: Given the position of the source (l, m), remove the corresponding circular source model visibilities from the sky model and then take the difference between the data and the sky model, so that the visibilities of the current source (with a reduced contamination from the others) are isolated. Apply faceting (Cornwell & Perley 1992): phase shift these visibilities in order to move the phase centre to the location of the source, by multiplying each visibility by the factor exp ( − 2πi(ul + vm)), and average them in a coarse grid (facet). This way, we reduce the field of view to a small patch around the source, with the advantage of reducing the number of visibilities used for the model fitting and therefore accelerating the computation. On the other hand, this procedure limits the maximum wavelength of the Fourier mode that can be measured because of the finite spacing of the facet uv points. Use the source visibilities for measuring the corresponding source ellipticity as in Section 2.1. Use the estimated ellipticity to improve model visibilities of the current source and remove them from the data. Repeat from step (i) until all sources are fitted. Note that in this algorithm the sky model is improved after each source fitting by replacing circular sources with the elliptical source that has been fitted. Moreover, by ordering the source extraction by decreasing flux, the source fitting is performed with a better approximation of the sky model for sources at low SNR. In the case of ‘bad measurements’, the corresponding sources are not removed in the first instance from the data and sky model visibilities, but they are re-fitted at the end of the procedure, when the ellipticities of all the other sources are measured and a better sky model is obtained. These unreliable fits are recognized by a too small standard deviation of the ellipticity likelihood to be realistic. This seems related to errors in the likelihood computation, when the cross-correlation function is not sufficiently smooth to be marginalized over the source position, possibly due to PSF sidelobes or too much noise in the data. Bad measurements are given weight zero in the shear computation. 3 SKA1 SIMULATIONS As in Paper I, we simulate SKA-MID 8-h observations of 1 deg2 at declination δ = −30 deg by using the SKA-MID3 Phase 1 antennas configuration provided in Heystek (2015). This integration time provides a complete circular coverage, i.e. without large gaps (because of the three telescope arms), and allow to reach a sensitivity of 10 μJy at 10σ. It would allow a targeted area of $$800 \deg ^2$$ to be observed with such sensitivity in 10 000 h in a forthcoming SKA1 radio weak lensing survey, sufficient for measuring cluster lensing. We choose the following conservative approximation of the frequency bandwidth: 950–1190 MHz, as proposed in Bonaldi et al. (2016). This seems to be the optimum observation frequency for a weak lensing survey with SKA1-MID, in case only 30 per cent of the full bandwidth of SKA Band 2 is usable (because of RFI problems, other surveys commensality, etc.). Visibilities are sampled every τacc = 60 s and we consider one large channel because smearing effects are not included on shorter time and bandwidth scales. The observations are simulated by using equations (1) and (6) and adding an uncorrelated Gaussian noise whose variance is given by Wrobel & Walker (1999):   \begin{eqnarray} \sigma ^2 = \frac{\mathrm{SEFD}^2}{2\eta ^2\Delta \nu \tau _\mathrm{acc}}, \end{eqnarray} (7)where Δν is the frequency channel bandwidth, SEFD = 400 Jy is the system equivalent flux density for SKA-MID dishes and η = 0.9 is the system efficiency (Braun 2014). For simplicity, we assume that the SKA1-MID core is composed of only SKA dish antennas even if part of it actually contains 64 MeerKAT dishes. SF galaxy populations are generated according to distributions estimated from archival data of VLA radio surveys at 1.4 GHz. As described in Rivi et al. (2016a), we estimated distributions of the flux S and the scalelength α of sources modelled by an exponential profile by the analysis of faint sources (order of tens μJy) catalogue of the SWIRE field survey (Owen & Morrison 2008). This is the radio catalogue with the largest number of SF galaxies, being related to the deepest survey so far in terms of the radio source density (source flux cut ∼ 15 μJy), and it contains source size measurements from the imaged data. The flux distribution is fitted by a power law: p(S) ∝ S−1.34. The scalelength is obtained by the relation with the measured full width at half-maximum:4 FWHM = 2αln (2), and its distribution is fitted dependently on source flux by a lognormal function with mean μ = ln (αmed) and variance σ ∼ 0.3, where αmed is given by equation (5). The variance value is suitably chosen in the middle of a range well-representing scalelength distributions for different flux values. The modulus e of the intrinsic ellipticities is generated according to a function proposed by Miller et al. (2013):   \begin{eqnarray} P(e) = \frac{Ne\left(1-\exp \left[ \frac{e-e_\mathrm{max}}{c}\right]\right)}{(1+e)(e^2+e_0^2)^{1/2}}. \end{eqnarray} (8)The parameter values, c = 0.2298 and e0 = 0.0732 are obtained from fitting the VLA-COSMOS field data. Although this survey is less deep than the SWIRE, but still detecting μJy sources (flux cut ∼75 μJy), we rely on a recent re-analysis of the L-band radio visibility data where the level of systematics in the measurement of the galaxy position angle is significantly reduced (Tunbridge et al. 2016). In fact, the previous analysis (Schinnerer et al. 2007) as well as the one of the VLA-SWIRE were mainly focused on faint source counts. The normalization factor is N = 2.595. We generate galaxy populations with flux densities ranging between 10 and 200 μJy. According to our flux distribution we obtain a source number density of 2.7 gal arcmin−2. To be consistent, we adopt this source number density in our simulation, although more accurate modelling from recent radio continuum surveys suggest that a higher source number density should be detected at such a flux density cut in real observations (Condon et al. 2012; Mancuso et al. 2017). This is the expected source number density for the proposed 2-yr SKA1 radio weak lensing survey covering 5000 $$\deg ^2$$ (Brown et al. 2015). 4 GALAXY SHAPE MEASUREMENT In this section, first we select the facet size to be used in the source extraction by simulating visibilities of individual sources. Then, we simulate populations of galaxies located simultaneously in the field of view in order to show the efficacy of our source extraction algorithm as a function of the source number density. 4.1 Facet size The facet size is affected by the weighting scheme (Briggs, Schwab & Sramek 1999) adopted in the gridding phase. For example, natural weighting optimises the sensitivity for detecting weak sources by emphasizing the data from short baselines. In this case, a relatively large facet size is expected even for covering a single galaxy because a large contribution to the signal is from long wavelength modes which must be adequately sampled by small facet cells. In effect, the source in the image domain turns to be convolved with a large natural-weighted PSF with a broad low-level plateau (see left-hand panel of Fig. 1). Uniform weighting will require instead a much smaller facet size because it emphasizes data from long baselines, where most of the source shape signal is contained. This is reflected by the small uniform-weighted PSF (see right-hand panel of Fig. 1). On the other side, the weighting scheme used in the faceting procedure shouldn't affect the model fitting, provided that the measurement uncertainties are propagated correctly in the likelihood computation (see equation 21 of Paper I for the natural case) and model visibilities are consistent with the observed data. This may not be true for measurement methods in the image domain, as shown in Tunbridge et al. (2016). Figure 1. View largeDownload slide Image of the PSF corresponding to our SKA1 uv coverage Left: natural weighting. The PSF has a broad, low-level plateau because uv points have short spacings close to zero, as they tend to spend more time per unit area near the uv origin. Right: uniform weighting. Figure 1. View largeDownload slide Image of the PSF corresponding to our SKA1 uv coverage Left: natural weighting. The PSF has a broad, low-level plateau because uv points have short spacings close to zero, as they tend to spend more time per unit area near the uv origin. Right: uniform weighting. Since we are interested in the detection of faint sources for radio weak lensing, we adopt a natural weighting scheme. To minimize the number of sources falling in the same facet, we define a facet size dependent on source flux, as it is related to the size of the source. We split the flux total range of the simulated population, i.e. 10–200 μJy, in seven bins as shown in Table 1. Facet uv point coordinates have to be re-computed only once per bin as the model fitting is performed according to source flux order. We chose larger bins at large fluxes because the sizes of such sources increase more gradually with the flux compared to the ones with low flux (see equation 5). To estimate the facet size for each bin, we simulate raw visibilities of a single galaxy in the primary beam in order to avoid any source contamination effects, and vary the noise added to the visibilities (to have SNR ≥ 100) in order to see the effect of the source size only. We measure the galaxy ellipticity after averaging visibilities in the facet. The best-fitting slope5 for the ellipticity measurements of 1000 galaxies is computed for different facet sizes and source flux ranging between the flux bin bounds. We select the facet size when a fixed best-fitting slope threshold of about 0.97 is reached, as listed in Table 1. Table 1. Facet sizes dependent on source flux range and corresponding best-fitting slopes for 1000 sources. S (μJy)  No. cells  Best-fitting slope      e1  e2  200–150  600  0.9774 ± 0.0025  0.9675 ± 0.0025  150–100  550  0.9795 ± 0.0030  0.9664 ± 0.0030  100–80  500  0.9797 ± 0.0032  0.9660 ± 0.0032  80–60  460  0.9774 ± 0.0029  0.9614 ± 0.0029  60–40  420  0.9760 ± 0.0031  0.9631 ± 0.0031  40–20  350  0.9756 ± 0.0028  0.9557 ± 0.0029  20–10  280  0.9765 ± 0.0030  0.9510 ± 0.0030  S (μJy)  No. cells  Best-fitting slope      e1  e2  200–150  600  0.9774 ± 0.0025  0.9675 ± 0.0025  150–100  550  0.9795 ± 0.0030  0.9664 ± 0.0030  100–80  500  0.9797 ± 0.0032  0.9660 ± 0.0032  80–60  460  0.9774 ± 0.0029  0.9614 ± 0.0029  60–40  420  0.9760 ± 0.0031  0.9631 ± 0.0031  40–20  350  0.9756 ± 0.0028  0.9557 ± 0.0029  20–10  280  0.9765 ± 0.0030  0.9510 ± 0.0030  View Large Note that the fitting for the first ellipticity component is better than the second one because of a slightly anisotropy of the PSF as discussed in Section 6. The selected facet sizes are consistent with the relation between the uv grid cell size Δu (in units of wavelengths) and the related field of view (in radians): ψ = 1/Δu (Briggs et al. 1999). We also note that for small sources (low flux), it is actually better to use smaller facets even in the case of a single source in the field of view. This is shown in Fig. 2 where the difference between binned measurements and true values of the first ellipticity component of 1000 galaxies, with realistic flux distribution in the range 10–200 μJy, are plotted both for the case where the facet size is constant and equal to 600 (left-hand panel) and when the facet has a variable size dependent on the source flux (right-hand panel). Similar plots are obtained for the second component. In the latter case we obtained 25 bad measurements (see Section 2.2) and the best-fitting slopes of the two ellipticity components are 0.9552 ± 0.0057 and 0.9426 ± 0.0061, respectively, whereas for the case of 600 × 600 facet the best-fitting slopes for the same galaxy population and noise are 0.9306 ± 0.0054 and 0.9135 ± 0.0056 and the number of bad measurements is three times larger. Figure 2. View largeDownload slide Binned measurements (minus true values) of the ellipticity first component of 1000 individual galaxies randomly located in the field of view and with flux ranging between 10 and 200 μJy. Left: facet size 600 × 600 for all sources. Right: facet size dependent on source flux. Figure 2. View largeDownload slide Binned measurements (minus true values) of the ellipticity first component of 1000 individual galaxies randomly located in the field of view and with flux ranging between 10 and 200 μJy. Left: facet size 600 × 600 for all sources. Right: facet size dependent on source flux. These results are due to the fact that we do not model exactly the primary beam because the model visibilities are directly sampled on the uv facet points. This means that in the image domain the sidelobes of the source model are not suppressed by any apodisation, whereas the gridding of the original uv coverage causes the full image to be apodised by a broad 2D sinc function which has the effect in the data of suppressing background sources that are a long way from the phase centre and the distant sidelobes from the primary source. The grid sampling causes the resulting image domain facet to become a small, but aliased version of the apodised image. The aliased model is an incorrect description of the apodised and aliased data, and the discrepancy will get worse for smaller facets and at large distances from the source. Fig. 2 shows that a suitable facet size dependent on the source flux/size may be a trade-off between these two effects. We could improve the model by applying the same gridding operations as in the data (sampling on the original uv coverage and then averaging in the facet), but this will add a large amount of computational time. Our results show that the adopted model approximation is acceptable, provided that the facet sizes are sufficiently large to not affect the significant sidelobes in the image domain. Otherwise we expect the discrepancy between data and model to become severe and the biases may become less robust and hence less calibratable. 4.2 Dependence on source number density We estimate the efficacy of the source extraction method by measuring the slope of the best-fitting line of 104 galaxy ellipticity measurements as a function of the source number density. For each measurement, we simulate sources located simultaneously in the field of view according to a uniform distribution. Results are plotted in Fig. 3 and show reliable fits, independent of the source density up to 2.8 gal arcmin−2. In this case, the best-fitting slopes of the correlation for each ellipticity component are 0.9365 ± 0.0017 and 0.9262 ± 0.0017 respectively and the number of bad measurements is about 1 per cent. At higher densities galaxy ellipticity measurement starts to deteriorate, as residuals of nearby galaxies affect the model fitting, but may still be good enough for shear measurement because of the improved statistics (as shown in Section 5). Shape measurements of galaxies may be improved by a joint fitting within facets by applying the Hamiltonian Monte Carlo technique (Neal 2011). RadioLensfit results used as starting points should reduce the burn-in phase and accelerate convergence. Since the number of sources in the facet will be relatively small this approach becomes more feasible and preliminary results with this method show a better accuracy in the galaxy ellipticity fitting, although requiring a large computational time (Rivi et al., in preparation). Figure 3. View largeDownload slide Best-fitting slope for both galaxy ellipticity components as a function of the source number density. Figure 3. View largeDownload slide Best-fitting slope for both galaxy ellipticity components as a function of the source number density. Using a single channel, the serial version of RadioLensfit running on an Intel Xeon E5-2650 takes on average about 10 s gal−1 computing time for the model fitting. As discussed in Paper I, the shared memory parallelization with OpenMP allows us to exploit all the computational resources when the amount of memory for the source model fitting requires the usage of the full CPU. Its implementation has been optimized by distributing to each thread the likelihood computation and marginalization over the position parameters for different scalelength values of the model. It doesn't scale linearly with the number of threads because the likelihood mariginalization over the scalelength parameter is not parallel and there is an overhead for the creation and destruction of OpenMP threads at each iteration of the likelihood maximization and sampling. This version running on all the eight cores of the CPU takes on average about 2.4 s gal−1. 5 SHEAR Following Paper I, for shear measurement we simulate galaxy populations as described in Section 3 in a field of view of 1 $$\deg ^2$$. We generate populations free of shape noise (Massey et al. 2007; Nakajima & Bernstein 2007): for each ellipticity modulus, 10 equally spaced galaxy orientations are generated so that the corresponding ellipticity values are distributed uniformly on a circle, and galaxies whose ellipticity is on the same ring are given the same size and flux. We generate sources randomly located according to a uniform distribution. All measurements are performed with facet size dependent on the source flux as defined in Table 1. We measure the reduced shear $$\boldsymbol {g}=g_1 + \mathrm{i}g_2$$ as the weighted average of the galaxies’ ellipticities using weights that approximate the inverse-variance of each ellipticity measurement. Error bars are given by the standard deviation of the shear values estimated from 1000 bootstrap resamples. Fig. 4 shows shear measurements as a function of the source number density up to 8 gal armin−2, when no input shear is applied, i.e. g1 = g2 = 0. At high densities, the larger number of sources compensates the less accurate galaxy shape fitting (Fig. 3), still producing shear values consistent with the results obtained at the SKA1 source density corresponding to a population of about 104 galaxies. Figure 4. View largeDownload slide Shear components estimated from a galaxy population in 1 $$\deg ^2$$ as a function of the source number density for input g = 0. Figure 4. View largeDownload slide Shear components estimated from a galaxy population in 1 $$\deg ^2$$ as a function of the source number density for input g = 0. For this population, we measure the shear also for input reduced shear values with amplitude g = 0.04 and eight different orientations. The input shear $$\boldsymbol {g}$$ action on the intrinsic galaxy ellipticity $$\boldsymbol {e}^s$$ is simulated following Seitz & Schneider (1997):   \begin{eqnarray} \boldsymbol {e} = \frac{\boldsymbol {e}^s + \boldsymbol {g}}{1+\boldsymbol {g}^*\boldsymbol {e}^s}, \end{eqnarray} (9)where $$\boldsymbol {g}^*$$ is the shear complex conjugate. We compare results with the optimal case where each galaxy is at the phase centre and the only one contained in the field of view, considering the same galaxy population. Results are plotted in Fig. 5, both for SNR ≥ 10 and SNR ≥ 25, where measurements from individual source visibilities are green crosses and the ones from the same population but with all sources simultaneously in the primary beam are black crosses. Clearly error bars (cross arms) are larger at SNR ≥ 10 as the galaxy population is dominated by lower flux sources. Figure 5. View largeDownload slide Comparison of shear measurements: input values are blue points, measured values from single sources at the phase centre are green crosses, measured values from sources simultaneously in the f.o.v are black crosses. Figure 5. View largeDownload slide Comparison of shear measurements: input values are blue points, measured values from single sources at the phase centre are green crosses, measured values from sources simultaneously in the f.o.v are black crosses. The measured shear bias, defined as   \begin{eqnarray} g_i^m - g_i = m_i g_i + c_i, \qquad i=1,2, \end{eqnarray} (10)is shown in Table 2. At SNR ≥ 10, the multiplicative biases mi for the two shear components are, respectively, 1.6 and 1.4 times the ideal case of a single source in the field of view, while additive bias ci are almost consistent. Selecting galaxies with SNR ≥ 25 the population reduces to 5810 sources (i.e. 1.6 gal arcmin−2). This should affect the bias uncertainty only, as the bias on the shear measurement should not depend on the number of sources. At this regime, the m values reduce by a factor two instead of three as in the single source case. This is due to the source signal contamination by residuals of nearby galaxies, which is a new contribution to the shear bias that seems to have an effect on the multiplicative terms only. It may be mitigated by refining ellipticity measurements by joint fitting within larger facets, as explained in Section 4.2. Note that ‘neighbour bias’ affects optical lensing measurements also (see Miller et al. 2013; Jarvis et al. 2016; Mandelbaum et al. 2017; Zuntz et al. 2017; Samuroff et al. 2018). Table 2. Shear bias components estimated from a realistic population of ∼104 galaxies randomly distributed in 1 $$\deg ^2$$, corresponding to a source number density of 2.8 gal arcmin−2.     m1  c1  m2  c2  SNR ≥ 10  Single source  − 0.0904 ± 0.0186  0.00655 ± 0.00050  − 0.1297 ± 0.0171  0.00632 ± 0.00046    Multiple sources  − 0.1428 ± 0.0274  0.00872 ± 0.00073  − 0.1864 ± 0.0256  0.00578 ± 0.00067  SNR ≥ 25  Single source  − 0.0352 ± 0.0127  − 0.00006 ± 0.00034  − 0.0506 ± 0.0132  − 0.00070 ± 0.00035    Multiple sources  − 0.0677 ± 0.0242  0.00089 ± 0.00064  − 0.0992 ± 0.0230  − 0.00112 ± 0.00060      m1  c1  m2  c2  SNR ≥ 10  Single source  − 0.0904 ± 0.0186  0.00655 ± 0.00050  − 0.1297 ± 0.0171  0.00632 ± 0.00046    Multiple sources  − 0.1428 ± 0.0274  0.00872 ± 0.00073  − 0.1864 ± 0.0256  0.00578 ± 0.00067  SNR ≥ 25  Single source  − 0.0352 ± 0.0127  − 0.00006 ± 0.00034  − 0.0506 ± 0.0132  − 0.00070 ± 0.00035    Multiple sources  − 0.0677 ± 0.0242  0.00089 ± 0.00064  − 0.0992 ± 0.0230  − 0.00112 ± 0.00060  View Large As discussed in Paper I, the noise bias values exceed SKA1 survey requirements6 except for the additive component at SNR ≥ 25. However they are comparable to the ones obtained from optical surveys using lensfit (Fenech Conti et al. 2017) and im3shape (Zuntz et al. 2017), where a shear calibration correction reduced the multiplicative bias to well below the percent level. Standard approaches in the optical domain derive such calibration by inferring the bias from simulated data matching the observations (Bruderer et al. 2016) or parametrizing the bias as a function of the observed galaxy properties (Kuijken et al. 2015; Jarvis et al. 2016). Recently, a self-calibration approach, implemented in the Metacalibration method (Huff & Mandelbaum 2017; Sheldon & Huff 2017) and used in the analysis of the Dark Energy Survey7 (DES) (Zuntz et al. 2017), proved to be the most efficient, being able to recover the input shear in realistic simulations to better than a part in a thousand. It also isotropises the PSF to remove any additive bias. The key idea of the method is to compute the shear estimator response for a shape measurement directly from observed data perturbed with a small known shear. This way all the features present in real data are already included, which are instead extremely difficult to model accurately in external simulations, and it can be applied to any shear measurement method based on averages of galaxy shapes. A similar approach may then be implemented quite easily in the interferometer data analysis, with the advantage that for the additive bias at radio wavelengths we know the PSF much better than at optical wavelengths and we can make the PSF isotropic directly by weighting the uv-plane (as discussed in Section 6). 6 ADDITIVE BIAS DEPENDENCE ON THE IMAGE-PLANE PSF SHAPE It is well known from weak lensing optical surveys that shear additive bias is dependent on the PSF shape (Miller et al. 2013). For radio interferometers, the PSF is deterministically defined by the uv coverage of the telescope (i.e. antennas locations and pointing direction). For example, if we increase the antenna pointing declination to a larger zenith distance the PSF shape becomes compressed along the y-axis. We measure the additive noise bias for different pointing declinations at various SNR values. Starting from our uv coverage (corresponding to the same declination as the observatory latitude), we simulate the effect on the uv points when the phase centre declination is increased by an angle ϕ. A plot of the image-plane PSF ellipticity components as functions of the zenith distance is given in Fig. 6. The R-squared size of the PSF slightly increases from 14.83 arcsec2 to 15.54 arcsec2. These values are computed from the quadrupole moments of the image domain as follows (Schneider 2006):   \begin{eqnarray} \boldsymbol {e} = \frac{Q_{xx}-Q_{yy}+2\mathrm{i}Q_{xy}}{Q_{xx}+Q_{yy}+2(Q_{xx}Q_{yy}-Q^2_{xy})^{\frac{1}{2}}}, \end{eqnarray} (11)  \begin{eqnarray} R = Q_{xx}+Q_{yy}. \end{eqnarray} (12) Figure 6. View largeDownload slide PSF ellipticity components versus the zenith distance at the starting of the observation. Figure 6. View largeDownload slide PSF ellipticity components versus the zenith distance at the starting of the observation. We simulate individual source visibilities, to avoid nearby source contamination effects, and assume a constant maximum facet size 1000 × 1000 to ensure that the galaxy convolved with the PSF is contained in the facet even when the PSF becomes highly anisotropic. We observe that the additive bias is dependent on source size. In fact measurements at the same SNR obtained by lowering the noise instead of increasing the source flux cut, produce larger bias values, meaning that the additive bias worsens when source sizes decreases. This is consistent with the analysis presented in Massey et al. (2013). Noise bias causes a correlation between measured shear and PSF ellipticity even when we correct for the PSF in the model fitting. This becomes noticeable at low SNR, where the first ellipticity component increases significantly towards larger negative values, as the PSF becomes more compressed along the y direction (see left-hand panel of Fig. 7). At large SNR the additive bias almost disappears independently of the PSF shape (see right-hand panel of Fig. 7). This is in good agreement with what we should expect. Because of our long integration time the PSF is not isotropic even when the declination equals the observatory latitude, as in the baseline simulation, besides the fact that the distribution of the SKA-MID baselines on the ground are not circularly symmetric. Therefore, a large additive bias is still measured at small zenith distances. The PSF anisotropy may be reduced by combining snapshots obtained over a range of hour angles, however this may not be sufficient to reach an additive bias acceptable for weak lensing surveys. To further reduce the noise bias at low SNR, we also need to weight the uv plane to ensure that the PSF is more isotropic. A standard technique in radio imaging to improve PSF shape is to use tapering functions (Briggs et al. 1999) to define uv points weights, although a more specific weighting scheme may be required. Moreover, as for the multiplicative bias, we can calibrate the additive shear bias with simulations. This is more feasible than in optical surveys because the PSF is deterministic at radio wavelengths. However, any such calibration would be strongly dependent on distributions of source properties, so isotropising the PSF is a much better option. Figure 7. View largeDownload slide Shear additive bias as a function of the zenith distance at different lower signal to noise. Facet size fixed at 1000 × 1000. Figure 7. View largeDownload slide Shear additive bias as a function of the zenith distance at different lower signal to noise. Facet size fixed at 1000 × 1000. Note that when using variable size faceting, the facet size must be dependent not only on the source size but also on the PSF shape. In fact, as the PSF becomes anisotropic the facet size may become too small relative to the size of the source convolved with the PSF, modifying the effective shape of the source. For example, the left-hand panel of Fig. 8 shows what happens at SNR ≥ 10 when we maintain the same flux dependent facet size (Table 1) for all pointing declinations: at large zenith distances the shortest baselines occupy smaller uv frequencies and therefore can measure wavelengths longer than the limit imposed by the small facet size used to extract the majority of the galaxy population. If we increase the facet size according to the source declination, e.g. 50 cells per side every 10 deg declination increment, we obtain consistent results with the case of one single large facet (see right-hand panel of Fig. 8). Figure 8. View largeDownload slide Shear additive bias at SNR ≥ 10 as a function of the zenith distance. Left: facet size dependent on source flux according to Table 1 for all declinations. Right: facet size increased with source declination produce similar results as for a constant large facet size. Figure 8. View largeDownload slide Shear additive bias at SNR ≥ 10 as a function of the zenith distance. Left: facet size dependent on source flux according to Table 1 for all declinations. Right: facet size increased with source declination produce similar results as for a constant large facet size. 7 CONCLUSIONS We have extended the presentation of the RadioLensfit method, introduced in Paper I for the simple case of individual galaxies located at the phase centre, to the real case where many galaxies are randomly located in the field of view. This has been done by isolating the visibilities of each source and shifting the phase centre to the source position so that a coarse grid can still be used to reduce nearby galaxy residuals contamination and accelerate ellipticity measurement computation. Source extraction has been performed by removing apart from the data the simulated visibilities of the sky model, but the source of interest, given the positions and fluxes of all sources in the field of view from the image, and down-weighting what remains of nearby source-contamination by averaging visibilities in a coarse grid (facet). For gridding, we adopted a natural weighting to maximize sensitivity and estimated the smallest facet size dependent on source flux thresholds in order to minimize the number of nearby galaxies included in the facet. We tested the source extraction procedure, simulating visibilities of SF galaxy populations observed by SKA1-MID in the first 30 per cent of frequency Band 2. We adopted flux and scalelength parameters distributions estimated from the VLA SWIRE catalogue and used the lensfit ellipticity prior with coefficients fitted from a new version of the VLA COSMOS catalogue optimized on shape measurements. We showed the efficacy of our source extraction algorithm as a function of the source number density, obtaining a reliable galaxy ellipticity fitting for the density expected from the current proposal of the SKA1 radio weak lensing survey. Shear measurements from 8-h observation of 1 deg2 show that the bias due to the extraction procedure mainly affects the multiplicative bias as no significant change has been observed for the additive bias when comparing with the bias obtained for the ideal case of a single source at the phase centre at a time. This bias may be mitigated by a second step in the galaxy ellipticity measurement, where a joint fitting within the facets is performed with HMC, starting from the values obtained with RadioLensfit. However, multiplicative noise bias calibration is also required as for the optical domain. We finally observed that because of our uv coverage the PSF is slightly anisotropic even if pointing close to the zenith, therefore we obtain a large additive bias (on average about 0.0068 ± 0.0006 at SNR ≥ 10). Although a suitable choice of the integration time (split and distributed along a longer period of time) may reduce the PSF anisotropy, a uv weighting scheme may still be required to satisfy weak lensing requirements. It should be optimized to avoid any significant reduction of the signal to noise. ACKNOWLEDGEMENTS We thank Tom Kitching for useful discussions and Ben Tunbridge for providing the distribution parameters of our ellipticity prior, obtained by fitting VLA-COSMOS data as for the prior function presented in his paper. We also thank Sphesihle Makhathini for the support in the simulation of the SKA1-MID uv coverages. MR acknowledges the support of the Science and Technology Facilities Council (STFC) via an SKA grant. LM acknowledges STFC grant ST/N000919/1. Footnotes 1 https://www.skatelescope.org 2 http://www.e-merlin.ac.uk/legacy/projects/superclass.html 3 SKA-MID latitude is −30°49΄48.00΄S. 4 The FWHM is derived from the Gaussian model fit of the source after PSF deconvolution. 5 Consistently with Paper I, we refer to the ellipticity best-fitting slope instead of the multiplicative bias when measuring galaxy shapes. This terminology is used to clearly distinguish it from the shear multiplicative bias, which is obtained from the best fit of shear measurements (each being the weighted average of galaxy ellipticities). 6 For a 2-yr SKA1-MID weak lensing survey over 5000 $$\deg ^2$$ and zmed = 1.0, the requirements for cosmological parameters measurements to be dominated by statistical rather than systematic errors are: multiplicative bias m < 0.0067, additive bias c < 0.00082 (Brown et al. 2015). They are derived using the rules provided in Amara & Refregier (2008). 7 https://www.darkenergysurvey.org REFERENCES Amara A., Refregier A., 2008, MNRAS , 391, 228 https://doi.org/10.1111/j.1365-2966.2008.13880.x CrossRef Search ADS   Becker R., White R., Helfand D. J., 1995, ApJ , 450, 559 https://doi.org/10.1086/176166 CrossRef Search ADS   Blain A., 2002, ApjL , 570, L51 https://doi.org/10.1086/341103 CrossRef Search ADS   Bonaldi A., Harrison I., Camera S., Brown M., 2016, MNRAS , 463, 3686 https://doi.org/10.1093/mnras/stw2104 CrossRef Search ADS   Braun R., 2014, SKA Memo SKA-TEL.SCI-SKO-SRQ-001  Bridle A., Schwab F., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser. Vol. 180, Synthesis Imaging in Radio Astronomy II . Astron. Soc. Pac., San Francisco, p. 371 Briggs D., Schwab F., Sramek R., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser. Vol. 180, Synthesis Imaging in Radio Astronomy II . Astron. Soc. Pac., San Francisco, p. 129 Brown M., Battye R., 2011, MNRAS , 410, 2057 Brown M. et al.  , 2015, Advancing Astrophysics with the Square Kilometre Array, Vol. 2, Proc. Sci. , SISSA, Trieste, p. 1365 Bruderer C., Chang C., Refregier A., Amara A., Bergé J., Gamper L., 2016, ApJ , 817, 25 https://doi.org/10.3847/0004-637X/817/1/25 CrossRef Search ADS   Camera S., Harrison I., Bonaldi A., Brown M., 2017, MNRAS , 464, 4747 https://doi.org/10.1093/mnras/stw2688 CrossRef Search ADS   Chang T., Refregier A., 2002, ApJ , 570, 447 https://doi.org/10.1086/339496 CrossRef Search ADS   Chang T., Refregier A., Helfand D., 2004, ApJ , 617, 794 https://doi.org/10.1086/425491 CrossRef Search ADS   Condon J. et al.  , 2012, AJ , 758, 23 https://doi.org/10.1088/0004-637X/758/1/23 CrossRef Search ADS   Cornwell T. J., Perley R. A., 1992, A&A , 261 Demetroullas C., Brown M., 2016, MNRAS , 456, 3100 https://doi.org/10.1093/mnras/stv2876 CrossRef Search ADS   Fenech Conti I., Herbonnet R., Hoekstra H., Merten J., Miller L., Viola M., 2017, MNRAS , 467, 1627 Feroz F., Marshall P., Hobson M., 2008, preprint (arXiv:0810.0781) Harrison I., Brown M., 2015, SKA ECP150007, preprint (arXiv:1507.06639) Harrison I., Camera S., Zuntz J., Brown M., 2016, MNRAS , 463, 3674 https://doi.org/10.1093/mnras/stw2082 CrossRef Search ADS   Heystek L., 2015, SKA Memo SKA-TEL-INSA-0000537  Huff E., Mandelbaum R., 2017, ApJ , preprint (arXiv:1702.02600) Jarvis M. et al.  , 2016, MNRAS , 460, 2245 https://doi.org/10.1093/mnras/stw990 CrossRef Search ADS   Kuijken K. et al.  , 2015, MNRAS , 454, 3500 https://doi.org/10.1093/mnras/stv2140 CrossRef Search ADS   Mancuso C. et al.  , 2017, AJ , 842, 95 https://doi.org/10.3847/1538-4357/aa745d CrossRef Search ADS   Mandelbaum R. et al.  , 2017, MNRAS , preprint (arXiv:1710.00885) Massey R. et al.  , 2007, MNRAS , 376, 13 https://doi.org/10.1111/j.1365-2966.2006.11315.x CrossRef Search ADS   Massey R. et al.  , 2013, MNRAS , 429, 661 https://doi.org/10.1093/mnras/sts371 CrossRef Search ADS   Melchior P., Viola M., 2012, MNRAS , 424, 2757 https://doi.org/10.1111/j.1365-2966.2012.21381.x CrossRef Search ADS   Melchior P., Böhnert A., Lombardi M., Bartelmann M., 2010, A&A , 510, A75 https://doi.org/10.1051/0004-6361/200912785 CrossRef Search ADS   Miller L. et al.  , 2013, MNRAS , 429, 2858 https://doi.org/10.1093/mnras/sts454 CrossRef Search ADS   Morales M., 2006, ApJ , 650, L21 https://doi.org/10.1086/508614 CrossRef Search ADS   Nakajima R., Bernstein G., 2007, ApJ , 133, 1763 https://doi.org/10.1086/511957 CrossRef Search ADS   Neal R. M., 2011, in Brooks S., Gelman A., Jones G., Meng X.-L., eds, MCMC using Hamiltonian dynamics. Astronomical Journal of the American Astronomical Society , Chapman & Hall/CRC Press Owen F., Morrison G., 2008, AJ , 136, 1889 https://doi.org/10.1088/0004-6256/136/5/1889 CrossRef Search ADS   Patel P., Bacon D. J., Beswick R. J., Muxlow T. W. B., Hoyle B., 2010, MNRAS , 401, 2572 https://doi.org/10.1111/j.1365-2966.2009.15836.x CrossRef Search ADS   Patel P. et al.  , 2015, Advancing Astrophysics with the Square Kilometre Array, Vol. 2, Proc. Sci. , SISSA, Trieste, p. 1427 Refregier A., Kacprzak T., Amara A., Bridle S., Rowe B., 2012, MNRAS , 425, 1951 https://doi.org/10.1111/j.1365-2966.2012.21483.x CrossRef Search ADS   Rivi M., Miller L., Makhathini S., Abdalla F., 2016a, EXTRA-RADSUR2015, Proc. Sci. , SISSA, Trieste Rivi M., Miller L., Makhathini S., Abdalla F. B., 2016b, MNRAS , 463, 1881 (Paper I) https://doi.org/10.1093/mnras/stw2041 CrossRef Search ADS   Samuroff S. et al.  , 2018, MNRAS , 475, 4524 https://doi.org/10.1086/516587 preprint (arXiv:1708.01534) Schinnerer E. et al.  , 2007, ApJS , 172, 46 https://doi.org/10.1086/516587 CrossRef Search ADS   Schneider P., 2006, in Meylan G., Jetzer P., North P., Schneider P., Kochanek C. S., Wambsganss J., eds, Saas-Fee Advanced Course 33: Gravitational Lensing: Strong, Weak and Micro. Springer-Verlag, Berlin, p. 269 Seitz C., Schneider P., 1997, A&A , 318, 687 Sheldon E. S., Huff E. M., 2017, ApJ , 841, 24 https://doi.org/10.3847/1538-4357/aa704b CrossRef Search ADS   Smirnov O. M., 2011, A&A , 527, A106 https://doi.org/10.1051/0004-6361/201016082 CrossRef Search ADS   Tunbridge B., Harrison I., Brown M., 2016, MNRAS , 463, 3339 https://doi.org/10.1093/mnras/stw2224 CrossRef Search ADS   Wrobel J., Walker R., 1999, in Taylor G., Carilli C., Perley R., eds, ASP Conf. Ser., Vol. 180, Synthesis Imaging in Radio Astronomy II , Astron. Soc. Pac., San Francisco, p. 171 Zuntz J. et al.  , 2017, MNRAS , preprint (arXiv:1708.01533) © 2018 The Author(s) Published by Oxford University Press on behalf of the Royal Astronomical Society

### Journal

Monthly Notices of the Royal Astronomical SocietyOxford University Press

Published: May 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations