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 GYREpygyre
— Python module for readingsummary.h5
gyre.in
— GYRE input file used to createsummary.h5
spb.mesa
— stellar model file used to createsummary.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:
[1]:
# 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.
[2]:
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:
[3]:
# 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:
[4]:
# 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()
[4]:
<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