304.2. Galaxy cluster red sequence galaxies#

304_2_cluster_red_sequence

304.2. Galaxy cluster red sequence galaxies¶

For the Rubin Science Platform at data.lsst.cloud.
Data Release: Data Preview 2
Container Size: Large
LSST Science Pipelines version: 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: How to identify the red sequence of a galaxy cluster.

LSST data products: deep_coadd images, Object table

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

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).

Find red sequence galaxies of PSZ2 G309.43-72.86 using color-magnitude diagrams (CMDs) and color-color plots (CCs). Visualize those galaxies in the coadded images.

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 packages for bulter, TAP, and display.

In [1]:

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
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 PSZ2 G309.43-72.86 Brightest Central Galaxy (BCG) is known at RA=10.2082deg, DEC=-44.1307deg.

Provide the center coordinates and the range in degree for making a cone search for the catalog.

In [5]:

ra_bcg = 10.2082
dec_bcg = -44.1307

range_deg = 0.1

Make a cone search for the catalog. This takes a few seconds.

In [6]:

query = "SELECT objectId, coord_ra, coord_dec, refExtendedness, " \
        "deblend_peak_center_x, deblend_peak_center_y, " \
        "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), " \
        "CIRCLE('ICRS', %f, %f, %f)) = 1 "%(ra_bcg, dec_bcg, range_deg)

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: ', len(obj_cat))

Job phase is COMPLETED

Number of objects:  10777

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

In [7]:

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

Use the butler to get the deep coadd image in r band for visualization.

In [8]:

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

deep_coadd = butler.get(dataset_refs[0])

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

3. Plot data¶

3.1. Select extended objects¶

Select the objects that do not have flags and negative fluxes to ensure the quality, and use extendedness to get galaxies. 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 CModel photometry (Bosch et al. 2018) for galaxies. Cut at SNR>3 for visualization and data quality.

Consider three bands r, i, z, and consider CMDs in colors r-i and i-z, where RS is clear to see. Finding RS in multiple colors instead of a single color improves the quality of RS galaxy selection (Fu et al. 2024).

In [9]:

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

sel = obj_cat["refExtendedness"] == 1

for band in band_list:
    sel &= (obj_cat[f"{band}_cModel_flag"] == False)
    sel &= obj_cat[f'{band}_cModel_SNR'] > 3.

obj_cat = obj_cat[sel]

3.2. Plot magnitude histograms¶

Visualize the histogram of magnitudes in each band and find the peak. Also, plot the magnitude vs. signal-to-noise-ratio (SNR) for comparision. At the magnitude histogram peak, the SNR is about 5 to 20.

In [10]:

for band in band_list:

    fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")

    mag_lim_min, mag_lim_max = 15, 28
    bins = np.linspace(mag_lim_min, mag_lim_max, 41)
    mid = 0.5 * (bins[1:] + bins[:-1])

    count, _, _ = ax.hist(obj_cat[f"{band}_cModelMag"], bins=bins,
                          histtype="step")
    mid_max = mid[np.argmax(count)]
    print(band, "peak mag: ", mid_max)
    ax.axvline(mid_max, c="r", ls="--", label="peak")

    ax.set_yscale("log")
    ax.set_xlabel("Magnitude")
    ax.set_ylabel("Count")
    ax.set_xlim([mag_lim_min, mag_lim_max])

    ax2 = ax.twinx()
    ax2.scatter(obj_cat[f"{band}_cModelMag"], obj_cat[f"{band}_cModel_SNR"],
                s=1, c="k")
    ax2.axhline(5, c="g", ls=":", label="SNR=5")

    ax2.set_yscale("log")
    ax2.set_ylabel("SNR")

    plt.title(f"{band} band")

    ax.legend(loc="lower left")
    ax2.legend(loc="upper left")

r peak mag:  24.9125
i peak mag:  24.5875
z peak mag:  24.5875

No description has been provided for this image

Figure 1: Distribution of magnitudes, and relation with the signal-to-noise-ratio (SNR).

3.3. Plot color-magnitude diagrams¶

It can be seen that the RS locates at r-i=0.6 and i-z=0.3, respecitively. Plot lines for RS, and select bright galaxies (< 23 mag) near RS.

In [11]:

color_dict = {
    "ri": 0.6,
    "iz": 0.3,
}

eps = 0.1
sel_rs = np.array([True] * len(obj_cat))
mag_max = 23

for ind in range(len(band_list)-1):

    band1 = band_list[ind]
    band2 = band_list[ind+1]

    fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")

    ax.scatter(obj_cat[f"{band1}_cModelMag"],
               obj_cat[f"{band1}_cModelMag"] - obj_cat[f"{band2}_cModelMag"],
               marker=".", s=0.3)

    ax.set_xlim([18, 25])
    ax.set_ylim([-2, 2])

    ax.set_xlabel(band1)
    ax.set_ylabel(f"{band1}-{band2}")

    color = color_dict[f"{band1}{band2}"]
    ax.axhline(color - eps, ls="--", c="k", alpha=0.3)
    ax.axhline(color + eps, ls="--", c="k", alpha=0.3)

    sel_rs &= np.abs(obj_cat[f"{band1}_cModelMag"] - obj_cat[f"{band2}_cModelMag"] - color) < eps
    sel_rs &= obj_cat[f"{band1}_cModelMag"] < mag_max

Figure 2: CMDs and RS in r-i and i-z.

3.4. Plot color-color diagrams¶

Next, plot galaxies in the color-color space, and highlight the RS galaxies.

In [12]:

ind = 0
band1 = band_list[ind]
band2 = band_list[ind+1]
band3 = band_list[ind+2]

fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")

ax.scatter(obj_cat[f"{band1}_cModelMag"] - obj_cat[f"{band2}_cModelMag"],
           obj_cat[f"{band2}_cModelMag"] - obj_cat[f"{band3}_cModelMag"],
           marker=".", s=0.5)

ax.scatter(obj_cat[f"{band1}_cModelMag"][sel_rs] - obj_cat[f"{band2}_cModelMag"][sel_rs],
           obj_cat[f"{band2}_cModelMag"][sel_rs] - obj_cat[f"{band3}_cModelMag"][sel_rs],
           marker=".", s=0.3, label="RS")

ax.axvline(color_dict[f"{band1}{band2}"], ls="--", c="k", alpha=0.3)
ax.axhline(color_dict[f"{band2}{band3}"], ls="--", c="k", alpha=0.3)

ax.legend()

ax.set_xlim([-1, 2])
ax.set_ylim([-1, 2])

ax.set_xlabel(f"{band1}-{band2}")
ax.set_ylabel(f"{band2}-{band3}")

Out[12]:

Text(0, 0.5, 'i-z')

Figure 3: RS in r-i vs i-z color-color space.

4. Red sequence galaxies¶

4.1. Spatial distribution¶

Plot the 2D distribution of those RS galaxies and use color to show the magnitude. The galaxies are concentrated at the cluster center, and there are some clumps in the distribution.

In [13]:

fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")
im = ax.scatter(obj_cat["coord_ra"][sel_rs],
                obj_cat["coord_dec"][sel_rs],
                s=2.,
                c=obj_cat["r_cModelMag"][sel_rs],
                cmap="viridis",
                )
ax.plot([ra_bcg], [dec_bcg], 'r*')
ax.invert_xaxis()
ax.set_xlabel("RA [deg]")
ax.set_ylabel("DEC [deg]")
cbar = fig.colorbar(im, ax=ax)
cbar.set_label("r [mag]")

Figure 4: 2D RS galaxies distribution.

4.2. Radial density profile¶

Plot the number density (number per unit area) of RS galaxies as a function of radial distance. The density drops when the radial distance gets larger.

In [14]:

coord_obj = SkyCoord(obj_cat["coord_ra"],
                     obj_cat["coord_dec"])
coord_bcg = SkyCoord(ra_bcg * u.deg, dec_bcg * u.deg)
sep = coord_bcg.separation(coord_obj)
radial_distance = sep.arcmin[sel_rs]

bin_num = 6
bins = np.linspace(0, range_deg*60, bin_num)
mid = 0.5 * (bins[1:] + bins[:-1])
area = (bins[1:]**2 - bins[:-1]**2) * np.pi

count, _ = np.histogram(radial_distance, bins=bins)

fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")
ax.scatter(mid, count/area)
ax.set_xlabel("Radial distance [arcmin]")
ax.set_ylabel("Density [arcmin^-2]")

Out[14]:

Text(0, 0.5, 'Density [arcmin^-2]')

Figure 5: Radial number density profile of RS galaxies.

4.3. Magnitude histogram¶

Plot the magnitude histogram of RS galaxies in each band. This is similar to the luminosity function, because galaxy member galaxies are almost at the same distance to the observer.

In [15]:

fig, ax = plt.subplots(figsize=(4, 3), layout="constrained")

bins = np.linspace(15, mag_max, 21)

for band in band_list:
    ax.hist(obj_cat[f"{band}_cModelMag"][sel_rs], bins=bins,
            histtype="step", label=band)

ax.set_yscale("log")
ax.set_xlabel("Magnitude")
ax.set_ylabel("Count")
ax.legend()

Out[15]:

<matplotlib.legend.Legend at 0x7eeea9db3610>

Figure 6: Magnitude histogram of RS galaxies.

4.4. Objects on the image¶

Plot galaxies on the coadd image.

First, display the image in the Firefly tab and set the mask plane to be transparent.

In [16]:

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

Overplot objects. Leave some blank region near the edge.

In [17]:

bbox = deep_coadd.bbox

sel_inside = obj_cat["x"][sel_rs] < (bbox.stop.x-100)
sel_inside &= obj_cat["y"][sel_rs] < (bbox.stop.y-100)

for i in obj_cat[sel_rs][sel_inside]:
    afw_display.dot("o", i["x"], i["y"], size=20, ctype="red")

The deep coadd image which is displayed in the Firefly tab does not cover the entire cluster, as PSZ2 G309.43-72.86 is very large.

Figure 7: A screenshot of the Firefly display window as it should appear for the deep_coadd image that overlaps the center of PSZ2 G309.43-72.86, with red sequence members marked with red circles.