304.3. Galaxy cluster weak lensing with HSM

304.3. Galaxy cluster weak lensing with HSM¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: large
LSST Science Pipelines version: Release r30.0.11
Last verified to run: 2026-09-24
Repository: github.com/lsst/tutorial-notebooks
DOI: 10.11578/rubin/dc.20250909.20

Learning objective: Evaluate the weak lensing effect of galaxy cluster PSZ2 G309.43-72.86 on background galaxies.

LSST data products: deep_coadd images, Object table

Packages: lsst.daf.butler, lsst.rsp, lsst.afw.display, lsst.geom.

Credit: Originally developed by the Rubin Community Science team. Please consider acknowledging them if this notebook is used for the preparation of journal articles, software releases, or other notebooks.

Get Support: Everyone is encouraged to ask questions or raise issues in the Support Category of the Rubin Community Forum. Rubin staff will respond to all questions posted there.

1. Introduction¶

PSZ2 G309.43-72.86 (Léget 2026) is a known rich cluster at redshift=0.35 covered by the Data Preview 2 (DP2). This notebook provides a demonstration of detecting the lensing signal with PSZ2 G309.43-72.86 and DP2 data.

Galaxy clusters are the largest gravitionally bound objects in the Universe. According to General Relativity, massive objects warp spacetime, and light follows the curvature of spacetime. Thus, massive objects act like lenses that bend the path of the light emitted from distant sources. This effect is called "gravitational lensing".

The large mass of a galaxy cluster coherently distorts the images of backgroud galaxies, and this distortion occurs over a large sky area around the cluster. The lensing effect on the shape of a single background galaxy far from the cluster center is small, but the lensing signal can be detected by averaging the shape of many background galaxies. This is because the original shapes of the galaxies are randomly oriented. This effect is called weak gravitational lensing (WL), which shows up in the statistics of a large sample of galaxy shapes. More detailed introduction to cluster WL can be found in review papers (e.g., Bartelmann & Schneider 2001; Umetsu 2020).

Note that WL also happens between galaxies and when the light passes through the large-scale structure of the Universe (cosmic shear). Compared to galaxy-galaxy lensing and cosmic shear, the cluster WL signal is generally about 1 magnitude higher (at the level of 10 times).

Related tutorials: The 100-level tutorials on using the butler and Table Access Protocol (TAP). The 200-level tutorials on object tables and deep coadds.

1.1. Import packages¶

Import general packages (numpy, matplotlib, astropy), the LSST Science Pipelines packages for the Butler, TAP, display, and geometry.

In [1]:

import numpy as np
from scipy.stats import binned_statistic
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord

from lsst.daf.butler import Butler
from lsst.rsp import RSPDiscovery
import lsst.afw.display as afwDisplay

1.2. Define parameters and functions¶

Set up Firefly for displaying images.

In [2]:

afwDisplay.setDefaultBackend("firefly")
afw_display = afwDisplay.Display(frame=1)

Get the butler via repo and collection.

In [3]:

butler = Butler("dp2", collections="dp2")

Get the TAP service.

In [4]:

discovery = RSPDiscovery("dp2")
service = discovery.get_tap_client()

2. Load data¶

The Brightest Cluster Galaxy (BCG) for PSZ2 G309.43-72.86 has coordinates of Right Ascension (RA), DEClination (DEC) = 10.2082deg, -44.1307deg

In [5]:

ra_bcg = 10.2082
dec_bcg = -44.1307

2.1. Choose the i-band¶

Using the i or r band shapes is common in weak lensing studies (Mandelbaum et al. 2018). The main reason is the balance between seeing and sky brightness for shape measurement (Fu et al. 2022). Another reason is that Differential Chromatic Refraction (DCR) is weaker in redder bands, which reduces the Point Spread Function (PSF) elongation along the zenith and makes the galaxy's pre-PSF shape easier to measure (Plazas et al. 2012). The pre-PSF shape means the galaxy shape before affected by the atmosphere.

Choose to load data obtained in the i-band filter for this lensing analysis.

In [6]:

band_shape = 'i'

2.2. Object catalog¶

Use the TAP service to make a cone search of the Object catalog for objects detected in the deep_coadd images within 0.4 degrees of the BCG.

Define the query to return shape parameters in the selected band, photometry in the riz bands, and other useful columns. Compute the signal-to-noise ratio (SNR) in the query since it is not included in the schema (otherwise, to save time, it is not recommended to do computation in the query).

In [7]:

radius_deg = 0.4

query = "SELECT objectId, coord_ra, coord_dec, refExtendedness, " \
        "deblend_peak_center_x, deblend_peak_center_y, " \
        f"{band_shape}_blendedness, " \
        f"{band_shape}_hsmShapeRegauss_e1, " \
        f"{band_shape}_hsmShapeRegauss_e2, " \
        f"{band_shape}_hsmShapeRegauss_sigma, " \
        f"{band_shape}_hsmShapeRegauss_flag, " \
        "r_cModelMag, i_cModelMag, z_cModelMag, " \
        "r_cModelMagErr, i_cModelMagErr, z_cModelMagErr, " \
        "r_cModelFlux/r_cModelFluxErr AS r_cModel_SNR, " \
        "i_cModelFlux/i_cModelFluxErr AS i_cModel_SNR, " \
        "z_cModelFlux/z_cModelFluxErr AS z_cModel_SNR, " \
        "r_cModel_flag, i_cModel_flag, z_cModel_flag " \
        "FROM dp2.Object " \
        "WHERE CONTAINS(POINT('ICRS', coord_ra, coord_dec), " \
        f"CIRCLE('ICRS', {ra_bcg}, {dec_bcg}, {radius_deg})) = 1 "

Run the query job. Check if the query job is completed and turn the result into a table. Print the length of the table (number of objects).

In [8]:

job = service.submit_job(query)
job.run()
job.wait(phases=['COMPLETED', 'ERROR'])
print('Job phase is', job.phase)
if job.phase == 'ERROR':
    job.raise_if_error()
assert job.phase == 'COMPLETED'
obj_cat = job.fetch_result().to_table()
print('Number of objects is', len(obj_cat))

Job phase is COMPLETED

Number of objects is 178612

Change the column names for simpler access later. These are coordinates of the peak after source detection.

In [9]:

obj_cat.rename_column('deblend_peak_center_x', 'x')
obj_cat.rename_column('deblend_peak_center_y', 'y')

2.3. Deep coadd image¶

Use the butler to get the deep_coadd images in the band_shape filter that contains the BCG of PSZ2 G309.43-72.86. Print the dataIds of the images and pick the first one.

In [10]:

dataset_refs = butler.query_datasets("deep_coadd",
                                     where="band.name = band AND \
                                     patch.region OVERLAPS POINT(ra, dec)",
                                     bind={"band": band_shape,
                                           "ra": ra_bcg,
                                           "dec": dec_bcg})
for ref in dataset_refs:
    print('dataId:', ref.dataId)

deep_coadd = butler.get(dataset_refs[0])

dataId: {band: 'i', skymap: 'lsst_cells_v2', tract: 2877, patch: 34}

Display the image in the Firefly tab and set the mask to be transparent.

In [11]:

afw_display.image(deep_coadd)
afw_display.setMaskTransparency(100)

Mark the location of the BCG with a big orange circle.

In [12]:

wcs = deep_coadd.sky_projection
coord = SkyCoord(ra_bcg, dec_bcg, frame='icrs', unit='deg')
pixels = wcs.sky_to_pixel(coord)
afw_display.dot("o", pixels.x, pixels.y, size=100, ctype="orange")

3. Identify candidate background galaxies¶

Select Objects that are likely to be extended (i.e., not point-like; not stars), and then select galaxies that are likely to be in the "background" (i.e., behind the cluster) based on their color.

Clusters generally show a "red sequence" (RS) in the "color-magnitude diagram" (CMD) due to evolution -- those red galaxies are the oldest and reddest in the cluster (Kodama and Arimoto 1997; Gladders and Yee 2000). Galaxies that are redder than the RS (i.e., are likely at higher redshifts and more redshifted) are candidate background galaxies. In practice, galaxy photometric redshifts (photo-z) are used to identify background galaxies. For this tutorial, only extendedness and color are used.

3.1. Select galaxies¶

From the returned Objects select the subset that are likely to be extended objects (galaxies, not stars) using the refExtendedness column (0 for point-like, 1 for extended). This is the extendedness in the reference band. The extendness is either 0 or 1, and it is based on the difference between the PSF magnitude and the Composite-Model (CModel) magnitude (Abazajian et al. 2004, Bosch et al. 2018). Note, the refExtendedness parameter has not yet been characterized for purity or completeness. See also the 311-series tutorial on extendedness and the documentation on catalog flag columns.

Use the riz CModel photometry (Bosch et al. 2018), which is appropriate for measurements of galaxies (extended objects), and retain only objects with SNR>5 to select well-detected objects. Select the objects that do not have flag parameters *_cModel_flag set to 1, to ensure the quality.

In [13]:

band_list = ["r", "i", "z"]

sel = obj_cat["refExtendedness"] == 1
for band in band_list:
    sel &= obj_cat[f"{band}_cModel_flag"] == 0
    sel &= obj_cat[f"{band}_cModel_SNR"] > 5.
obj_cat = obj_cat[sel]

print('Number of objects is', len(obj_cat))

Number of objects is 78908

3.2. Color-magnitude diagrams¶

Create CMDs of galaxies using the three bands $riz$, and two colors r-i and i-z with which the RS will be clearly visible.

Finding the RS in multiple colors instead of a single color improves the determination of background galaxies (Fu et al. 2024).

In [14]:

plt_dict = {0: ['r', 'i'], 1: ['i', 'z']}
fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=False, sharex=False)
for i in [0, 1]:
    band1 = plt_dict[i][0]
    band2 = plt_dict[i][1]
    ax[i].scatter(obj_cat[band1+"_cModelMag"],
                  obj_cat[band1+"_cModelMag"] - obj_cat[band2+"_cModelMag"],
                  marker=".", s=0.3, alpha=0.5, color='k')
    ax[i].set_xlabel(band1 + ' [AB mag]')
    ax[i].set_ylabel(band1 + '-' + band2 + ' [AB mag]')
    ax[i].set_xlim([17, 25])
    ax[i].set_ylim([-0.5, 1.5])
plt.tight_layout()
plt.show()

Figure 1: The color-magnitude diagram for all galaxies near PSZ2 G309.43-72.86.

3.3. Cluster red-sequence galaxies and background galaxies¶

Figure 1 shows that the RS is located approximately at $r-i=0.6$ and $i-z=0.3$ mag, respectively, when the galaxies have different brightness ($r$ or $i$ band magnitude). That concentration represents the RS.

Add lines marking the RS to the CMD using an RS width of 0.2 mag, and select bright RS galaxies (sel_rs, $<20$ mag).

Select background galaxies (sel_bg) as those redder than the RS and sufficiently faint ($>18$ mag). Allow some gap from the the RS because of magnitude measurement uncertainties.

In [15]:

color_dict = {"ri": 0.6, "iz": 0.3}
rs_half_width = 0.1
sel_rs = np.array([True] * len(obj_cat))
sel_bg = np.array([True] * len(obj_cat))
bg_mag_min = 18
rs_mag_max = 20

fig, ax = plt.subplots(1, 2, figsize=(8, 3), sharey=False, sharex=False)
for i in [0, 1]:
    band1 = plt_dict[i][0]
    band2 = plt_dict[i][1]
    ax[i].scatter(obj_cat[band1+"_cModelMag"],
                  obj_cat[band1+"_cModelMag"] - obj_cat[band2+"_cModelMag"],
                  marker=".", s=0.3, alpha=0.5, color='k', label='All galaxies')

    color = color_dict[f"{band1}{band2}"]
    ax[i].axhline(color - rs_half_width, ls="--", color="grey", alpha=1)
    ax[i].axhline(color + rs_half_width, ls="--", color="grey", alpha=1)

    sel_bg &= (obj_cat[f"{band1}_cModelMag"]
               - obj_cat[f"{band2}_cModelMag"]) > (color + rs_half_width)
    sel_rs &= np.abs((obj_cat[f"{band1}_cModelMag"]
                      - obj_cat[f"{band2}_cModelMag"]) - color) < rs_half_width

    ax[i].set_xlabel(band1 + ' [AB mag]')
    ax[i].set_ylabel(band1 + '-' + band2 + ' [AB mag]')
    ax[i].set_xlim([17, 25])
    ax[i].set_ylim([-0.5, 1.5])

for band in band_list:
    sel_bg &= obj_cat[f"{band}_cModelMag"] > bg_mag_min
    sel_rs &= obj_cat[f"{band}_cModelMag"] < rs_mag_max

for i in [0, 1]:
    band1 = plt_dict[i][0]
    band2 = plt_dict[i][1]
    ax[i].scatter(obj_cat[band1 + "_cModelMag"][sel_bg],
                  obj_cat[band1 + "_cModelMag"][sel_bg] - obj_cat[band2 + "_cModelMag"][sel_bg],
                  marker=".", s=0.3, alpha=1, color='blue', label='Background galaxies')
    ax[i].scatter(obj_cat[band1 + "_cModelMag"][sel_rs],
                  obj_cat[band1 + "_cModelMag"][sel_rs] - obj_cat[band2 + "_cModelMag"][sel_rs],
                  marker=".", s=0.3, alpha=1, color='red', label='RS galaxies')
    ax[i].legend()

plt.tight_layout()
plt.show()

Figure 2: This figure shows the CMD as in Figure 1, but with the RS marked with dashed lines and the background galaxies colored blue. Red markers show RS galaxies. Black markers are the same as Figure 1 (all galaxies).

3.4. Spatial distribution of galaxies¶

Plot the 2D distribution of the RS and background galaxies. As expected, the RS galaxies are clustered and the background galaxies have a nearly uniform distribution.

In [16]:

fig, ax = plt.subplots(figsize=(4, 4))
ax.scatter(obj_cat["coord_ra"][sel_rs],
           obj_cat["coord_dec"][sel_rs],
           s=3., color='red')
ax.scatter(obj_cat["coord_ra"][sel_bg],
           obj_cat["coord_dec"][sel_bg],
           s=2., alpha=0.2, color='blue')
ax.invert_xaxis()
ax.set_xlabel("RA [deg]")
ax.set_ylabel("Dec [deg]")
plt.show()

Figure 3: The 2D distribution of RS (red) and background (blue) galaxies. The southern concentration of RS galaxies comes from the ACOS 67 galaxy cluster.

3.5. Mark galaxies on the deep coadd¶

For background (RS) galaxies that overlap the deep_coadd image displayed in Firefly, mark them on the image with a blue (red) circle. Skip galaxies that are too close to the boundary.

First, get the bbox for the boundary.

In [17]:

bbox = deep_coadd.bbox

Mark background galaxies (skip the ones within 100 pixels of the boundary for visualization). The background galaxies are smaller and fainter visually.

In [18]:

sel_inside = obj_cat["x"][sel_bg] < (bbox.stop.x-100)
sel_inside &= obj_cat["y"][sel_bg] < (bbox.stop.y-100)
sel_inside &= obj_cat["x"][sel_bg] > (bbox.start.x+100)
sel_inside &= obj_cat["y"][sel_bg] > (bbox.start.y+100)
with afw_display.Buffering():
    for obj in obj_cat[sel_bg][sel_inside]:
        afw_display.dot("o", obj["x"], obj["y"], size=20, ctype="blue")

Mark RS galaxies (skip the ones within 100 pixels of the boundary for visualization). Those bright large cluster galaxies are the RS galaxies.

In [19]:

sel_inside = obj_cat["x"][sel_rs] < (bbox.stop.x-100)
sel_inside &= obj_cat["y"][sel_rs] < (bbox.stop.y-100)
sel_inside &= obj_cat["x"][sel_rs] > (bbox.start.x+100)
sel_inside &= obj_cat["y"][sel_rs] > (bbox.start.y+100)
with afw_display.Buffering():
    for obj in obj_cat[sel_rs][sel_inside]:
        afw_display.dot("o", obj["x"], obj["y"], size=20, ctype="red")

4. Lensing analysis¶

WL slightly distorts galaxy shapes, which are described by ellipticity. The ellipticity is constructed from the second moments of the object's 2D flux distribution, corrected for PSF effects. The Hirata-Seljak-Mandelbaum (HSM) shape measurement is available in the catalog. It uses a Regaussianization process to correct for the PSF effect that blurs the galaxy images. The details of HSM are described in the work of Hirata and Seljak (2003) and Mandelbaum et al. (2005).

For HSM shapes, the mean ellipticity divided by 2 approximates the lensing shear. Between the measured shear and the true shear there is a small bias, which is caused by galaxy shape dispersion, measurement noise, pixelization, and other effects. Usually, a shear calibration corrects for this bias. However, accurately determining shear calibration parameters requires further image simulation, which is beyond the scope of this notebook. Thus, for demonstration, skip shear calibration but focus on studying the mean galaxy shape in this notebook.

4.1. Select galaxy samples¶

*_hsmShapeRegauss_e1,2 are the ellipiticity components corrected for PSF. The total ellipiticity is $\sqrt{e_1^2+e_2^2}$. Require the measured ellipticity to be within 2. This removes galaxies with very large measured ellipticities that are unphysical, but allows galaxies with measured ellipticities slightly larger than 1 caused by noise (Mandelbaum et al. 2018).

Perform cuts on the background galaxy sample to obtain high-quality objects for lensing analysis (Mandelbaum et al. 2018). Here, *_hsmShapeRegauss_sigma is the catalog estimate of the shape measurement uncertainty due to pixel noise. *_hsmShapeRegauss_flag is the flag set for any failure with the Regaussianziation shapes. *_blendedness measures how much the flux is affected by neighbors. The reference gives more details of the parameters and thresholds.

In [20]:

sel_bg &= np.sqrt(obj_cat[f'{band_shape}_hsmShapeRegauss_e1']**2
                  + obj_cat[f'{band_shape}_hsmShapeRegauss_e2']**2) < 2
sel_bg &= obj_cat[f'{band_shape}_hsmShapeRegauss_sigma'] <= 0.4
sel_bg &= obj_cat[f'{band_shape}_hsmShapeRegauss_flag'] == 0
sel_bg &= obj_cat[f'{band_shape}_blendedness'] < 0.42

Get the background galaxy sample.

In [21]:

obj_cat_bg = obj_cat[sel_bg]

4.2. Ellipticity explanation and calculation¶

Compute the tangential and cross ellipticities of background galaxies with respect to the cluster center. The lensing effect stretchs the galaxy image tangetially to the mass center (i.e., the lensing shear) if the mass is axisymmetric, while there is no effect on the cross direction. Therefore, the mean tangential ellipticity of background galaxies will be generally positive and the mean cross ellipticity will be close to zero.

In the following formulae, $e_\textrm{T}$ is the tangential ellipticity, $e_\textrm{X}$ is the cross ellipticity, $\varphi$ is the position angle towards the cluster center (top North and left East), $e_1$ and $e_2$ are the ellipticity components along the x-axis and along the direction 45 degrees counterclockwise from the x-axis. Note, in the formulae the position angle starts from the negative RA direction (West) and increases counterclockwise, but in astropy the position angle starts from North and increases counterclockwise. Thus, add $\pi/2$ to the position angle computed by astropy.

$e_\textrm{T}= - e_1 \cos(2 \varphi)- e_2 \sin(2 \varphi)$

$e_\textrm{X}= e_1 \sin(2 \varphi) - e_2 \cos(2 \varphi)$

Here, the tangential direction is perpendicular to the line connecting the galaxy and the cluster center, while the cross direction is 45 degrees counterclockwise from the tangential direction.

Set the coordinates.

In [22]:

center_ra = ra_bcg
center_dec = dec_bcg
ra = obj_cat_bg["coord_ra"]
dec = obj_cat_bg["coord_dec"]

coord0 = SkyCoord(center_ra, center_dec, frame='icrs', unit='deg')
coord1 = SkyCoord(ra, dec, frame='icrs', unit='deg')

Compute the position angle.

In [23]:

position_angle = coord0.position_angle(coord1).rad + np.pi/2.

Compute ellipticity components.

In [24]:

e1 = obj_cat_bg[f'{band_shape}_hsmShapeRegauss_e1']
e2 = obj_cat_bg[f'{band_shape}_hsmShapeRegauss_e2']

e_t = - e1 * np.cos(2.*position_angle) - e2 * np.sin(2.*position_angle)
e_x = + e1 * np.sin(2.*position_angle) - e2 * np.cos(2.*position_angle)

Also, compute the radial distance in arcminutes.

In [25]:

r = coord0.separation(coord1).arcmin

4.3. Physical interpretation¶

Examine the distributions of the tangential and cross ellipticities. They both show Gaussian-like distributions (one peak with similar spread and nearly symmetric), but the tangential ellipticity distribution has a mean slightly above zero, while the cross ellipticity distribution has a mean close to zero (Kaiser 1995; Umetsu 2020).

In the code the underscore _ for plt.hist corresponds to bin edges, which is equivalent to the variable bins, and thus skip it here. The _ for plt.legend corresponds to a legend variable and is not used later.

In [26]:

bins = np.linspace(-1.2, 1.2, 21)
mid = 0.5 * (bins[1:] + bins[:-1])

plt.figure()

n_e_t, _, im = plt.hist(e_t, bins=bins, histtype='step',
                        density=True, label='e_t')
n_e_x, _, im = plt.hist(e_x, bins=bins, histtype='step',
                        density=True, ls='--', label='e_x')

plt.axvline(np.mean(e_t), color='C0', label='mean[e_t]')
plt.axvline(np.mean(e_x), color='C1', ls='--', label='mean[e_x]')

plt.xlabel('Ellipticity')
plt.ylabel('Count')

_ = plt.legend()

Figure 4: Distributions of the tangential ellipticities and the cross ellipticities of the selected background galaxies.

The difference between the two distributions shows a noisy dipole-like feature -- this is caused by weak lensing, which slightly distorts the galaxy shapes along the tangential direction and leads to a small offset of the distribution peak, but the spreads of the two distributions are similar. Then, the difference between those two distributions is close to a "dipole" (zero->negative->positive->zero). See also an example in the work of Dell'Antonio et al. (2020).

In [27]:

plt.figure()
plt.step(mid, n_e_t - n_e_x, where='mid')
plt.axvline(0, color='k', ls=':')
plt.xlabel('Ellipticity')
plt.ylabel('Count difference')

Out[27]:

Text(0, 0.5, 'Count difference')

Figure 5: Difference between the distributions of the tangential ellipticities and the cross ellipticities.

Skip the cluster central region where the data is noisy. Bin the data by radial distance to get the mean tangential and cross ellipticities. The result shows mean shapes as a function of radial distance (cluster profiles). Note, these are mean ellipticity profiles rather than shear profiles. The observed ellipticities of those background galaxies (rather than their intrinsic shapes) have been distorted by lensing. Shear calibration is not performed in this notebook for demonstration, as mentioned earlier, and thus lensing shear is not computed in this notebook.

In [28]:

r_max = np.max(r)
print('r_max [arcmin]:', r_max)

bins = np.logspace(np.log10(4), np.log10(r_max), 6)
mid = 0.5 * (bins[1:] + bins[:-1])

e_t_mean, _, _ = binned_statistic(r, e_t, bins=bins)
e_x_mean, _, _ = binned_statistic(r, e_x, bins=bins)

e_t_std, _, _ = binned_statistic(r, e_t, bins=bins,
                                 statistic='std')
e_x_std, _, _ = binned_statistic(r, e_x, bins=bins,
                                 statistic='std')

count, _, _ = binned_statistic(r, e_t, bins=bins,
                               statistic='count')

r_max [arcmin]: 23.99955839135898

Estimate the error bar by the standard deviation divided by the square root of the number of data points ($std/\sqrt{N}$). Here, the shape noise (which can be as large as the ellipticity itself) is usually about 10 times bigger than other sources of error -- for example, the large-scale structure (LSS) can affect the measurement of the cluster lensing signal, because the mass in LSS can cause lensing effects as well. Therefore, using the statistical error is sufficient in this case, assuming each galaxy's shape is an independent measurement.

In [29]:

e_t_err = e_t_std / count ** 0.5
e_x_err = e_x_std / count ** 0.5

When plotting the cross component result, shift the data point slightly (by a factor of 1.05) for better visualization.

In [30]:

plt.figure()
plt.errorbar(mid, e_t_mean, e_t_err, fmt='o', label='e_t')
plt.errorbar(mid*1.05, e_x_mean, e_x_err, fmt='x', label='e_x')
plt.axhline(0, color='k', ls=':')
plt.gca().set_xscale('log')
plt.xlabel('Cluster-centric radial distance [arcmin]')
plt.ylabel('Mean shape <e_i>')
_ = plt.legend()

Figure 6: Mean shape profile with error bars estimated by statistical uncertainty.

Figure 6 shows that this cluster produces a clear lensing signal (positive mean tangential ellipticity) and the signal is stronger (the mean tangential ellipticity is larger) near the cluster center, which is good. Additionally, Figure 6 shows that the mean cross ellipticity is just fluctuating around zero, indicating no significant systematics (as expected). The error bar can be further reduced with deeper observations (more background galaxies) or by combining the data from a large sample of clusters.

5. Exercises for the learner¶

• Repeat the notebook but use r as the selected band (instead of the i band in this tutorial) for the lensing analysis, and test whether the final plot looks like the i band result (they should be generally consistent because lensing effect is independent of color).

Note: It is encouraged that researchers join the LSST Dark Energy Science Collaboration (DESC) to carry out cosmological studies using the Rubin data.