Retrievals¶

This section shows how to setup atmospheric retrievals with Pyrat Bay.

Setting up
Retrieval parameters
Full examples

Setting up¶

Multinest¶

Since version 2.0, Pyrat Bay enables atmospheric retrievals using the nested sampling algorithm [Skilling2004] [Skilling2006], via the MultiNest implementation [Feroz2009] [Buchner2014]. This is the recommended retrieval algorithm.

make sure to install multinest and MPI on your machine. This can be quite specific for each machine, so I cannot help much there. Here and here are some installation guides that may help.

Then install their Python wrappers, e.g., with these commands:

pip install pymultinest
pip install mpi4py

Configuration file¶

These are the requirements for the configuration file of a retrieval run:

General

# Pyrat Bay run mode: [tli atmosphere spectrum opacity retrieval radeq]
runmode = retrieval

# Output log and spectrum file names:
logfile = ret_wasp18b_all/WASP18b_eclipse_all.log

# Observing geometry [transit eclipse f_lambda]
rt_path = transit

First, we need to set runmode = retrieval in the configuration file. Second, set a folder (and optionally, a path) to store the output files. Lastly, specify the observing geometry.

Note

logfile can point to a folder that does not exist yet, which will be created when launching the run.

Observations

# The observations
dunits = ppm
obsfile = inputs/obs_wasp18b_eclipse_all.dat

A retrieval run must include the obsfile argument to define the observed dataset (a transit, eclipse, or direct-imaging spectrum) which will constrain the model parameters. TBD: link to obsfile explanation.

The dunits sets the units of the data outputs (i.e., figures and screen outputs).

Parameters

# Name      value   pmin    pmax    step  prior  prior_sigma
retrieval_params =
    T_iso       2850.0  500.0  3000.0    1.0
    M_planet    0.27      0.1     0.4    1.0  0.266  0.033
    log_p_ref   -1.0     -9.0     2.0    1.0
    [M/H]        1.5     -2.0     2.5    1.0
    C/O          0.5     -2.0     2.5    1.0

Define retrieval_params to set the free parameters to retrieve. The full list of possible parameters is shown in the Retrieval parameters section.

Sampler

# Retrieval setup:
sampler = multinest
nlive = 1000
resume = True
post_processing = True

# Retrieval temperature boundaries:
tlow  =  800
thigh = 2500

Finally, set the retrieval algorithm and other specific configurations. For multinest, we need to define the number of live points nlive.

The resume option indicates whether to resume sampling from a previous run, or start from scratch.

Launch / multi-processing¶

Retrievals runs in single-CPU mode can be started from the command-line as:

pbay -c wasp39b_retrieval_transit_jwst.cfg

Normally, we will want to use parallel computing to speed up our runs. Multinest runs can be parallelized via MPI. This is the command to do so, e.g., with 128 parallel CPUs:

mpirun -n 128 pbay -c wasp39b_retrieval_transit_jwst.cfg

Retrieval parameters¶

This is the list of all available free parameters (retrieval_params argument) to be fit in a retrieval. Note that there are requirements to enable some of them.

Type	Parameter	Model requirement	Comments
Temperature	`T_iso`	`tmodel = isothermal`	see Isothermal TP section

	`log_p1`	`tmodel = madhu`	see Madhu TP section
	`log_p2`
	`log_p3`
	`a1`
	`a2`
	`T0`

	`log_kappa'`	`tmodel = guillot`	see Guillot TP section
	`log_gamma1`
	`log_gamma2`
	`alpha`
	`T_irr`
	`T_int`

Composition	`log_X`	`log_X` in `vmr_vars`	Constant VMR profile. `X` is a species name
	`slope_X`	`slant_X` in `vmr_vars`	`X` is a species name
	`log_VMR0_X`
	`log_p0_X`
	`min_log_X`
	`max_log_X`

	`[M/H]`	`chemistry = tea`	metallicity scale factor
	`[X/H]`		`X` is an element name
	`X/Y`		`X` and `Y` are element names
	`log_X`	`chemistry = tea` and `log_X` in `vmr_vars`	Constant VMR embedded in equilibrium atmosphere. `X` is a species name

Isotopic ratios	`iso_X`	`X` in `isotope_ratios`	\( \log_{10}({f}) \), where \( f \) is the isotopic fraction for isotope `X`

Clouds	`log_p_cl`	`deck` in `clouds`	pressure at top of opaque cloud deck. \( \log_{10}(p/{\rm bar}) \) units
	`log_k_ray`	`lecavelier` in `rayleigh`
	`alpha_ray`	`lecavelier` in `rayleigh`
	`f_patchy`

Solo	`R_planet`		radius at `log_p_ref`
	`log_p_ref`		pressure at `R_planet`. \( \log_{10}(p/{\rm bar}) \) units
	`M_planet`
	`f_dilution`		Emission dilution factor as in TBD
	`T_eff`		stellar effective temperature
	`rv_shift`		radial-velocity offset in km s^-1

Data	`offset_X`	`offset_X` in `offset_inst`	shift depth to data points containing `X` in name. Units as given in `dunits` argument

	`scale_X`	`scale_X` in `uncert_scaling`	scale uncertainty of data points containing `X` in name
	`quadrature_X`	`quadrature_X` in `uncert_scaling`	add noise in quadrature to data points containing `X` in name. Units as given in `dunits` argument

The retrieval_params argument defines the free parameters that can be retrieved. See for example the extract below:

#   Name       value      min     max   step  prior  prior_sigma
retrieval_params =
    T_iso       2850    500.0  3000.0    1.0
    M_planet    0.27      0.1     0.4    1.0  0.266  0.033
    log_p_ref   -1.0     -9.0     2.0    1.0
    [M/H]        1.5     -2.0     3.0    1.0
    C/O          0.5      0.0     2.0    1.0

In the retrieval_params variable, users configure one free parameter per row, where each input field defines (from left to right):

the free-parameter name (see Table above)
an initial value (useful for forward-model calculations or when keeping a fixed value)
the minimum and maximum values to be explored for the parameter
a step which defines whether the parameter is free or remains fixed
(optionally) a Gaussian prior estimation and uncertainty

Priors¶

By default, parameters priors are assumed to be uniform between the minimum and maximum boundaries defined in retrieval_params. To impose Gaussian priors to a free parameter, users can define the prior value and uncertainties in the fields after step. In the example above, we impossed a Gaussian prior on the planetary mass of \(M = 0.266 \pm 0.033\) \(M_{\rm jup}\), while keeping unfiform priors for all other parameters.

Similarly, users can impose uneven priors by defining both prior uncertainties. For example, to impose a prior mass estimation of the form \(M = 0.266^{+0.3}_{-0.4}\) \(M_{\rm jup}\):

#   Name       value      min     max   step  prior  prior_lo  prior_hi
retrieval_params =
    T_iso       2850    500.0  3000.0    1.0
    M_planet    0.27      0.1     0.4    1.0  0.266     0.040     0.030
    log_p_ref   -1.0     -9.0     2.0    1.0
    [M/H]        1.5     -2.0     3.0    1.0
    C/O          0.5      0.0     2.0    1.0

Fixed parameters¶

Users can keep a parameter fixed by setting the step field to zero. In such case the parameter will adopt the provided value during a run. For example, to explore metallicities while keeping a solar C/O ratio, we would set:

#   Name       value      min     max   step  prior  prior_sigma
retrieval_params =
    T_iso       2850    500.0  3000.0    1.0
    M_planet    0.27      0.1     0.4    1.0  0.266  0.033
    log_p_ref   -1.0     -9.0     2.0    1.0
    [M/H]        0.0     -2.0     3.0    1.0
    C/O          0.59     0.0     2.0    0.0

Shared parameters¶

Users can have multiple parameters sharing the value of a single free parameter by setting step to a negative value, which is the index of the other parameter to share.

For example, if we want to explore abundances for an atmosphere with two parameters, one for the alkali metals (Na and K) and one for everything else, we can define a sodium free parameter [Na/H] (the fifth row in retrieval_params, below) and a potassium free parameter [K/H] that shares the metallicity of sodium (with index negative 5):

#   Name       value      min     max   step  prior  prior_sigma
retrieval_params =
    T_iso       2850    500.0  3000.0    1.0
    M_planet    0.27      0.1     0.4    1.0  0.266  0.033
    log_p_ref   -1.0     -9.0     2.0    1.0
    [M/H]        0.0     -2.0     3.0    1.0
    [Na/H]       0.0     -2.0     3.0    1.0
    [K/H]        0.0     -2.0     3.0   -5.0

Full examples¶

Here are a couple of examples to reproduce retrieval analyses from peer-reviewed articles: