305.2. Stellar variability characterization

305.2. Stellar variability characterization¶

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

Learning objective: Learn a way of identifying candidate variable stars in DP2.

LSST data products: DiaObject, ForcedSourceOnDiaObject, CcdVisit

Packages: matplotlib, numpy, Astropy, astroquery, lsst.rsp, lsst.utils.plotting

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¶

This notebook demonstrates how to use the lightcurve statistics in the DiaObject table to identify candidate variable stars. The candidates are then run through a Lomb-Scargle periodogram code to identify likely periodicity, after which phased lightcurves are plotted.

Related tutorials: There are 200-level tutorials on difference_images as well as the DiaSource, DiaObject and ForcedSourceOnDiaObject catalogs. There is another notebook in this 305 series that presents how to make a light curve of a known variable star. The 207-level tutorials explain timeseries quantities that are available in the DiaObject table.

1.1. Import packages¶

Import numpy, a fundamental package for scientific computing with arrays in Python (numpy.org), and matplotlib, a comprehensive library for data visualization (matplotlib.org; matplotlib gallery). From Astropy (astropy.org) import the units and SkyCoord modules and the Lomb-Scargle periodogram tools.

From the lsst package, import the RSPDiscovery module for accessing the Table Access Protocol (TAP) service, and some plotting utilities.

In [1]:

import matplotlib.pyplot as plt
import numpy as np
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.table import vstack
from astropy.timeseries import LombScargle, LombScargleMultiband

from lsst.rsp import RSPDiscovery
from lsst.utils.plotting import (get_multiband_plot_colors,
                                 get_multiband_plot_symbols)

1.2. Define parameters and functions¶

Create an instance of the TAP service.

In [2]:

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

Define filter names, plot markers, linestyles, and colors for plotting

In [3]:

filter_colors = get_multiband_plot_colors()
filter_names = filter_colors.keys()
filter_symbols = get_multiband_plot_symbols()

2. Find candidate variable stars¶

The DiaObject table contains many lightcurve summary statistics, which are derived from the contents of the DiaSource table. See the schema for DiaObject for a list of columns. Details of these statistical quantities summarizing the properties of each DiaObject's lightcurve are given in the 207-series tutorials on "Timeseries values" and "Timeseries distributions."

The goal of this tutorial is to identify variable stars in our own Milky Way galaxy. Conduct a search in the relatively low Galactic latitude Rubin_SV-225_-40 field from the DP2 dataset. Not only will this field have relatively high stellar density, but it is also well sampled, with a total of ~1000 visits in DP2.

2.1. Get a sample of DiaObjects¶

To identify candidates such as the periodic variable star from the first notebook in the 305-series, one should select stars that are varying in brightness significantly. Do so by using [f]_psfFluxSigma, which measures the standard deviation of the PSF fluxes from the difference images.

For this tutorial, focus on the r-band measurements. Note that the quantities selected below are all based on the difference-image point source flux values, except for r_scienceFlux*, which are based on the calibrated visit_image.

Query the DiaObject table in a 1-degree radius about a position in the low-latitude field, extracting stars with direct-image signal-to-noise (S/N) greater than 5, and with at least 100 DiaSource measurements:

In [4]:

ra = 226.7
dec = -40.3
search_radius = 1.0

In [5]:

query = "SELECT ra, dec, diaObjectId, "\
        "r_psfFluxMean, nDiaSources, r_scienceFluxMean, r_scienceFluxMeanErr, "\
        "r_psfFluxNdata, r_psfFluxSigma, "\
        "scisql_nanojanskyToAbMag(r_scienceFluxMean) as rmag, "\
        "scisql_nanojanskyToAbMag(g_scienceFluxMean) as gmag "\
        "FROM dp2.DiaObject "\
        "WHERE CONTAINS (POINT('ICRS', ra, dec), "\
        "CIRCLE('ICRS', " + str(ra) + ", " + str(dec) + ", " + str(search_radius) + ")) = 1 "\
        "AND r_scienceFluxMean/r_scienceFluxMeanErr > 5 "\
        "AND nDiaSources > 100"
print(query)

SELECT ra, dec, diaObjectId, r_psfFluxMean, nDiaSources, r_scienceFluxMean, r_scienceFluxMeanErr, r_psfFluxNdata, r_psfFluxSigma, scisql_nanojanskyToAbMag(r_scienceFluxMean) as rmag, scisql_nanojanskyToAbMag(g_scienceFluxMean) as gmag FROM dp2.DiaObject WHERE CONTAINS (POINT('ICRS', ra, dec), CIRCLE('ICRS', 226.7, -40.3, 1.0)) = 1 AND r_scienceFluxMean/r_scienceFluxMeanErr > 5 AND nDiaSources > 100

Run the TAP search and fetch the results in table form.

In [6]:

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'

DiaObjsFull = job.fetch_result().to_table()

Job phase is COMPLETED

In [7]:

len(DiaObjsFull)

Out[7]:

37799

2.2. Lightcurve statistics for DiaObjects¶

Plot the distribution of r_psfFluxSigma/r_scienceFluxMean (i.e., scale the sigma of the flux by the mean value from the direct images) as a function of r magnitude.

In [8]:

fig, ax = plt.subplots(1, 1, figsize=(6, 4))

plt.plot(DiaObjsFull['rmag'],
         DiaObjsFull['r_psfFluxSigma']/DiaObjsFull['r_scienceFluxMean'],
         '.', ms=2, alpha=0.5, color=filter_colors['r'])
plt.xlabel('r magnitude')
plt.ylabel('r_psfFluxSigma/mean')
plt.xlim(16.8, 21.3)
plt.ylim(-0.02, 0.45)
plt.minorticks_on()
plt.show()

Figure 1: Lightcurve statistics of all candidates from the TAP search, with r_psfFluxSigma scaled by dividing by the mean science image flux for each object.

There is a well-defined locus of points in the above plot, where most objects have r_psfFluxSigma (normalized by the mean flux) near zero, implying that they are non-variable.

2.3. Filter for variable star candidates¶

Filter the DiaObject results (defined as DiaObjsFull) with the following conditions:

• nDiaSources > 300: Candidate has more than 300 DiaSource detections (total over all bands); applied in TAP query.
• r_scienceFluxMean S/N > 5: Candidate has signal-to-noise larger than 5 in the visit_images; applied in TAP query.
• r_psfFluxSigma/r_scienceFluxMean > 0.05: Candidate has flux variations with sigma larger than 5 percent.
• 18 < r < 19.5: Candidate has mean r magnitude between 18 to 19.5 (i.e., faint enough to avoid saturation in any images as the star varies in brightness, while also avoiding faint sources with large measurement uncertainties).
• -0.2 < (g-r) < 0.4: Candidate has a g-r color within the instability strip (i.e., RR Lyrae or delta Scuti-type variables).

Apply these cuts to the DiaObjsFull table, extracting a list of candidates as pick_all.

In [9]:

sigcut = 0.05
ndiasourcescut = 300

magrange = (DiaObjsFull['rmag'] > 18.0) & (DiaObjsFull['rmag'] < 19.5)
pick_color = ((DiaObjsFull['gmag']-DiaObjsFull['rmag']) < 0.4) &\
             ((DiaObjsFull['gmag']-DiaObjsFull['rmag']) > -0.2)
pick_sig = (DiaObjsFull['r_psfFluxSigma']/DiaObjsFull['r_scienceFluxMean'] > sigcut)
pick_ndia = (DiaObjsFull['nDiaSources'] > ndiasourcescut)

pick_all = np.where(magrange & pick_sig & pick_color & pick_ndia)
len(pick_all[0])

Out[9]:

4

Plot the distribution of psfFluxSigma again, but this time marking the cuts applied and the pick_all sample of candidate variables. Plot them on a log scale to more clearly show differences.

In [10]:

fig, ax = plt.subplots(1, 1, figsize=(6, 4))

plt.hlines(sigcut, 16.8, 21.3, color='gray', linestyle='--')
plt.vlines([18.0, 19.5], 5e-2, 8e-1, color='gray', linestyle=':')
plt.plot(DiaObjsFull['rmag'],
         DiaObjsFull['r_psfFluxSigma']/DiaObjsFull['r_scienceFluxMean'],
         '.', ms=2, alpha=0.5, color=filter_colors['r'])
plt.plot(DiaObjsFull['rmag'][pick_all],
         DiaObjsFull['r_psfFluxSigma'][pick_all]/DiaObjsFull['r_scienceFluxMean'][pick_all],
         'k*', ms=5)
plt.xlabel('r magnitude')
plt.ylabel('r_psfFluxSigma/mean')
plt.yscale('log')
plt.xlim(16.8, 21.3)
plt.ylim(8e-4, 8e-1)
plt.minorticks_on()
plt.show()

Figure 2: The same quantities as in Figure 1, but with the selected candidate variables highlighted as black stars, demonstrating that the candidates are within the thresholds used for selection.

The candidates are all within the thresholds used for selection.

3. Plot lightcurves of variable star candidates¶

3.1. Plot the raw, as-measured lightcurves¶

Query the ForcedSourceOnDiaObject table for all measurements of each object, based on their diaObjectIds. Then loop over the selected objects, plotting the lightcurve for each.

Option to view the results table:

In [11]:

# DiaObjsFull[pick_all]

In [12]:

id_list = list(DiaObjsFull[pick_all]['diaObjectId'])
id_string = "(" + ", ".join(str(value) for value in id_list) + ")"

results = service.search("SELECT fsodo.diaObjectId, "
                         "fsodo.visit, fsodo.band, "
                         "fsodo.psfFlux, fsodo.psfFluxErr, "
                         "fsodo.psfDiffFlux, fsodo.psfDiffFluxErr, "
                         "fsodo.psfDiffFlux_flag, fsodo.diff_PixelFlags_nodataCenter, "
                         "fsodo.pixelFlags_saturatedCenter, fsodo.invalidPsfFlag, "
                         "fsodo.psfFlux_flag, fsodo.pixelFlags_bad, "
                         "vis.expMidptMJD "
                         "FROM dp2.ForcedSourceOnDiaObject as fsodo "
                         "JOIN dp2.Visit as vis "
                         "ON (vis.visit = fsodo.visit)"
                         "WHERE pixelFlags_bad = 0 "
                         "AND psfDiffFlux_flag = 0 AND diff_PixelFlags_nodataCenter = 0 "
                         "AND pixelFlags_saturatedCenter = 0 AND invalidPsfFlag = 0 "
                         "AND psfFlux_flag = 0 AND pixelFlags_bad = 0 "
                         "AND diaObjectId IN "+id_string)

results = results.to_table()

In [13]:

filters = ['g', 'r', 'i']

for id in id_list:
    print(f"diaObjectId: {id}")

    print(f"number of measurements: {len(results[(results['diaObjectId'] == id)])}")

    fig = plt.figure(figsize=(6, 4))

    for f, filt in enumerate(filters):
        fx = np.where((results['band'] == filt) & (results['diaObjectId'] == id))[0]
        plt.plot(results['expMidptMJD'][fx], results['psfDiffFlux'][fx],
                 filter_symbols[filt], ms=10, mew=0, alpha=0.5,
                 color=filter_colors[filt])

    plt.title(str(id))
    plt.xlabel('MJD (days)')
    plt.ylabel('psfDiffFlux (nJy)')

    plt.tight_layout()
    plt.show()

diaObjectId: 743572294793691373
number of measurements: 537

diaObjectId: 744863671200514193
number of measurements: 1075

diaObjectId: 743582465276248313
number of measurements: 873

diaObjectId: 743580953447760076
number of measurements: 584

Figure 3: Lightcurves of all selected candidates. Each panel shows the measured difference-image flux (psfDiffFlux) in nanoJanskies vs. observation time in MJD days.

3.2. Plot phased lightcurves¶

The lightcurves above show that these objects clearly vary in brightness, but it is unclear whether there is a pattern (for example a periodic signal) in that variability.

Loop over the candidates, extracting the frequency of the highest peak in the Lomb-Scargle results, and placing each measurement into the correct phase of the star's variability period.

In [14]:

periods = {}

for id in id_list:
    print(f"diaObjectId: {id}")

    res = results[(results['diaObjectId'] == id)]
    print(f"number of measurements: {len(res)}")

    minfreq = 1 / (0.95*u.d)
    maxfreq = 1 / (0.05*u.d)

    obj_mjd_days = np.array(res['expMidptMJD']) * u.day
    obj_fluxes = np.array(res['psfFlux'])
    obj_flux_errs = np.array(res['psfFluxErr'])
    ls_model = LombScargleMultiband(
        obj_mjd_days,
        obj_fluxes,
        res['band'],
        obj_flux_errs
    )
    obj_frequency, obj_power = ls_model.autopower(
        minimum_frequency=minfreq,
        maximum_frequency=maxfreq,
    )
    max_power = np.argmax(obj_power)
    obj_freq = obj_frequency[max_power]
    obj_period = 1.0 / obj_freq
    print(f"period: {obj_period: .8F}")
    periods[id] = obj_period

    t0 = 0.0
    obj_mjd_norm = (obj_mjd_days.value - t0) / obj_period.value
    obj_phase = np.mod(obj_mjd_norm, 1.0)

    mag = (obj_fluxes*u.nJy).to(u.ABmag).value

    fig = plt.figure(figsize=(6, 4))

    for f, filt in enumerate(filter_names):
        fx = np.where(res['band'] == filt)[0]

        plt.plot(obj_phase[fx], mag[fx],
                 filter_symbols[filt], color=filter_colors[filt],
                 ms=8, mew=0, alpha=0.5, label=filt)

    plt.ylim(np.nanpercentile(mag, 98) + 0.4, np.nanpercentile(mag, 2) - 0.4)
    plt.xlabel('phase')
    plt.ylabel('mag')
    plt.title(str(id)+", phased")
    plt.legend()
    plt.minorticks_on()
    plt.tight_layout()
    plt.show()

    print('\n')

diaObjectId: 743572294793691373
number of measurements: 537

period:  0.26624586 d

diaObjectId: 744863671200514193
number of measurements: 1075

period:  0.63283593 d

diaObjectId: 743582465276248313
number of measurements: 873

period:  0.34091785 d

diaObjectId: 743580953447760076
number of measurements: 584

period:  0.15581064 d

Figure 4: As in Figure 3, the figure shows lightcurves of all selected candidates. In this plot, the vertical axis shows PSF magnitudes from forced measurements on the direct images. The horizontal axis shows phase, which runs from 0 to 1, placing each observed point into the proper phase of the variable's estimated period.

3.3. Refine the period search¶

It is clear that the stars above exhibit some periodic variable behavior, but the phased lightcurves don't look as "clean" as one would like. This is likely because the period estimates are not quite correct. Re-run the Lomb-Scargle algorithm for each star, changing the following two things:

1. Zoom in on a narrower range of periods to search, using the period from our initial estimates to define the new range, and
2. Increase the "samples_per_peak" in the Lomb-Scargle search to more finely sample the distribution of frequencies.

Wrap the steps in a function that can be called for each star.

In [15]:

def period_search(results, my_id, period0, pick_band='r'):
    min_period = period0 - 0.2
    max_period = period0 + 0.2

    if min_period < 0:
        min_period = 0.05

    min_period = min_period * u.day
    max_period = max_period * u.day

    min_freq = 1.0 / max_period
    max_freq = 1.0 / min_period

    res = results[(results['band'] == pick_band) & (results['diaObjectId'] == my_id)]

    obj_mjd_days_band = np.array(res['expMidptMJD']) * u.day
    obj_fluxes_band = np.array(res['psfFlux'])
    obj_frequency, obj_power =\
        LombScargle(obj_mjd_days_band, obj_fluxes_band).autopower(minimum_frequency=min_freq,
                                                                  maximum_frequency=max_freq,
                                                                  samples_per_peak=2000)

    peakbin = np.argmax(obj_power)
    mean_peak_freq = obj_frequency[peakbin].value

    fig, ax = plt.subplots(1, 2, figsize=(6, 3))

    plt.sca(ax[0])
    plt.plot(obj_frequency, obj_power)
    plt.vlines(mean_peak_freq, 0, 0.3, linestyle='--', color='red')
    plt.minorticks_on()
    plt.xlabel('frequency (1/d)')
    plt.ylabel('power')

    plt.sca(ax[1])
    plt.plot(1 / obj_frequency, obj_power)
    plt.vlines(1/mean_peak_freq, 0, 0.3, linestyle='--', color='red')
    plt.minorticks_on()
    plt.xlabel('period (d)')
    plt.ylabel('power')

    fig.suptitle(my_id)
    plt.tight_layout()
    plt.show()

    max_power = np.argmax(obj_power)
    obj_freq = obj_frequency[max_power]
    obj_period = 1.0 / obj_freq

    t0 = 0.0
    obj_mjd_days = np.array(results['expMidptMJD']) * u.day
    obj_fluxes = np.array(results['psfFlux'])

    obj_mjd_norm = (obj_mjd_days.value - t0) / obj_period.value
    obj_phase = np.mod(obj_mjd_norm, 1.0)

    mag = (obj_fluxes*u.nJy).to(u.ABmag).value

    fig = plt.figure(figsize=(6.5, 3))

    for band in filter_names:
        results_band = (results['band'] == band) & (results['diaObjectId'] == my_id)

        plt.plot(obj_phase[results_band], mag[results_band],
                 filter_symbols[band], color=filter_colors[band],
                 ms=8, mew=0, alpha=0.5, label=band)

    plt.gca().invert_yaxis()
    plt.legend(ncol=4, loc="lower right")
    plt.title(f"Period: {obj_period: .8F}")
    plt.xlabel('phase')
    plt.ylabel('magnitude')
    plt.ylim(np.nanpercentile(mag, 98)+0.3, np.nanpercentile(mag, 2)-0.3)
    plt.minorticks_on()
    plt.show()

Execute the period refinement function to improve the precision of each candidate and review the resulting phased light curves.

In [16]:

for my_id in id_list:
    period_search(results, my_id, periods[my_id].value, pick_band='r')

Figure 5: As in Figure 4, the lower panel of each figure shows a lightcurve for each selected candidate, in magnitude vs. phase. The upper plot for each candidate shows Lomb-Scargle "power" vs. frequency on the left, and power vs. period on the right. The selected "best" period is marked with a vertical red line in each panel. The lightcurves are phased with this period, and the period is labeled at the top of the lightcurve.

3.4. Check objects in Simbad¶

Query Simbad using the ra, dec position of each object to see if they have already been classified as variable objects. Start by compiling the results in a dict, then combine them all into an Astropy table.

In [17]:

simbad = Simbad()
simbad.add_votable_fields("otype", "allfluxes")

simbad_dict = {}

for obj in DiaObjsFull[pick_all]:
    simbad_dict[obj['diaObjectId']] = simbad.query_region(SkyCoord(obj['ra'], obj['dec'],
                                                          unit=(u.deg, u.deg),
                                                          frame='icrs'), radius=2 * u.arcsec)

In [18]:

for i, key in enumerate(simbad_dict.keys()):
    if i == 0:
        tab = simbad_dict[key]
    else:
        tab = vstack(list(simbad_dict.values()), metadata_conflicts='silent')

Print out the stacked table.

In [19]:

tab

Out[19]:

Table length=4

main_id	ra	dec	coo_err_maj	coo_err_min	coo_err_angle	coo_wavelength	coo_bibcode	otype	B	F150W	F200W	F444W	g	G	H	i	I	J	K	r	R	u	U	V	z
	deg	deg	mas	mas	deg
object	float64	float64	float32	float32	int16	str1	object	object	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64	float64
Gaia DR3 6004904483618003200	226.5106485661879	-41.01792212063167	0.1095	0.1045	90	O	2020yCat.1350....0G	EB*	--	--	--	--	--	17.95594024658203	--	--	--	--	--	--	--	--	--	--	--
CRTS J150748.5-394118	226.95255925441373	-39.688372268556655	0.1613	0.1253	90	O	2020yCat.1350....0G	RR*	--	--	--	--	--	18.49850082397461	--	--	--	--	--	--	--	--	--	--	--
Gaia DR3 6005347144427786624	226.72443472038665	-40.30404863746611	0.1364	0.1078	90	O	2020yCat.1350....0G	RR*	--	--	--	--	--	18.265031814575195	--	--	--	--	--	--	--	--	--	--	--
Gaia DR3 6005276599582740096	227.0514067069929	-40.73327023308945	0.2008	0.1637	90	O	2020yCat.1350....0G	EB*	--	--	--	--	--	18.944852828979492	--	--	--	--	--	--	--	--	--	--	--

The "otype" column shows that these candidate variables are all known eclipsing binaries ("EB*") or RR Lyrae ("RR*"). Some refinement may be necessary to get their periods correct, but they are all obviously periodic according to their phased lightcurves.

In [ ]: