In-Class Assignment 6#

A Simple Type II Supernova Lightcurve#

import numpy as np
import matplotlib.pyplot as plt

A simple model of a supernova lightcurve is just to add up the energy of the photons that result from the decay:

\[{}^{56}\mathrm{Ni} \rightarrow {}^{56}\mathrm{Co} \rightarrow {}^{56}\mathrm{Fe}\]

The evolution equations for the number densities appear as:

\[\begin{split} \frac{dn_\mathrm{Ni}}{dt} = -\lambda_\mathrm{Ni} n_\mathrm{Ni} \\ \frac{dn_\mathrm{Co}}{dt} = +\lambda_\mathrm{Ni} n_\mathrm{Ni} - \lambda_\mathrm{Co} n_\mathrm{Co} \end{split}\]

where \(\lambda_\mathrm{Ni}\) is the decay constant for \({}^{56}\mathrm{Ni}\) and \(\lambda_\mathrm{Co}\) is the decay constant for \({}^{56}\mathrm{Co}\).

The initial conditions are:

\[\begin{split} n_\mathrm{Ni}(t = 0) = n_0 \\ n_\mathrm{Co}(t = 0) = 0 \end{split}\]

We can integrate the Ni equation directly, giving:

\[n_\mathrm{Ni} = n_0 e^{-\lambda_\mathrm{Ni} t} \equiv X_{\rm{Ni}}\]

then the Co equation becomes:

\[\frac{dn_\mathrm{Co}}{dt} = \lambda_\mathrm{Ni} n_0 e^{-\lambda_\mathrm{Ni} t} - \lambda_\mathrm{Co} n_\mathrm{Co}\]

This can be integrated by applying an integrating factor:

\[\frac{d}{dt} \left ( e^{\lambda_\mathrm{Co} t} n_\mathrm{Co} \right ) = e^{\lambda_\mathrm{Co} t} \frac{dn_\mathrm{Co}}{dt} + \lambda_\mathrm{Co} e^{\lambda_\mathrm{Co} t} n_\mathrm{Co} = e^{\lambda_\mathrm{Co} t} \left [ \frac{dn_\mathrm{Co}}{dt} + \lambda_\mathrm{Co} n_\mathrm{Co} \right ]\]

Now the term in the \([ \ldots ]\) is just \(\lambda_\mathrm{Ni} n_\mathrm{Ni}\), so this becomes:

\[\frac{d}{dt} \left ( e^{\lambda_\mathrm{Co} t} n_\mathrm{Co} \right ) = e^{\lambda_\mathrm{Co} t} \lambda_\mathrm{Ni} n_0 e^{-\lambda_\mathrm{Ni} t}\]

This gives:

\[n_\mathrm{Co} = \frac{\lambda_\mathrm{Ni}}{\lambda_\mathrm{Co} - \lambda_\mathrm{Ni}} n_0 \left ( e^{-\lambda_\mathrm{Ni} t} - e^{-\lambda_\mathrm{Co} t} \right ) \equiv X_{\rm{Co}}\]

With these number densities, we can get the total amount of Ni or Co at any point in time:

\[N_\mathrm{Ni} = n_\mathrm{Ni} V = \frac{X_\mathrm{Ni} \rho}{m_\mathrm{Ni}} V = \frac{X_\mathrm{Ni} M_\mathrm{exp}}{m_\mathrm{Ni}}\]

and

\[N_\mathrm{Co} = \frac{X_\mathrm{Co} M_\mathrm{exp}}{m_\mathrm{Co}}\]

This gives:

\[ L(t) = -Q_\mathrm{Ni} (-\lambda_\mathrm{Ni} N_\mathrm{Ni} ) - Q_\mathrm{Co} (-\lambda_\mathrm{Co} N_\mathrm{Co} ) \]
\[ L(t) \approx \left [ Q_\mathrm{Ni} \lambda_\mathrm{Ni} X_\mathrm{Ni} + Q_\mathrm{Co} \lambda_\mathrm{Co} X_\mathrm{Co} \right ] \frac{M_\mathrm{exp}}{m_{56}} \]

where \(m_{56}\) is just \(56 m_u\), with \(m_u\) the atomic mass unit (this approximates the mass of Ni and Co).

Implementation#

some physical constants

seconds_per_day = 24*3600
erg_per_MeV = 1.6e-12*1.e6
m_u = 1.66e-24     # CGS

solar_mass = 2.e33          # CGS
solar_luminosity = 4.e33    # CGS

the core LightCurve class

# uncomment to begin
"""
def LightCurve(M_exp, X_Ni_0=1.0):

    # note: see Nadyozhin 1994, ApJ, 92, 527 for Q values of
    # just photons.  They suggest Q_Ni = 1.75 MeV, Q_Co = 3.73 MeV

    thalf_Ni = 6.1 * seconds_per_day      # half-life of Ni-56
    lambda_Ni = np.log(2) / thalf_Ni
    Q_Ni = 1.75 * erg_per_MeV
    m_Ni = 56 * m_u

    thalf_Co = 77.1*seconds_per_day     # half-life of Co-56
    lambda_Co = np.log(2) / thalf_Co
    Q_Co = 3.73 * erg_per_MeV
    m_Co = 56 * m_u

    def X_Ni(t):
        return ##

    def X_Co(t):
        return ##

    def L(t):
        return ##
    
"""
'\ndef LightCurve(M_exp, X_Ni_0=1.0):\n\n    # note: see Nadyozhin 1994, ApJ, 92, 527 for Q values of\n    # just photons.  They suggest Q_Ni = 1.75 MeV, Q_Co = 3.73 MeV\n\n    thalf_Ni = 6.1 * seconds_per_day      # half-life of Ni-56\n    lambda_Ni = np.log(2) / thalf_Ni\n    Q_Ni = 1.75 * erg_per_MeV\n    m_Ni = 56 * m_u\n\n    thalf_Co = 77.1*seconds_per_day     # half-life of Co-56\n    lambda_Co = np.log(2) / thalf_Co\n    Q_Co = 3.73 * erg_per_MeV\n    m_Co = 56 * m_u\n\n    def X_Ni(t):\n        return ##\n\n    def X_Co(t):\n        return ##\n\n    def L(t):\n        return ##\n    \n'

Lightcurve from \(1~M_\odot\) of \({}^{56}\mathrm{Ni}\)#

# initial condition
M_exp = 1.0*solar_mass

"""
lc = LightCurve(##
t = np.logspace(np.log10(1.0), np.log10(3.e7), 100)

fig, ax = plt.subplots()
ax.plot(t / seconds_per_day, ##
"""
'\nlc = LightCurve(##\nt = np.logspace(np.log10(1.0), np.log10(3.e7), 100)\n\nfig, ax = plt.subplots()\nax.plot(t / seconds_per_day, ##\n'

Different masses of \({}^{56}\mathrm{Ni}\)#

t = np.logspace(np.log10(1.0), np.log10(3.e7), 100)

## create the same plot for M_Ni = [0.4,0.6,0.8,1.0,1.2]Msun

Describe the trend you observe from the plot above. What is the result of changing the total Nickel mass?