Including clouds
Powerlaw cloud models
Let’s initialize an atmosphere in the usual way, see “Getting Started”:
[1]:
import numpy as np
import matplotlib.pyplot as plt
from petitRADTRANS.radtrans import Radtrans
from petitRADTRANS import physical_constants as cst
from petitRADTRANS.physics import temperature_profile_function_guillot_global
atmosphere = Radtrans(pressures = np.logspace(-6,2,100),
line_species = ['H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'],
rayleigh_species = ['H2', 'He'],
gas_continuum_contributors = ['H2-H2', 'H2-He'],
wavelength_boundaries = [0.3, 15])
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--H2/H2--H2-NatAbund/H2--H2-NatAbund__BoRi.R831_0.6-250mu.ciatable.petitRADTRANS.h5'... Done.
Loading CIA opacities for H2-He from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--He/H2--He-NatAbund/H2--He-NatAbund__BoRi.DeltaWavenumber2_0.5-500mu.ciatable.petitRADTRANS.h5'... Done.
Successfully loaded all CIA opacities
Successfully loaded all opacities
Units in petitRADTRANS: all units inside petitRADTRANS are in cgs. However, when interfacing with the code, you are expected to provide pressures in bars (more intuitive), and the mean molecular mass in atomic mass units. They will be converted to cgs units within the code.
Next we set up the atmospheric parameters, with a parameter selection identical to the “Getting Started” case.
[2]:
planet_radius = 1.0*cst.r_jup_mean
reference_gravity = 10**3.5
reference_pressure = 0.01
pressures_bar = atmosphere.pressures*1e-6 # cgs to bar
kappa_IR = 0.01
gamma = 0.4
T_int = 200.
T_equ = 1500.
temperatures = temperature_profile_function_guillot_global(pressures_bar, kappa_IR, gamma, reference_gravity, T_int, T_equ)
mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperatures)
mass_fractions['He'] = 0.24 * np.ones_like(temperatures)
mass_fractions['H2O'] = 0.001 * np.ones_like(temperatures)
mass_fractions['CO-NatAbund'] = 0.01 * np.ones_like(temperatures)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperatures)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperatures)
mass_fractions['Na'] = 0.00001 * np.ones_like(temperatures)
mass_fractions['K'] = 0.000001 * np.ones_like(temperatures)
mean_molar_masses = 2.33 * np.ones_like(temperatures)
Abundances in petitRADTRANS: remember that abundances in pRT are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). One 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 molar mass, and \(n_i\) is the VMR of species \(i\). This is implemented in petitRADTRANS.chemistry.utils.mass_fractions2volume_mixing_ratios() and petitRADTRANS.chemistry.utils.volume_mixing_ratios2mass_fractions().
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]:
# Clear
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure)
clear = transit_radii/cst.r_jup_mean
First, \(\gamma = -4\) (Rayleigh-like):
[4]:
power_law_opacity_350nm = 0.01
power_law_opacity_coefficient = -4.
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
power_law_opacity_350nm=power_law_opacity_350nm,
power_law_opacity_coefficient=power_law_opacity_coefficient)
m4 = transit_radii/cst.r_jup_mean
Second, \(\gamma = -2\) (weaker scattering power law):
[5]:
power_law_opacity_350nm = 0.01
power_law_opacity_coefficient = -2.
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
power_law_opacity_350nm=power_law_opacity_350nm,
power_law_opacity_coefficient=power_law_opacity_coefficient)
m2 = transit_radii/cst.r_jup_mean
Third, \(\gamma = 0\) (flat opacity):
[6]:
power_law_opacity_350nm = 0.01
power_law_opacity_coefficient = 0.
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
power_law_opacity_350nm=power_law_opacity_350nm,
power_law_opacity_coefficient=power_law_opacity_coefficient)
m0 = transit_radii/cst.r_jup_mean
Fourth, the exoctic case of, \(\gamma = 1\) (positive opacity slope):
[7]:
power_law_opacity_350nm = 0.01
power_law_opacity_coefficient = 1.
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
power_law_opacity_350nm=power_law_opacity_350nm,
power_law_opacity_coefficient=power_law_opacity_coefficient)
p1 = transit_radii/cst.r_jup_mean
# Make plot
fig, ax = plt.subplots(figsize = (10,6))
ax.plot(wavelengths*1e4, clear, label = 'Clear')
ax.plot(wavelengths*1e4,
m4,
label = r'Powerlaw cloud, $\gamma = -4$')
ax.plot(wavelengths*1e4,
m2,
label = r'Powerlaw cloud, $\gamma = -2$')
ax.plot(wavelengths*1e4,
m0,
label = r'Powerlaw cloud, $\gamma = 0$')
ax.plot(wavelengths*1e4,
p1,
label = r'Powerlaw cloud, $\gamma = 1$')
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
ax.legend(loc = 'best')
plt.show()
Gray cloud deck and/or scaled Rayleigh 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
Do not use this gray cloud deck implementation if you calculate pRT emission spectra in scattering mode. We increase the opacity to a VERY high value, :math:`10^{99}` cm^2/g, at pressures > cloud top pressure. This can mess with the matrix inversion method used during the flux calculation. Please usethe treatment explained hereto define your own custom gray opacity in this case.
[8]:
fig, ax = plt.subplots(figsize = (10,6))
# Clear
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure)
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = 'Clear')
# Gray cloud deck at 0.01 bar
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
opaque_cloud_top_pressure = 0.01)
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = 'Gray cloud deck at 0.01 bar')
# Haze (10 x gas Rayleigh scattering)
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
haze_factor = 10.)
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = 'Rayleigh haze')
# Haze + cloud deck
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
opaque_cloud_top_pressure = 0.01,
haze_factor = 10.)
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = 'Rayleigh haze + cloud deck')
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
ax.legend(loc = 'best')
Defining arbitrary opacity functions (can be used to parameterize clouds)
What happens if you have a favorite cloud parameterization, but its implementation is not offered in pRT? Fear not, there is an option for you! When running the method calculate_transit_radii() (or calculate_flux()) you can hand it any arbitrary function that returns a scattering or absorption opacity as a function of wavelength and pressure to pRT. This is then going to be used in that particular calculate_transit_radii() or calculate_flux() run. An example can be found below,
where we implement the following cloud model:
and \(0\) if \(P\geq P_{\rm base}\). Here, \(\lambda\) and \(P\) are the wavelength in micron and pressure in bar, respectively, and \(\kappa_0\), \(\lambda_0\), \(p\), \(P_{\rm base}\) and \(f_{\rm sed}\) are free parameters, thought to parameterize the opacity (in \(\rm cm^2/g\)) at the reference wavelength \(\lambda_0\), the power law dependence of the opacity with wavelengh at \(\lambda \gg \lambda_0\), the cloud base pressure in bar, and the cloud scale height descrease factor, respectively.
Here an example implementation of this cloud model:
[9]:
def cloud_opas(kappa_0,
lambda_0,
p,
P_base,
f_sed):
# in micron, in bar
def give_opacity(wlen, press):
retVal = np.zeros((len(wlen), len(press)))
for i_p in range(len(press)):
if press[i_p] < P_base:
retVal[:, i_p] = kappa_0 / (1. + (wlen/lambda_0)**p) * (press[i_p]/P_base)**f_sed
else:
retVal[:, i_p] = 0.
return retVal
return give_opacity
We can now hand this function to calculate_transit_radii() or (calculate_flux()). If we want to use it as an absorption opacity, we have to use the additional_absorption_opacities_function keyword. If we want to use it as a scattering opacity, we have to use the additional_scattering_opacities_function keyword. Note that in the transmission spectrum case the results will not differ between the two options, because both scattering and absorption cross-sections simply contribute to
the total extinction opacity.
Note that the function you hand to pRT has to take in wavelengths in microns and pressure in bar as input, and return the opacities as a matrix of shape number of wavelengths :math:`times` number of pressure layers, identical to how it is implemented above.
Below some examples:
[10]:
fig, ax = plt.subplots(figsize = (10,6))
# Clear
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure)
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = 'Clear')
# External cloud absorption function, test 1
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
additional_absorption_opacities_function=\
cloud_opas(1e-2, 1., 1., 0.01, 2))
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = r'User specified cloud function ' \
r'$(\kappa_0=0.01 \ {\rm cm^2/g}$, $\lambda_0=1 \ {\rm \mu m}$, $p=1$, ' \
r'$P_{\rm base} = 0.01 \ {\rm bar}$ and $f_{\rm sed}=2$)')
# External cloud absorption function, test 2
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
additional_absorption_opacities_function=\
cloud_opas(1e1, 0.5, 4, 0.001, 1))
ax.plot(wavelengths*1e4,
transit_radii/cst.r_jup_mean,
label = r'User specified cloud function ' \
r'$(\kappa_0=10 \ {\rm cm^2/g}$, $\lambda_0=0.5 \ {\rm \mu m}$, $p=4$, ' \
r'$P_{\rm base} = 0.001 \ {\rm bar}$ and $f_{\rm sed}=1$)')
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
plt.legend(loc = 'best')
plt.show()
If you use this option in correlated-k (``c-k``) mode, please only parameterize opacities that vary slowly enough with wavelength. Molecular and atomic line opacities must not be added in this way, for example, for these the correlated-k opacity treatment is required. Clouds (even if they have an absorption feature) and other pseudo-continuum opacities such as collision induced absorption (CIA) could be handled well with the treatment above, however.
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 take 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 solid (s), crystalline, irregularly shaped particles (DHS method). The mode identifier for the \(\rm Mg_2SiO_4\) opacity therefore is 'Mg2SiO4(s)_crystalline__DHS'. Note that for some species only crystalline or amorphous cross-sections are available.
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.
[11]:
atmosphere = Radtrans(pressures = np.logspace(-6,2,100),
line_species = ['H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'],
rayleigh_species = ['H2', 'He'],
gas_continuum_contributors = ['H2-H2', 'H2-He'],
cloud_species = ['Mg2SiO4(s)_crystalline__DHS'],
wavelength_boundaries = [0.3, 15])
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--H2/H2--H2-NatAbund/H2--H2-NatAbund__BoRi.R831_0.6-250mu.ciatable.petitRADTRANS.h5'... Done.
Loading CIA opacities for H2-He from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--He/H2--He-NatAbund/H2--He-NatAbund__BoRi.DeltaWavenumber2_0.5-500mu.ciatable.petitRADTRANS.h5'... Done.
Successfully loaded all CIA opacities
Loading opacities of cloud species 'Mg2SiO4(s)_crystalline__DHS' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/clouds/Mg2SiO4(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000__DHS.R39_0.1-250mu.cotable.petitRADTRANS.h5'(crystalline_000, using DHS scattering)... Done.
Successfully loaded all clouds opacities
Successfully loaded all opacities
Units in petitRADTRANS: remember that all units inside petitRADTRANS are in cgs. However, when interfacing with the code, you are expected to provide pressures in bars (more intuitive). Pressures will be converted to cgs units within the code.
Next we will setup the amtosphere, with a parameter selection almost identical to the “Getting Started” case. We also add the mass fraction of the cloud species we are interested in now.
[12]:
mass_fractions = {}
mass_fractions['H2'] = 0.74 * np.ones_like(temperatures)
mass_fractions['He'] = 0.24 * np.ones_like(temperatures)
mass_fractions['H2O'] = 0.001 * np.ones_like(temperatures)
mass_fractions['CO-NatAbund'] = 0.01 * np.ones_like(temperatures)
mass_fractions['CO2'] = 0.00001 * np.ones_like(temperatures)
mass_fractions['CH4'] = 0.000001 * np.ones_like(temperatures)
mass_fractions['Na'] = 0.00001 * np.ones_like(temperatures)
mass_fractions['K'] = 0.000001 * np.ones_like(temperatures)
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = 0.0000005 * np.ones_like(temperatures)
Abundances in petitRADTRANS: remember that abundances in pRT are in units of mass fractions, not number fractions (aka volume mixing ratio, VMR). One 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 molar mass of a molecule/atom/ion/… of species \(i\), \(\mu\) is the atmospheric mean molar mass, and \(n_i\) is the VMR of species \(i\). This is implemented in petitRADTRANS.chemistry.utils.mass_fractions2volume_mixing_ratios() and petitRADTRANS.chemistry.utils.volume_mixing_ratios2mass_fractions().
Setting the cloud particle size
We have to define a few additional parameters for the clouds, to start the calculation. Since we here assume a log-normal paricle size distribution, we give the mean particle size and width of the distribution:
[13]:
cloud_particles_mean_radii = {}
cloud_particles_mean_radii['Mg2SiO4(s)_crystalline__DHS'] = 0.00005*np.ones_like(temperatures)
# I.e. a 0.5-micron particle size (0.00005 cm)
# Here it is vertically constant, but it can vary!
cloud_particle_radius_distribution_std = 1.05
# a value of 1.0 would be a delta function, so we assume a very narrow distribtion here.
Now, let’s calculate a clear and cloudy spectrum to compare:
[14]:
fig, ax = plt.subplots(figsize = (10,6))
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
cloud_particles_mean_radii=cloud_particles_mean_radii,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'cloudy', zorder = 2)
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = np.zeros_like(temperatures)
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
cloud_particles_mean_radii=cloud_particles_mean_radii,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'clear', zorder = 1)
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
ax.legend(loc='best')
plt.show()
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
Below we will not predefine the particle size. Instead we will calculate it with the “Eddysed” approach presented in Ackerman & Marley (2001). For this, one specifies an eddy diffusion coefficient (\(K_{zz}\), with units of cm\(^2\)/s), which describes the strength of vertical atmospheric mixing. The cloud particle size is then controlled by \(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. Here we use an atmospheric structure almost identical to the one above, with the difference that we also add amorphous and spherical iron cloud particles ('Fe(s)_amorphous__Mie').
[15]:
atmosphere = Radtrans(pressures = np.logspace(-6,2,100),
line_species = ['H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'],
rayleigh_species = ['H2', 'He'],
gas_continuum_contributors = ['H2-H2', 'H2-He'],
cloud_species = ['Mg2SiO4(s)_crystalline__DHS',
'Fe(s)_amorphous__Mie'],
wavelength_boundaries = [0.3, 15])
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--H2/H2--H2-NatAbund/H2--H2-NatAbund__BoRi.R831_0.6-250mu.ciatable.petitRADTRANS.h5'... Done.
Loading CIA opacities for H2-He from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--He/H2--He-NatAbund/H2--He-NatAbund__BoRi.DeltaWavenumber2_0.5-500mu.ciatable.petitRADTRANS.h5'... Done.
Successfully loaded all CIA opacities
Loading opacities of cloud species 'Mg2SiO4(s)_crystalline__DHS' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/clouds/Mg2SiO4(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000__DHS.R39_0.1-250mu.cotable.petitRADTRANS.h5'(crystalline_000, using DHS scattering)... Done.
Loading opacities of cloud species 'Fe(s)_amorphous__Mie' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/clouds/Fe(s)_amorphous/Fe-NatAbund(s)_amorphous/Fe-NatAbund(s)_amorphous__Mie.R39_0.1-250mu.cotable.petitRADTRANS.h5'(amorphous, using Mie scattering)... Done.
Successfully loaded all clouds opacities
Successfully loaded all opacities
Here we define the EddySed parameters. Below we use different \(f_{\rm sed}\) for the two different cloud species:
cloud_f_sed = {}
cloud_f_sed['Mg2SiO4(s)_crystalline__DHS'] = 2.
cloud_f_sed['Fe(s)_amorphous__Mie'] = 3.
If all the clouds should use the same \(f_{\rm sed}\) we can just set:
cloud_f_sed = 2.
[16]:
eddy_diffusion_coefficients = np.ones_like(temperatures)*1e1**7.5
cloud_f_sed = {}
cloud_f_sed['Mg2SiO4(s)_crystalline__DHS'] = 2.
cloud_f_sed['Fe(s)_amorphous__Mie'] = 3.
# We still assume a log-normal particle size distribution:
cloud_particle_radius_distribution_std = 1.05
Again, let’s calculate a clear and cloudy spectrum to compare:
[17]:
fig, ax = plt.subplots(figsize = (10,6))
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = 0.0000005 * np.ones_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = 0.0000005 * np.ones_like(temperatures)
wavelengths, transit_radii, additional_output = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'cloudy', zorder = 2)
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = np.zeros_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = np.zeros_like(temperatures)
wavelengths, transit_radii, _ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'clear', zorder = 1)
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
ax.legend(loc='best')
plt.show()
The resulting mean particle sizes can be accessed like this:
[18]:
fig, ax = plt.subplots(figsize = (10,6))
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylim([1e2,1e-6])
ax.set_ylabel('P (bar)')
ax.set_xlabel('Average particle size (microns)')
for i_cloud in range(len(atmosphere.cloud_species)):
ax.plot(additional_output['cloud_particles_mean_radii'][:,i_cloud]/1e-4,
atmosphere.pressures/1e6,
label = atmosphere.cloud_species[i_cloud])
plt.legend()
plt.show()
Particle sizes: the opacity of particles smaller than 1 nm and larger than 10 cm is set to zero.
Hansen Distribution Clouds
While the EddySed model used here typically uses a log-normal distribution for the cloud particle radii, it is also possible to use the Hansen (1971) distribution. To do this, you’ll need to set the cloud_particles_radius_distribution parameter in calculate_flux() or calculate_transmission() to 'hansen', and set the effective width of the distribution with the cloud_hansen_b parameter. Alternately, you
can specify both the distribution width and center using the cloud_hansen_b and cloud_hansen_a arguments. Note that the value of cloud_hansen_b will be different from that of cloud_particle_radius_distribution_std to give a similar distribution, because the effective width is weighted by the particle area. See Hansen (1971) for more details.
Here, we’ll calculate the mean particle size as a function of altitude, so cloud_hansen_a is not used. We will also look at the particle sizes again.
[19]:
fig, ax = plt.subplots(figsize = (10,6))
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = 0.0000005 * np.ones_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = 0.0000005 * np.ones_like(temperatures)
eddy_diffusion_coefficients = np.ones_like(temperatures)*1e1**7.5
cloud_f_sed = {}
cloud_f_sed['Mg2SiO4(s)_crystalline__DHS'] = 2.
cloud_f_sed['Fe(s)_amorphous__Mie'] = 3.
cloud_particle_radius_distribution_std = 1.05
cloud_hansen_b = 0.01
wavelengths, transit_radii, additional_output = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std,
cloud_particles_radius_distribution = 'hansen',
cloud_hansen_b = cloud_hansen_b)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'cloudy', zorder = 2)
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = np.zeros_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = np.zeros_like(temperatures)
wavelengths, transit_radii, __ = atmosphere.calculate_transit_radii(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
planet_radius=planet_radius,
reference_pressure=reference_pressure,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std,
cloud_particles_radius_distribution = 'hansen',
cloud_hansen_b = cloud_hansen_b)
ax.plot(wavelengths, transit_radii/cst.r_jup_mean, label = 'clear', zorder = 1)
ax.set_xscale('log')
ax.set_xlabel('Wavelength (microns)')
ax.set_ylabel(r'Transit radius ($\rm R_{Jup}$)')
ax.legend(loc='best')
plt.show()
fig, ax = plt.subplots(figsize = (10,6))
ax.set_yscale('log')
ax.set_xscale('log')
ax.set_ylim([1e2,1e-6])
ax.set_ylabel('P (bar)')
ax.set_xlabel('Average particle size (microns)')
for i_cloud in range(len(atmosphere.cloud_species)):
ax.plot(additional_output['cloud_particles_mean_radii'][:,i_cloud]/1e-4,
atmosphere.pressures/1e6,
label = atmosphere.cloud_species[i_cloud])
plt.legend()
plt.show()
Particle sizes: the opacity of particles smaller than 1 nm and larger than 10 cm 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 with scattering_in_emission = True, (see “Scattering for Emission Spectra” for an example on how to do this in detail).
[20]:
fig, ax = plt.subplots(figsize = (10,6))
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = np.zeros_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = np.zeros_like(temperatures)
wavelengths, flux, _ = atmosphere.calculate_flux(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths*1e4, flux, label = 'clear, no scattering')
mass_fractions['Mg2SiO4(s)_crystalline__DHS'] = 0.0000005 * np.ones_like(temperatures)
mass_fractions['Fe(s)_amorphous__Mie'] = 0.0000005 * np.ones_like(temperatures)
wavelengths, flux, _ = atmosphere.calculate_flux(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths*1e4, flux, label = 'cloudy, no scattering')
atmosphere_scat = Radtrans(pressures = np.logspace(-6,2,100),
line_species = ['H2O',
'CO-NatAbund',
'CH4',
'CO2',
'Na',
'K'],
rayleigh_species = ['H2', 'He'],
gas_continuum_contributors = ['H2-H2', 'H2-He'],
cloud_species = ['Mg2SiO4(s)_crystalline__DHS',
'Fe(s)_amorphous__Mie'],
wavelength_boundaries = [0.3, 15],
scattering_in_emission=True)
wavelengths, flux, _ = atmosphere_scat.calculate_flux(temperatures=temperatures,
mass_fractions=mass_fractions,
mean_molar_masses=mean_molar_masses,
reference_gravity=reference_gravity,
eddy_diffusion_coefficients=eddy_diffusion_coefficients,
cloud_f_sed=cloud_f_sed,
cloud_particle_radius_distribution_std = cloud_particle_radius_distribution_std)
ax.plot(wavelengths*1e4, flux, label = 'cloudy, scattering')
ax.set_xscale('log')
ax.set_xlabel('Wavelength [microns]')
ax.set_ylabel(r'Planet flux, $F_{\lambda}$ [erg cm$^{-2}$ s$^{-1}$ cm$^{-1}$]')
ax.legend()
plt.show()
Loading Radtrans opacities...
Loading line opacities of species 'H2O' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/H2O/1H2-16O/1H2-16O__POKAZATEL.R1000_0.3-50mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO-NatAbund' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO/C-O-NatAbund/C-O-NatAbund__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CH4' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CH4/12C-1H4/12C-1H4__HITEMP.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'CO2' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/CO2/12C-16O2/12C-16O2__UCL-4000.R1000_0.3-50mu.ktable.petitRADTRANS.h5.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'Na' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/Na/23Na/23Na__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Loading line opacities of species 'K' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/lines/correlated_k/K/39K/39K__Allard.R1000_0.1-250mu.ktable.petitRADTRANS.h5'... Done.
Successfully loaded all line opacities
Loading CIA opacities for H2-H2 from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--H2/H2--H2-NatAbund/H2--H2-NatAbund__BoRi.R831_0.6-250mu.ciatable.petitRADTRANS.h5'... Done.
Loading CIA opacities for H2-He from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/collision_induced_absorptions/H2--He/H2--He-NatAbund/H2--He-NatAbund__BoRi.DeltaWavenumber2_0.5-500mu.ciatable.petitRADTRANS.h5'... Done.
Successfully loaded all CIA opacities
Loading opacities of cloud species 'Mg2SiO4(s)_crystalline__DHS' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/clouds/Mg2SiO4(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000/Mg2-Si-O4-NatAbund(s)_crystalline_000__DHS.R39_0.1-250mu.cotable.petitRADTRANS.h5'(crystalline_000, using DHS scattering)... Done.
Loading opacities of cloud species 'Fe(s)_amorphous__Mie' from file '/Users/molliere/Desktop/input_data_v3/input_data/opacities/continuum/clouds/Fe(s)_amorphous/Fe-NatAbund(s)_amorphous/Fe-NatAbund(s)_amorphous__Mie.R39_0.1-250mu.cotable.petitRADTRANS.h5'(amorphous, using Mie scattering)... Done.
Successfully loaded all clouds opacities
Successfully loaded all opacities
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).
Standard flux units: before pRT3 flux was accessed as atmosphere.flux after running atmosphere.calc_flux(), which contained flux as \(F_\nu\), so in units of erg cm\(^{-2}\) s\(^{-1}\) Hz\(^{-1}\). pRT’s calculate_flux() method now returns wavelength and flux as \(F_\lambda\) in its standard setting, so flux in erg cm\(^{-2}\) s\(^{-1}\) cm\(^{-1}\). To return frequencies and \(F_\nu\), instead of wavelengths and \(F_\lambda\),
please set the keyword frequencies_to_wavelengths=False when calling calculate_flux().