Note

Download full notebook here

Getting Set up to Run the Code

In this tutorial you will learn:

1. How to pick condensates to run in a cloud model

2. What is \(f_{sed}\) and what number is right for my model?

3. What are the chemical limitations in the code

4. How to compute initial Mie scattering grid

If you are already familiar with these things, I recommend moving to the next tutorial

from bokeh.io import output_notebook
from bokeh.plotting import show, figure
from bokeh.palettes import Colorblind
output_notebook()
import numpy as np
import pandas as pd
import astropy.units as u

#here is virga
import virga.justdoit as jdi

Loading BokehJS ...

We will be pretty thorough in this first tutorial. In the very first step, we will set up our calculation run. This includes:

1. Defining what gas condensates we want to include in the calculation.

2. Defining the sedimentation efficiency \(f_{sed}\).

3. Defining metallicity and atmospheric mean molecular weight

4. Defining grid of particle radii for the Mie code

Let’s briefly go through these one-by-one.

How to Pick Gas Condensates

We highly recommend looking at your climate profile first to give an informed decision to the cloud model. First check out what current condensates are available

#you can print the available condensate species in the code
jdi.available()

['Al2O3',
 'CH4',
 'CaAl12O19',
 'CaTiO3',
 'Cr',
 'Fe',
 'H2O',
 'KCl',
 'Mg2SiO4',
 'MgSiO3',
 'MnS',
 'NH3',
 'Na2S',
 'SiO2',
 'TiO2',
 'ZnS']

Then you can let virga give you a recommendation with some visualization aid

pressure = np.logspace(-5,3,30) #simple isotherml PT profile (kelvin)
temperature = np.zeros(30)+1300
metallicity = 1 #atmospheric metallicity relative to Solar
mean_molecular_weight = 2.2 # atmospheric mean molecular weight

#get virga recommendation for which gases to run
recommended = jdi.recommend_gas(pressure, temperature, metallicity,mean_molecular_weight,
                #Turn on plotting
                 plot=True)
#print the results
print(recommended)

['Cr', 'Fe', 'Mg2SiO4', 'MgSiO3', 'MnS', 'SiO2']

How to Compute Temperature Condensation Curves

jdi.recommend is comparing your pressure-temperature curve against condensation temperature curves for each species. If you are interested in computing your own you can do so!

cond_t_fig = figure(y_axis_type='log',y_range=[1e2,1e-3],height=400, width=400,
                  y_axis_label='Pressure(bars)',x_axis_label='Temperature (K)',
                   title='Condensation Temperatures')
metallicity = 1
mean_molecular_weight = 2.2
for gas_name, c in zip(['Cr','MgSiO3','MnS'],Colorblind[8]): #case sensitive names
    #grab p,t from eddysed
    p,t = jdi.condensation_t(gas_name, metallicity, mean_molecular_weight)
    #plot it up
    cond_t_fig.line(t,p, legend_label=gas_name,color=c,line_width=3)

show(cond_t_fig)

How to pick \(f_{sed}\)

This parameter will set the vertical extent of our cloud. The general equation is:

condensate concentration (z) \(\propto\) concentration at cloud base \(exp(-f_{sed}z/H)\)

Therefore, high values of \(f_{sed}\) deplete the cloud quickly, low values of \(f_{sed}\) create large vertically thick clouds. For context, it is thought that Hot Jupiters have fractional :math:`f_{sed}` values around 0.1, while Jupiter clouds can be parameterized with much larger values of 3-6.

Here is a little schematic to show the general behavior of \(f_{sed}\)

#km, arbitrary planet scale height
H = 27

#define the cloud base
cloud_base = 0.8 #bars

#and some arbitrary height
z = np.linspace(0,300,100)#km
pressure_grid=100*np.exp(-z/H)#bars

#setup figure
toy_model = figure(x_axis_type='log',x_range=[1e-3,1.1],y_range=[1,1e-2],height=400, width=400,
                  y_axis_label='Pressure(bars)',x_axis_label='Arbitrary Fake Cloud Concentration Units',
                  title='F_sed Toy Model')

#loop through some fseds and plot them up
for fsed,c in zip([10,6,3,1,0.1],['red','orange','green','blue','purple']):
    cloud=np.zeros(100)+1e-13 #start with empty array

    #exponential equation I mentioned above
    fake_cloud = np.exp(-fsed*z[pressure_grid<cloud_base]/H) #fill with exponential fake cloud
    cloud[pressure_grid<cloud_base] = fake_cloud/np.max(fake_cloud) #fill cloud with fake cloud and normalize
    toy_model.line(cloud,pressure_grid,color=c,legend_label='F_sed='+str(fsed))
    toy_model.legend.location='bottom_left'
show(toy_model)

Chemical Limitations of the Code (e.g. M/H)

Although we are trying to make virga as general as possible. There are limitations to the metallicity and mean molecular weights that users enter.

There are two values that are hard-coded for H2 atmospheres

1. The diameter of the atmospheric molecule (We assume d = 2.827e-8 for H2)

2. The depth of the Lennard-Jones potential well for the atmosphere (we use 59.7 K, the H2 value)

You may ask… Seems REAL easy to make these an input.. What gives!?

The reason we fix this is because the chemical concentrations (taken from Morley et al 2012) cannot be trivially scaled and applied to non-H2 very high metallicity atmospheres.

To see what gas properties are currently in the code, you can explore gas_properties

atmo_mean_molecular_weight = 2.2
metallicity = 1 #non-log, relative to solar

#retrieve :
#condensate molecular weight, gas mass mixing ratio, density of the condensate
gas_mw, gas_mmr, rho_p = jdi.gas_properties.H2O(atmo_mean_molecular_weight, metallicity )

Some of our functions have a metallicity dependence (such as ZnS)

#of course there are other gases which do contain metallicity dependencies
#e.g. ZnS
gas_mw, gas_mmr, rho_p = jdi.gas_properties.ZnS(atmo_mean_molecular_weight, 100*metallicity )

Creating the Mie scattering Database

The wavelength dependence of the code relies on computing the scattering and absorbing coefficients from Mie theory.

Our python code uses PyMieScatt, which has extensive documentation on the how the calculations are done assuming that all particles are spherical and homogeneous.

What you need to know is that before running virga you must first create a database of these parameters given a particles complex refractive index (\(m = n+ik\)), which are basically just a condensate species optical properties. These optical properties are highly sought after properties, which are measured in the lab.

Complex refractive indices

\(n\) and \(k\) are a function of wavelength. Mie scattering and absorbing coefficients are a function of wavelength and particle size.

You can learn a lot about a particle itself by exploring the complex refractive indices. Let’s look at both.

virga comes with a folder with files that end in .refrind. We’ve fixed the grid of these refractive indices. However, note that the wavelength grid needs to be fine enough to resolve any spectral features you are interested in.

#directory where you have put `refrind`
refrind_dir = '/data/virga'
wave, n, k  = jdi.get_refrind('Na2S', refrind_dir)

Lets reproduce Figure 2 in Morley et al 2012.

The single most useful thing to remember about refractive indices is: imaginary part (also called the extinction coefficient) relates to the absorbing power of the species.

from bokeh.layouts import row
k_fig = figure(height = 200,width=300 ,y_axis_label = 'Imaginary Part',x_axis_type='log',y_axis_type='log',x_axis_label='Wave (um)')
n_fig = figure(height = 200, width=300,y_axis_label = 'Real Part',x_axis_type='log',x_axis_label='Wave (um)')
n_fig.line(wave,n)
k_fig.line(wave,k)
show(row([n_fig, k_fig]))

Pre-set wavelength Grid

Note that we have supplied \(n\) and \(k\) on a pre-set wavelength grid. Some people have inquired about increasing the resolution of the grid. Doing so requires more robust laboratory data. Ask your local theorist what might be available for your problem.

Running the Mie Database Code

The function below will create a set of files called e.g. ZnS.Mieff in whatever directory you specify. So beware of file overwrites.

Just for example, we will run a single gas with a 10 radius grid. BUT IRL USE AT LEAST THE DEFAULT nradii=40 and for all the gases

import os
out_dir = os.getcwd()

#JUST FOR THE EXAMPLE, let's run one species on a small radii grid
gas_name = 'H2O'
qext, qscat, g_qscat, radii,wave = jdi.calc_mie_db([gas_name],
                                refrind_dir, out_dir, rmin = 1e-5, nradii = 10)

#Your IRL run, which will take some time to compute
#qext, qscat, cos_qscat, radii = ed.calc_mie_db(ed.avaialable(), refrind_dir, out_dir)

Computing H2O
/home/nbatalh1/codes/virga/docs/notebooks/H2O.mieff

Analyzing Mie Parameters

For reflected light calculations, the most common/useful plots to look at are

• Asymmetry (g_qscat/qscat) : will show how forward/backscattering the gas is

• Qscat / Qext (qscat/qext) : will show you how bright the condensate is expected to be

from bokeh.palettes import viridis
from bokeh.models import LogColorMapper, ColorBar, LogTicker

colors = viridis(len(radii))
qfig = figure(width=300, height=300,
              x_axis_type='log',y_axis_label ='Asymmetry',
              x_axis_label='Wavelength(um)',title=gas_name)

wfig = figure(width=400, height=300,
              x_axis_type='log',y_axis_label ='Qscat/Qext',
              x_axis_label='Wavelength(um)',title=gas_name)

for i in range(len(radii)):
    qfig.line(wave, g_qscat[:,i,0]/qscat[:,i,0], color=colors[i])
    wfig.line(wave, qscat[:,i,0]/qext[:,i,0], color=colors[i])

#color bar
color_mapper = LogColorMapper(palette="Viridis256", low=min(radii*1e4), high=max(radii*1e4))
color_bar = ColorBar(color_mapper=color_mapper, ticker=LogTicker(),
                     label_standoff=12, border_line_color=None, location=(0,0),title="Radius (um)")
wfig.add_layout(color_bar, 'right')

show(row(qfig, wfig))

You can also quickly grab your files after you’ve computed them.

qext, qscat, g_qscat, nwave,radii,wave = jdi.get_mie(gas_name,out_dir)

qext.shape

(196, 10)