Note

Download full notebook here

Diving into code approximations

Like any code, there are several approximations that go into computing intesity from various phase functions. In reflected light models, some of these approximations drastically change the output spectra.

Here we go through each approximation to show what the effects are of each. Also, to show how to toggle between various inputs.

import warnings
warnings.filterwarnings('ignore')
#picaso
from picaso import justdoit as jdi
from picaso import justplotit as jpi
from bokeh.plotting import show, figure
from bokeh.layouts import column
from bokeh.palettes import Colorblind8
from bokeh.io import output_notebook
import astropy.units as u
output_notebook()

Loading BokehJS ...

Using the approx key

opa = jdi.opannection()
cloud_free = jdi.inputs()

Notice that all the approx keys have predefined inputs. These are our recommendations for how to run the code. But, users should always be weary of these and test their sensitivity to your results.

What approximations exist?

We can easily print out all our options for what approximations exist in the code

print('Options for Direct Scattring Phase: ', jdi.single_phase_options())

print('Options for Multiple Scattring Phase: ', jdi.multi_phase_options())

print('Options for Raman Scattring: ', jdi.raman_options())

Can also set functional form of forward/back scattering in approx['TTHG_params']
Options for Direct Scattring Phase:  ['cahoy', 'OTHG', 'TTHG', 'TTHG_ray']
Can also set delta_eddington=True/False in approx['delta_eddington']
Options for Multiple Scattring Phase:  ['N=2', 'N=1']
Options for Raman Scattring:  ['oklopcic', 'pollack', 'none']

Set inputs normally

cloud_free.phase_angle(0) #phase in radians

cloud_free.gravity(gravity=25, gravity_unit=u.Unit('m/(s**2)' ))

#set star
cloud_free.star(opa, 6000, 0.0122, 4.437)

#set atmosphere comp and temp
cloud_free.atmosphere(filename=jdi.jupiter_pt(),delim_whitespace=True)

#make a copy to have a separate cloud input dict
from copy import deepcopy
cloudy=deepcopy(cloud_free)
cloudy.clouds( filename=jdi.jupiter_cld(),delim_whitespace=True)

Direct Scattering Approximation

The derivation documentation has a full description of these direct scattering approximations. Briefly I’ll describe them here:

At the center of each is the One Term Henyey-Greenstein Phase Function (OTHG) and the Two Term HG Phase Function (TTHG).

We also know that planet atmospheres have high degrees of Rayleigh scattering. Cahoy+2010 developed a methodology for incorporating Rayleigh into the direct scattering component.

A more rrobust way of dealing with Rayleigh is to directly fold it’s phase function into the TTHG phase function (TTHG_Ray).

We’ll run each case with and without a cloud so you can see what happens in both regimes

#let's make two different figures for this
fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)
fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)

#define our options
options = jdi.single_phase_options()
colors = Colorblind8[0:len(options)]

#loop through all approximations
for approx, c in zip(options, colors):
    #set approximations
    cloud_free.approx(single_phase = approx)
    cloudy.approx(single_phase = approx)
    df = cloud_free.spectrum(opa)
    wno_nc, alb_nc = df['wavenumber'] , df['albedo']
    df = cloudy.spectrum(opa)
    wno_c, alb_c = df['wavenumber'] , df['albedo']
    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=approx, color=c, line_width=3)
    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)
jpi.plot_format(fig_cloudy)
jpi.plot_format(fig_no_cloud)
show(column(fig_no_cloud,fig_cloudy ))

Can also set functional form of forward/back scattering in approx['TTHG_params']

Multiple Scattering Approximations

Again, derivation documentation has a full description of these multiple scattering approximations.

To complete the multiple scattering integration over all diffuse angles, we have to use some mathematical tricks. Legendre Polynomials are often used to complete this integration to varying degrees. For Solar System/Exoplanet papers, we often stop the expansion at either N=1 or N=2. Our standard will be to run with N=2, but below we show how to run each.

fig = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300,
             x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)

options = jdi.multi_phase_options()
colors = Colorblind8[0:len(options)*2]

for approx, c1,c2 in zip(options, colors[0:2], colors[2:]):
    cloud_free.approx(multi_phase= approx)
    cloudy.approx(multi_phase = approx)
    df = cloud_free.spectrum(opa)
    wno_nc, alb_nc = df['wavenumber'] , df['albedo']
    df = cloudy.spectrum(opa)
    wno_c, alb_c = df['wavenumber'] , df['albedo']
    fig.line(1e4/wno_nc, alb_nc, color=c1, line_width=3)
    fig.line(1e4/wno_c, alb_c,  color=c2, line_width=3)
jpi.plot_format(fig)
show(fig)

Can also set delta_eddington=True/False in approx['delta_eddington']

Raman Scattering Approximations

We all know the importance of Rayleigh scattering in planetary atmospheres. Raman scattering also has important implications for our spectra (these features have been observed in Solar System planets). In particular, at short wavelengths, Raman scattering imprints molecular features from the star on the planetary spectrum.

The most complete analysis of all Raman approximations is in Sromosvky+2005. From these, we use the Pollack Approximation that was used in Cahoy+2010 and others.

We include the original Pollack methodology, but also include a modified version with Oklopcic et al 2018 cross sections and updated methodology to include effects of stellar spectrum.

fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)
fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)

options = jdi.raman_options()
colors = Colorblind8[0:len(options)]

for approx, c in zip(options, colors):
    cloud_free.approx(raman = approx)
    cloudy.approx(raman = approx)
    df = cloud_free.spectrum(opa)
    wno_nc, alb_nc = df['wavenumber'] , df['albedo']
    df = cloudy.spectrum(opa)
    wno_c, alb_c = df['wavenumber'] , df['albedo']
    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=approx, color=c, line_width=3)
    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)

jpi.plot_format(fig_cloudy)
jpi.plot_format(fig_no_cloud)
show(column(fig_no_cloud,fig_cloudy ))

The Effect of Stellar Spectrum

With the updated Raman scattering approximation, you will notice imprints of the stellar spectrum in the planet reflected light spectrum. Take a look below.

fig_cloudy = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                    ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)
fig_no_cloud = figure(x_range=[0.3,1],y_range=[0,1.2], width=500, height=300
                      ,x_axis_label='Wavelength [μm]', y_axis_label='Geometric Albedo',)

#lets play around with different stellar spectra Teff
stellar_teff = [6000,4000,3500]

colors = Colorblind8[0:len(options)]

cloud_free.approx(raman = 'oklopcic')
cloudy.approx(raman = 'oklopcic')

for approx, c in zip(stellar_teff, colors):
    cloud_free.star(opa, approx, 0.0122, 4.437)
    cloudy.star(opa, approx ,0.0122, 4.437)
    df = cloud_free.spectrum(opa)
    wno_nc, alb_nc = df['wavenumber'] , df['albedo']
    df = cloudy.spectrum(opa)
    wno_c, alb_c = df['wavenumber'] , df['albedo']
    fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label=str(approx), color=c, line_width=3)
    fig_cloudy.line(1e4/wno_c, alb_c,  color=c, line_width=3)

cloud_free.approx(raman = 'pollack')
cloudy.approx(raman = 'pollack')

df = cloud_free.spectrum(opa)
wno_nc, alb_nc = df['wavenumber'] , df['albedo']
df = cloudy.spectrum(opa)
wno_c, alb_c = df['wavenumber'] , df['albedo']
fig_no_cloud.line(1e4/wno_nc, alb_nc, legend_label='Pollack', color='black', line_width=2, line_dash='dashed')
fig_cloudy.line(1e4/wno_c, alb_c,  color='black', line_width=2, line_dash='dashed')


jpi.plot_format(fig_cloudy)
jpi.plot_format(fig_no_cloud)
show(column(fig_no_cloud,fig_cloudy ))