<jupyter_start><jupyter_text>SN IaThis notebook explores how to normalize the SN Ia DTD given by Maoz and Graur 2017.<jupyter_code>import tabulation
import numpy as np
import betterplotlib as bpl
from scipy import integrate
from astropy import units as u
from astropy import cosmology

bpl.presentation_style()
bpl.presentation_style()
lifetimes = tabulation.Lifetimes("Raiteri_96")<jupyter_output><empty_output><jupyter_text>There are options in exactly how to parametrize the DTD. We can do$$ r(t) \propto t^{-1.13} $$ where we just use the raw age, or we can scale the age to the age of the 8 solar mass star.$$ r(t) \propto \left( \frac{t}{\tau(8)} \right)^{-1.13} $$These options matter since the lifetimes vary with metallicity, changing $\tau(8)$. We want to choose one of these options. The first is the same at late times for all metallicities, while the second has the same peak height for all metallicities.<jupyter_code>new_art_s = 1.13

def z_dependent_rate_scaled(age, z):
    age_8 = lifetimes.lifetime(8.0, z)
    if age < age_8:
        return 0
    else:
        return (age / age_8) ** (-new_art_s)
    
def z_dependent_rate_raw(age, z):
    age_8 = lifetimes.lifetime(8.0, z)
    if age < age_8:
        return 0
    else:
        return age ** (-new_art_s)
z_test = [1E-6, 1E-5, 1E-4, 1E-3, 1E-2, 0.02, 0.04]
z_plot = [1E-4, 1E-3, 2E-2]<jupyter_output><empty_output><jupyter_text>Plot them to see these differences visually.<jupyter_code>fig, axs = bpl.subplots(ncols=2, figsize=[15, 6])
ages = np.logspace(7, 10, 1000)
for z in z_plot:
    rate_raw = [z_dependent_rate_raw(a, z) for a in ages]
    rate_scaled = [z_dependent_rate_scaled(a, z) for a in ages]
    
    axs[0].plot(ages, rate_scaled, label=z)
    axs[1].plot(ages, rate_raw, label=z)

for ax in axs:
    ax.add_labels("Time [yr]", "SN Ia rate [SN year$^{-1}$ M$^{-1}_\odot$]")
    ax.set_yscale("log")
    ax.set_xscale("log")
    ax.legend()
    
axs[0].set_title("Scaled")
axs[1].set_title("Raw");<jupyter_output>/Users/gillenb/miniconda3/envs/general/lib/python3.5/site-packages/matplotlib/mathtext.py:839: MathTextWarning: Substituting with a symbol from Computer Modern.
  MathTextWarning)
/Users/gillenb/miniconda3/envs/general/lib/python3.5/site-packages/matplotlib/figure.py:2022: UserWarning: This figure includes Axes that are not compatible with tight_layout, so results might be incorrect.
  warnings.warn("This figure includes Axes that are not compatible "<jupyter_text>Test how they integrate. For better normalization, we'd like them to be pretty consistent at all metallicities.<jupyter_code>for z in z_test:
    integral = integrate.quad(z_dependent_rate_scaled, 0, 14E9, args=[z])[0]
    age_8 = lifetimes.lifetime(8, z)
    print("{:7f} {:.3E} {:.2E}".format(z, age_8, integral))<jupyter_output>0.000001 2.850E+07 1.21E+08
0.000010 2.850E+07 1.21E+08
0.000100 3.076E+07 1.30E+08
0.001000 4.172E+07 1.70E+08
0.010000 4.078E+07 1.67E+08
0.020000 3.798E+07 1.57E+08
0.040000 3.594E+07 1.49E+08<jupyter_text>Here metallicities with shorter lifetimes contribute less over roughly a Hubble time.<jupyter_code>for z in z_test:
    integral = integrate.quad(z_dependent_rate_raw, 0, 14E9, args=[z])[0]
    age_8 = lifetimes.lifetime(8, z)
    print("{:7f} {:.3E} {:.2E}".format(z, age_8, integral))<jupyter_output>0.000001 2.850E+07 4.57E-01
0.000010 2.850E+07 4.57E-01
0.000100 3.076E+07 4.49E-01
0.001000 4.172E+07 4.17E-01
0.010000 4.078E+07 4.19E-01
0.020000 3.798E+07 4.27E-01
0.040000 3.594E+07 4.32E-01<jupyter_text>Here metallicities with shorter lifetimes contribute more over a Hubble time. This seems more justified, since more white dwarfs will be created (since more stars die), therefore more SN Ia is reasonable.The raw version is clearly more stable. This is because the early times don't contribute much. Late times are much more important (since we integrate over time, not log time as shown in the plots). Having equal values at late times is more important to keeping the rates equal.The raw version seems to be better.  Determining normalizationWe want the total number of SN Ia to match the normalization given in Maoz and Graur 2017, Table 3 for field galaxies: $N/M_{\star} = 1.6^{+0.3}_{-0.3} \times 10^{-3}$ over a Hubble time. Since the normalization will depend on stellar mass, we have to pick a metallicity at which to normalize. We choose solar.However, I want to be general so that the user in ART can specifiy a different $N$ and have it work easily. I'll instead normalize to 1 so that we can simply multiply by the desired $N$ to get it to work. I'll check this later in the file.We need to integrate over a Hubble time. What we choose for that might or might not matter. I'll test to see what variations happen with different choices for the Hubble time.<jupyter_code>for t_h in [13E9, 14E9, 15E9]:
    age_8 = lifetimes.lifetime(8, z)
    print(integrate.quad(z_dependent_rate_raw, age_8, t_h, args=[z])[0])<jupyter_output>0.42876138150908166
0.4323337970416205
0.4356288425889357<jupyter_text>So the exact choice doesn't really matter. I'll use the cosmology given in the paper (Hui's first paper).<jupyter_code>omega_matter = 0.304
hubble_const = 0.681
omega_lambda = 1.0 - omega_matter
hubble_const_units = 100 * hubble_const * u.km / (u.second * u.Mpc)

sim_cosmo = cosmology.LambdaCDM(H0=hubble_const_units, 
                                Om0=omega_matter, 
                                Ode0=omega_lambda)
t_hubble = sim_cosmo.age(0)
print(t_hubble)
t_hubble = t_hubble.to(u.year).value<jupyter_output>13.791485359241204 Gyr<jupyter_text>Then we can get the normalizaton. We want it to integrate to the N given above$$ N = \int_{\tau(8 M_\odot)}^{t_H} A f(t) dt $$$$ A = \frac{N}{\int_{\tau(8 M_\odot)}^{t_H} f(t) dt} $$<jupyter_code>true_N = 1.6E-3
norm_N = 1.0
age_8_solar = lifetimes.lifetime(8, 0.02)
raw_integral = integrate.quad(z_dependent_rate_raw, age_8_solar, t_hubble, args=[0.02])[0]

normalization = norm_N / raw_integral<jupyter_output><empty_output><jupyter_text>So here is the final functional form, with the parameters printed below. This still needs to be multiplied by the number of SN Ia to get the real rate, however.<jupyter_code>def base_art_rate(age, z):
    age_8 = lifetimes.lifetime(8.0, z)
    if age < age_8:
        return 0
    else:
        return normalization * age ** (-new_art_s)
    
def final_art_rate(age, z, N_SNIa):
    return N_SNIa * base_art_rate(age, z)
    
print("A = {:.10f}".format(normalization))
print("s = {:.2f}".format(new_art_s))<jupyter_output>A = 2.3480851917
s = 1.13<jupyter_text>As a check, integrate to infinity for all metallicities, to see how much it varies.<jupyter_code>for z in z_test:
    age_8 = lifetimes.lifetime(8, z)
    integral = integrate.quad(base_art_rate, age_8, t_hubble, args=[z])[0]
    print("{:5f} {:.4E}".format(z, integral))<jupyter_output>0.000001 1.0711E+00
0.000010 1.0711E+00
0.000100 1.0519E+00
0.001000 9.7735E-01
0.010000 9.8284E-01
0.020000 1.0000E+00
0.040000 1.0135E+00<jupyter_text>These values are pretty close. They are within the roughly 10% error in the rate in the Maoz and Graur paper, so I think this is fine.<jupyter_code>fig, ax = bpl.subplots()
ages = np.logspace(7, 10, 1000)
for z in z_plot:
    rate_final = [final_art_rate(a, z, true_N) for a in ages]
    ax.plot(ages, rate_final, label=z)
    
ax.add_labels("Time [yr]", "SN Ia rate [SN year$^{-1}$ M$^{-1}_\odot$]")
ax.set_yscale("log")
ax.set_xscale("log")
ax.legend()<jupyter_output><empty_output><jupyter_text>Checking variations in the normalizationIn ART there will be a parameter that controls the number of SN Ia per SSP. If this is varied I want to be able to vary the prescription. I am 100% sure I can just scale the normalization, since it's linear, but want numerical proof. I'll just check this at solar metallicity, since that is where the normalization is done.<jupyter_code>z = 0.02
for N in [1.3E-3, 1E-3, 1E-2, 1.787982E-2]:
    age_8 = lifetimes.lifetime(8, z)
    integral= integrate.quad(final_art_rate, age_8, t_hubble, args=(z, N))[0]
    print("{:.6E}  {:.6E}".format(N, integral))<jupyter_output>1.300000E-03  1.300000E-03
1.000000E-03  1.000000E-03
1.000000E-02  1.000000E-02
1.787982E-02  1.787982E-02
{'fraction_non_alphanumeric': 0.1158188660920238, 'fraction_numerical': 0.08765327182226539, 'mean_word_length': 4.301158301158301, 'pattern_counts': {'":': 0, '<': 44, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '281223', 'license': 'MIT', 'n_tokens_mistral': 3516, 'n_tokens_neox': 3019, 'n_words': 920, 'path': 'notebooks/sn_Ia.ipynb', 'repo_name': 'gillenbrown/Tabulation', 'size': '8237'}