Light Curves for an SPB Star

This Cookbook recipe demonstrates how to synthesize light curves for a slowly pulsating B-type (SPB) star, using the MSG Python interace. The light curves are based on surface perturbation data for a \(5\,\Msun\) model, produced by the GYRE stellar oscillation code.

Preparation

Download spb-light-curve.tar.bz2 and unpack it in your working directory. This archive contains the following items:

• summary.h5 — surface perturbation data produced by GYRE

• pygyre — Python module for reading summary.h5

• pg-*.h5 — photgrid files for each passband [1]

• gyre.in — GYRE input file used to create summary.h5

• spb.mesa — stellar model file used to create summary.h5 [2]

Note that gyre.in and spb.mesa aren’t needed for the recipe; they’re included for those who wish to recreate summary.h5, or are curious about its provenance.

Initialization

To start, import modules and set other parameters:

# Import modules

import numpy as np
import matplotlib.pyplot as plt
import scipy.special as ss

import astropy.constants as con

import pymsg
import pygyre

# Set plot parameters

%matplotlib inline
plt.rcParams.update({'font.size': 16})

Define a Synthesis Function

Next, define a function that synthesizes a light curve given the following arguments:

• phi — oscillation phase array (radians)

• i — observer inclination angle (radians)

• photgrid_name — name of photgrid file

This function follows the semi-analytical approach described in Section 2 of Townsend (2003); the comments in the function definition highlight the specific equations from that paper.

def synthesize_light_curve(phi, i, photgrid_name):

    # Create the PhotGrid object

    photgrid = pymsg.PhotGrid(photgrid_name)

    # Read the GYRE output file into a table

    tbl = pygyre.read_output('summary.h5')

    # Select the table row corresponding to the l=m=2, g_30 mode

    mask = (tbl['l'] == 2) & (tbl['m'] == 2) & (tbl['n_pg'] == -30)

    row = tbl[mask][0]

    # Extract stellar parameters

    M_star = row['M_star']
    R_star = row['R_star']
    L_star = row['L_star']

    Teff = (L_star/(4*np.pi*con.sigma_sb.cgs.value*R_star**2))**0.25
    g = con.G.cgs.value*M_star/R_star**2

    # Extract surface perturbation parameters

    l = row['l']
    m = row['m']

    omega = row['omega']

    dR_R = row['xi_r_ref']      # Delta_R (T03, eqn. 4)
    dL_L = row['lag_L_ref']
    dT_T = 0.25*(dL_L - 2*dR_R) # Delta_T (T03, eqn. 5)
    dg_g = -(2 + omega**2)*dR_R # Delta_g (T03, eqn. 6)

    # Evaluate intensity moments & their partial derivatives (T03, eqn. 15)

    x = {'Teff': Teff, 'log(g)': np.log10(g)}

    I_0 = photgrid.D_moment(x, 0)
    I_x = photgrid.D_moment(x, l)

    dI_x_dlnT = photgrid.D_moment(x, l, deriv={'Teff': True}) * Teff
    dI_x_dlng = photgrid.D_moment(x, l, deriv={'log(g)': True}) * np.log(np.exp(1.))

    # Set up differential flux functions (T03, eqns. 12-14)

    R_lm = (2+l)*(1-l) * I_x/I_0 * ss.sph_harm(m, l, 0., i)
    T_lm = dI_x_dlnT/I_0 * ss.sph_harm(m, l, 0., i)
    G_lm = dI_x_dlng/I_0 * ss.sph_harm(m, l, 0., i)

    # Construct the light curve (T03, eqn. 11)

    dF_F = ((dR_R*R_lm + dT_T*T_lm + dg_g*G_lm) * np.exp(1j*phi)).real

    # Return it

    return dF_F

Evaluate Light Curves

Using this function, evaluate light curves in each passband:

# Set up phase array

phi = np.linspace(0, 4*np.pi, 1001)

# Set observer inclination

i = np.radians(75)

# Evaluate light curves

dF_F_Kepler = synthesize_light_curve(phi, i, 'pg-Kepler.h5')
dF_F_TESS = synthesize_light_curve(phi, i, 'pg-TESS.h5')
dF_F_Gaia_B = synthesize_light_curve(phi, i, 'pg-Gaia-B.h5')
dF_F_Gaia_R = synthesize_light_curve(phi, i, 'pg-Gaia-R.h5')
dF_F_JWST = synthesize_light_curve(phi, i, 'pg-JWST.h5')

Note that the pg-JWST.h5 file is based on the JWST NIRCam F150W filter — chosen to demonstrate what the light curve will look like in the near infrared, but with the recognition that it’s very unlikely JWST will ever observe an SPB star!

Plot the Light Curves

As a final step, plot the light curve in each filter:

# Plot light curves

plt.figure()

plt.plot(phi/(2*np.pi), dF_F_Gaia_B, label='Gaia B', color='blue')
plt.plot(phi/(2*np.pi), dF_F_Kepler, label='Kepler', color='green')
plt.plot(phi/(2*np.pi), dF_F_TESS, label='TESS', color='yellow')
plt.plot(phi/(2*np.pi), dF_F_Gaia_R, label='Gaia R', color='orange')
plt.plot(phi/(2*np.pi), dF_F_JWST, label='JWST', color='red')

plt.xlabel('Phase (cycles)')
plt.ylabel('Delta F/F')
plt.legend()

<matplotlib.legend.Legend at 0x7f8a86fb86a0>

A clear trend can be seen in the light curves: from red wavelengths to blue (i.e., JWST → Gaia R → TESS → Kepler → Gaia B), the amplitude becomes larger and maximum/minimum light occur at later phases. The amplitude trend arises because the visible and NIR fall on the Rayleigh-Jeans tail of the star’s spectral energy distribution, and therefore the flux variations due to temperature perturbations vary as \(\lambda^{-4}\). The phase trend is a consequence of non-adiabatic effects that cause (for this particular mode) the radius perturbations to lag the temperature perturbations by \(\sim 45^{\circ}\).

Note that the overall scaling of the light curves is not physically meaningful; as with any linear oscillation calculation, the normalization of perturbations is arbitrary. In the case of GYRE, this normalization is chosen to yield a mode inertia of \(M R^{2}\).

Footnotes

[1]

Generated by applying the specgrid_to_photgrid tool to an extended-wavelength version of the demo grid.

[2]

Generated using the MESA stellar evolution code.