Chempy

we will now introduce the Chempy function which will calculate the chemical evolution of a one-zone open box model



In [1]:

    
%pylab inline









    



Populating the interactive namespace from numpy and matplotlib



In [2]:

    
# loading the default parameters

from Chempy.parameter import ModelParameters
a = ModelParameters()

Loading all the input

solar abundances
SFR
infall
initial abundances and inflowing abundances



In [3]:

    
# Initialising sfr, infall, elements to trace, solar abundances

from Chempy.wrapper import initialise_stuff
basic_solar, basic_sfr, basic_infall = initialise_stuff(a)
elements_to_trace = a.elements_to_trace

Elemental abundances at start

We need to define the abundances of:

The ISM at beginning
The corona gas at beginning
The cosmic inflow into the corona for all times. For all we chose primordial here.



In [4]:

    
# Setting the abundance fractions at the beginning to primordial

from Chempy.infall import INFALL, PRIMORDIAL_INFALL
basic_primordial = PRIMORDIAL_INFALL(list(elements_to_trace),np.copy(basic_solar.table))
basic_primordial.primordial()
basic_primordial.fractions









    Out[4]:





array([ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.76,  0.24,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  0.  ,
        0.  ,  0.  ,  0.  ,  0.  ,  0.  ])

Initialising the element evolution matrix

We now feed everything into the abundance matrix and check its entries



In [5]:

    
# Initialising the ISM instance

from Chempy.time_integration import ABUNDANCE_MATRIX
cube = ABUNDANCE_MATRIX(np.copy(basic_sfr.t),np.copy(basic_sfr.sfr),np.copy(basic_infall.infall),list(elements_to_trace),list(basic_primordial.symbols),list(basic_primordial.fractions),float(a.gas_at_start),list(basic_primordial.symbols),list(basic_primordial.fractions),float(a.gas_reservoir_mass_factor),float(a.outflow_feedback_fraction),bool(a.check_processes),float(a.starformation_efficiency),float(a.gas_power), float(a.sfr_factor_for_cosmic_accretion), list(basic_primordial.symbols), list(basic_primordial.fractions))


# All the entries of the ISM instance

print(list(cube.cube.dtype.names))


# Helium at start

print('Primordial ratio of H to He: ',cube.cube['H'][0]/cube.cube['He'][0])
print('Helium over time: ',cube.cube['He'])









    



['sfr', 'infall', 'time', 'feedback', 'mass_in_remnants', 'stars', 'gas', 'Z', 'alpha', 'sn1a', 'sn2', 'pn', 'bh', 'hn', 'Al', 'Ar', 'B', 'Be', 'C', 'Ca', 'Cl', 'Co', 'Cr', 'Cu', 'F', 'Fe', 'Ga', 'Ge', 'H', 'He', 'K', 'Li', 'Mg', 'Mn', 'N', 'Na', 'Ne', 'Ni', 'O', 'P', 'S', 'Sc', 'Si', 'Ti', 'V', 'Zn']
Primordial ratio of H to He:  3.16666666667
Helium over time:  [  9.26869415e-05   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00
   0.00000000e+00   0.00000000e+00   0.00000000e+00   0.00000000e+00]

Time integration

With the advance_one_step method we can evolve the matrix in time, given that we provide the feedback from each steps previous SSP.



In [6]:

    
# Now we run the time integration

from Chempy.wrapper import SSP_wrap
basic_ssp = SSP_wrap(a)
for i in range(len(basic_sfr.t)-1):
    j = len(basic_sfr.t)-i
    
    # The metallicity needs to be passed for the yields to be calculated as well as the initial elemental abundances

    
    element_fractions = []
    for item in elements_to_trace:
        element_fractions.append(float(np.copy(cube.cube[item][max(i-1,0)]/cube.cube['gas'][max(i-1,0)])))## gas element fractions from one time step before	   
    metallicity = float(cube.cube['Z'][i])
    time_steps = np.copy(basic_sfr.t[:j])
    basic_ssp.calculate_feedback(float(metallicity), list(elements_to_trace), list(element_fractions), np.copy(time_steps))
    cube.advance_one_step(i+1,np.copy(basic_ssp.table),np.copy(basic_ssp.sn2_table),np.copy(basic_ssp.agb_table),np.copy(basic_ssp.sn1a_table),np.copy(basic_ssp.bh_table))
print(cube.cube['He'])









    



[  9.26869415e-05   1.40862570e-02   2.46524030e-02   3.23386825e-02
   3.76707786e-02   4.11093438e-02   4.30465899e-02   4.38022068e-02
   4.36425361e-02   4.27879025e-02   4.14303451e-02   3.97149234e-02
   3.77562120e-02   3.56443278e-02   3.34470737e-02   3.12219166e-02
   2.90100138e-02   2.68477097e-02   2.47550032e-02   2.27529674e-02
   2.08520421e-02   1.90581357e-02   1.73784234e-02   1.58118906e-02
   1.43563052e-02   1.30109258e-02   1.17715493e-02   1.06334198e-02]

Making abundances from element fractions

The cube stores everything in elemental fractions, we use a tool to convert these to abundances scaled to solar:



In [7]:

    
# Turning the fractions into dex values (normalised to solar [X/H])

from Chempy.making_abundances import mass_fraction_to_abundances
abundances,elements,numbers = mass_fraction_to_abundances(np.copy(cube.cube),np.copy(basic_solar.table))
print(abundances['He'])









    



[-0.02951008 -0.02942689 -0.0221913  -0.01535444 -0.00917462 -0.00349022
  0.00185113  0.00685321  0.01154271  0.01595232  0.0203323   0.02472094
  0.02911161  0.03348354  0.03775661  0.04195023  0.04600286  0.04999718
  0.05383017  0.05758016  0.06122583  0.06471174  0.06816704  0.07150743
  0.0746963   0.07781998  0.08085494  0.08381115]



In [8]:

    
## Alpha enhancement over time

plot(cube.cube['time'][1:],abundances['O'][1:]-abundances['Fe'][1:])
plt.xlabel('time in Gyr')
plt.ylabel('[O/Fe]')









    Out[8]:





Text(0,0.5,'[O/Fe]')



In [9]:

    
# [X/Fe] vs. [Fe/H]

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], label = 'O')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], label = 'Mn')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], label = 'N')
plt.xlabel('[Fe/H]')
plt.ylabel('[X/Fe]')
plt.legend()









    Out[9]:





<matplotlib.legend.Legend at 0x7f610176ff60>

Likelihood calculation

There are a few build-in functions (actually representing the observational constraints from the Chempy paper) which return a likelihood. One of those is called sol_norm and compares the proto-solar abundances with the Chempy ISM abundances 4.5 Gyr ago.



In [10]:

    
# Here we load a likelihood test for the solar abundances
# This is how it looks for the prior parameters with the default yield set

from Chempy.data_to_test import sol_norm
probabilities, abundance_list, element_names = sol_norm(True,a.name_string,np.copy(abundances),np.copy(cube.cube),elements_to_trace,a.element_names,np.copy(basic_solar.table),a.number_of_models_overplotted,a.produce_mock_data,a.use_mock_data,a.error_inflation)









    





<matplotlib.figure.Figure at 0x7f610183a8d0>

Net vs. total yield

Now we will change a little detail in the time-integration. Instead of letting unprocessed material that is expelled from the stars ('unprocessed_mass_in_winds' in the yield tables) being composed of the stellar birth material, which would be consistent (and is what I call 'net' yield), we now use solar-scaled material which only has the same metallicity as the stellar birth material (This is what happens if yield tables are giving the total yield including the unprocessed material, which means that the author usually uses solar-scaled material which is then expelled by the star, but might not even be produced by it). Therefore we see a difference in the likelihood which is better for the total yields case (-180.05 vs -198.30). We see the difference especially well in K and Ti.



In [11]:

    
cube = ABUNDANCE_MATRIX(np.copy(basic_sfr.t),np.copy(basic_sfr.sfr),np.copy(basic_infall.infall),list(elements_to_trace),list(basic_primordial.symbols),list(basic_primordial.fractions),float(a.gas_at_start),list(basic_primordial.symbols),list(basic_primordial.fractions),float(a.gas_reservoir_mass_factor),float(a.outflow_feedback_fraction),bool(a.check_processes),float(a.starformation_efficiency),float(a.gas_power), float(a.sfr_factor_for_cosmic_accretion), list(basic_primordial.symbols), list(basic_primordial.fractions))
basic_ssp = SSP_wrap(a)
for i in range(len(basic_sfr.t)-1):
    j = len(basic_sfr.t)-i
    metallicity = float(cube.cube['Z'][i])
    
    # Instead of using the ISM composition we use solar scaled material
    
    solar_scaled_material = PRIMORDIAL_INFALL(list(elements_to_trace),np.copy(basic_solar.table))
    solar_scaled_material.solar(np.log10(metallicity/basic_solar.z))
    time_steps = np.copy(basic_sfr.t[:j])
    basic_ssp.calculate_feedback(float(metallicity), list(elements_to_trace), list(solar_scaled_material.fractions), np.copy(time_steps))
    cube.advance_one_step(i+1,np.copy(basic_ssp.table),np.copy(basic_ssp.sn2_table),np.copy(basic_ssp.agb_table),np.copy(basic_ssp.sn1a_table),np.copy(basic_ssp.bh_table))
abundances,elements,numbers = mass_fraction_to_abundances(np.copy(cube.cube),np.copy(basic_solar.table))


# We do the solar abundance test again and see that the likelihood improves

probabilities, abundance_list, element_names = sol_norm(True,a.name_string,np.copy(abundances),np.copy(cube.cube),elements_to_trace,a.element_names,np.copy(basic_solar.table),a.number_of_models_overplotted,a.produce_mock_data,a.use_mock_data,a.error_inflation)









    





<matplotlib.figure.Figure at 0x7f61281d8ba8>

Making chemical evolution modelling fast and flexible

Now we have all ingredients at hand. We use a wrapper function were we only need to pass the ModelParameters.



In [12]:

    
# This is a convenience function 

from Chempy.wrapper import Chempy
a = ModelParameters()
cube, abundances = Chempy(a)
plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], label = 'O')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], label = 'Mn')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], label = 'N')
plt.xlabel('[Fe/H]')
plt.ylabel('[X/Fe]')
plt.legend()









    Out[12]:





<matplotlib.legend.Legend at 0x7f60f86f8b38>

IMF effect

now we can easily check the effect of the IMF on the chemical evolution



In [13]:

    
# prior IMF

a = ModelParameters()
a.imf_parameter= (0.69, 0.079,-2.29)
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], label = 'O', color = 'b')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], label = 'Mn', color = 'orange')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], label = 'N', color = 'g')


# top-heavy IMF

a = ModelParameters()
a.imf_parameter = (0.69, 0.079,-2.09)
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = ':')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = ':')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'g', linestyle = ':')

# bottom-heavy IMF

a = ModelParameters()
a.imf_parameter = (0.69, 0.079,-2.49)
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = '--')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = '--')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'g', linestyle = '--')


plt.xlabel('[Fe/H]')
plt.ylabel('[X/Fe]')
plt.title('IMF effect: top-heavy as dotted line, bottom-heavy as dashed line')
plt.legend()









    Out[13]:





<matplotlib.legend.Legend at 0x7f60f8f32978>

SFR effect

We can do the same for the peak of the SFR etc...



In [14]:

    
# Prior SFR

a = ModelParameters()
a.sfr_scale = 3.5
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], label = 'O', color = 'b')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], label = 'Mn', color = 'orange')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], label = 'N', color = 'g')


# Early peak in the SFR

a = ModelParameters()
a.sfr_scale = 1.5
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = ':')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = ':')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'green', linestyle = ':')


# late peak in the SFR

a = ModelParameters()
a.sfr_scale = 6.5
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = '--')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = '--')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'green', linestyle = '--')


plt.xlabel('[Fe/H]')
plt.ylabel('[X/Fe]')
plt.title('SFR effect: early peak as dotted line, late peak as dashed line')
plt.legend()









    Out[14]:





<matplotlib.legend.Legend at 0x7f60f8eb3710>

Time resolution

The time steps are equidistant and the resolution is flexible. Even with coarse 0.5Gyr resolution the results are quite good, saving a lot of computational time. Here we test different time resolution of 0.5, 0.1 and 0.025 Gyr.

All results converge after metallicity increases above -1. The shorter time sampling allows more massive stars to explode first which generally have alpha enhanced yields, therefore the [O/Fe] is higher in the beginning.



In [15]:

    
## 0.5 Gyr resolution

a = ModelParameters()
a.time_steps = 28 # default
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], label = 'O', color = 'b')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], label = 'Mn', color = 'orange')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], label = 'N', color = 'g')


# 0.1 Gyr resolution

a = ModelParameters()
a.time_steps = 136
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = ':')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = ':')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'green', linestyle = ':')


# 25 Myr resolution

a = ModelParameters()
a.time_steps = 541
cube, abundances = Chempy(a)

plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], color = 'b', linestyle = '--')
plot(abundances['Fe'][1:],abundances['Mn'][1:]-abundances['Fe'][1:], color = 'orange', linestyle = '--')
plot(abundances['Fe'][1:],abundances['N'][1:]-abundances['Fe'][1:], color = 'green', linestyle = '--')


plt.xlabel('[Fe/H]')
plt.ylabel('[X/Fe]')
plt.title('Time resolution effect: 0.5 solid, 0.1 dotted, 0.025Gyr dashed line')
plt.legend()









    Out[15]:





<matplotlib.legend.Legend at 0x7f60f86c3eb8>

A note on chemical evolution tracks and 'by eye' fit

Sometimes Astronomers like to show that their chemical evolution track runs through some stellar abundance data points. But if we want the computer to steer our result fit we need to know the selection function of the stars that we try to match and we need to take our star formation history into account (Maybe there are almost no stars formed after 8Gyr).

We assume that we have an unbiased sample of red clump stars
We reproduce its selection function by multiplying their age-distribution for a flat SFR with the SFR. (for the age distribution I have included a cut from a mock catalogue according to Just&Rybizki 2016 but you could also use the analytic formula from Bovy+2014)
Then we sample some synthetic stars (with observational errors) along the chemical evolutionary track



In [21]:

    
# Default model parameters
from Chempy import localpath

a = ModelParameters()
a.check_processes = True
cube, abundances = Chempy(a)


# Red clump age distribution

selection = np.load(localpath + "input/selection/red_clump_new.npy")
time_selection = np.load(localpath + "input/selection/time_red_clump_new.npy")

plt.plot(time_selection,selection)
plt.xlabel('Age in Gyr')
plt.title('Age distribution of Red clump stars')
plt.show()


# We need to put the age distribution on the same time-steps as our model

selection = np.interp(cube.cube['time'], time_selection[::-1], selection)
plt.plot(cube.cube['time'],selection)
plt.xlabel('time in Gyr')
plt.title('Normalisation for a population of Red clump stars')
plt.show()


# Comparing to the SFR

plt.plot(cube.cube['time'],cube.cube['sfr'])
plt.xlabel('time in Gyr')
plt.title('SFR')
plt.show()


# Convolution of SFR and Red clump age distribution

weight = cube.cube['sfr']*selection
plt.plot(cube.cube['time'],weight)
plt.xlabel('time in Gyr')
plt.title('Weight to reproduce red clump stellar sample')
plt.show()



In [22]:

    
# Here we sample 1000 stars with this age-distribution

from Chempy.data_to_test import sample_stars
sample_size = 1000
x,y = sample_stars(cube.cube['sfr'][1:],selection[1:],abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:],float(basic_solar.table['error'][np.where(basic_solar.table['Symbol']=='Fe')]),float(basic_solar.table['error'][np.where(basic_solar.table['Symbol']=='O')]),int(sample_size))
plt.plot(x,y,"g.", alpha = 0.2, label = '(%d) synthesized red clum stars' %(int(sample_size)))
plt.plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:], 'r', label = 'evolutionary track')
plt.xlabel('[Fe/H]')
plt.ylabel('[O/Fe]')
plt.title("Sampling from SFH and red clump age distribution")
plt.legend(bbox_to_anchor = [1,1.5])
plt.show()



In [23]:

    
# And we plot the 2d histogramm where we see that our model prediction for a red clump population


plt.hist2d(x,y,20)
plt.plot(abundances['Fe'][1:],abundances['O'][1:]-abundances['Fe'][1:],'r')
plt.xlabel('[Fe/H]')
plt.ylabel('[O/Fe]')
plt.title("Sampling from SFH and red clump age distribution")
plt.show









    Out[23]:





<function matplotlib.pyplot.show>

This PDF can then be compared to real data to get a realistic likelihood.

The nucleosynthetic feedback per element

With the plot_processes routine we can plot the total feedback of each element and the fractional contribution from each nucleosynthetic feedback for a specific Chempy run.



In [24]:

    
# Loading the routine and plotting the process contribution into the current folder
# Total enrichment mass in gray to the right, single process fractional contribution to the left

from Chempy.data_to_test import plot_processes
plot_processes(True,a.name_string,cube.sn2_cube,cube.sn1a_cube,cube.agb_cube,a.element_names,np.copy(cube),a.number_of_models_overplotted)









    Out[24]:





[0]






    





<matplotlib.figure.Figure at 0x7f60f8e57470>



In [ ]: