{ "cells": [ { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "Exospheric metals\n", "===============\n", "\n", "Some extremely hot planets may be in a state of hydrodynamic escape, in which the outflow of H is so intense, that it can drag up heavier species, like C and O to the upper atmosphere of the planet. In this notebook, we will use `p-winds` to quantify the amount of C and O in the exosphere of the hot Jupiter HD 209458 b using the modules `carbon` and `oxygen` (version 1.4.1 and onwards).\n", "\n", "As always, let's start by importing the necessary packages." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import matplotlib.pylab as pylab\n", "import astropy.units as u\n", "import astropy.constants as c\n", "from p_winds import tools, parker, hydrogen, helium, carbon, oxygen, lines, transit\n", "\n", "pylab.rcParams['figure.figsize'] = 9.0,6.5\n", "pylab.rcParams['font.size'] = 18" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We are going to replicate the [quickstart example for HD 209458 b](https://p-winds.readthedocs.io/en/latest/quickstart.html) 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](https://ui.adsabs.harvard.edu/abs/2016A%26A...586A..75S/abstract). We will also assume: \n", "* The atmosphere is mostly made up of H and He with number fractions $0.9$ and $0.1$, respectively\n", "* C and O are trace elements with solar abundance based on [Asplund et al. 2009](https://ui.adsabs.harvard.edu/abs/2009ARA%26A..47..481A/abstract). `p-winds` uses these solar values by default, but they can be set by the user if preferred.\n", "* The H and He nuclei are fully neutral near the planet's surface (this is going to be self-consistently calculated later). In the case of C, we assume they are fully singly-ionized near the surface.\n", "\n", "We will also need to know other parameters, namely: the stellar mass and radius, and the semi-major axis of the orbit." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# HD 209458 b planetary parameters, measured\n", "R_pl = 1.39 # Planetary radius in Jupiter radii\n", "M_pl = 0.73 # Planetary mass in Jupiter masses\n", "impact_parameter = 0.499 # Transit impact parameter\n", "a_pl = 0.04634 # Orbital semi-major axis in astronomical units\n", "\n", "# HD 209458 stellar parameters\n", "R_star = 1.20 # Stellar radius in solar radii\n", "M_star = 1.07 # Stellar mass in solar masses\n", "\n", "# A few assumptions about the planet's atmosphere\n", "m_dot = 10 ** 10.27 # Total atmospheric escape rate in g / s\n", "T_0 = 9100 # Wind temperature in K\n", "h_fraction = 0.90 # H number fraction\n", "he_fraction = 1 - h_fraction # He number fraction\n", "he_h_fraction = he_fraction / h_fraction\n", "mean_f_ion = 0.0 # Mean ionization fraction (will be self-consistently calculated later)\n", "mu_0 = (1 + 4 * he_h_fraction) / (1 + he_h_fraction + mean_f_ion) \n", "# mu_0 is the constant mean molecular weight (assumed for now, will be updated later)\n", "\n", "# The trace abundances of C and O\n", "c_abundance = 8.43 # Asplund et al. 2009\n", "c_fraction = 10 ** (c_abundance - 12.00)\n", "o_abundance = 8.69 # Asplund et al. 2009\n", "o_fraction = 10 ** (o_abundance - 12.00)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "units = {'wavelength': u.angstrom, 'flux': u.erg / u.s / u.cm ** 2 / u.angstrom}\n", "spectrum = tools.make_spectrum_from_file('../../data/solar_spectrum_scaled_lambda.dat',\n", " units)\n", "plt.loglog(spectrum['wavelength'], spectrum['flux_lambda'])\n", "plt.ylim(1E-5, 1E4)\n", "plt.xlabel(r'Wavelength (${\\rm \\AA}$)')\n", "plt.ylabel(r'Flux density (erg s$^{-1}$ cm$^{-2}$ ${\\rm \\AA}^{-1}$)')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can calculate the distribution of ionized/neutral hydrogen and the structure of the upper atmosphere." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "initial_f_ion = 0.0\n", "r = np.logspace(0, np.log10(20), 100) # Radial distance profile in unit of planetary radii\n", "\n", "f_r, mu_bar = hydrogen.ion_fraction(r, R_pl, T_0, h_fraction, \n", " m_dot, M_pl, mu_0, star_mass=M_star,\n", " semimajor_axis=a_pl,\n", " spectrum_at_planet=spectrum, exact_phi=True,\n", " initial_f_ion=initial_f_ion, relax_solution=True,\n", " return_mu=True)\n", "\n", "f_ion = f_r\n", "f_neutral = 1 - f_r\n", "\n", "vs = parker.sound_speed(T_0, mu_bar) # Speed of sound (km/s, assumed to be constant)\n", "rs = parker.radius_sonic_point_tidal(M_pl, vs, M_star, a_pl) # Radius at the sonic point (jupiterRad)\n", "rhos = parker.density_sonic_point(m_dot, rs, vs) # Density at the sonic point (g/cm^3)\n", "\n", "r_array = r * R_pl / rs\n", "v_array, rho_array = parker.structure_tidal(r_array, vs, rs, M_pl, M_star, a_pl)\n", "\n", "# Convenience arrays for the plots\n", "r_plot = r_array * rs / R_pl\n", "v_plot = v_array * vs\n", "rho_plot = rho_array * rhos\n", "\n", "# Finally plot the structure of the upper atmosphere\n", "# The circles mark the velocity and density at the sonic point\n", "ax1 = plt.subplot()\n", "ax2 = ax1.twinx()\n", "ax1.semilogy(r_plot, v_plot, color='C0')\n", "ax1.plot(rs / R_pl, vs, marker='o', markeredgecolor='w', color='C0')\n", "ax2.semilogy(r_plot, rho_plot, color='C1')\n", "ax2.plot(rs / R_pl, rhos, marker='o', markeredgecolor='w', color='C1')\n", "ax1.set_xlabel(r'Radius (R$_{\\rm pl}$)')\n", "ax1.set_ylabel(r'Velocity (km s$^{-1}$)', color='C0')\n", "ax2.set_ylabel(r'Density (g cm$^{-3}$)', color='C1')\n", "ax1.set_xlim(1, 10)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will also need the neutral and ion fractions of helium." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate the helium ion fraction\n", "f_he_plus = helium.ion_fraction(\n", " r, v_array, rho_array, f_ion,\n", " R_pl, T_0, h_fraction, vs, rs, rhos, spectrum,\n", " initial_f_he_ion=0.0, relax_solution=True)\n", "\n", "# Hydrogen atom mass\n", "m_h = c.m_p.to(u.g).value\n", "\n", "# Number density of helium nuclei \n", "he_fraction = 1 - h_fraction\n", "n_he = (rho_array * rhos * he_fraction / (h_fraction + 4 * he_fraction) / m_h)\n", "n_he_ion = f_he_plus * n_he" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With all this setup done, we will proceed to calculate the number densities of neutral, singly-ionized and doubly-ionized C. We will assume that the C nuclei are all singly-ionized near the surface (so `initial_f_c_ion` is `[1.0, 0.0]`)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_c_ii, f_c_iii = carbon.ion_fraction(radius_profile=r,\n", " velocity=v_array,\n", " density=rho_array,\n", " hydrogen_ion_fraction=f_ion,\n", " helium_ion_fraction=f_he_plus,\n", " planet_radius=R_pl,\n", " temperature=T_0,\n", " h_fraction=h_fraction,\n", " c_fraction=c_fraction,\n", " speed_sonic_point=vs,\n", " radius_sonic_point=rs,\n", " density_sonic_point=rhos,\n", " spectrum_at_planet=spectrum,\n", " initial_f_c_ion=np.array([1.0, 0.0]),\n", " method='Radau',\n", " relax_solution=True)\n", "\n", "# Number density of carbon nuclei \n", "n_c = (rho_array * rhos * c_fraction / (h_fraction + 4 * he_fraction + 12 * c_fraction) / m_h)\n", "\n", "n_c_i = (1 - f_c_ii - f_c_iii) * n_c\n", "n_c_ii = f_c_ii * n_c\n", "n_c_iii = f_c_iii * n_c\n", "\n", "plt.semilogy(r, n_c_i, color='C0', label='C I')\n", "plt.semilogy(r, n_c_ii, color='C1', label='C II')\n", "plt.semilogy(r, n_c_iii, color='C2', label='C III')\n", "plt.xlabel(r'Radius (R$_\\mathrm{pl}$)')\n", "plt.ylabel('Number density (cm$^{-3}$)')\n", "plt.xlim(1, 10)\n", "plt.ylim(1E-3, 1E7)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In the plot above we see that most of the C nuclei in the upper atmosphere of HD 209458 b are singly-ionized, with only a small fraction of them in neutral state by a factor of 3 orders of magnitude. Furthermore, doubly-ionized C are not present.\n", "\n", "Next, we do the same exercise for oxygen. We will assume that all O nuclei are neutral near the surface at first." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_o_ii = oxygen.ion_fraction(radius_profile=r,\n", " velocity=v_array,\n", " density=rho_array,\n", " hydrogen_ion_fraction=f_ion,\n", " helium_ion_fraction=f_he_plus,\n", " planet_radius=R_pl,\n", " temperature=T_0,\n", " h_fraction=h_fraction,\n", " o_fraction=o_fraction,\n", " speed_sonic_point=vs,\n", " radius_sonic_point=rs,\n", " density_sonic_point=rhos,\n", " spectrum_at_planet=spectrum,\n", " initial_f_o_ion=0.0,\n", " relax_solution=True)\n", "\n", "# Number density of oxygen nuclei \n", "n_o = (rho_array * rhos * o_fraction / (h_fraction + 4 * he_fraction + 16 * o_fraction) / m_h)\n", "\n", "n_o_i = (1 - f_o_ii) * n_o\n", "n_o_ii = f_o_ii * n_o\n", "\n", "plt.semilogy(r, n_o_i, color='C0', label='O I')\n", "plt.semilogy(r, n_o_ii, color='C1', label='O II')\n", "plt.xlabel(r'Radius (R$_\\mathrm{pl}$)')\n", "plt.ylabel('Number density (cm$^{-3}$)')\n", "plt.xlim(1, 10)\n", "plt.ylim(1E-3, 1E7)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And there we see that most of the oxygen nuclei are singly-ionized, except in the innermost upper atmosphere layers, where neutral oxygen dominates.\n", "\n", "For the next part of this tutorial, we will estimate the in-transit absorption profiles for the relevant wavelengths in real observations. These exospheric carbon and oxygen lines are located in the ultraviolet, which are accessible only with the *Hubble Space Telescope*.\n", "\n", "The tricky part of the C II triplet is that the line at 133.4 nm arises from the ground state, and the other two at 133.5 nm are a doublet arising from the first excited state. The number density of C II we calculated above assumes that all nuclei are in the ground state, and we need to estimate how many of them are in the excited state as well.\n", "\n", "There are a few ways of calculating the population of C II nuclei. I used the [`ChiantiPy` code](https://github.com/chianti-atomic/ChiantiPy/), and estimated that, for the exosphere of HD 209458 b (temperature 9100 K and density of electrons $\\sim 6 \\times 10^6$ cm$^{-3}$), we have roughly 33% of C II ions in the ground state, and 66% in the excited state." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Set up the ray tracing. We will use a coarse 100-px grid size,\n", "# but we use supersampling to avoid hard pixel edges.\n", "# We convert everything to SI units because they make our lives\n", "# much easier.\n", "R_pl_physical = R_pl * 71492000 # Planet radius in m\n", "r_SI = r * R_pl_physical # Array of altitudes in m\n", "v_SI = v_array * vs * 1000 # Velocity of the outflow in m / s\n", "n_c_ii_SI = n_c_ii * 1E6 # Volumetric densities in 1 / m ** 3\n", "planet_to_star_ratio = 0.12086\n", "\n", "flux_map, t_depth, r_from_planet = transit.draw_transit(\n", " planet_to_star_ratio, \n", " planet_physical_radius=R_pl_physical, \n", " impact_parameter=impact_parameter, \n", " phase=0.0,\n", " supersampling=10,\n", " grid_size=100)\n", "\n", "# Retrieve the properties of the C II lines; they were hard-coded\n", "# using the tabulated values of the NIST database\n", "# wX = central wavelength, fX = oscillator strength, a_ij = Einstein coefficient\n", "w0, w1, w2, f0, f1, f2, a_ij_0, a_ij_1, a_ij_2 = lines.c_ii_properties()\n", "\n", "m_C = 12 * 1.67262192369e-27 # Carbon atomic mass in kg\n", "wl = np.linspace(1332.9, 1337.1, 200) * 1E-10 # Wavelengths in Angstrom\n", "\n", "method = 'average'\n", "\n", "spectrum_0 = transit.radiative_transfer_2d(flux_map, r_from_planet, # Ground state\n", " r_SI, n_c_ii_SI * 0.33, v_SI, w0, f0, a_ij_0,\n", " wl, T_0, m_C, wind_broadening_method=method) \n", "spectrum_1 = transit.radiative_transfer_2d(flux_map, r_from_planet, # Excited state\n", " r_SI, n_c_ii_SI * 0.66, v_SI, w1, f1, a_ij_1,\n", " wl, T_0, m_C, wind_broadening_method=method)\n", "spectrum_2 = transit.radiative_transfer_2d(flux_map, r_from_planet, # Excited state\n", " r_SI, n_c_ii_SI * 0.66, v_SI, w2, f2, a_ij_2,\n", " wl, T_0, m_C, wind_broadening_method=method)\n", "\n", "plt.plot(wl * 1E10, spectrum_0, ls='--', label='C II, ground state')\n", "plt.plot(wl * 1E10, spectrum_1, ls='--', label='C II, excited state')\n", "plt.plot(wl * 1E10, spectrum_2, ls='--', label='C II, excited state')\n", "plt.legend()\n", "plt.xlabel(r'Wavelength in air (${\\rm \\AA}$)')\n", "plt.ylabel('Normalized flux')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now for the oxygen lines. Similarly, the first line is a resonant line, and the other two are doubles arising from two different excited states. \n", "\n", "From `ChiantiPy`, the ratios for the ground state, the first and second excited states are 0.85, 0.10 and 0.05." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n_o_i_SI = n_o_i * 1E6 # Volumetric densities in 1 / m ** 3\n", "\n", "# Retrieve the properties of the O I triplet; they were hard-coded\n", "# using the tabulated values of the NIST database\n", "# wX = central wavelength, fX = oscillator strength, a_ij = Einstein coefficient\n", "w0, w1, w2, f0, f1, f2, a_ij_0, a_ij_1, a_ij_2 = lines.o_i_properties()\n", "\n", "m_O = 16 * 1.67262192369e-27 # Oxygen atomic mass in kg\n", "wl = np.linspace(1300.9, 1307.1, 200) * 1E-10 # Wavelengths in Angstrom\n", "\n", "method = 'average'\n", "\n", "spectrum_0 = transit.radiative_transfer_2d(flux_map, r_from_planet, # Ground state\n", " r_SI, n_o_i_SI * 0.85, v_SI, w0, f0, a_ij_0,\n", " wl, T_0, m_O, wind_broadening_method=method)\n", "spectrum_1 = transit.radiative_transfer_2d(flux_map, r_from_planet, # First excited state\n", " r_SI, n_o_i_SI * 0.10, v_SI, w1, f1, a_ij_1,\n", " wl, T_0, m_O, wind_broadening_method=method)\n", "spectrum_2 = transit.radiative_transfer_2d(flux_map, r_from_planet, # Second excited state\n", " r_SI, n_o_i_SI * 0.05, v_SI, w2, f2, a_ij_2,\n", " wl, T_0, m_O, wind_broadening_method=method)\n", "\n", "plt.plot(wl * 1E10, spectrum_0, ls='--', label='O I, ground state')\n", "plt.plot(wl * 1E10, spectrum_1, ls='--', label='O I, first excited state')\n", "plt.plot(wl * 1E10, spectrum_2, ls='--', label='O I, second excited state')\n", "plt.legend()\n", "plt.xlabel(r'Wavelength in air (${\\rm \\AA}$)')\n", "plt.ylabel('Normalized flux')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What we can conclude from this exercise is that, if we assume a Parker-wind like escape for HD 209458 b, solar abundances of C and O, and a solar-like high-energy spectrum for the host star, then the excess in-transit signature of neutral O will be about 0.5% in the core of the strongest O I line in the UV. For C II, a much promising signature stronger than 7.5% will be present in the strongest line of the transmission spectrum." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 1 }