Simple Harmonic Oscillator Celerite Kernels for Stellar Granulation and Oscillations

This is the documentation for shocksgo. The goal of shocksgo is to generate light curves of stars accounting for the effects of granulation, super-granulation and p-mode oscillations.

You can view the source code and/or contribute to shocksgo via GitHub.

Documentation

Installing shocksgo

You can install shocksgo with pip:

pip install shocksgo

You can install shocksgo from the source code by doing the following:

git clone https://github.com/bmorris3/shocksgo.git
cd shocksgo
python setup.py install

shocksgo requires python >=3.5, numpy, celerite, and astropy.

If you have any trouble installing shocksgo, feel free to post an issue on GitHub.

Getting Started

This tutorial will show some examples of how to make solar and stellar light curves using shocksgo.

Generating a Solar Light Curve

To generate a sample of ten hours of solar fluxes at 60 second cadence, we can use generate_solar_fluxes:

import matplotlib.pyplot as plt
import astropy.units as u

from shocksgo import generate_solar_fluxes

times, fluxes, kernel = generate_solar_fluxes(duration=10*u.hour, cadence=60*u.s)

plt.plot(times.to(u.hour), 1e6 * fluxes)
plt.gca().set(xlabel='Time [hours]', ylabel='Relative Flux [ppm]')
plt.show()

(Source code)

We can check that the power spectrum of the fluxes that we’ve generated reproduce the solar power spectrum:

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u

from shocksgo import generate_solar_fluxes, power_spectrum

times, fluxes, kernel = generate_solar_fluxes(duration=100*u.day, cadence=60*u.s)

freq, power = power_spectrum(fluxes, d=60)

plt.loglog(freq * 1e6, power, ',', label='Samples')
plt.loglog(freq * 1e6, 1e6 * kernel.get_psd(2*np.pi*freq)/(2*np.pi), alpha=0.7, label='Kernel')
plt.ylim([1e-5, 1e3])
plt.xlim([1e-2, 1e4])
plt.gca().set(xlabel='Frequency [$\mu$Hz]', ylabel='Power [ppm$^2$/$\mu$Hz]')
plt.show()

(Source code)

Zooming into the p-mode oscillations, we can see the peaks are reproduced:

(Source code)

Generating a Stellar Light Curve

To generate a sample of steller fluxes at 60 second cadence, we can use generate_stellar_fluxes:

import matplotlib.pyplot as plt
import astropy.units as u
from astropy.constants import M_sun, L_sun, R_sun

from shocksgo import generate_stellar_fluxes

# Stellar properties
M = 0.9 * M_sun
T_eff = 5340 * u.K
L = 0.56 * L_sun
R = 0.7 * R_sun

times, fluxes, kernel = generate_stellar_fluxes(duration=100*u.day, M=M, T_eff=T_eff, R=R, L=L, cadence=60*u.s)

plt.plot(times.to(u.day), 1e6 * fluxes)
plt.gca().set(xlabel='Time [days]', ylabel='Relative Flux [ppm]', title='G9V star')
plt.show()

(Source code)

We can see the shift in the p-mode oscillations relative to the solar ones above if we plot the power spectrum:

import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.constants import M_sun, L_sun, R_sun

from shocksgo import generate_stellar_fluxes, power_spectrum

# Stellar properties
M = 0.9 * M_sun
T_eff = 5340 * u.K
L = 0.56 * L_sun
R = 0.876 * R_sun

times, fluxes, kernel = generate_stellar_fluxes(duration=10*u.day, M=M, T_eff=T_eff, R=R, L=L, cadence=60*u.s)

freq, power = power_spectrum(fluxes, d=60)

plt.semilogy(freq * 1e6, power, ',', label='Samples')
plt.semilogy(freq * 1e6, 1e6 * kernel.get_psd(2*np.pi*freq)/(2*np.pi), alpha=0.7, label='Kernel')
plt.ylim([1e-5, 1e-1])
plt.xlim([2500, 5000])
plt.gca().set(xlabel='Frequency [$\mu$Hz]', ylabel='Power [ppm$^2$/$\mu$Hz]')
plt.show()

(Source code)

Custom Frequencies

Suppose you have a list of model p-mode frequencies, and you would like to generate a light curve from your custom list of frequencies (without scaling from the solar values). You can do so using a different set of keyword arguments in the generate_stellar_fluxes function, like so:

import numpy as np
import matplotlib.pyplot as plt
import astropy.units as u
from astropy.constants import R_sun, M_sun, L_sun

from shocksgo import generate_stellar_fluxes

M = 1*M_sun
L = 1*L_sun
T_eff = 5777 * u.K
R = 1*R_sun

freqs = np.linspace(2000, 4000, 10)  # in microHertz
log_amps = np.exp(-0.5 * (freqs - 3000)**2 / 1000**2) - 32
log_lifetimes = np.ones_like(freqs) * 7
duration = 2 * u.min

times, fluxes, kernel = generate_stellar_fluxes(duration, M, T_eff, R, L,
                                                frequencies=freqs,
                                                log_amplitudes=log_amps,
                                                log_mode_lifetimes=log_lifetimes)

test_freqs = np.logspace(-2, np.log10(freqs.max()), 1e6)
plt.loglog(test_freqs, 1e6/(2*np.pi) * kernel.get_psd(2*np.pi*test_freqs*1e-6))
plt.gca().set(xlabel='Frequency [$\mu$Hz]', ylabel='Power [ppm$^2$/$\mu$Hz]')
plt.show()

(Source code)

shocksgo Documentation

This is the documentation for shocksgo.

Reference/API

shocksgo Package

Functions

generate_solar_fluxes(duration[, cadence])	Generate an array of fluxes with zero mean which mimic the power spectrum of the SOHO/VIRGO SPM observations.
generate_stellar_fluxes(duration, M, T_eff, R, L)	Generate an array of fluxes with zero mean which mimic the power spectrum of the SOHO/VIRGO SPM observations, scaled for a star with a given mass, effective temperature, luminosity and radius.
interpolate_missing_data(times, fluxes[, …])	Assuming times are uniformly spaced with missing cadences, fill in the missing cadences with linear interpolation.
power_spectrum(fluxes[, d])	Compute the power spectrum of fluxes in units of [ppm^2 / microHz].
test(**kwargs)	Run the tests for the package.

Classes

UnsupportedPythonError

Class Inheritance Diagram

Overview

Methods

We compute these light curves efficiently by taking advantage of celerite, a fast Gaussian process regression package, which we use to approximate solar and stellar power spectrum densities with sums of simple harmonic oscillator (SHO) kernels of the form:

\[S(\omega) = \sqrt{\frac{2}{\pi}} \frac{S_0\,\omega_0^4} {(\omega^2-{\omega_0}^2)^2 + {\omega_0}^2\,\omega^2/Q^2}\]

where \(\omega = 2\pi f\) is the angular frequency. We use one SHO kernel term for super/meso-granulation, another for ordinary granulation, and about 50 terms for the comb of p-mode peaks.

Scaling relations for p-modes

For computation of stellar p-mode oscillation frequencies, we use the scaling relations found in Huber et al. (2012) and references therein (e.g. Kjeldsen & Bedding 1995 ), namely Equation 4:

\[\nu_\textrm{max} \propto M R^{-2} T_{\rm eff}^{-1/2},\]

and Equation 3

\[\Delta \nu_\textrm{max} \propto M^{1/2} R^{-3/2}.\]

The amplitude scaling of the p-mode oscillations is given by Equation 9 of Huber et al. (2011):

\[A \propto \frac{L^s}{M^t T_\textrm{eff}^{r-1} c(T_\textrm{eff})}\]

where \(r = 2\), \(s = 0.886\), \(t = 1.89\) and

\[c(T_\textrm{eff}) = \left( \frac{T_\textrm{eff}}{5934 \textrm{K}} \right)^{0.8}.\]

Scaling relations for granulation

For computation of large and small scale stellar surface granulation frequencies, we use the scaling relation found in Kallinger et al. (2014):

\[\tau_\textrm{eff} \propto \nu^{-0.89}_\textrm{max},\]

where \(\tau_\textrm{eff}\) is the characteristic granulation timescale. The amplitudes of granulation scale as

\[a \propto \nu^{-2}_\textrm{max}.\]

(Kjeldsen & Bedding, 2011).