Mixed modes

Mixed modes are usually difficult to deal with in peakbagging because their disposition differs from the regular pattern that can be found e.g. for solar-type stars on the main-sequence. This tutorial therefore presents the apollinaire tools specifically dedicated to \(\ell=1\) mixed mode analysis.

import apollinaire as apn

apn.__version__

'1.5.0'

We consider the case of KIC8684542, a red giant star where evidence of a magnetic signature in the mixed mode splittings was found by Li et al. (2022).

A first look at the data

We start by loading the data and taking a visual look at the power spectral density (PSD).

t, lc = apn.timeseries.load_light_curve (star=8684542)
dt = np.median (np.diff (t*86400))
freq, psd = apn.psd.series_to_psd (lc, dt)

fig, (ax1, ax2) = plt.subplots (1, 2, figsize=(10,4))

ax1.plot (freq*1e6, psd*1e-6, color="darkgrey")

ax1.set_xlabel (r"$\nu$ ($\mu$Hz)")
ax1.set_ylabel (r"PSD (ppm$^2 / \mu$Hz)")
ax1.set_xscale ("log")
ax1.set_yscale ("log")

cond = (freq>120e-6)&(freq<230e-6)
ax2.plot (freq[cond]*1e6, psd[cond]*1e-6, color="darkgrey")
ax2.set_xlabel (r"$\nu$ ($\mu$Hz)")

fig.tight_layout ()

Fitting the \(\ell=0,2\) acoustic pattern

The global seismic parameters are available in the literature, so we can directly use them here. The first step of our analysis consist in fitting the pattern of pure acoustic modes \(\ell=0,2\).

teff = 4796 # From the TASOC database
dnu = 13.484 #From Li 2022
numax = 176.8 #From Li 2022

mcmcDir = "mixed_modes_tutorial"
if not os.path.exists (mcmcDir) :
    os.mkdir (mcmcDir)

apn.peakbagging.stellar_framework(freq*1e6, psd*1e-6,
                                  dnu=dnu,
                                  numax=numax,
                                  back_mle=True,
                                  teff=teff,
                                  low_cut=1,
                                  format_cornerplot='png',
                                  parallelise=True,
                                  fit_l1=False,
                                  hide_l1=True,
                                  n_width=10,
                                  progress=True,
                                  mcmcDir=mcmcDir,
                                  filename_back=os.path.join (mcmcDir, 'background.png'),
                                  filemcmc_back=os.path.join (mcmcDir, 'mcmc_sampler_background.h5'),
                                  filename_pattern=os.path.join (mcmcDir, 'pattern.png'),
                                  filemcmc_pattern=os.path.join (mcmcDir, 'mcmc_sampler_pattern.h5'),
                                  filename_peakbagging=os.path.join (mcmcDir, 'summary_peakbagging.png'),
                                  a2z_file=os.path.join (mcmcDir, 'modes_param.a2z'),
                                  n_order=5,
                                  n_order_peakbagging=5)

Searching for mixed mode frequencies

Now that the reference pattern surrounding the mixed modes is identified, we can attempt at detecting the power excess corresponding to them. We follow the idea of Corsaro (2019) and Corsaro et al. (2020) to identify these power excesses thanks to a high-threshold nested sampling exploration. In order to identify a maximum of mixed mode frequency, the nested sampling exploration is deliberately left non-converged.

back = np.loadtxt (os.path.join (mcmcDir, 'background.dat'))
pkb = np.loadtxt (os.path.join (mcmcDir, 'modes_param.pkb'))

l1_freq, l1_np, samplers, fig = apn.peakbagging.search_mixed_mode_frequencies (
                                   pkb, freq*1e6, psd*1e-6, back,
                                   parallelise=True,
                                   verbose=True,
                                   bounds_log_a=(1,10),
                                   burnt_in=2000,
                                   bins=60,
                                   count_threshold=75,
                                   neighbour_threshold=3,
                                   n_width=1,
                                   linewidth=0.1,
                                   filename=os.path.join (mcmcDir, 'mixed_mode_search.dat'),
                                   dpi=300,
                                   figsize=(10,6),
                                   seed=2488557248,
                                   ignore_existing=True)

Fitting the mixed modes

Now that we have a good guess of the mixed mode frequencies, we can attempt at fitting their properties more in details.

(labels,
 mode_indexes,
 p_orders,
 median,
 e_param,
 E_param) = apn.peakbagging.explore_distribution_mixed_modes (l1_freq, pkb,
                                                              freq*1e6, psd*1e-6, back,
                                                              l1_np=l1_np,
                                                              fit_splittings=True,
                                                              fmin=None, fmax=None,
                                                              individual_amplitudes=True,
                                                              width_init=0.1,
                                                              max_split=0.5,
                                                              hide_l02=False,
                                                              n_width=1,
                                                              f_window=2.5,
                                                              filename=os.path.join (mcmcDir, 'mixed_mode_splittings.h5'),
                                                              save_chain=True,
                                                              existing_chains='sample',
                                                              parallelise=True,
                                                              instr="geometric",
                                                              progress=True,
                                                              nwalkers=256,
                                                              nsteps=0,
                                                              thin=20,
                                                              discard=750,
                                                              show_corner=True,
                                                              format_cornerplot='png',
                                                              bins=30,
                                                              verbose=True,
                                                              plot_datapoints=False,
                                                              dpi=50)

Highlighting the splitting asymmetries

From the fit we performed, we can now compute the asymmetry of the mode splittings

cond = labels=="freq"
(nu_zonal,
 e_nu_zonal,
 E_nu_zonal) = median[cond], e_param[cond], E_param[cond]
cond = labels=="split_prograde"
(split_prograde,
 e_split_prograde,
 E_split_prograde) = median[cond], e_param[cond], E_param[cond]
cond = labels=="split_retrograde"
(split_retrograde,
 e_split_retrograde,
 E_split_retrograde) = median[cond], e_param[cond], E_param[cond]

delta_asym = (- split_retrograde + split_prograde)*1e3
e_delta_asym = np.sqrt (e_split_prograde**2 + e_split_retrograde**2)*1e3
E_delta_asym = np.sqrt (E_split_prograde**2 + E_split_retrograde**2)*1e3

fig, ax = plt.subplots (1, 1, figsize=(5,4))

# Showing only component with a clear asymmetry detection
cond = delta_asym - e_delta_asym > 0

ax.errorbar (nu_zonal[cond], delta_asym[cond],
             xerr=[e_nu_zonal[cond], E_nu_zonal[cond]],
             yerr=[e_delta_asym[cond], E_delta_asym[cond]],
            capsize=3, color="blue", marker="o", ls="--")

ax.set_xlabel (r"Frequency ($\mu$Hz)")
ax.set_ylabel (r"$\delta_{\rm asym}$ (nHz)")