import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import solve_for_masses as em
import mass_loss as ms
from numpy.random import normal
%matplotlib inline
from importlib import reload
reload(ms)
<module 'mass_loss' from '/Users/bcamposestrada/Library/CloudStorage/OneDrive-UniversityofCopenhagen/EvapMass/mass_loss.py'>
As an example we work with the famous Kepler-36 system, which contains a rocky and gaseous planet. Planet parameters are taken from Carter et al. (2012, Science 337 556).
To start with we need to make some basic choices. First we need to assume the composition of the solid cores. This uses the Fortney et al. (2007, ApJ 659 1661) mass-radius relation. They can either be iron-rock or water-rock (both the iron fraction and ice fraction cannot be non-zero). Xiron=1/3 will consider an "Earth-like" composition of 1/3 iron and 2/3 silicate rock.
Xiron = 1./3.
Xice = 0.
As discussed in the paper (Owen & Campos Estrada, 2020), one needs to choose a Kelvin-Helmholtz timescale for the H/He atmospheres at which to do the comparision (i.e. at what age is mass-loss important), the answer is very weakly dependent on this value, here we pick the standard value of 100 Myr.
Tkh_Myr=100.
Now we input the planetary and stellar parameters and their errors. Any errors in the Carter et al. (2012) paper that are not symmetric we just make symmetric crudely using addition in quadrature.
K36s_M = 1.071 # stellar mass, solar masses
K36s_Mer = 0.043 # stellar mass error, solar masses
K36s_Teff = 5911 # Stellar effective temperature, K
K36s_Teffer = 66 # stellar effective temperature error, K
K36s_R = 1.626 # stellar radius, solar radii - note we should really use the MS radii but this is and example.
K36s_Rer=0.019 # stellar radius error, solar radii
K36s_age=6800 # stellar age, Myr
K36s_age_er=1000. #stellar age error, Myr
#radius in earth unit and period in days + errors for planet b
K36b_R = 1.486
K36b_Rer=0.035
K36b_P = 13.83989
K36b_Per = np.sqrt(0.00082**2.+0.00060**2.) # make period errors symmetric
#radius in earth unit and period in days + errors for planet c
K36c_R = 3.679
K36c_Rer=0.054
K36c_P = 16.23855
K36c_Per = np.sqrt(0.00038**2.+0.00054**2.)
A note on the efficiency parameter. The code contains three options for the efficiency parameter in the mass-loss models.
# set the efficiency option
eff_option = 3
Now do the calculation for the minimum mass of Kepler-36c to be consistent with photoevaporation assuming the mean values for all parameters.
The location of the radius valley can be chosen. The default is set to 1.8 Rearth, but this is generally only true for sun-like stars. One should change the location of the valley accordingly (for example for M-dwarfs the valley is generally at smaller radii - see van Eylen et al. 2021).
valley_loc=1.8
# evaluate minimum mass for the mean value and print it
system = em.psystem('Kepler36')
system.add_planet('36b',K36b_R,K36b_Rer,K36b_P,K36b_Per)
system.add_planet('36c',K36c_R,K36c_Rer,K36c_P,K36c_Per)
system.star.mass=K36s_M
system.star.radius=K36s_R
system.star.Teff=K36s_Teff
system.star.age = K36s_age
system.update_planet_info()
system.above_or_below_valley()
system.mass_rocky(Xiron,Xice)
ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)
Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)
Mmedian = np.copy(Mout)
print(Mmedian,'Mearth')
8.412384425520381 Mearth
However, given the stellar and planetary parameters contain errors these must be included, therefore we assume the errors are gaussian and randonly sample them. Here we use 1000 samples, but a value should be chosen such that the upper-limit is converged.
Mout_error= np.zeros(1000)
flag_out = np.zeros(1000,dtype=np.int8)
for i in range(np.size(Mout_error)):
K36b_R_use = normal(K36b_R,K36b_Rer,1)
K36b_P_use = normal(K36b_P,K36b_Per,1)
K36c_R_use = normal(K36c_R,K36c_Rer,1)
K36c_P_use = normal(K36c_P,K36c_Per,1)
system = em.psystem('Kepler36_%d' %i)
system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)
system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)
system.star.mass=normal(K36s_M,K36s_Mer,1)
system.star.radius=normal(K36s_R,K36s_Rer,1)
system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)
system.star.age = normal(K36s_age,K36s_age_er,1)
system.update_planet_info()
system.above_or_below_valley()
system.mass_rocky(Xiron,Xice)
ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)
Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)
Mout_error[i] = Mout
flag_out[i]=flag
Now plot the resulting distribution and calculate the 95% upper-limit, the 95% upperlimit is plotted as a dashed line, the actual measured mass is plotted as the dotted line. This indicates (as known from previous work - Lopez & Fortney 2013, Owen & Morton 2016) that Kepler-36 is consistent with the photoevaporation model. The value to quote, as we are concerned with the minimum mass to be consistent is something like the 95% upper-limit, not the mean or median value.
plt.hist(Mout_error,bins=30,density=True,histtype='step',lw=2)
#note if density=True this simply means are under the histogram is 1, not the height of the bars sums to 1!
plt.plot([8.08,8.08],[0.,1.],':',color='b', lw=1.5) #Must set this to the meassured mass of the planet
plt.plot([(np.percentile(Mout_error,5.)),(np.percentile(Mout_error,5.))],[0.,1.],'--',color='k', lw=1.5)
plt.ylim((0.,1.))
plt.ylabel('Probability density',fontsize=12)
plt.xlabel(r'Minimum Core Mass [M$_\oplus$]',fontsize=12)
print('The 95% upper limit to be consistent with photoevaporation is')
print((np.percentile(Mout_error,5.)),'Mearth')
The 95% upper limit to be consistent with photoevaporation is 6.930341420699428 Mearth
Now we will run for the two other efficiency options to show the effect on the result
eff_option = 1
Mout_error_1= np.zeros(1000)
flag_out_1= np.zeros(1000,dtype=np.int8)
for i in range(np.size(Mout_error)):
K36b_R_use = normal(K36b_R,K36b_Rer,1)
K36b_P_use = normal(K36b_P,K36b_Per,1)
K36c_R_use = normal(K36c_R,K36c_Rer,1)
K36c_P_use = normal(K36c_P,K36c_Per,1)
system = em.psystem('Kepler36_%d' %i)
system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)
system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)
system.star.mass=normal(K36s_M,K36s_Mer,1)
system.star.radius=normal(K36s_R,K36s_Rer,1)
system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)
system.star.age = normal(K36s_age,K36s_age_er,1)
system.update_planet_info()
system.above_or_below_valley()
system.mass_rocky(Xiron,Xice)
ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)
Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)
Mout_error_1[i] = Mout
flag_out_1[i]=flag
eff_option = 2
Mout_error_2= np.zeros(1000)
flag_out_2= np.zeros(1000,dtype=np.int8)
for i in range(np.size(Mout_error)):
K36b_R_use = normal(K36b_R,K36b_Rer,1)
K36b_P_use = normal(K36b_P,K36b_Per,1)
K36c_R_use = normal(K36c_R,K36c_Rer,1)
K36c_P_use = normal(K36c_P,K36c_Per,1)
system = em.psystem('Kepler36_%d' %i)
system.add_planet('36b',K36b_R_use,K36b_Rer,K36b_P_use,K36b_Per)
system.add_planet('36c',K36c_R_use,K36c_Rer,K36c_P_use,K36c_Per)
system.star.mass=normal(K36s_M,K36s_Mer,1)
system.star.radius=normal(K36s_R,K36s_Rer,1)
system.star.Teff=normal(K36s_Teff,K36s_Teffer,1)
system.star.age = normal(K36s_age,K36s_age_er,1)
system.update_planet_info()
system.above_or_below_valley()
system.mass_rocky(Xiron,Xice)
ms.find_hardest_rocky(system,Tkh_Myr,Xiron,Xice,eff_option)
Mout, flag = ms.min_mass_gaseous(system.planets[0],system.planets[1],Tkh_Myr,Xiron,Xice,system.star.age,eff_option)
Mout_error_2[i] = Mout
flag_out_2[i]=flag
## Now plot
plt.hist(Mout_error,bins=30,density=True,histtype='step',lw=2,label='Owen & Jackson (2012) efficiency')
plt.hist(Mout_error_1,bins=30,density=True,histtype='step',lw=2,label='Constant Efficiency')
plt.hist(Mout_error_2,bins=30,density=True,histtype='step',lw=2,label='Escape velocity scaling')
plt.ylim((0.,1.))
plt.ylabel('Probability Density',fontsize=12)
plt.xlabel(r'Minimum Core Mass [M$_\oplus$]',fontsize=12)
plt.legend(loc=0)
<matplotlib.legend.Legend at 0x7fba44259fd0>
We see for planets like Kepler-36c, which contain a large amount of H/He today (relative to other sub-neptunes ~10% compared to 1-2%) the choice of efficiency parameter matters. At young ages these inflated planets with massive H/He envelopes are large. At this point advective energy losses become important the efficiency of the outflow drops with falling escape velocity.