Tidal effects

In this notebook we explore the Roche lobe effects in one-dimensional Parker wind models for exospheres of hot exoplanets. These effects were notably discussed in Erkaev et al. (2007), and were implemented in p-winds version 1.3.0.

The process of including the Roche lobe (or tidal) effects in p-winds models is simple, and it involves changing how the radius at the sonic point and outflow structure are calculated.

We start by importing all the necessary packages.

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

pylab.rcParams['figure.figsize'] = 9.0,6.5
pylab.rcParams['font.size'] = 18

/home/docs/checkouts/readthedocs.org/user_builds/p-winds/envs/dev/lib/python3.11/site-packages/p_winds/tools.py:24: UserWarning: Environment variable PWINDS_REFSPEC_DIR is not set.
  warn("Environment variable PWINDS_REFSPEC_DIR is not set.")

We are going to replicate the quickstart example for HD 209458 b and include the tidal effects. We will assume that the planet has an isothermal upper atmosphere with temperature of \(9\,100\) K and a total mass loss rate of \(2 \times 10^{10}\) g s\(^{-1}\) based on the results from Salz et al. 2016. We will also assume: * The atmosphere is made up of only H and He * The H number fraction is \(0.9\) * Initially a fully neutral atmosphere (this is going to be self-consistently calculated later).

We will also need to know other parameters, namely: the stellar mass and radius, and the semi-major axis of the orbit.

# HD 209458 b planetary parameters, measured
R_pl = 1.39  # Planetary radius in Jupiter radii
M_pl = 0.73  # Planetary mass in Jupiter masses
impact_parameter = 0.499  # Transit impact parameter
a_pl = 0.04634  # Orbital semi-major axis in astronomical units

# HD 209458 stellar parameters
R_star = 1.20  # Stellar radius in solar radii
M_star = 1.07  # Stellar mass in solar masses

# A few assumptions about the planet's atmosphere
m_dot = 10 ** 10.27  # Total atmospheric escape rate in g / s
T_0 = 9100  # Wind temperature in K
h_fraction = 0.90  # H number fraction
he_fraction = 1 - h_fraction  # He number fraction
he_h_fraction = he_fraction / h_fraction
mean_f_ion = 0.0  # Mean ionization fraction (will be self-consistently calculated later)
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)

Next, we retrieve the high-energy spectrum of the host star with fluxes at the planet. For this example, we use the solar spectrum for convenience.

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()

Now we can calculate the distribution of ionized/neutral hydrogen. And here is where things differ a bit to include the tidal effects. We need to specify what is the stellar mass and the semi-major axis of the planet.

initial_f_ion = 0.0
r = np.logspace(0, np.log10(20), 100)  # Radial distance profile in unit of planetary radii

f_r, mu_bar = hydrogen.ion_fraction(r, R_pl, T_0, h_fraction,
                            m_dot, M_pl, mu_0, star_mass=M_star,
                            semimajor_axis=a_pl,
                            spectrum_at_planet=spectrum, exact_phi=True,
                            initial_f_ion=initial_f_ion, relax_solution=True,
                            return_mu=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('Number fraction')
plt.xlim(1, 10)
plt.ylim(0, 1)
plt.legend()
plt.show()

Instead of using the functions parker.radius_sonic_point() and parker.structure(), we use parker.radius_sonic_point_tidal() and parker.structure_tidal(). Their inputs are a bit different, but the outputs are the same. As before, the velocities and densities calculated by parker.structure_tidal() are measured in units of sound speed and density at the sonic point, respectively.

vs = parker.sound_speed(T_0, mu_bar)  # Speed of sound (km/s, assumed to be constant)
rs = parker.radius_sonic_point_tidal(M_pl, vs, M_star, a_pl)  # Radius at the sonic point (jupiterRad)
rhos = parker.density_sonic_point(m_dot, rs, vs)  # Density at the sonic point (g/cm^3)

r_array = r * R_pl / rs
v_array, rho_array = parker.structure_tidal(r_array, vs, rs, M_pl, M_star, a_pl)

# 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()

Let’s compare the structure with and without tidal effects to see how they differ. We will plot the profiles with the effects implemented in full lines, and those without in dashed lines.

f_r_nte, mu_bar_nte = hydrogen.ion_fraction(r, R_pl, T_0, h_fraction,
                            m_dot, M_pl, mu_0,
                            spectrum_at_planet=spectrum, exact_phi=True,
                            initial_f_ion=initial_f_ion, relax_solution=True,
                            return_mu=True)

vs_nte = parker.sound_speed(T_0, mu_bar_nte)  # Speed of sound (km/s, assumed to be constant)
rs_nte = parker.radius_sonic_point(M_pl, vs_nte)  # Radius at the sonic point (jupiterRad)
rhos_nte = parker.density_sonic_point(m_dot, rs_nte, vs_nte)  # Density at the sonic point (g/cm^3)

r_array_nte = r * R_pl / rs_nte
v_array_nte, rho_array_nte = parker.structure(r_array_nte)

# Convenience arrays for the plots
v_plot_nte = v_array_nte * vs_nte
rho_plot_nte = rho_array_nte * rhos_nte

# 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()

# With tidal effects
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')

# No tidal effects
ax1.semilogy(r_plot, v_plot_nte, color='C0', ls='--')
ax1.plot(rs_nte / R_pl, vs_nte, marker='o', markeredgecolor='w', color='C0')
ax2.semilogy(r_plot, rho_plot_nte, color='C1', ls='--')
ax2.plot(rs_nte / R_pl, rhos_nte, 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()

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 inner layer of the atmosphere.

# In the initial state, the fraction of singlet and triplet helium
# are, respectively, 1.0 and 0.0
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_fraction, vs, rs, rhos, spectrum,
    initial_state=initial_state, 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
he_fraction = 1 - h_fraction
n_he = (rho_array * rhos * he_fraction / (h_fraction + 4 * he_fraction) / 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()

As you can see, the resulting structure and He number density profiles are significantly different than the ones from the quickstart example, so for HD 209458 b the tidal effects are very important. You can find more details about these effects in Erkaev et al. (2007), Murray-Clay et al. (2009) and in Appendix E in Vissapragada et al. (2022).