The bright star named KIC 10963064 is a solar-like, main sequence star that was observed by NASA's Kepler spacecraft between 2009 and 2013. Because the object is fairly bright (9th magnitude), and because Kepler observed it in its high-cadence mode (1 minute exposures), the Kepler data for this object is of extremely high quality.
The quality of the data is so good, in fact, that it can be used to study the minor brightness fluctuations in the star that are caused by sound waves propagating through the star's atmosphere. The study of such oscillations can reveal the mass, size, and internal structure of a star using a technique that astronomers call asteroseismology.
This tutorial will demonstrate how the Lightkurve Python package can be used to carry out a quick-look asteroseismic analysis of the star in five easy steps:
We'll start by searching Kepler's public data archive hosted by MAST for files which contain the short-cadence light curves recorded by Kepler:
import lightkurve as lk
search = lk.search_lightcurvefile('KIC10963065', cadence='short', mission='Kepler')
print(search)
SearchResult containing 27 data products. # observation target_name productFilename distance --- ----------------- ------------- ------------------------------------ -------- 0 Kepler Quarter 2 kplr010963065 kplr010963065-2009259162342_slc.fits 0.0 1 Kepler Quarter 5 kplr010963065 kplr010963065-2010111051353_slc.fits 0.0 2 Kepler Quarter 5 kplr010963065 kplr010963065-2010140023957_slc.fits 0.0 3 Kepler Quarter 5 kplr010963065 kplr010963065-2010174090439_slc.fits 0.0 4 Kepler Quarter 6 kplr010963065 kplr010963065-2010203174610_slc.fits 0.0 5 Kepler Quarter 6 kplr010963065 kplr010963065-2010234115140_slc.fits 0.0 6 Kepler Quarter 6 kplr010963065 kplr010963065-2010265121752_slc.fits 0.0 7 Kepler Quarter 7 kplr010963065 kplr010963065-2010296114515_slc.fits 0.0 8 Kepler Quarter 7 kplr010963065 kplr010963065-2010326094124_slc.fits 0.0 9 Kepler Quarter 7 kplr010963065 kplr010963065-2010355172524_slc.fits 0.0 ... ... ... ... ... 16 Kepler Quarter 13 kplr010963065 kplr010963065-2012121044856_slc.fits 0.0 17 Kepler Quarter 13 kplr010963065 kplr010963065-2012151031540_slc.fits 0.0 18 Kepler Quarter 13 kplr010963065 kplr010963065-2012179063303_slc.fits 0.0 19 Kepler Quarter 14 kplr010963065 kplr010963065-2012211050319_slc.fits 0.0 20 Kepler Quarter 14 kplr010963065 kplr010963065-2012242122129_slc.fits 0.0 21 Kepler Quarter 14 kplr010963065 kplr010963065-2012277125453_slc.fits 0.0 22 Kepler Quarter 15 kplr010963065 kplr010963065-2012310112549_slc.fits 0.0 23 Kepler Quarter 15 kplr010963065 kplr010963065-2012341132017_slc.fits 0.0 24 Kepler Quarter 15 kplr010963065 kplr010963065-2013011073258_slc.fits 0.0 25 Kepler Quarter 17 kplr010963065 kplr010963065-2013121191144_slc.fits 0.0 26 Kepler Quarter 17 kplr010963065 kplr010963065-2013131215648_slc.fits 0.0 Length = 27 rows
It looks like the target was observed in Kepler's Quarters 2 through 17. We'll only download the data for Quarters 5 through 7 to make this tutorial faster to run, but everything we'll do below should work with the full data set as well!
files = search[1:10].download_all()
files
LightCurveFileCollection of 9 objects: KIC 10963065 (9 KeplerLightCurveFiles) Quarters: 5,5,5,6,6,6,7,7,7
The data is now available to us as a LightCurveFileCollection
object, which groups a number of light curve files.
Lightkurve provides a convenience method to extract the recommended PDCSAP_FLUX
light curve from each file, and stitch them together into a single LightCurve
object:
lc = files.PDCSAP_FLUX.stitch()
Next, let's inspect our time series by plotting the data:
lc.plot();
The light curve looks a bit messy due to the presence of outliers. This can be explained by a variety of detector effects (e.g. cosmic rays). Discussing these outliers is the topic of a different tutorial. For now, let's tidy the light curve by removing the outliers and NaN values:
lc = lc.remove_outliers().remove_nans()
lc.plot();
That looks better!
Next, we'll carry out a Lomb Scargle Periodogram analysis on the light curve.
Note that astronomers who study solar-like oscillations tend to prefer a periodogram to use units of "power spectral density" (psd), which simply means that the X-axis will show frequency units in $\mu Hz$ and the Y-axis will show units in power density instead of amplitude. (Hint: give normalization = 'amplitude'
a try to see a different type of normalization.)
pg = lc.to_periodogram(method='lombscargle', normalization='psd')
pg.plot(scale='log');
The figure shown above is often called a power spectrum, and the bump near $2000\, \mu \rm Hz$ is called a power excess.
Let's take a closer look at the area near the power excess. For this, we can use the minimum_frequency
and maximum_frequency
keywords. This truncates the frequency range evaluated when performing the Lomb Scargle transformation.
(Hint: you can also set minimum_period
and maximum_period
if you prefer using periods instead of frequencies!)
pg = lc.to_periodogram(method='lombscargle', normalization='psd',
minimum_frequency=1000, maximum_frequency=3000)
To get a better look at the data, we'll also plot a smoothed periodogram over the top using the periodogram's smooth()
method.
ax = pg.plot()
pg.smooth(method='boxkernel', filter_width=1.).plot(ax=ax, label='Smoothed', c='red', lw=2);
An important measurement in asteroseismology is the central frequency of the power peaks, which is often denoted as $\nu_{\rm max}$ (pronounced "nu max").
A very easy (but somewhat inaccurate) estimators for this quantity is the frequency of the highest peak in the periodogram. For convenience, you can easily obtain the highest peak in a periodogram object using its frequency_at_max_power
property:
pg.frequency_at_max_power
The highest peak is located near $2244 \mu Hz$.
You can also use the show_properties()
method to inspect this and all other properties of the periodogram object:
pg.show_properties()
lightkurve.Periodogram properties: Attribute Description Units ---------------------- -------------- ------- nterms 1 targetid 10963065 default_view frequency label KIC 10963065 ls_method fast frequency_at_max_power 2244.0321 uHz max_power 0.0 1 / uHz nyquist 8496.4268 uHz period_at_max_power 0.0004 1 / uHz frequency array (47707,) uHz period array (47707,) 1 / uHz power array (47707,) 1 / uHz meta <class 'dict'>
The value of pg.frequency_at_max_power
is not suitable for a proper asteroseismic analysis however, because the properties of a star scale with the center of the envelope that surrounds the peaks in the power spectrum, rather than the frequency of the maximum peak. We can see that pg.frequency_at_max_power
is offset from the center of power:
ax = pg.plot()
ax.axvline(pg.frequency_at_max_power.value, lw=5, ls='dashed');
We will improve our estimate of $\nu_{\rm max}$ next.
The first step towards estimating $\nu_{\rm max}$ and $\Delta\nu$ is to remove the background from the periodogram, so that we obtain a Signal-to-Noise (SNR) spectrum. This will prevent our analysis from being confused by low-frequency noise.
Lightkurve provides a fairly primitive flattening tool for quicklook analyses. More advanced pipelines exist which e.g. offer the user to fit the granulation background in a more careful way.
The flatten()
method returns a SNRPeriodogram
. This behaves in the same way as a regular LombScarglePeriodogram
, except the power is now a unitless signal-to-noise ratio.
snr = pg.flatten()
snr.plot();
From this periodogram, we can create a Seismology
helper object which will enable us to extract interesting quantities from the periodogram:
seis = snr.to_seismology()
seis
Seismology(ID: KIC 10963065) - no values have been computed so far.
You can see that seismology is letting us know no values have been estimated so far.
The Seismology
class provides an estimate_numax()
method which uses an Autocorrelation Function (following the techniques in Viani et al. 2019) to find where the power excess is highest in a shape approximately close to what we would expect for a solar-like oscillator's modes.
numax = seis.estimate_numax()
numax
It is not wise to trust the value returned by this method blindly. Fortunately, we have the option to verify the accuracy of the estimate using the diagnose_numax()
method, which returns a set of diagnostic plots:
seis.diagnose_numax();
These diagnostic plots suggest that the 2D autocorrelation method returned a reasonable estimate for $\nu_{\rm max}$. At time of writing, Lightkurve does not yet return uncertainties on these values. This may change in the future!
The next step is finding $\Delta\nu$. This also uses an autocorrelation method, but increases the width of the window so it covers the whole mode envelope.
Passing in the $\nu_{\rm max}$ autocorrelates a region approximately equivalent to the spread of the seismic mode envelope. Because this is where the p-modes reside, the resulting function should show neatly spaced peaks at a frequency lag equal to $\Delta\nu$. Seismology stores all its calculated values. You can pass a numax into this function, but ideally we'll want to have run .estimate_numax()
first.
We estimate the $\Delta\nu$ first, and then have a look at the diagnostics:
deltanu = seis.estimate_deltanu()
ax = seis.diagnose_deltanu();
That looks sensible. More advanced mode-fitting pipelines will likely be able to do slightly better, but this is a good first-look result.
The real hallmark of a good $\Delta\nu$ estimate is that it produces a nice echelle diagram. All periodogram objects have access to the plot_echelle()
method, which only requires a $\Delta\nu$. If we pass it a $\nu_{\rm max}$ too, it just shows the area around the seismic mode envelope we care about!
seis.plot_echelle();
If you prefer to smooth your echelle diagrams a little, we can do so using the smooth_filter_width
kwarg, which calls seismology.periodogram.smooth(filter_width=smooth_filter_width)
internally. We can also play around with the plotting keywords a little, to make it more visible.
seis.plot_echelle(smooth_filter_width=3., scale='log', cmap='viridis');
Thanks to the high signal-to-noise of these data, you can see the radial, dipole, and quadrupole modes clearly in the spectrum, including slight curvature in the mode frequencies.
If we have a look at the Seismology module now, we'll notice it's holding on to our estimated $\nu_{\rm max}$ and $\Delta\nu$ values! We can recover them as attributes of seismology:
print(seis)
seis.deltanu
Seismology(ID: KIC 10963065) - computed values: * numax: 2175.00 uHz (method: ACF2D) * deltanu: 103.05 uHz (method: ACF2D)
Finally, let's use these asteroseismic quantities to estimate the stellar parameters of the target star.
Lightkurve provides easy access to basis asteroseismic scaling relationships. This comes with a few caveats, e.g. uncertainties are not supported at this time, and Lightkurve is not not equipped (yet) to deal with corrections to the scaling relations that are required to get accurate stellar parameters for more evolved stars.
Using the values reported in the Asteroseismic Modelling Portal (AMP) reported in Mathur et al. 2012, we know the temperature of our target to be $T_{\rm eff} = 6046\, K$, and we have pre-exiting literature values available:
Teff = 6046
trueradius = 1.20
truemass = 1.03
truelogg = 4.29
Here is how we can approach these values using Lightkurve's built-in scaling relations:
mass = seis.estimate_mass(Teff)
mass
We can do the same again for radius and surface gravity.
radius = seis.estimate_radius(Teff)
radius
logg = seis.estimate_logg(Teff)
logg
The literature values for these are:
print('Radius : {} vs {:.2f} Rsol'.format(trueradius, radius.value))
print('Mass : {} vs {:.2f} Msol'.format(truemass, mass.value))
print('logg : {} vs {:.2f} dex'.format(truelogg, logg.value))
Radius : 1.2 vs 1.24 Rsol Mass : 1.03 vs 1.10 Msol logg : 4.29 vs 4.30 dex
That's a good match!
This ends this tutorial on Lightkurve's suite of asteroseismology tools. We hope it will make quick-look asteroseismic analyses more accessible!