Retrievals: Emission Spectra Retrieval
Written by Evert Nasedkin. Please cite pRT’s retrieval package (Nasedkin et al. 2024) in addition to pRT (Mollière et al. 2019) if you make use of the retrieval package for your work.
This retrieval is based on the forward model used in Mollière et al. (2020) for HR8799e, and shows a more realistic example of how to set up a retrieval. Emission spectra retrievals, particularly when multiple datasets are included and scattering is taken into account, can take a lot of computational resources to run. Therefore, this notebook outlines the setup for such a retrieval, but we advise only running this retrieval on a cluster (it can take multiple days on hundreds of cores).
For this example, the module chemistry.pre_calculated_chemistry is used to solve for disequilibrium chemistry (actually equilibrium chemistry with a simple quenching treatment), and the chemistry.clouds module is used for condensation. Note that you can import other models from petitRADTRANS/retrieval/models.py, see the API
documentation. This file is included in the pRT package folder. Alternatively you can access them through git (just click the link). Of course you can also write your own model. For this we recommend using existing models as template. Remember to respect the input and output format of the model functions described in “Basic Retrieval Tutorial”
The model here uses a simple adaptive mesh refinement (AMR) algorithm to improve the pressure resolution around the location of the cloud bases.
Getting started
Please make sure to have worked through the“Basic Retrieval Tutorial”before looking at the material below.
In this tutorial, we will outline the process of setting up a RetrievalConfig object, which is the class used to set up a pRT retrieval. The basic process is always to set up the configuration, and then pass it to the Retrieval class to run the retrieval using, for example, pyMultiNest. Like mentioned in the “Basic Retrieval Tutorial” several standard plotting outputs will also be produced by the retrieval class. Most of the classes and functions used in this tutorial have more advanced features than what will be explained here, so it’s highly recommended to take a look at the code and API documentation. There should be enough flexibility built in to cover most typical retrieval studies, but if you have feature requests please get in touch, or open an issue on gitlab.
[1]:
# Let's start by importing everything we need.
import os
# To not have numpy start parallelizing on its own
os.environ["OMP_NUM_THREADS"] = "1"
import numpy as np
import matplotlib.pyplot as plt
import petitRADTRANS as prt
from petitRADTRANS import physical_constants as cst
# Import the class used to set up the retrieval.
from petitRADTRANS.retrieval import Retrieval,RetrievalConfig
# Import Prior functions, if necessary.
from petitRADTRANS.retrieval.utils import gaussian_prior
# Import atmospheric model function
from petitRADTRANS.retrieval.models import emission_model_diseq
[2]:
##########################
# Define the pRT run setup
##########################
retrieval_config = RetrievalConfig(retrieval_name="HR8799e_example", # Give a useful name for your retrieval
run_mode="retrieve", # 'retrieve' to run, or 'evaluate' to make plots
amr=True, # Adaptive mesh refinement, slower if True
scattering_in_emission=True) # Add scattering for emission spectra clouds
For this example we include the GRAVITY data as published in Mollière et al. (2020). To reproduce the published results, please also include the archival SPHERE and GPI data from Zurlo et al. (2015) and Greenbaum et al. (2016).
[ ]:
##################
# Read in Data
##################
# Here we import petitRADTRANS to find the path of the example files on your machine.
# In general this is not required, you just put the files in the folder that you are running
# Your script in, for example.
import petitRADTRANS # need to get the name for the example data
path_to_data = ""
path_to_data = petitRADTRANS.__file__.split('__init__.py')[0] # Default location for the example
retrieval_config.add_data('GRAVITY',
f"{path_to_data}retrieval/examples/emission/observations/HR8799e_Spectra.fits",
data_resolution = 500,
model_resolution = 1000,
model_generating_function = emission_model_diseq)
# Note that in Mollière et al. (2020), additional data sets from SPHERE and GPI were used
Photometric data
If we want to add photometry, we can do that as well! The photometry file should have the format:
# Name, lower wavelength edge [um], upper wavelength edge [um], flux density [W/m2/micron], flux error [W/m2/micron]
You are required to provide a model function for calculating the spectrum, as with spectral data, but also a photometric transformation function, which is used to convert the model spectrum into synthetic photometry. This would typically make use of the transmission function for a particular filter. We recommend the use of the species package (https://species.readthedocs.io/), in particular the SyntheticPhotometry module to provide these functions. If no function is provided, the
RetrievalConfig will attempt to import species to use this module, using the name provided as the filter name.
If you are using transmission spectra, your photometric transformation function should model the difference between the clear and occulted stellar spectrum, returning the difference in (planet radius/stellar radius)^2.
[ ]:
retrieval_config.add_photometry(path_to_data + 'retrieval/examples/emission/observations/hr8799e_photometry.txt',
emission_model_diseq,
model_resolution = 40)
Parameters and Priors
Here we add all of the parameters used in the retrieval.
[ ]:
#################################################
# Add parameters, and priors for free parameters.
#################################################
# This run uses the model of Molliere (2020) for HR8799e
# The lambda function provide uniform priors.
# Distance to the planet in cm
retrieval_config.add_parameter(name = 'D_pl', free = False, value = 41.2925*nc.pc)
# Log of the surface gravity in cgs units.
retrieval_config.add_parameter('log_g',True,
transform_prior_cube_coordinate = \
lambda x : 2.+3.5*x)
# Planet radius in cm
retrieval_config.add_parameter('planet_radius', True, \
transform_prior_cube_coordinate = \
lambda x : ( 0.7+1.3*x)*cst.r_jup_mean)
# Temperature in Kelvin
retrieval_config.add_parameter('T_int', True, \
transform_prior_cube_coordinate = \
lambda x : 300.+2000.*x,
value=0.0)
# Spline temperature structure parameters. T1 < T2 < T3
# As these priors depend on each other, they are implemented in the model function
retrieval_config.add_parameter('T3', True, \
transform_prior_cube_coordinate = \
lambda x : x,
value=0.0)
retrieval_config.add_parameter('T2', True, \
transform_prior_cube_coordinate = \
lambda x : x,
value=0.0)
retrieval_config.add_parameter('T1', True, \
transform_prior_cube_coordinate = \
lambda x : x)
# Optical depth model
# power law index in tau = delta * press_cgs**alpha
retrieval_config.add_parameter('alpha', True, \
transform_prior_cube_coordinate = \
lambda x :1.0+x)
# proportionality factor in tau = delta * press_cgs**alpha
retrieval_config.add_parameter('log_delta', True, \
transform_prior_cube_coordinate = \
lambda x : x)
# Chemistry
# A 'free retrieval' would have each line species as a parameter
# Using a (dis)equilibrium model, we only supply bulk parameters.
# Carbon quench pressure
retrieval_config.add_parameter('log_pquench',True,\
transform_prior_cube_coordinate = \
lambda x : -6.0+9.0*x)
# Metallicity [Fe/H]
retrieval_config.add_parameter('Fe/H',True,\
transform_prior_cube_coordinate = \
lambda x : -1.5+3.0*x)
# C/O ratio
retrieval_config.add_parameter('C/O',True,\
transform_prior_cube_coordinate = \
lambda x : 0.1+1.5*x)
# Clouds
# Based on an Ackermann-Marley (2001) cloud model
# Width of particle size distribution
retrieval_config.add_parameter('sigma_lnorm', True,
transform_prior_cube_coordinate = \
lambda x : 1.05 + 1.95*x)
# Vertical mixing parameters
retrieval_config.add_parameter('log_kzz',True,\
transform_prior_cube_coordinate = \
lambda x : 5.0+8.*x)
# Sedimentation parameter
retrieval_config.add_parameter('fsed',True,\
transform_prior_cube_coordinate = \
lambda x : 1.0 + 10.0*x)
[ ]:
#######################################################
# Define opacity species to be included
#######################################################
retrieval_config.set_rayleigh_species(['H2', 'He'])
retrieval_config.set_continuum_opacities(['H2-H2', 'H2-He'])
retrieval_config.set_line_species([
'H2O__POKAZATEL',
'CO-NatAbund',
'CH4',
'CO2',
'HCN',
'FeH',
'H2S',
'NH3',
'PH3',
'Na',
'K',
'TiO',
'VO',
'SiO'],
eq = True
)
retrieval_config.add_cloud_species('Fe(s)_crystalline__DHS', eq = True, abund_lim = (-3.5,1.0))
retrieval_config.add_cloud_species('MgSiO3(s)_crystalline__DHS', eq = True, abund_lim = (-3.5,1.0))
[ ]:
retrieval = Retrieval(retrieval_config,
output_dir = "",
sample_spec = False)
[ ]:
# Before we run the retrieval, let's set up plotting.
##################################################################
# Define what to put into corner plot if run_mode == 'evaluate'
##################################################################
retrieval_config.parameters['planet_radius'].plot_in_corner = True
retrieval_config.parameters['planet_radius'].corner_label = r'$R_{\rm P}$ ($\rm R_{Jup}$)'
retrieval_config.parameters['planet_radius'].corner_transform = lambda x : x/nc.r_jup_mean
retrieval_config.parameters['log_g'].plot_in_corner = True
retrieval_config.parameters['log_g'].corner_ranges = [2., 5.]
retrieval_config.parameters['log_g'].corner_label = "log g"
retrieval_config.parameters['fsed'].plot_in_corner = True
retrieval_config.parameters['log_kzz'].plot_in_corner = True
retrieval_config.parameters['log_kzz'].corner_label = "log Kzz"
retrieval_config.parameters['C/O'].plot_in_corner = True
retrieval_config.parameters['Fe/H'].plot_in_corner = True
retrieval_config.parameters['log_pquench'].plot_in_corner = True
retrieval_config.parameters['log_pquench'].corner_label = "log pquench"
for spec in retrieval_config.cloud_species:
cname = spec.split('_')[0]
retrieval_config.parameters['log_X_cb_'+cname].plot_in_corner = True
retrieval_config.parameters['log_X_cb_'+cname].corner_label = cname
##################################################################
# Define axis properties of spectral plot if run_mode == 'evaluate'
##################################################################
retrieval_config.plot_kwargs["spec_xlabel"] = 'Wavelength [micron]'
retrieval_config.plot_kwargs["spec_ylabel"] = "Flux [W/m2/micron]"
retrieval_config.plot_kwargs["y_axis_scaling"] = 1.0
retrieval_config.plot_kwargs["xscale"] = 'log'
retrieval_config.plot_kwargs["yscale"] = 'linear'
retrieval_config.plot_kwargs["resolution"] = 100.# Maximum resolution, will rebin the data
retrieval_config.plot_kwargs["nsample"] = 100 # If we want a plot with many sampled spectra
##################################################################
# Define from which observation object to take P-T
# in evaluation mode (if run_mode == 'evaluate'),
# add PT-envelope plotting options
##################################################################
retrieval_config.plot_kwargs["take_PTs_from"] = 'GRAVITY'
retrieval_config.plot_kwargs["temp_limits"] = [150, 3000]
retrieval_config.plot_kwargs["press_limits"] = [1e2, 1e-5]
[ ]:
retrieval = Retrieval(retrieval_config,
output_dir = "",
sample_spec = False)
retrieval.run(n_live_points = 2000,
sampling_efficiency=0.05,
const_efficiency_mode=True,
resume = True)
Once the retrieval is complete, the easiest way to generate standard output plots is to use the plot_all function.
[ ]:
retrieval.plot_all(contribution = True)
Contact
If you need any additional help, don’t hesitate to contact Evert Nasedkin.