Interface with ATES
ATES is a hydrodynamic escape code written in Fortran by Andrea Caldiroli and published in Caldiroli et al. 2021. It is a more sophisticated formulation of the problem than the isothermal Parker wind assumption used in the original p-winds code. In this notebook, we show you how to use the fluid module as an interface with the ATES code with an API that is familiar
to p-winds users. Then, we show how to use the ATES outputs to calculate a transmission spectrum using the radiative transfer tools of p-winds.
[1]:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import astropy.constants as c
import astropy.units as u
import pickle
from astropy.convolution import convolve
from scipy.optimize import minimize
from p_winds import parker, hydrogen, helium, transit, lines, fluid, tools
pylab.rcParams['figure.figsize'] = 9.0,6.5
pylab.rcParams['font.size'] = 18
/home/docs/checkouts/readthedocs.org/user_builds/p-winds/envs/latest/lib/python3.11/site-packages/p_winds/tools.py:25: UserWarning: Environment variable PWINDS_REFSPEC_DIR is not set.
warn("Environment variable PWINDS_REFSPEC_DIR is not set.")
/home/docs/checkouts/readthedocs.org/user_builds/p-winds/envs/latest/lib/python3.11/site-packages/p_winds/fluid.py:33: UserWarning: Environment variable ATES_DIR is not set.
warn("Environment variable ATES_DIR is not set.")
/home/docs/checkouts/readthedocs.org/user_builds/p-winds/envs/latest/lib/python3.11/site-packages/p_winds/fluid.py:39: UserWarning: No Fortran compiler found in your PATH.
warn("No Fortran compiler found in your PATH.")
Let’s start with the observation of the He triplet transmission spectrum of HAT-P-11 b using the CARMENES spectrograph. This data is openly available in the DACE platform. But we will retrieve it from a public Gist for convenience.
[2]:
# The observed transmission spectrum
data_url = 'https://gist.githubusercontent.com/ladsantos/a8433928e384819a3632adc469bed803/raw/a584e6e83073d1ad3444248624927838588f22e4/HAT-P-11_b_He.dat'
# We skip 2 rows instead of 1 to have an odd number of rows and allow a fast convolution later
He_spec = np.loadtxt(data_url, skiprows=2)
wl_obs = He_spec[:, 0] # Angstrom
f_obs = 1 - He_spec[:, 1] * 0.01 # Normalized flux
u_obs = He_spec[:, 2] * 0.01 # Flux uncertainty
# Convert in-vacuum wavelengths to in-air
s = 1E4 / np.mean(wl_obs)
n = 1 + 0.0000834254 + 0.02406147 / (130 - s ** 2) + 0.00015998 / (38.9 - s ** 2)
wl_obs /= n
# We will also need to know the instrumental profile that
# widens spectral lines. We take the width from Allart et al. (2018),
# the paper describing the HAT-P-11 b data.
def gaussian(x, mu=0.0, sigma=1.0):
return 1 / sigma / (2 * np.pi) ** 0.5 * np.exp(-0.5 * (x - mu) ** 2 / sigma ** 2)
instrumental_profile_width_v = 3.7 # Instrumental profile FWHM in km / s (assumed Gaussian)
sigma_wl = instrumental_profile_width_v / (2 * (2 * np.log(2)) ** 0.5) / \
c.c.to(u.km / u.s).value * np.mean(wl_obs) # Same unit as wl_obs
instrumental_profile = gaussian(wl_obs, np.mean(wl_obs), sigma=sigma_wl)
plt.errorbar(wl_obs, f_obs, yerr=u_obs)
plt.xlabel(r'Wavelength (${\rm \AA}$)')
plt.ylabel('Normalized flux')
plt.show()
Now we set up the simulation. This is very similar to the quick start example and the advanced tutorial.
[3]:
# Set up the simulation
# Fixed parameters of HAT-P-11 b (not to be sampled)
R_pl = 0.389 # Planetary radius (Jupiter radii)
M_pl = 0.09 # Planetary mass (Jupiter masses)
a_pl = 0.05254 # Semi-major axis (au)
pl_teq = 847 # Planet's equilibrium temperature (K)
M_star = 0.809 # Stellar host mass (Solar masses)
R_star = 0.683 # Stellar host radius (Solar radii)
a_star = ((a_pl * u.au) / (R_star * u.solRad)).decompose().value
planet_to_star_ratio = 0.057989
impact_parameter = 0.132
# mean_f_ion = 0.90 # Initially assumed, but the model relaxes for it
# mu_0 = (1 + 4 * he_h_fraction) / (1 + he_h_fraction + mean_f_ion)
# # mu_0 is the constant mean molecular weight (assumed for now, will be updated later)
# Physical constants
m_h = c.m_p.to(u.g).value # Hydrogen atom mass in g
m_He = 4 * 1.67262192369e-27 # Helium atomic mass in kg
k_B = 1.380649e-23 # Boltzmann's constant in kg / (m / s) ** 2 / K
# Free parameters: ATES has several free parameters, which you can consult in the documentation
# For now, our only free parameter will be the H number fraction
h_fraction = 0.90 # H number fraction, chosen solar for this example
he_fraction = 1 - h_fraction # He number fraction
he_h_fraction = he_fraction / h_fraction
# Altitudes samples (this can be a very important setting)
# r = np.logspace(0, np.log10(20), 100)
# First guesses of fractions (not to be fit, but necessary for the calculation)
# initial_f_ion = 0.0 # Fraction of ionized hydrogen
# initial_f_he = np.array([1.0, 0.0]) # Fraction of singlet, triplet helium
# Model settings
# relax_solution = True # This will iteratively relax the solutions until convergence
# exact_phi = True # Exact calculation of H photoionization
sample_phases = np.linspace(-0.50, 0.50, 5) # Phases that we will average to obtain the final spectrum
# The phases -0.5 and +0.5 correspond to the times of first and fourth transit contact
w0, w1, w2, f0, f1, f2, a_ij = lines.he_3_properties()
w_array = np.array([w0, w1, w2]) # Central wavelengths of the triplet
f_array = np.array([f0, f1, f2]) # Oscillator strengths of the triplet
a_array = np.array([a_ij, a_ij, a_ij]) # This is the same for all lines in then triplet
n_samples = len(sample_phases)
transit_grid_size = 100 # Also very important to constrain computation time
supersampling = 5 # This is used to improve the hard pixel edges in the ray tracing
The full spectrum of HAT-P-11 until 2600 Å is not known. But we can use a proxy for which we do have a full spectrum: HD 40307. It has a similar size, spectral type, effective temperature, and surface gravity as HAT-P-11. We take the spectrum from the MUSCLES database. There is a convenience function in tools that calculates the spectrum arriving at a planet based on the MUSCLES SEDs.
[4]:
data_url = 'https://gist.githubusercontent.com/ladsantos/5fc31bfbb49812ae3aa321f34cc30106/raw/be774e96c80504a9b43d95b4bfc2264e417c7b5c/HAT-P-11b_full_spec.dat'
spec = np.loadtxt(data_url, skiprows=1)
host_spectrum = {'wavelength': spec[:, 0], 'flux_lambda': spec[:, 1],
'wavelength_unit': u.angstrom,
'flux_unit': u.erg / u.s / u.cm ** 2 / u.angstrom}
plt.loglog(host_spectrum['wavelength'], host_spectrum['flux_lambda'])
plt.xlabel(r'Wavelength (${\rm \AA}$)')
plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')
plt.show()
Next, we can jump into modeling the planet’s upper atmosphere. ATES is significantly more computationally expensive than p-winds, so this online tutorial will not run ATES and instead has the lines below commented out. If you are running this notebook in your own machine, you need to uncomment the lines in the cell below.
The momentum threshold parameter sets when the model is converged and its stopping condition. For this first run, we use momentum_threshold=0.05.
[5]:
# Calculate the model: if you are running this notebook in your machine, uncomment the following lines
# ates_results = fluid.ates_model(
# planet_radius=R_pl,
# planet_mass=M_pl,
# planet_equilibrium_temperature=pl_teq,
# semi_major_axis=a_pl,
# stellar_mass=M_star,
# spectrum_at_planet=host_spectrum,
# he_h_ratio=he_h_fraction,
# # The next three parameters are important for this first run of the ATES model.
# # Check the documentation and Caldiroli+2021 for explanations on these.
# load_ic='False',
# reconstruction='PLM',
# momentum_threshold=0.05
# )
Alright, so for the second run we load the previous model as initial condition using load_ic=True and use reconstruction='WENO3' and set the momentum threshold to 0.02.
[6]:
# Note: if you are running this notebook in your machine, uncomment the following lines
# ates_results = fluid.ates_model(
# planet_radius=R_pl,
# planet_mass=M_pl,
# planet_equilibrium_temperature=pl_teq,
# semi_major_axis=a_pl,
# stellar_mass=M_star,
# spectrum_at_planet=host_spectrum,
# he_h_ratio=he_h_fraction,
# # The next three parameters are changed in the second run of the ATES model
# load_ic=True,
# reconstruction='WENO3',
# momentum_threshold=0.02
# )
Since it’s expensive to run an ATES model, you can also save it using pickle. To do that, uncomment the following command. For this online notebook, we shall leave it commented out since we didn’t run ATES here.
[7]:
# Save results
# with open('hat_p_11b_ates_result.pkl', 'wb') as f:
# pickle.dump(ates_results, f)
To read previously saved results, run the following command.
[8]:
# Open results
with open('hat_p_11b_ates_result.pkl', 'rb') as f:
ates_results = pickle.load(f)
The ATES hydrodynamic escape model is done! The results are contained within a dict object, which is trivial to parse. Please consult the documentation to see what comes inside it.
[9]:
log_m_dot_0 = ates_results['log_m_dot'] # Log10 of sub-stellar mass loss rate (g / s)
m_dot = 10 ** log_m_dot_0 # Sub-stellar mass loss rate (g / s)
r = ates_results['r']
v_ates = ates_results['velocity']
temp_ates = ates_results['temperature']
# Some of other profiles that we can analyze
f_h_ates = ates_results['n_h_ii'] / ates_results['n_h_i'] # Fraction of ionized H
n_he_i_ates = ates_results['n_he_i'] # Number density of neutral He
n_he_ii_ates = ates_results['n_he_ii'] # Number density of singly-ionized He
n_he_23s_ates = ates_results['n_he_23s'] # Number density of neutral He in the metastable 2^3S state -- this is what we observe in exoplanets
plt.plot(r, temp_ates)
plt.xlabel(r'Radial distance (R$_{\rm pl}$)')
plt.ylabel(r'Temperature (K)')
[9]:
Text(0, 0.5, 'Temperature (K)')
With that, we are ready to calculate the transmission spectrum resulting from the ATES model. This is similar to the quickstar and advanced tutorial notebooks.
[10]:
# The transmission spectrum model
def transmission_model(wavelength_array, v_wind, n_he_3_distribution, log_T, v_array, v_turb=False, v_rot=0.0):
# Set up the transit configuration. We use SI units to avoid too many
# headaches with unit conversion
R_pl_physical = R_pl * 71492000 # Planet radius in m
r_SI = r * R_pl_physical # Array of altitudes in m
v_SI = v_array * 1000 # Velocity of the outflow in m / s
n_he_3_SI = n_he_3_distribution * 1E6 # Volumetric densities in 1 / m ** 3
# Set up the ray tracing
f_maps = []
t_depths = []
r_maps = []
for i in range(n_samples):
flux_map, transit_depth, r_map = transit.draw_transit(
planet_to_star_ratio,
impact_parameter=impact_parameter,
supersampling=supersampling,
phase=sample_phases[i],
planet_physical_radius=R_pl_physical,
grid_size=transit_grid_size
)
f_maps.append(flux_map)
t_depths.append(transit_depth)
r_maps.append(r_map)
# Do the radiative transfer
spectra = []
for i in range(n_samples):
spec = transit.radiative_transfer_2d(f_maps[i], r_maps[i],
r_SI, n_he_3_SI, v_SI, w_array, f_array, a_array,
wavelength_array, 10 ** log_T, m_He, bulk_los_velocity=v_wind,
wind_broadening_method='formal', turbulence_broadening=v_turb, v_rotation=v_rot)
# We add the transit depth because ground-based observations
# lose the continuum information and they are not sensitive to
# the loss of light by the opaque disk of the planet, only
# by the atmosphere
spectra.append(spec + t_depths[i])
spectra = np.array(spectra)
# Finally we take the mean of the spectra we calculated for each phase
spectrum = np.mean(spectra, axis=0)
return spectrum
log_T_ates = np.log10(ates_results['temperature'])
v_wind_0 = -2E3 # Line-of-sight wind velocity (m / s)
# Here we divide wl_obs by 1E10 to convert angstrom to m
t_spectrum = transmission_model(wl_obs / 1E10, v_wind_0, n_he_23s_ates, log_T_ates, v_ates)
plt.errorbar(wl_obs, f_obs, yerr=u_obs)
plt.plot(wl_obs, t_spectrum, color='k', lw=2)
plt.xlabel(r'Wavelength (${\rm \AA}$)')
plt.ylabel('Normalized flux')
plt.show()
As you can see, the ATES model overpredicts the metastable He absorption compared to the observation, by a large factor! One solution is to increase the H fraction to values closer to unity, such as 0.99, instead of assuming solar value. We leave this change as an exercise to the reader.