Including clouds

Powerlaw cloud models

Let’s initialize an atmosphere in the usual way, see “Getting Started”:

[1]:
import numpy as np
from petitRADTRANS import Radtrans

atmosphere = Radtrans(line_species = ['H2O', 'CO_all_iso', 'CH4', 'CO2', 'Na', 'K'], \
      rayleigh_species = ['H2', 'He'], \
      continuum_opacities = ['H2-H2', 'H2-He'], \
      wlen_bords_micron = [0.3, 15])

pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
 Done.

Units in petitRADTRANS: remember that all units in petitRADTRANS are in cgs, except for pressure, which is in bars, and the mean molecular weight (MMW), which is in units of atomic mass units.

Next we set up the atmospheric parameters, with a parameter selection identical to the “Getting Started” case.

[2]:
from petitRADTRANS import nat_cst as nc
R_pl = 1.838*nc.r_jup_mean
gravity = 1e1**2.45
P0 = 0.01

kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ)

mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperature)
mass_fractions['He'] = 0.24 * np.ones_like(temperature)
mass_fractions['H2O'] = 0.001 * np.ones_like(temperature)
mass_fractions['CO_all_iso'] = 0.01 * np.ones_like(temperature)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperature)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Na'] = 0.00001 * np.ones_like(temperature)
mass_fractions['K'] = 0.000001 * np.ones_like(temperature)

MMW = 2.33 * np.ones_like(temperature)

Abundances in petitRADTRANS: remember that abundances in pRT are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). You can convert between mass fractions and VMRs by using

\begin{equation} X_i = \frac{\mu_i}{\mu}n_i, \end{equation}

where \(X_i\) is the mass fraction of species \(i\), \(\mu_i\) the mass of a single molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molecular weight, and \(n_i\) is the VMR of species \(i\).

Now, let’s calculate cloudy spectra! The first available mode is the so-called power law cloud. In this case the following opacity is added to the scattering cross-section:

\begin{equation} \kappa = \kappa_0\left(\frac{\lambda}{\lambda_0}\right)^\gamma \end{equation}

Where \(\kappa_0\) is the opacity in units of cm\(^2\)/g, at \(\lambda_0=0.35 \ {\rm \mu m}\). The power law index \(\gamma\) fixes the wavelength dependence, Rayleigh-like scattering would be obtained for \(\gamma=-4\). Hence \(\kappa_0\) and \(\gamma\) are free parameters.

We now calculate a clear atmosphere, and then four power law cloudy models, with \(\kappa_0 = 0.01\) cm\(^2\)/g, and four different \(\gamma\) values.

Clear model:

[3]:
import pylab as plt
plt.rcParams['figure.figsize'] = (10, 6)

# Clear
atmosphere.calc_transm(temperature, mass_fractions, \
                       gravity, MMW, R_pl=R_pl, P0_bar=P0)
clear = atmosphere.transm_rad/nc.r_jup_mean

First, \(\gamma = -4\) (Rayleigh-like):

[4]:
kappa_zero = 0.01
gamma_scat = -4.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m4 = atmosphere.transm_rad/nc.r_jup_mean

Second, \(\gamma = -2\) (weaker scattering power law):

[5]:
kappa_zero = 0.01
gamma_scat = -2.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m2 = atmosphere.transm_rad/nc.r_jup_mean

Third, \(\gamma = 0\) (flat opacity):

[6]:
kappa_zero = 0.01
gamma_scat = 0.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

m0 = atmosphere.transm_rad/nc.r_jup_mean

Fourth, the exoctic case of, \(\gamma = 1\) (positive opacity slope):

[7]:
kappa_zero = 0.01
gamma_scat = 1.

atmosphere.calc_transm(temperature, mass_fractions, \
                    gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                    kappa_zero = kappa_zero, gamma_scat = gamma_scat)

p1 = atmosphere.transm_rad/nc.r_jup_mean

# Make plot

plt.plot(nc.c/atmosphere.freq/1e-4, clear, label = 'Clear')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m4, \
         label = r'Powerlaw cloud, $\gamma = -4$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m2, \
         label = r'Powerlaw cloud, $\gamma = -2$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         m0, \
         label = r'Powerlaw cloud, $\gamma = 0$')
plt.plot(nc.c/atmosphere.freq/1e-4, \
         p1,\
         label = r'Powerlaw cloud, $\gamma = 1$')
plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
plt.clf()
../../_images/content_notebooks_clouds_19_0.png
<Figure size 720x432 with 0 Axes>

Gray cloud deck and/or scaled Rayhleigh scattering

In addition to power law clouds, petitRADTRANS can also simply add a gray cloud deck, or scale the Rayleigh scattering opacities of the gas by a factor specified by the user. For these examples we will reuse the Radtrans object from the demonstration of the power law cloud above, there’s no need to intialize it again.

Let’s calculate three spectra: - A clear model - A cloudy model with a gray cloud deck at 0.01 bar - A haze-like model, mimicked by scaling the Rayleigh scattering of the gas by a factor 10 - A combination of the two above cloud models: gray cloud deck plus haze

[8]:
import pylab as plt
plt.rcParams['figure.figsize'] = (10, 6)

# Clear
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Clear')

# Gray cloud deck at 0.01 bar
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      Pcloud = 0.01)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Gray cloud deck at 0.01 bar')

# Haze (10 x gas Rayleigh scattering)
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      haze_factor = 10)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Rayleigh haze')

# Haze + cloud deck
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, R_pl=R_pl, P0_bar=P0, \
                      haze_factor = 10, Pcloud = 0.01)
plt.plot(nc.c/atmosphere.freq/1e-4, \
         atmosphere.transm_rad/nc.r_jup_mean, label = 'Rayleigh haze + cloud deck')

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
plt.clf()
../../_images/content_notebooks_clouds_23_0.png
<Figure size 720x432 with 0 Axes>

Condensate clouds from real optical constants

Let’s calculate some spectra using opacities derived from optical constants of various materials. In this example we will start of with forsterite, that is \(\rm Mg_2SiO_4\). For the list of available cloud species see “Available opacity species”.

Choose your favorite particle setup: petitRADTRANS offers multiple opacity versions for every given condensate, assuming either spherical (Mie scattering) or irreagularly shaped (DHS method = distribution of hollow spheres) particles. Moreover, opacities assuming a crystalline and/or amorphous internal structure can be used, where available.

We set up the atmosphere like before, this time loading the opacity of crystalline (c), irregularly shaped particles (DHS method, abbreviated d). The mode identifier for the \(\rm Mg_2SiO_4\) opacity therefore is Mg2SiO4(c)_cd. (c) stands for condensed. These are all possible modes, note that for some species only crystalline or amorphous cross-sections are available.

Opacity code Description
Species(c)_cm crystalline Mie opacity of Species
Species(c)_cd crystalline DHS opacity of Species
Species(c)_am amorphous Mie opacity of Species
Species(c)_ad amorphous DHS opacity of Species

But now, let’s start! We simply give one additional cloud_species list to the Radtrans object, containing the name of the opacity species we want to use.

[9]:
import numpy as np
from petitRADTRANS import Radtrans

atmosphere = Radtrans(line_species = ['H2O', 'CO_all_iso', 'CH4', 'CO2', 'Na', 'K'], \
      cloud_species = ['Mg2SiO4(c)_cd'], \
      rayleigh_species = ['H2', 'He'], \
      continuum_opacities = ['H2-H2', 'H2-He'], \
      wlen_bords_micron = [0.3, 15])

pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
 Done.

Units in petitRADTRANS: remember that all units in petitRADTRANS are in cgs, except for pressure, which is in bars, and the mean molecular weight (MMW), which is in units of atomic mass units.

Next we set up the atmospheric parameters, with a parameter selection almost identical to the “Getting Started” case. We now also add the mass fraction of the cloud species we are interested in. Note that the mass fractions only contains the name Mg2SiO4(c) now, the _cd suffix has to be omitted.

[10]:
from petitRADTRANS import nat_cst as nc
R_pl = 1.838*nc.r_jup_mean
gravity = 1e1**2.45
P0 = 0.01

kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperature = nc.guillot_global(pressures, kappa_IR, gamma, gravity, T_int, T_equ)

mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperature)
mass_fractions['He'] = 0.24 * np.ones_like(temperature)
mass_fractions['H2O'] = 0.001 * np.ones_like(temperature)
mass_fractions['CO_all_iso'] = 0.01 * np.ones_like(temperature)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperature)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Na'] = 0.00001 * np.ones_like(temperature)
mass_fractions['K'] = 0.000001 * np.ones_like(temperature)
mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

MMW = 2.33 * np.ones_like(temperature)

Abundances in petitRADTRANS: remember that abundances in petitCODE are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). You can convert between mass fractions and VMRs by using

\begin{equation} X_i = \frac{\mu_i}{\mu}n_i, \end{equation}

where \(X_i\) is the mass fraction of species \(i\), \(\mu_i\) the mass of a single molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molecular weight, and \(n_i\) is the VMR of species \(i\).

Setting the cloud particle size

We have to define a few additional parameters for the clouds, to start the calculation. We can give the mean particle size and width of the log-normal particle distribution:

[11]:
radius = {}
radius['Mg2SiO4(c)'] = 0.00005*np.ones_like(temperature) # I.e. a 0.5-micron particle size (0.00005 cm)

sigma_lnorm = 1.05

Now, let’s calculate a clear and cloudy spectrum to compare:

[12]:
atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       radius = radius, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'cloudy', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       radius = radius, sigma_lnorm = sigma_lnorm)
plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'clear', zorder = 1)

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc='best')
plt.show()
plt.clf()
../../_images/content_notebooks_clouds_39_0.png
<Figure size 720x432 with 0 Axes>

Here one sees that there is a lot of additional absorption in the optical and near-IR. Also note the silicate (Si-O stretching mode) becoming visible at 10 micron.

Calculating the particle size

Alternatively, one can specify an eddy diffusion parameter (\(K_{zz}\), with units of cm\(^2\)/s), which expresses how strongly atmospheric mixing is. The cloud particle size is then controlled by specifying \(K_{zz}\), and the unitless settling parameter \(f_{\rm sed}\), which expresses the particles’ mass averaged settling velocity, when compared to the local atmospheric mixing speed. Also, a width for the log-normal particle size distribution must be specified:

[13]:
Kzz = np.ones_like(temperature)*1e1**7.5
fsed = 2.
sigma_lnorm = 1.05

Again, let’s calculate a clear and cloudy spectrum to compare:

[14]:
mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'cloudy', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_transm(temperature, mass_fractions, gravity, MMW, \
                       R_pl=R_pl, P0_bar=P0, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.transm_rad/nc.r_jup_mean, label = 'clear', zorder = 1)

plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc='best')
plt.show()
plt.clf()
../../_images/content_notebooks_clouds_45_0.png
<Figure size 720x432 with 0 Axes>

The resulting mean particle sizes can be accessed like this:

[15]:
plt.yscale('log')
plt.xscale('log')

plt.ylim([1e2,1e-6])

plt.ylabel('P (bar)')
plt.xlabel('Average particle size (microns)')

plt.plot(atmosphere.r_g[:,atmosphere.cloud_species.index('Mg2SiO4(c)')]/1e-4, pressures)
plt.show()
plt.clf()
../../_images/content_notebooks_clouds_47_0.png
<Figure size 720x432 with 0 Axes>
Minimum particle size: the opacity of particles smaller than 1 nm is set to zero.

Emission spectra

We can use the Radtrans object created above to calculate emission spectra as well. Here we will also generate an instance of pRT with scattering turned on, (see “Scattering for Emission Spectra” for an example on how to do this in detail).

[16]:
mass_fractions['Mg2SiO4(c)'] = 0.0000005 * np.ones_like(temperature)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, \
         color = 'black', label = 'cloudy, no scattering', zorder = 1)

# Load scattering version of pRT
atmosphere = Radtrans(line_species = ['H2O', 'CO_all_iso', 'CH4', 'CO2', 'Na', 'K'], \
      cloud_species = ['Mg2SiO4(c)_cd'], \
      rayleigh_species = ['H2', 'He'], \
      continuum_opacities = ['H2-H2', 'H2-He'], \
      wlen_bords_micron = [0.3, 15], \
      do_scat_emis = True)
pressures = np.logspace(-6, 2, 100)
atmosphere.setup_opa_structure(pressures)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm, \
                       add_cloud_scat_as_abs = True)
plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, \
         label = 'cloudy, including scattering', zorder = 2)

mass_fractions['Mg2SiO4(c)'] = np.zeros_like(temperature)

atmosphere.calc_flux(temperature, mass_fractions, gravity, MMW, \
                       Kzz = Kzz, fsed=fsed, sigma_lnorm = sigma_lnorm)

plt.plot(nc.c/atmosphere.freq/1e-4, atmosphere.flux/1e-6, '-', \
         color = 'red', label = 'clear', zorder = 0)

plt.legend(loc='best')
plt.xscale('log')
plt.xlabel('Wavelength (microns)')
plt.ylabel(r'Planet flux $F_\nu$ (10$^{-6}$ erg cm$^{-2}$ s$^{-1}$ Hz$^{-1}$)')
plt.show()
plt.clf()

  Read CIA opacities for H2-H2...
  Read CIA opacities for H2-He...
 Done.

../../_images/content_notebooks_clouds_51_1.png
<Figure size 720x432 with 0 Axes>

Here we plotted the clear spectrum, neglecting the cloud opacity, the cloudy spectrum only considering the absorption of the cloud particles, and the cloudy spectrum with scattering, using the scattering version of pRT.

Scattering and petitRADTRANS: remember that scattering is included for emission spectra in petitRADTRANS only if requested specifically when generating the Radtrans object, as it increases the runtime (see “Scattering for Emission Spectra” for an example how to do this).