In this notebook you will learn the basics of using p-winds
to model the upper atmosphere (up to many planetary radii) of a H/He-dominated planet, and go all the way to calculating the transit signal in the Helium triplet at 1.083$\mu$m based on an atmospheric escape model.
p-winds
is largely based on the theoretical framework of Oklopčić & Hirata (2018) and Lampón et al. (2020), which themselves based their work on the stellar wind model of Parker (1958). Some of the radiative transfer calculations were based on Koskinen et al. (2010) and Allan & Vidotto (2019).
Notice: p-winds
is not suitable for simulations of the lower atmosphere (above $\sim 10^{-7}$ bar). Also, this page is a Jupyter notebook, and you can download it here.
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.pylab as pylab
import astropy.units as u
import astropy.constants as c
from p_winds import tools, parker, hydrogen, helium, transit, microphysics
# Uncomment the next line if you have a MacBook with retina screen
# %config InlineBackend.figure_format = 'retina'
pylab.rcParams['figure.figsize'] = 9.0,6.5
pylab.rcParams['font.size'] = 18
Let's start by setting the planetary and stellar parameters of HD 209458 b. We will assume that our planet has an isothermal upper atmosphere with temperature of $9 \times 10^3$ K and a total mass loss rate of $8 \times 10^{10}$ g s$^{-1}$ based on the results from Salz et al. 2016. We will also assume a H vs. He fraction (in number of atoms) of $0.9$, and an average ion fraction of $0.7$.
# HD 209458 b
R_pl = 1.39 # Planetary radius in Jupiter radii
M_pl = 0.73 # Planetary mass in Jupiter masses
m_dot = 8E10 # Total atmospheric escape rate in g / s
T_0 = 9000 # Wind temperature in K
h_he = 0.9 # H to He fraction
mean_f_ion = 0.7 # Mean ionization fraction
The next step is to calculate the structure of the planetary atmosphere in terms of densities and velocities.
Note: p-winds
requires the use of astropy.Quantity
for some of its inputs to avoid errors due to unit conversion. Furthermore, some of the quantities that the code calculates are in "convenience units" to avoid numerical overflows (e.g., velocities and densities calculated by parker.structure()
are measured in units of sound speed and density at the sonic point, respectively).
vs = parker.sound_speed(T_0, h_he, mean_f_ion) # Speed of sound (km/s, assumed to be constant)
rs = parker.radius_sonic_point(M_pl, vs) # Radius at the sonic point (jupiterRad)
rhos = parker.density_sonic_point(m_dot, rs, vs) # Density at the sonic point (g/cm^3)
# The `parker.structure()` function requires us to pass values of radius in units of
# radius at the sonic point (`rs`). So, first we build an `r_array` from 1 to 15
# planetary radii, than change its unit to `rs`
r_array = np.linspace(1, 15, 500) * R_pl / rs
v_array, rho_array = parker.structure(r_array)
# Convenience arrays for the plots
r_plot = r_array * rs / R_pl
v_plot = v_array * vs
rho_plot = rho_array * rhos
# Finally plot the structure of the upper atmosphere
# The circles mark the velocity and density at the sonic point
ax1 = plt.subplot()
ax2 = ax1.twinx()
ax1.semilogy(r_plot, v_plot, color='C0')
ax1.plot(rs / R_pl, vs, marker='o', markeredgecolor='w', color='C0')
ax2.semilogy(r_plot, rho_plot, color='C1')
ax2.plot(rs / R_pl, rhos, marker='o', markeredgecolor='w', color='C1')
ax1.set_xlabel(r'Radius (R$_{\rm pl}$)')
ax1.set_ylabel(r'Velocity (km s$^{-1}$)', color='C0')
ax2.set_ylabel(r'Density (g cm$^{-3}$)', color='C1')
ax1.set_xlim(1, 10)
plt.show()
The next step is to calculate the steady-state distribution of H (neutral or ionized) in the atmosphere. To do that, first we need to retrieve the high-energy spectrum of the host star with fluxes at the planet. For convenience, there is a text file in the data
folder of the p-winds
package containing the spectrum arriving at HD 209458 b (HD209458b_spectrum_lambda.dat
). This spectrum was retrieved from the X-exoplanets database, the unit of energy was changed from photons to erg, and the flux scaled to the semi-major axis of HD 209458 b. But, for now, we shall instead use the solar spectrum, because it covers all the wavelengths important for the helium steady-state.
There is a convenience method in tools.make_spectrum_from_file()
that takes text files as input for the spectrum and transforms it into a dict
that can be used as input for our calculations.
Note: We will see below that it is also possible to make calculations without the need of input spectra, but we will need monochromatic fluxes in the X-rays + extreme ultraviolet (0-911 Å and 0-504 Å) and far- to middle-ultraviolet (911-2600 Å).
units = {'wavelength': u.angstrom, 'flux': u.erg / u.s / u.cm ** 2 / u.angstrom}
spectrum = tools.make_spectrum_from_file('../../data/solar_spectrum_scaled_lambda.dat',
units)
plt.loglog(spectrum['wavelength'], spectrum['flux_lambda'])
plt.ylim(1E-5, 1E4)
plt.xlabel(r'Wavelength (${\rm \AA}$)')
plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\rm \AA}^{-1}$)')
plt.show()
Finally, we can calculate the distribution of ionized/neutral hydrogen. This involves calculating the differential equation 13 in Oklopčić & Hirata (2018). To achieve this, we start from an initial state at the innermost layer of the atmosphere and utilize the hydrogen.ion_fraction()
function. It takes as input many of the parameters we already set above.
One thing that you may want to change is the initial_f_ion
of the integration, which is an optional parameter. The initial state is the y0
of the differential equation to be solved. This array has the initial value of f_ion
(ionization fraction) at the surface of the planet. The standard value for this parameter is 0.0
, i.e., completely neutral near the surface.
# We compute `f_ion` from 1 to 15 planetary radii
r = np.linspace(1.0, 20, 500)
initial_f_ion = 0.0
f_r = hydrogen.ion_fraction(r, R_pl, T_0, h_he,
m_dot, M_pl, mean_f_ion,
spectrum_at_planet=spectrum,
initial_f_ion=initial_f_ion, relax_solution=True)
And now we plot the result.
f_ion = f_r
f_neutral = 1 - f_r
plt.plot(r, f_neutral, color='C0', label='f$_\mathrm{neutral}$')
plt.plot(r, f_ion, color='C1', label='f$_\mathrm{ion}$')
plt.xlabel(r'Radius (R$_\mathrm{pl}$)')
plt.ylabel('Number fraction')
plt.xlim(1, 10)
plt.ylim(0, 1)
plt.legend()
plt.show()
In the next step, we calculate the population of helium in singlet, triplet and ionized states using helium.population_fraction()
. Similar to the hydrogen module, we integrate starting from the outer layer of the atmosphere. Also, in this module we change the absolute and relative numerical tolerances of the integrator so it doesn't lose precision. This is necessary because, in some cases, the fractions of singlet and triplet state can differ by a few orders of magnitude, and since everything is solved at the same time, the solver can easily lose precision.
atol = 1E-8 # Absolute numerical tolerance for the solver
rtol = 1E-8 # Relative numerical tolerance
# In the initial state, the fraction of singlet and triplet helium
# is 1E-6, and the optical depths are null
initial_state = np.array([1.0, 0.0])
f_he_1, f_he_3 = helium.population_fraction(
r, v_array, rho_array, f_ion,
R_pl, T_0, h_he, vs, rs, rhos, spectrum,
initial_state=initial_state, atol=atol, rtol=rtol, relax_solution=True)
Finally, we plot the number densities of ionized helium, and helium in singlet and triplet states.
# Hydrogen atom mass
m_h = c.m_p.to(u.g).value
# Number density of helium nuclei
n_he = (rho_array * rhos * (1 - h_he) / (1 + 4 * (1 - h_he)) / m_h)
n_he_1 = f_he_1 * n_he
n_he_3 = f_he_3 * n_he
n_he_ion = (1 - f_he_1 - f_he_3) * n_he
plt.semilogy(r, n_he_1, color='C0', label='He singlet')
plt.semilogy(r, n_he_3, color='C1', label='He triplet')
plt.semilogy(r, n_he_ion, color='C2', label='He ionized')
plt.xlabel(r'Radius (R$_\mathrm{pl}$)')
plt.ylabel('Number density (cm$^{-3}$)')
plt.xlim(1, 10)
plt.ylim(1E-2, 1E10)
plt.legend()
plt.show()
Now that we have the 1-D profile of volumetric densities of the He triplet, we can calculate the spectroscopic transit signal using simplified ray tracing and radiative transfer. For now we will assume an impact factor of zero (so the planet is in the dead center of the star) and also a phase of zero (so the planet is at mid-transit phase). The phases of $-0.5$ and $+0.5$ correspond to the ingress and egress.
# First convert everything to SI units because they make our lives
# much easier.
R_pl_physical = R_pl * 71492000 # Planet radius in m
r_SI = r * R_pl_physical # Array of altitudes in m
n_he_3_SI = n_he_3 * 1E6 # Volumetric densities in 1 / m ** 3
planet_to_star_ratio = 0.12086
# Set up the simplified ray tracing. We will use a coarse 101-px grid size,
# we do not need more resolution than that. Higher resolutions can slow
# down the calculation significantly.
flux_map, density_map = transit.draw_transit(
planet_to_star_ratio,
impact_parameter=0.0, # We set the impact parameter to transit center
phase=0.0, # Phase is also set to transit center,
# check documentation for the allowed ranges
density_profile=n_he_3_SI,
profile_radius=r_SI,
planet_physical_radius=R_pl_physical,
grid_size=101
)
# And now we plot it just to check how the transit looks
plt.imshow(flux_map / density_map, origin='lower')
plt.show()
Next, we calculate the radiative transfer to estimate the spectroscopic signal in the He triplet. Since it is a triplet, we need to know the properties of the transitions. There is a handy function in the module microphysics
that returns these properties.
# Retrieve the properties of the triplet; they were hard-coded
# using the tabulated values of the NIST database
# wX = central wavelength, fX = oscillator strength, a_ij = Einstein coefficient
w0, w1, w2, f0, f1, f2, a_ij = microphysics.he_3_properties()
m_He = 4 * 1.67262192369e-27 # Helium atomic mass in kg
wl = np.linspace(1.0827, 1.0832, 1000) * 1E-6 # Wavelengths in m
v_wind = -2E3 # Line-of-sight wind velocity in m / s
# First, let's do the radiative transfer for each line of the triplet
# separately. Check the documentation to understand what are the
# input parameters, as there are many of them.
spectrum_0 = transit.radiative_transfer(flux_map, density_map,
wl, w0, f0, a_ij, T_0, m_He, v_wind)
spectrum_1 = transit.radiative_transfer(flux_map, density_map,
wl, w1, f1, a_ij, T_0, m_He, v_wind)
spectrum_2 = transit.radiative_transfer(flux_map, density_map,
wl, w2, f2, a_ij, T_0, m_He, v_wind)
# Finally let's calculate the combined spectrum of all lines in the triplet
# To do that, we combine all the line properties in their respective arrays
w_array = np.array([w0, w1, w2])
f_array = np.array([f0, f1, f2])
a_array = np.array([a_ij, a_ij, a_ij]) # This is the same for all lines in then triplet
spectrum = transit.radiative_transfer(flux_map, density_map,
wl, w_array, f_array, a_array, T_0, m_He, v_wind)
plt.plot(wl * 1E6, spectrum_0, ls='--')
plt.plot(wl * 1E6, spectrum_1, ls='--')
plt.plot(wl * 1E6, spectrum_2, ls='--')
plt.plot(wl * 1E6, spectrum, color='k', lw=2)
plt.xlabel('Wavelength in air ($\mu$m)')
plt.ylabel('Normalized flux')
plt.show()
In case you don't have an input spectrum handy, you can also calculate ionization fractions of H and the He population using monochromatic fluxes at specific wavelength channels. The results will probably be a bit different from the ones calculated with a spectrum. Based on the solar spectrum arriving at HD 209458 b, a reasonable estimate of these fluxes is:
$f_{504} = 500$ erg/s/cm$^2$ (wavelength range to ionize He singlet, 0 - 504 Å)
$f_{911} = 1000$ erg/s/cm$^2$ (wavelength range to ionize H, 0 - 911 Å)
$f_\mathrm{uv} = 5 \times 10^6$ erg/s/cm$^2$ (wavelength range to ionize He triplet, 911 - 2593 Å)
f504 = 500 #* u.erg / u.s / u.cm ** 2
f911 = 1000 #* u.erg / u.s / u.cm ** 2
fuv = 5E6 #* u.erg / u.s / u.cm ** 2
# We compute `f_ion` from 1 to 15 planetary radii
r = np.linspace(1.0, 20, 500)
initial_f_ion = 0.0
f_r = hydrogen.ion_fraction(r, R_pl, T_0, h_he,
m_dot, M_pl, mean_f_ion,
flux_euv=f911,
initial_f_ion=initial_f_ion, relax_solution=True)
f_ion = f_r
f_neutral = 1 - f_r
plt.plot(r, f_neutral, color='C0', label='f$_\mathrm{neutral}$')
plt.plot(r, f_ion, color='C1', label='f$_\mathrm{ion}$')
plt.xlabel(r'Radius (R$_\mathrm{pl}$)')
plt.ylabel('H number fraction')
plt.xlim(1, 10)
plt.ylim(0, 1)
plt.legend()
plt.show()
Finally, we calculate the helium population with the monocromatic fluxes.
atol = 1E-8 # Absolute numerical tolerance for the solver
rtol = 1E-8 # Relative numerical tolerance
# In the initial state, the fraction of singlet and triplet helium is 1E-6, and the optical depths are null
initial_state = np.array([1.0, 0.0])
f_he_1, f_he_3 = helium.population_fraction(
r, v_array, rho_array, f_ion,
R_pl, T_0, h_he, vs, rs, rhos, flux_euv=f504, flux_fuv=fuv,
initial_state=initial_state, atol=atol, rtol=rtol, relax_solution=True)
# Number density of helium nuclei
n_he = (rho_array * rhos * (1 - h_he) / (1 + 4 * (1 - h_he)) / m_h)
n_he_1 = f_he_1 * n_he
n_he_3 = f_he_3 * n_he
n_he_ion = (1 - f_he_1 - f_he_3) * n_he
plt.semilogy(r, n_he_1, color='C0', label='He singlet')
plt.semilogy(r, n_he_3, color='C1', label='He triplet')
plt.semilogy(r, n_he_ion, color='C2', label='He ionized')
plt.xlabel(r'Radius (R$_\mathrm{pl}$)')
plt.ylabel('Number density (cm$^{-3}$)')
plt.xlim(1, 10)
plt.ylim(1E-2, 1E10)
plt.legend()
plt.show()