Asteroseismology is the study of stellar oscillations. To see them, we usually want to transfer from the time domain that light curves are in to the frequency domain. We can do that with methods such as Fourier Transforms or Lomb-Scargle Periodograms. lightkurve
has built in methods for working with data in the frequency domain, in the Periodogram
class.
Below we demonstrate some of the new functionalities in lightkurve
and how to use them for asteroseismology.
% load_ext autoreload
% autoreload 2
from lightkurve import KeplerTargetPixelFile
import matplotlib.pyplot as plt
import numpy as np
from lightkurve import log
import astropy.units as u
from tqdm import tqdm
log.setLevel('ERROR')
As an example, we can use a red giant star from Campaign 5. You can read more about asteroseismology with giants in K2 in Stello et al 2015 here.
First we'll use from_archive
to download the Target Pixel File for the target.
ID = 211416749
tpf = KeplerTargetPixelFile.from_archive(ID, campaign=5, cadence='short')
Be aware that this is a short cadence light curve, which means that this dataset is much larger than others we've used in previous tutorials. If you're rerunning this tutorial locally, it might take a few minutes to rerun the later steps.
tpf.flux.shape
(109343, 11, 11)
Let's plot the target below to see if the aperture is centered. We're plotting in 'log' scale here so that we can see the wings of the PSF.
tpf.plot(scale='log', aperture_mask=tpf.pipeline_mask)
<matplotlib.axes._subplots.AxesSubplot at 0x1c1f896438>
We can now create a KeplerLightCurve
using Simple Aperture Photometry with the to_lightcurve
method. We can then normalize, remove NaN values and remove outliers in one easy line. We can also use the new fill_gaps
method to fill in any gaps there are in our data using linear interpolation of the nearest neighbours. This creates an almost perfectly sampled time array with no gaps.
lc = tpf.to_lightcurve()
lc = lc.normalize().remove_nans().remove_outliers().fill_gaps()
We can plot our light curve easily with the plot
method.
lc.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x1c275fac50>
Let's take a look at the periodogram for this light curve. You can create a periodogram
object using lc.periodogram()
. This can be passed any array of frequencies that the user wants, but by default will create an array of frequency for you. You can also change the frequency units, depending on the range of frequencies you're looking for. In asteroseismology we usually use $\mu Hz$.
p = lc.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=10)
With this new periodogram object we can now plot the power spectrum of the data. Let's do that now!
ax = p.plot(c='k');
Unfortunately there is some ringing going on...there is a periodicity in the data due to the K2 motion. We can see below when we plot the data in "Period" space instead of "Frequency" space there is a significant periodicity at ~ 6 hours.
ax = p.plot(c='k', unit=u.hour, format='Period')
ax.text(6, 3e7, '6 Hour Roll Motion');
We can apply our own Self Flat Fielding (SFF) correction to the light curve to remove the K2 roll motion. Below we correct with SFF using the default correct
function, with windows=10
. (You can read more about SFF in our other tutorials.)
lc = tpf.to_lightcurve().normalize().remove_nans().remove_outliers()
clc = lc.correct(windows=10).remove_outliers().fill_gaps()
clc.plot()
<matplotlib.axes._subplots.AxesSubplot at 0x1c21dc89e8>
Our new, corrected lightcurve looks much flatter and has a much smaller CDPP. Let's try plotting the periodogram again.
p_clean = clc.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=10)
ax = p_clean.plot(c='k')
This is much better, although there is clearly some periodicity still left in the data. However, now we can see the oscillation modes of the target! Let's zoom in to the region of interest.
p_clean = clc.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=100)
ax = p_clean.plot(c='k')
We can even find out what effect our detrending parameters have on our target. Below we change the detrending windows from $w=2\; ... \; 10$ and plot the periodogram in each case. We can see different detrending parameters has an effect on each of the modes.
# Loop over several windows
for windows in tqdm(np.arange(2, 12, 2)):
# Create the light curve
lc = tpf.to_lightcurve().normalize().remove_nans().remove_outliers()
clc = lc.correct(windows=windows).remove_outliers().fill_gaps()
# Create the periodogram
p_clean = clc.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=100)
# Plot the periodogram
if windows == 2:
ax = p_clean.plot(alpha=0.4, label='Windows: {}'.format(windows))
else:
p_clean.plot(ax=ax, alpha=0.4, label='Windows: {}'.format(windows))
ax.legend()
100%|██████████| 5/5 [02:46<00:00, 33.38s/it]
<matplotlib.legend.Legend at 0x1c204118d0>
It looks like there is a significant effect on the oscillation modes as we vary the window size. When using K2 data for detecting these oscillations it is important to understand and mitigate the effects of detrending to ensure the modes have the highest possible signal to noise ratio.
We can also easily change the aperture size. If the aperture is too small, increasing the aperture size should allow us to include more of the target flux. If there is a contaminant nearby, decreasing the aperture may increase our singal to noise.
Below we create a new aperture for the target. In this case we've used an aperture where pixels have a value greater than the 70th percentile.
aper = np.nanmedian(tpf.flux, axis=0) > np.nanpercentile(np.nanmedian(tpf.flux, axis=0), 70)
/Users/ch/miniconda3/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:907: RuntimeWarning: All-NaN slice encountered result = np.apply_along_axis(_nanmedian1d, axis, a, overwrite_input) /Users/ch/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:1: RuntimeWarning: invalid value encountered in greater """Entry point for launching an IPython kernel.
We can plot up our new aperture against the pipeline aperture for comparison.
# Two subplots
fig, axs = plt.subplots(1, 2, figsize=(12, 4))
# Plot pipeline aperture mask
tpf.plot(axs[0], scale='log', aperture_mask=tpf.pipeline_mask)
axs[0].set_title('Pipeline Aperture')
# Plot larger aperture mask
tpf.plot(axs[1], scale='log', aperture_mask=aper)
axs[1].set_title('Larger Aperture')
Text(0.5,1,'Larger Aperture')
Now let's create our two light curves, one with the pipeline aperture and one with the new aperture.
tpf.to_lightcurve()
creates our SAP flux with the pipeline aperture as a default. To use a new aperture we simple use tpf.to_lightcurve(aperture_mask=NEW_APERTURE)
.
# Create the light curve with the pipelien aperture.
lc_pipeline = tpf.to_lightcurve().normalize().remove_nans().remove_outliers()
lc_pipeline = lc_pipeline.correct(windows=10).remove_outliers().fill_gaps()
# Create a light curve with a slightly larger aperture
lc_larger_aperture = tpf.to_lightcurve(aperture_mask=aper).normalize().remove_nans().remove_outliers()
lc_larger_aperture = lc_larger_aperture.correct(windows=10).remove_outliers().fill_gaps()
When we plot these two light curves we can see that the larger aperture is much less noisy.
#Plot the pipeline and large aperture light curves
ax = lc_pipeline.plot(label='Pipeline Aperture')
lc_larger_aperture.plot(ax=ax, label='Larger Aperture')
<matplotlib.axes._subplots.AxesSubplot at 0x1c213beb70>
Finally, when we plot the periodogram we can see we've increased the signal to noise ratio of our stellar oscillations.
# Create the periodograms
p_pipeline = lc_pipeline.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=100)
p_larger_aperture = lc_larger_aperture.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=100)
# Plot the periodograms
ax = p_pipeline.plot(c='k', alpha=0.5, label='Pipeline')
p_larger_aperture.plot(ax=ax, c='C3', alpha=0.5, label='Larger Aperture')
ax.legend();
By increasing our aperture size and including more pixels here we have increased the signal to noise of the oscillation modes of this red giant.
Finally, this target was actually observed twice by K2, once in April 2015 in Campaign 5 and once in December 2017 in Campaign 16. We can get both of these data sets and compare the results.
# Download the C16 TPF
tpf_c16 = KeplerTargetPixelFile.from_archive(ID, campaign=16, cadence='short')
# Create a new "large" aperture where flux is greater than the 70th percentile.
aper = np.nanmedian(tpf_c16.flux, axis=0) > np.nanpercentile(np.nanmedian(tpf_c16.flux, axis=0), 70)
# Create a new light curve for C16
lc_c16 = tpf_c16.to_lightcurve(aperture_mask='all').normalize().remove_nans().remove_outliers()
lc_c16 = lc_c16.correct(windows=10).remove_outliers().fill_gaps()
/Users/ch/miniconda3/lib/python3.6/site-packages/numpy/lib/nanfunctions.py:907: RuntimeWarning: All-NaN slice encountered result = np.apply_along_axis(_nanmedian1d, axis, a, overwrite_input) /Users/ch/miniconda3/lib/python3.6/site-packages/ipykernel_launcher.py:5: RuntimeWarning: invalid value encountered in greater """
# Create a periodogram for c16 data
p_c16 = lc_c16.periodogram(freq_unit=u.microHertz, max_frequency=400, min_frequency=100)
# Create subplots to plot into
ax = p_pipeline.plot(c='k', alpha=0.5, label='C5')
p_c16.plot(ax=ax, c='C3', alpha=0.5, label='C16')
ax.set_xlim(100, 400)
ax.set_ylim(0, 50000)
ax.legend()
<matplotlib.legend.Legend at 0x1c23cb19b0>
It looks like the two data sets provide similar modes, however the two campaigns have very different instrument systematics. To find the true answer, users should iterate over many different detrending parameters and aperture sizes, and combine datasets to increase signal to noise.