Home

Welcome to the iris Wiki

iris is a Python library to generate slab models of molecular emission in disks.

The modules in iris are optimized with jax, and set up to efficiently model the emission of multiple species simultaneously, as well as the effects of radial temperature gradients and Keplerian line broadening.

---

Table of Contents

1. Installation
2. Background
   2.1 Slab Models
   2.2 Keplerian Broadening
3. Crash Course
   3.1 Initial Settings
   3.2 Fetching HITRAN Data
   3.3 Making a Slab Model
   3.4 Compiling the Code
4. Examples
5. FAQ

---

Installation

Install the package via pip:

pip install iris-jwst

iris requires the following packages be installed beforehand: jax, pandas, astropy, astroquery

and Python >= 3.8

Background

Slab Models

With iris, the molecular emission lines are modeled using a slab technique, where the intensities are set by an average excitation temperature $(T_{ex})$, a column density ($N_{mol}$), and an emitting area $(A_{em})$. The population levels are assumed to be in Local Thermodynamic Equilibrium (LTE). That is, the excitation temperature and the gas kinetic temperature are assumed to be the same.

We assume the lines are Gaussians and have an intrinsic full-width-at-half-maximum $(\Delta V)$ which includes turbulent and thermal broadening. Note that $\Delta V$ is a free parameter set by the user.

iris generates an optical depth grid at a very high resolution ($\lambda / \delta \lambda \sim 10^6$) which is then used to calculate the opacity-weighted flux density. For each species, we estimate the total optical depth of lines $i=1,2,...,N$ as

$\tau (\lambda) = \sum_{i}^{N} \tau_{0,i} e^{(-(\lambda - \lambda_0)^2/2\sigma_{\lambda}^2)} $

where $\lambda_0$ is the rest wavelength of line $i$, $\sigma_{\lambda}$ is the line standard deviation $\left( \frac{\Delta V}{2\sqrt{2 \log 2}} \right)$, and the line-center optical depth is

$\tau_0 = \sqrt{\frac{\ln 2}{\pi}} \frac{A_{ul} N_{mol} \lambda_0^3}{4 \pi \Delta V} \left( x_l \frac{g_u}{g_l} - x_u \right)$.

Here $A_{ul}$ is the Einstein coefficient, $x_{l,u}$ the lower- and upper-level populations, and $g_{l,u}$ the lower- and upper-level degeneracies, respectively.

Setting a distance to the source of $d$, we can then compute the flux assuming emission from a Blackbody $B_{\nu}$ at temperature $T = T_{ex}$:

$F(\lambda) = \left(\frac{A_{em}}{d^2}\right) B_{\nu}(T) \left( 1 - e^{-\tau(\lambda)} \right)$.

The rightmost term above sets the final line profile, accounting for saturation effects at the line center:

The flux model is then post-processed by iris to a user-specified resolving power (R) and wavelength grid.

Keplerian Broadening

Optionally, we can go a step further and include the effects of Keplerian broadening on the line profiles. iris does this in 4 steps, which only add a couple milliseconds to the compute time:

1. Make a standard flux model $F_{\lambda}$ with turbulent+thermal broadening only.
2. Generate a flat, Keplerian disk model with inclination $i$, stellar mass $M_{*}$, innermost radius $r_{in}$, and outermost radius $r_{out}$.
3. Finely sample the line-of-sight velocities $v_{LOS} = v_{K} \sin(i) \cos(\phi)$ within the emitting region, where $v_{K} = (G M_{*}/r)^{0.5}$ is the Keplerian velocity and $\phi$ is the azimuthal angle. Make a normalized histogram ($k_{LOS}$) of $v_{LOS}$.
4. Convolve $F_{\lambda}$ with $k_{LOS}$.

The effects will be most noticeable for disks at high inclinations, very small inner radii, and/or with higher stellar masses.

Crash Course

Initial Settings

Whenever we work with iris, first we need to import jax in 64-bit mode. This is very important for the kind of calculations made by iris:

import jax 
jax.config.update("jax_enable_x64", True)

We can then import any other useful packages and iris itself:

import jax.numpy as jnp

import numpy as np
import matplotlib.pyplot as plt

import iris as iris
from iris import setup

Fetching HITRAN Data

The very first time we model a molecule, we need to download the corresponding HITRAN spectroscopic data. iris will do this for you. First we specify a path where we want the data to be downloaded, and then we use setup.setup_linelists to fetch and format the data:

path_to_moldata = './'  # path where we want to save the HITRAN data

# Get the data for H$_2$O and $^{13}$CO
setup.setup_linelists('H2O', 'H2O', 1, path_to_moldata)
setup.setup_linelists('CO2', 'CO2', 1, path_to_moldata)
setup.setup_linelists('13CO', 'CO', 2, path_to_moldata)

Note the signature of this function: the first argument is the isotopolog name, this can be whatever you want to name this species, e.g. '13CO' for $^{13} \text{CO}$. The second argument is the HITRAN main isotopolog name, this is the name for the main isotopolog in HITRAN, e.g. 'CO'. The third argument is the HITRAN isotopolog number, this is the ID number for the particular isotopolog of interest. For $^{13} \text{CO}$, the ID number is 2. You can easily find the ID numbers in the HITRAN website.

Making a Slab Model

1. Define a wavelength grid (in micron) to evaluate opacities. The spacing should correspond to an $R \sim 10^{5-6}$ or whatever is enough to sample the intrinsic line widths $\Delta V$ you want.

fine_wgrid = np.arange(11.5, 18.5, 1e-5)

Here we make a wavelength grid from 11.5 to 18.5 micron, spaced out by $10^{-5}$ micron intervals.

2. Define a wavelength to downsample the final flux model. This will usually be the wavelength grid of your data.

obs_wgrid = np.arange(12.0,18.0,0.002)

Be careful obs_wgrid is within fine_wgrid, otherwise it's impossible to downsample the model.

3. Define the resolving power of the instrument

R = 3200

iris will convolve the flux density model with a Gaussian of equivalent width to the resolving power provided.

4. Initialize the slab object. You only need to do this once, it will take a few seconds to parse the molecular data.

slab = iris.slab(molecules=['CO2', 'H2O'], wlow=11.0, whigh=19.0, path_to_moldata=path_to_moldata)

We have created a slab object to model $\text{CO}_2$ and $\text{H}_2\text{O}$, including data from 11 to 19 micron, reading from the folder defined above. Some important things to consider:

• iris will model the molecules in the order provided.
• The wlow and whigh min. and max. wavelengths should extend beyond the wavelength range spanned by fine_wgrid. Otherwise, you risk missing molecular lines. Still, the code will be faster when using a smaller wavelength range.

5. Set up the wavelength grid in the slab:

slab.setup_grid(fine_wgrid, obs_wgrid, R)

We call setup_grid any time we want to update the wavelength grid.

6. Define your disk parameters.

Distance to source:

distance = 120 # pc

The temperature, column density, emitting area, and $\Delta V$ parameters must each be (M, N) arrays: the model parameters are matrices with M rows and N columns, where M is the number of species and N is the number of components.

# Excitation temperatures for each molecule in K
T_ex =  np.array([np.array([800.0, 400.0, 200.0]),
                  np.array([800.0, 400.0, 200.0])])
# column densities in cm^-2
N_mol = np.array([np.array([1e16, 5e16, 1e16]),
                  np.array([5e18, 5e17, 1e17])])
# emitting areas in au^2
A_au =  np.array([np.array([0.1, 1.0, 10.0]),
                  np.array([0.1, 1.0, 10.0])])

# intrinsic line widths in km/s (line FWHM)
dV =    np.array([np.array([2.0, 2.0, 2.0]),
                  np.array([2.0, 2.0, 2.0])])

Here we are modeling 3 slabs for $\text{CO}_2$ and $\text{H}_2\text{O}$ each. Note that each slab has a different temperature (800, 400, 200 K), column density and area, but all have the same line width of 2 km/s. This is one way we can model radial temperature gradients in disks.

7. Set up the disk object. For now, let's ignore Keplerian broadening.

We call the slab.setup_disk module every time we wish to update a physical parameter.

slab.setup_disk(distance, T_ex=T_ex, N_mol=N_mol, A_au=A_au, dV=dV)

8. Generate the slab model.

We use the slab.simulate module to calculate the slab model.

slab.simulate()

If we wanted to include Keplerian effects, we would define an inclination, stellar mass, and innermost radius, and call slab.simulate_keplerian()

The full flux density model can then be accessed as flux.flux_model. The downsampled and convolved model is slab.downsampled_flux.

plt.figure(figsize=(17,4))

plt.plot(obs_wgrid, slab.downsampled_flux, color='orangered', lw=1)

plt.xlim(12, 18.0)

plt.tick_params(direction='in', top=True, right=True, length=5, labelsize=14)
plt.xlabel('Wavelength ($\mu$m)', fontsize=17)
plt.ylabel('Flux (Jy)', fontsize=17)

Compiling the Code

Thanks to jax, we can greatly speed-up performance by jit-compiling the code. In short, we want to define a function that calls slab.setup_disk, slab.setup_grid, and slab.simulate(). Then we compile using jax.jit():

def compiled_slab(distance, T_ex, N_mol, A_au, dV, fine_wgrid, wavelength, R):

    slab.setup_disk(distance, T_ex, N_mol, A_au, dV) # initialize object
    slab.setup_grid(fine_wgrid, wavelength, R) # set up wavelength and model parameters
    slab.simulate() # make model

    return slab.downsampled_flux, slab.flux_model

'''Just-in-time compilation'''
compiled_slab_jit = jax.jit(compiled_slab)
_ = compiled_slab_jit(distance, T_ex, N_mol, A_au, dV, fine_wgrid, wavelength, R)

Calling compiled_slab will be $\sim 100$ times faster after compiling. Ideally, this is the function you would call inside an inference loop when fitting a slab model to data.

Examples

Please see the Examples page for tutorials.

---

FAQ

1. Do I need a GPU?

Not necessarily. You can install the CPU-only version of jax to use iris. But the code will be substantially slower, since the computations cannot be vectorized that way. It will take about $1000 \times$ longer to generate a slab (~ a few seconds).

2. How do I install jax?

Please see: https://jax.readthedocs.io/en/latest/installation.html and follow the instructions for your particular GPU and operating system.

3. Can I run this on a MacBook Pro M(1-4)?

The short answer is no. iris is not set up to run on a personal laptop, and jax is not yet optimized to use the MacBook GPU. You can still install the CPU-only version of jax (See 1).

4. How much GPU memory do I need?

It depends. For a few molecules and temperature components, 16 GB should be enough. To model temperature gradients with tens of components or full retrievals of many species and isotopologs, you'll need access to a computing cluster with an A100 GPU with 40 GB or more. Services like Google Colab provide these for $1.50 / hr.