Tutorial 5: Using the LDTk limb darkening model

The TSModel transmission spectroscopy transit model in ExoIris is based on PyTransit’s RRModel and includes its ability to model stellar limb darkening using an arbitrary rotationally symmetric function. The LDTk limb darkening model (exoiris.LDTkLD) takes advantage of this and models stellar limb darkening using theoretical PHOENIX limb darkening profiles directly instead of representing it as simple analytical functions.

for all the passbands parameterised by three stellar parameters: the effective temperature, log g, and metallicity.

The LDTk limb darkening is parameterised by (\(T_\mathrm{Eff}\), \(\log g\), \(z\)) model calculates the limb darkening profiles automatically for all passbands defined by a set of filters for a set of stellar parameters. This is very useful for transmission spectroscopy since the number of limb darkening parameters stays constant (three) no matter how many passbands are modelled simultaneously.

Author: Hannu Parviainen Edited: 7 November 2024

%run ../setup_multiprocessing.py

%matplotlib inline

from multiprocessing import Pool
from xarray import load_dataset
from scipy.interpolate import splev, splrep
from numpy import array, geomspace, linspace
from matplotlib.pyplot import subplots, setp

from exoiris import ExoIris, TSData, LDTkLD
from ldtk import BoxcarFilter

Data preparation

def read_data(fname, name=""):
    with load_dataset(fname) as ds:
        return TSData(time=ds.time.values, wavelength=ds.wavelength.values, fluxes=ds.flux.values, errors=ds.error.values, name=name)

d1 = read_data('data/nirHiss_order_1.h5', "WASP-39b JWST NIRISS Order 1")
d1.remove_outliers()
d1.mask_transit(t0=2459783.5015, p=4.0552842, t14=0.13)
d1.normalize_to_poly()

d2 = read_data('data/nirHiss_order_2.h5', "WASP-39b JWST NIRISS Order 2")
d2.remove_outliers()
d2.crop_wavelength(d2.bbox_wl[0], 1.0)
d2.mask_transit(ephemeris=d1.ephemeris)
d2.normalize_to_poly()

TSData Name:'WASP-39b JWST NIRISS Order 2' [0.57 - 1.00] nwl=957 npt=518

db = d2.bin_wavelength(r = 100) + d1.bin_wavelength(r = 100)
fig = db.plot()
fig.tight_layout()

Setup the LDTk limb darkening model

The way LDTkLD works is familiar to anyone who has used LDTk before. You initialise the model by giving it the stellar effective temperature, surface gravity, and metallicity estimates as (mean, \(1\sigma\) uncertainty) tuples and a set of filters. You can also optionally tell the model which theoretical stellar spectrum dataset it should use (vis, vis-lowres, visir, or visir-lowres). For JWST transmission spectroscopy, one should choose either visir, or visir-lowres.

ldm = LDTkLD(db, teff=(5327, 139), logg=(4.38, 0.09), metal=(-0.01, 0.1), dataset='visir-lowres')

ts = ExoIris('05a_ldtk_limb_darkening', ldmodel=ldm, data=db, nk=50, nldc=10, nthreads=1)

Now, when we look at the model parameterisation, we notice the limb darkening is parameterised only by the three stellar parameters teff, logg, and metal. The priors are taken automatically from LDTkLD.

ts.ps[:10]

[  0 |G| rho            U(a = 0.1, b = 25.0)                     [    0.00 ..      inf],
   1 |G| tc             N(μ = 0.0, σ = 0.1)                      [    -inf ..      inf],
   2 |G| p              N(μ = 1.0, σ = 1e-05)                    [    0.00 ..      inf],
   3 |G| b              U(a = 0.0, b = 1.0)                      [    0.00 ..      inf],
   4 |G| secw           N(μ = 0.0, σ = 1e-05)                    [   -1.00 ..     1.00],
   5 |G| sesw           N(μ = 0.0, σ = 1e-05)                    [   -1.00 ..     1.00],
   6 |G| teff           N(μ = 5327.0, σ = 139.0)                 [    0.00 ..      inf],
   7 |G| logg           N(μ = 4.38, σ = 0.09)                    [    0.00 ..      inf],
   8 |G| metal          N(μ = -0.01, σ = 0.1)                    [    -inf ..      inf],
   9 |G| k_00.57318     U(a = 0.02, b = 0.2)                     [    0.00 ..      inf]]

ts.set_prior('tc', 'NP', 2459694.286, 0.003)
ts.set_prior('p', 'NP', 4.05487, 1e-5)

ts.set_radius_ratio_prior('UP', 0.14, 0.15)

Customise the radius ratio knot locations

ts.set_radius_ratio_knots(geomspace(ts.data[0].wavelength[5], ts.data[1].wavelength[-5], 50))
ts.add_radius_ratio_knots(linspace(0.768-0.015, 0.768+0.015, 5))
ts.set_radius_ratio_prior('UP', 0.14, 0.15)
ts.plot_setup();

ts.fit_white()
fig = ts.plot_white(figsize=(13,4))
fig.tight_layout()

ts.save(overwrite=True)

Set up multiprocessing

Now we’re almost ready for fitting and MCMC sampling, but since we’re using multiprocessing to parallelise the process, we first need to take some extra steps to make sure everything works the way it’s supposed to.

First, we need to define a log-posterior function that calls the ExoIris log-posterior method. This must be done so that Python can pickle the method for parallelisation (if you know a better way, please let me know). Next, we also create a multiprocessing pool that will be used by the global optimiser and the MCMC sampler.

def lnpostf(pv):
    return ts.lnposterior(pv)

pool = Pool(8)

Fit the transmission spectrum

ts.fit(niter=2500, npop=250, pool=pool, lnpost=lnpostf)

ts.plot_fit(result='fit', figsize=(13,10), height_ratios=(0.3, 0.7),
            res_args=dict(pmin=5, pmax=95),
            trs_args=dict(xscale='log', ylim=(2.0,2.25), xticks=[0.6, 0.88, 1.16, 1.44, 1.72, 2.0, 2.30, 2.8]));

We can immediately see from the residual plot that the current best-fit solution is bad. This is something we could have guessed already by looking at the log posterior distribution in the optimiser plot that shows a spread of tens of thousands. Let’s continue the fitting for another 25 iterations. Each successive ExoIris.fit call continues optimisation from the solution of the previous call, and it also plots the old and new posterior and parameter populations at the end of the fit to visualise how the population is changing.

So, it is clear at this stage that running the optimiser for 25 iterations is not doing much. Let’s continue the optimisation for 5000 iterations and see if this does the trick. This should take 2-3 min (or less if you initialise the pool with more processes) and the optimisation should finish before it reaches 5000 iterations (the progress bar turns red).

ts.fit(niter=5000, npop=250, pool=pool, lnpost=lnpostf)

Let’s plot everything once again. This time, the residuals look good (there’s no trace of a transit signal there), the transmission spectrum looks something else than just noise, and the limb darkening parameters also look good, so we’re ready to move to the final step: the MCMC sampling.

ts.plot_fit(result='fit', figsize=(13,7),
            res_args=dict(pmin=5, pmax=95),
            trs_args=dict(xscale='log', ylim=(2.0,2.25), xticks=[0.6, 0.88, 1.16, 1.44, 1.72, 2.0, 2.30, 2.8]))

Save the optimisation results

ts.save(overwrite=True)

MCMC sampling

Next comes the final part, obtaining a posterior sample using Markov Chain Monte Carlo (MCMC) sampling. This is done using the ExoIris.sample method, that starts with the parameter vector population from the global optimisation.

As a first step, let’s still inflate the limb darkening parameter prior widths by ten to make sure we’re not constraining them too much in the sampling phase. So, now we set uncertainty_multiplier=100 instead of uncertainty_multiplier=10 that we used in fitting.

ts.sample(1000, thin=100, repeats=3, pool=pool, lnpost=lnpostf)

df = ts._tsa.posterior_samples()

import ldtk

ldtk.__version__

Version('1.8.4')

ts.plot_fit(result='mcmc', figsize=(13,7),
            res_args=dict(pmin=5, pmax=95),
            trs_args=dict(xscale='log', ylim=(2.0,2.25), xticks=[0.6, 0.88, 1.16, 1.44, 1.72, 2.0, 2.30, 2.8]));

Save the results and export the transmission spectrum

Congratulations, you now have a low-resolution transmission spectrum for your observations! Let’s save the model one more time to also store the posterior samples inside the fits file.

ts.save(overwrite=True)

ts.transmission_spectrum

Table length=175

wavelength	radius_ratio	radius_ratio_e	area_ratio	area_ratio_e
um
float64	float64	float64	float64	float64
0.5731806007740325	0.14015972198590915	0.0002946211952941102	0.019644747667203616	8.260477409168684e-05
0.5789412098270379	0.14255500096953785	0.0001970778993912647	0.02032192830145446	5.6193077320264214e-05
0.5847597144484151	0.1442734705162264	0.00016348873700431724	0.02081483429480218	4.717799283285726e-05
0.5906366965031731	0.14540652332722578	0.00016201349387189487	0.021143057026111134	4.712005749875868e-05
0.59657274370421	0.14604012229579705	0.0001621752253652372	0.021327717320235594	4.737224591468161e-05
0.6025684496710865	0.14626237897238517	0.00015550270570965281	0.021392683502665752	4.5490326640701425e-05
0.6086244139893888	0.14615818417316961	0.0001463822205448086	0.021362214800803664	4.2790612092112884e-05
0.6147412422706893	0.14582020417274466	0.0001430773724558657	0.021263531944991892	4.1729795817659045e-05
0.6209195462131082	0.1453646369239378	0.00014944788735951639	0.02113087766802886	4.345358064343948e-05
...	...	...	...	...
2.611551468253669	0.14635091939362085	0.0005656118463229492	0.02141859160737751	0.00016556809356418434
2.637798216678329	0.14660253414539987	0.0005386993533271586	0.021492303017863576	0.0001579321469294182
2.6643087515193176	0.14683599772480388	0.0005736270827847371	0.021560810228489276	0.00016844435880324725
2.6910857238963963	0.14701690595236902	0.0006587785775378388	0.02161397063582576	0.00019371368161863418
2.718131811573747	0.14713266864626942	0.0007277630223983029	0.021648022182974812	0.00021414604848910907
2.7454497192277545	0.14710933590575853	0.0007336594024687242	0.02164115671076942	0.00021582529436975226
2.77304217871748	0.14688344691892896	0.0007366228927499031	0.02157474697903777	0.0002163720976768297
2.800911949357857	0.14645971011502953	0.0010294058779203182	0.021450446687030122	0.00030135243153311947
2.8215442681551877	0.1460142722140494	0.001582545896465745	0.021320167698805948	0.00046154002777712774

---

©2024 Hannu Parviainen