AGN Damped Random Walk Example

In this notebook we run a simulation of a damped random walk AGN variability model built-in with LightCurveLynx. The model is similiar to one used in ELAsTiCC simulations, it uses the standard disk spectrum as the SED “baseline” and perturbs it stochastically.

import numpy as np

from astropy.cosmology import Planck18
from lightcurvelynx.base_models import FunctionNode
from lightcurvelynx.math_nodes.np_random import NumpyRandomFunc
from lightcurvelynx.math_nodes.ra_dec_sampler import ObsTableRADECSampler
from lightcurvelynx.math_nodes.scipy_random import SamplePDF
from lightcurvelynx.models.agn import AGN
from lightcurvelynx.obstable.opsim import OpSim
from lightcurvelynx.simulate import simulate_lightcurves
from lightcurvelynx.survey_info import SurveyInfo
from lightcurvelynx.utils.plotting import plot_lightcurves

Initialize Survey-related Objects

First we will load the survey information, which consists of the observations (obstable), passband information, and noise model. In this example, we use Rubin OpSim information from predownloaded files in the data/ directory in the root directory. Note that nothing in this directory is saved to github, so the files might have to be downloaded initially.

For Rubin, a large number of OpSims can be found at https://s3df.slac.stanford.edu/data/rubin/sim-data/.

Since we are using a standard survey type (LSST), we can just pass the name to the SurveyInfo info constructor and have it load the default passbands and noise model.

# Load the OpSim data.
obstable = OpSim.from_url(
    "https://s3df.slac.stanford.edu/data/rubin/sim-data/sims_featureScheduler_runs5.3/baseline/baseline_v5.3.0_10yrs.db",
)
print(f"Loaded OpSim with {len(obstable)} rows.")

survey_info = SurveyInfo(obstable=obstable, survey_name="LSST")

Loaded OpSim with 1844571 rows.

Create Parameter Distributions

With LightCurveLynx model parameters may be both fixed numerical values or to be drawn from random distributions. Here we create few of AGN model parameters.

# Draw redshift from the uniform distribution
redshift = NumpyRandomFunc("uniform", low=0.1, high=1.0)

# Draw decimal logarithm of the black hole mass (in solar masses)
# from the uniform distribution
lg_bh_mass = NumpyRandomFunc("uniform", low=7.0, high=9.0)
# Transform to mass in solar masses
bh_mass = FunctionNode(
    lambda lg_mass: 10**lg_mass,
    lg_mass=lg_bh_mass,
    node_label="bh_mass",
)


# Eddington ratio for the accretion rate, e.g. accretion rate to critical value ratio
# First, we define a PDF function, for "blue" galaxies, see
# https://github.com/burke86/imbh_forecast/blob/master/var.ipynb
# https://ui.adsabs.harvard.edu/abs/2019ApJ...883..139S/abstract
# https://iopscience.iop.org/article/10.3847/1538-4357/aa803b/pdf
def edd_ratio_pdf(value):
    xi = 10**-1.65
    lambda_br = 10**-1.84
    delta1 = 0.471 - 0.7
    delta2 = 2.53
    min_lambda = 0.01
    max_lambda = 1.0
    value = np.asarray(value)
    fill_mask = (value >= min_lambda) & (value <= max_lambda)
    prob = np.zeros_like(value)
    prob[fill_mask] = xi / (
        (value[fill_mask] / lambda_br) ** delta1 + (value[fill_mask] / lambda_br) ** delta2
    )
    return prob


# We use a special random node to sample from this PDF
edd_ratio = SamplePDF(edd_ratio_pdf)

# We use a special distribution to draw positions from the survey coverage area
radec = ObsTableRADECSampler(
    obstable,
    radius=3.0,  # degrees
    node_label="ra_dec_sampler",
)

Create a Model and Run Simulations

It is time to create the model and simulate few light curves!

model = AGN(
    # Reference time is not important for a stochastic model
    t0=obstable.time_bounds()[0],  # reference date, the survey start
    blackhole_mass=bh_mass,  # black hole mass in solar masses
    edd_ratio=edd_ratio,  # accretion rate to critical accretion rate raio
    redshift=redshift,
    ra=radec.ra,
    dec=radec.dec,
    cosmology=Planck18,
)

# Make it reproducible
rng = np.random.default_rng(42)

df = simulate_lightcurves(
    model=model,
    num_samples=10,
    survey_info=survey_info,
    param_cols=[
        "bh_mass.lg_mass",
        "AGN_0.edd_ratio",
    ],
    rng=rng,
)
# params is too large to show =)
df.drop(columns=["params"])

Simulating: 100%|██████████| 10/10 [00:00<00:00, 10.49obj/s]

	id	ra	dec	nobs	t0	z	bh_mass_lg_mass	AGN_0_edd_ratio	lightcurve
0	0	139.858196	5.627499	668	61208.201382	0.275175	7.378943	0.040020	mjd filter ... obs_idx is_saturated 61218.951231 z ... 4952 False +667 rows ... ... ... ...
1	1	20.117732	-12.482954	712	61208.201382	0.520049	7.259843	0.017156	mjd filter ... obs_idx is_saturated 61240.423002 u ... 16365 False +711 rows ... ... ... ...
2	2	49.830898	-27.321914	706	61208.201382	0.139423	7.951410	0.014922	mjd filter ... obs_idx is_saturated 61301.392936 y ... 36659 False +705 rows ... ... ... ...
3	3	81.318890	-15.613943	711	61208.201382	0.238861	7.453819	0.028611	mjd filter ... obs_idx is_saturated 61327.292817 u ... 47433 False +710 rows ... ... ... ...
4	4	23.753189	-10.960924	630	61208.201382	0.714744	8.339628	0.012180	mjd filter ... obs_idx is_saturated 61240.42139 u ... 16362 False +629 rows ... ... ... ...
5	5	280.082832	-41.721618	925	61208.201382	0.770286	7.874304	0.013223	mjd filter ... obs_idx is_saturated 61209.279026 u ... 896 False +924 rows ... ... ... ...
6	6	197.537710	-33.589037	767	61208.201382	0.970759	8.665356	0.010108	mjd filter ... obs_idx is_saturated 61219.071402 y ... 5214 False +766 rows ... ... ... ...
7	7	295.161710	-67.694582	697	61208.201382	0.393243	8.400530	0.037724	mjd filter ... obs_idx is_saturated 61209.336866 u ... 1005 False +696 rows ... ... ... ...
8	8	11.187587	-73.916434	657	61208.201382	0.433414	7.624733	0.027525	mjd filter ... obs_idx is_saturated 61221.432763 y ... 7530 False +656 rows ... ... ... ...
9	9	162.854353	-1.282532	733	61208.201382	0.522600	8.664520	0.030152	mjd filter ... obs_idx is_saturated 61213.955359 r ... 2993 False +732 rows ... ... ... ...

10 rows x 9 columns

Plot few light curves

for idx in range(5):
    # Extract the row for this object.
    row = df.iloc[idx]
    lc = row["lightcurve"]

    plot_lightcurves(
        fluxes=lc["flux"],
        times=lc["mjd"],
        fluxerrs=lc["fluxerr"],
        filters=lc["filter"],
    )