40K antineutrino spectrum

Evaluation of electron antineutrino spectrum in beta-minus decay of 40K.

Units: mass and energy in keV

References listed at the end.

Ondřej Šrámek, Charles University, ondrej.sramek@gmail.com, 3 Oct 2017

In [1]:
# inputs
mass_k40 = 39.963998166 # u [NIST Atomic Weights]
mass_ca40 = 39.962590863 # u [NIST Atomic Weights]
MeV_per_u = 931.4940954 # [NIST http://physics.nist.gov/cgi-bin/cuu/Value?muc2mev]
Z = 19 # atomic number of K
A = 40 # mass number of 40K
dE = 1.0 # keV energy discretization
me = 510.9989461 # keV mass of electron [NIST http://physics.nist.gov/cgi-bin/cuu/Value?mec2mev]
alpha = 7.2973525664e-3 # fine-structure constant [NIST https://physics.nist.gov/cgi-bin/cuu/Value?alph]

Q-value from nuclide masses

In [2]:
Q = (mass_k40 - mass_ca40) * MeV_per_u * 1000 # keV Q-value
print ('Q =', Q, 'keV')
Q = 1310.894434939248 keV

Branching ratio of 40K beta- decay from Renne et al (2011)

In [3]:
lambda_beta = 4.9548e-10 # [Renne et al 2011]
lambda_ec = 5.757e-11 # [Renne et al 2011]
lambda1 = lambda_beta + lambda_ec
branch = lambda_beta/lambda1
print ('branching =', branch)
branching = 0.8959045294277189
In [4]:
import numpy as np
import scipy.special
from math import sqrt, pi, exp
import matplotlib.pyplot as plt
In [5]:
plt.rcParams['figure.figsize'] = (8,5)

Following Dye (2012), section 2

dye2012eqns3to5

In [6]:
E = np.arange(0.5*dE, Q, dE)
we = Q + me - E                   # Dye eqn. 3a
pe = np.sqrt(we**2 - me**2)       # Dye eqn. 3b
gamma = sqrt(1-alpha**2*(Z+1)**2) # Dye eqn. 5a
eta = alpha*(Z+1)*we/pe           # Dye eqn. 5b
dndE_dye = we * E**2 * pe**(2*gamma-1) * np.exp(pi*eta) * abs( scipy.special.gamma( gamma + eta*1j ) )**2

Remains to normalize the spectrum. Normalizing to branching ratio of Renne et al 2011 (=0.8959). Integrating by simple trapezoidal rule.

In [7]:
sum = np.sum(dndE_dye) * dE
dndE_dye = dndE_dye / sum * branch

Following Enomoto (2005), eqns. 2.5–2.7

enomoto2005phdeqns5to7

In [8]:
Emax = Q + me
Zprime = Z+1 # assuming Z is Z of daughter
gamma = sqrt(1-(alpha*Zprime)**2) # Enomoto eqn. 2.7
R = 0.426 * alpha * A**(1/3) # Enomoto eqn. 2.7
Ee = Emax - E # Enomoto eqn. 2.6
W = Ee/me # Enomoto eqn. 2.7
y = alpha * Zprime * W/np.sqrt(W**2-1) # Enomoto eqn. 2.7
F = 2*(1+gamma) * (2*np.sqrt(W**2-1)*R)**(2*(gamma-1)) * np.exp(pi*y) * abs( scipy.special.gamma( gamma + y*1j ) )**2 / abs( scipy.special.gamma( 2*gamma+1 ) )**2 # Enomoto eq. 2.7
dndE_eno = F * (Emax-Ee)**2 * np.sqrt(Ee**2-me**2) * Ee # Enomoto eqn. 2.5, up to a scalar factor

# normalize spectrum
sum = np.sum(dndE_eno) * dE
dndE_eno = dndE_eno / sum * branch

Difference between Dye and Enomoto

In [9]:
plt.plot(E, (dndE_dye-dndE_eno)*2/(dndE_dye+dndE_eno) )
plt.xlabel('E (keV)')
plt.ylabel('Relative difference in dn/dE')
plt.title('Relative difference Dye–Enomoto: 2*(D-E)/(D+E)')
plt.show()

That is, Dye and Enomoto formulas give identical spectrum.

Assuming Z in Enomoto's eqns. 2.7 is Z of daughter. The relative differences are at the 10-14 level at most, which is numerical noise.

Cross check against Enomoto's file

Plot spectrum and compare to Enomoto's spectrum (datafile AntineutrinoSpectrum_40K.knt from Enomoto's Tohoku website).

In [10]:
# read-in Enomoto's spectrum
x, y = np.loadtxt('AntineutrinoSpectrum_40K.knt', usecols=(0, 1), unpack=True)
dx = 1.

plt.semilogy(E, dndE_dye, color="C0") # evaluation after Dye
#plt.semilogy(E, dndE_eno, '--', color="C2") # evaluation after Enomoto
plt.semilogy(x, y, color="C1")    # Enomoto's file
plt.text(600, 8e-5, 'My evaluation', color="C0")
#plt.text(600, 4e-5, 'My evaluation after Enomoto', color="C2")
plt.text(600, 4e-5, 'Enomoto data file', color="C1")
plt.xlim(0,1400)
plt.ylim(1e-6,2e-3)
plt.xlabel('E (keV)')
plt.ylabel('dn/dE (1/keV/decay)')
plt.show()

Clearly we are not done yet!

Third unique forbidden transition correction

Enomoto (2006): "...the 3rd unique forbidden transition formula is used for 40K..."

[Note that the forbidden transition nature of 40K decay is not mentioned in Enomoto's PhD nor in Dye (2012). Fiorentini et al (2007) talk about it, see footnote 15 on page 126.]

From Table 15.1 of http://www.umich.edu/~ners311/CourseLibrary/bookchapter15.pdf, I use the following correction or a multiplicative shape factor $S_L(p_e,p_\nu)$, where $p_e$ and $p_\nu$ are the momenta of the electron and the antineutrino, for the third unique forbidden transition:

$$S_L=\frac{p_e^6 + 7p_e^4p_\nu^2 + 7p_e^2p_\nu^4 + p_\nu^6}{(m_ec)^6}$$

(I will say more about this at the end...)

In [11]:
corr = pe**6 + 7*pe**4*E**2 + 7*pe**2*E**4 + E**6 # cause in these units p_nu = E_nu
plt.semilogy(E, corr)
plt.xlabel('E (keV)')
plt.ylabel('$S^L(p,q)$')
plt.title('3rd unique forbidden transition shape factor (correction)')
plt.show()
In [12]:
dndE_dye2 = dndE_dye * corr

# normalize
sum = np.sum(dndE_dye2) * dE
dndE_dye2 = dndE_dye2 / sum * branch

plt.semilogy(E, dndE_dye2)
plt.semilogy(x, y, '--') # Enomoto's file
plt.text(600, 8e-5, 'My evaluation', color="C0")
plt.text(600, 4e-5, 'Enomoto data file', color="C1")
plt.xlim(0,1400)
plt.ylim(1e-6,2e-3)
plt.xlabel('E (keV)')
plt.ylabel('dn/dE (1/keV/decay)')
plt.show()

Relative differences from Enomoto's file

In [13]:
x2 = x[0:E.size]
y2 = y[0:E.size]

dd = (dndE_dye2-y2)/y2

plt.plot(E, dd)
plt.xlabel('E (keV)')
plt.ylabel('Relative difference in dn/dE')
plt.title('Relative difference (me-file)/file')
plt.show()

plt.plot(E, dd)
plt.xlabel('E (keV)')
plt.ylabel('Relative difference in dn/dE')
plt.title('Relative difference (me-file)/file, energy range 0–40 keV')
plt.xlim(0,40)
plt.show()

plt.plot(E, dd)
plt.xlabel('E (keV)')
plt.ylabel('Relative difference in dn/dE')
plt.title('Relative difference (me-file)/file, energy range >40 keV]')
plt.xlim(40,1400)
plt.ylim(-0.004,0.004)
plt.show()

I don't know where this discrepancy comes from. Severe at the first few keV.

Average antineutrino energy

In [14]:
print ("My evaluation:")
print ('    neutrinos per decay =', np.sum(dndE_dye2)*dE, '(this is just the normalization)')
print ('    average neutrino energy =', np.sum(dndE_dye2*E)*dE, 'keV')
print()

print ("Enomoto's file, before renormalization to my branching:")
print ('    neutrinos per decay =', np.sum(y2)*dx)
print ('    average neutrino energy =', np.sum(y2*x2)*dx, 'keV')
print()

# normalize Enomoto file to our branching
sum = np.sum(y2) * dx
y3 = y2 / sum * branch
print ("Enomoto's file, after renormalization to my branching:")
print ('    neutrinos per decay =', np.sum(y3)*dx)
print ('    average neutrino energy =', np.sum(y3*x2)*dx, 'keV')
My evaluation:
    neutrinos per decay = 0.8959045294277189 (this is just the normalization)
    average neutrino energy = 647.6090763091524 keV

Enomoto's file, before renormalization to my branching:
    neutrinos per decay = 0.89276544580554
    average neutrino energy = 645.3878266393695 keV

Enomoto's file, after renormalization to my branching:
    neutrinos per decay = 0.8959045294277189
    average neutrino energy = 647.6570972143852 keV

Ratio $\overline E_\nu/Q$

In [15]:
Enu_ave = np.sum(dndE_dye2*E)*dE
print ("Enu/Q =", Enu_ave/Q)
Enu/Q = 0.4940207686053418

Main results

  1. We can recover the 40K antineutrino spectrum shape of Enomoto.

  2. The change in branching results in change of mean antineutrino energy from 645.4 keV (Enomoto's file) to 647.6 keV.

  3. On average about one half of the $Q$ is carried away by the antineutrino.

Some notes on spectrum shape and forbidden transitions

[This will be brief and incomplete. More insight in relevant textbooks, e.g., Williams (1991). Other sources needed.]

In 1934 Fermi presented the first quantitative theory of β-decay. It validated Pauli's hypothesis and was the first serious extension of the concept of creating (massless) photons to creating (massive) particles.

Fermi's Golde Rule #2 gives the rate $\omega$ of transition from initial state $i$ to final state $f$

$$\omega = \frac{2\pi}{\hbar}\big|M_{if}\big|^2\frac{dN}{dE}$$

where $E$ is the total energy available in the center-of-mass of the final system and $N$ is the number of final states available with energy up to $E$; $M_{if}$ is the matrix element of the $i\rightarrow f$ interaction Hamiltonian.

Fermi then used these assumptions:

  1. $M_{if}=G_F/V$ where $G_F$ in the Fermi coupling constant and $V$ is the volume in which particle wave-functions are normalized.
  2. Effects of angular momentum consevation and parity non-conservation can be neglected. (That is, excludes forbidden transitions.)
  3. $Q$ is much smaller than the rest mass energies of the parent and daughter nuclei. That is, the daughter nucleus does not recoil, electron and neutrino momenta are uncorrelated.
  4. Neutrinos are mass-less.

Total energy $E$ is the sum of energy of the electron $E_e=m_ec^2+T_e$ and of the antineutrino $E_\nu=T_\nu$, where $T_e$ and $T_\nu$ are kinetic energies of the particles (neutrino mass is neglected, nuclei assumed stationary):

$$E = E_e + E_\nu = m_ec^2 + T_e + T_\nu = m_ec^2 + Q$$

For the massive electron, the relativistic relation between its energy $E_e$ and momentum $p_e$ is

$$E_e^2 = p_e^2c^2+m_e^2c^4,$$

and for the mass-less antineutrino it simplifies to

$$E_\nu = p_\nu c.$$

We observe that the electron momentum spectrum is continuous, therefore we want to write $d\omega/dp_e$, that is, we need to find

$$\frac{d^2N}{dp_edE} = ~?$$

It becomes an exercise on quantum-mechanical density of states. For a fixed $p_e$, the number of electron states $N_e$ having momentum $<p_e$ is

$$N_e = \frac{4\pi p_e^3V}{3(2\pi\hbar)^3},$$

(counting the number of discrete wavefunctions of wavelength greater than $2\pi\hbar/p_e$ that fit in a box of volume $V$) and the same can be written for antineutrino states, $N_\nu$ in terms of $p_\nu$. The electron and antineutrino states can be combined independently into a momentum state, therefore

$$N = N_eN_\nu = \frac{16\pi^2 p_e^3 p_\nu^3 V^2}{9(2\pi\hbar)^6}$$

Now we differentiate $N(p_e,p_\nu)$ with respect to $p_e$ and $E$ to find $d\omega/dp_e$ at a given $p_e$ (therefore $dE=dE_\nu=cdp_\nu$), to find

$$\frac{d^2N}{dp_edE} = \frac{16\pi^2 p_e^2 p_\nu^2 V^2}{(2\pi\hbar)^6c}.$$

Putting this together with the Golder Rule #2, we get the electron spectrum

$$\frac{d\omega}{dp_e} = \frac{G_F^2 p_e^2 p_\nu^2}{2\pi^3\hbar^7c},$$

which can be written in terms of the total energy $E$ and electron energy $E_e$ (using the energy–momentum relations above) as

$$\frac{d\omega}{dE_e} = \frac{G_F^2}{2\pi^3\hbar^7c^6} \sqrt{E_e^2-m_e^2c^4} E_e (E-E_e)^2.$$

Another factor $\big|M\big|$ may appear in the equation which accounts for the contribution of the initial and final nuclear state, thus relaxing the first Fermi's assumption (of $M_{if}=G_F/V$).

Coulombic interaction of daughter nucleus with emitted electron introduces a multiplicative factor (function), so called Fermi Function, $F(Z',p_e)$, which depends on the charge $Z'$ of the daughter nucleus and the electron momentum $p_e$.

Allowed transitions are those where the electron and antineutrino carry away a total of either one or zero units of angular momentum (spin triplet state = Gamow–Teller transitions, or spin singlet state = Fermi transitions). Forbidden transitions are those where the change in nuclear spin is greater than 0 or 1, or where nuclear parity changes from parent to daughter. The first forbidden transitions have 0, 1 or 2 units change of nuclear spin and a change in nuclear parity. The second forbidden transitions have 1, 2 or 3 units change of nuclear spin and no change in nuclear parity. The 40K β-decay is a third forbidden transition (-4 units change in nuclear spin and a change in parity):

40K beta-decay scheme

In the case of forbidden transitions, the angular momentum barrier depresses the rate and the spectrum shape changes. An additional multiplicative factor is present, the shape factor $S_L(p_e,p_\nu)$, where $p_e$ and $p_\nu$ are the momenta of the electron and the antineutrino. The shape factor can be calculated from relativistic quantum mechanics in the special case of unique transitions. Unique transitions are those where the angular momentum vector and the two lepton spins are all aligned. For the first three unique forbidden transition we have:

bielajew_bookchapter15_tab15.1

$$S_3(p_e,p_\nu) = \frac{p_e^6 + 7p_e^4p_\nu^2 + 7p_e^2p_\nu^4 + p_\nu^6}{(m_ec)^6}$$

Overall, we get

$$\frac{d\omega}{dE_e} = \frac{G_F^2\big|M\big|^2}{2\pi^3\hbar^7c^6} S_L(p_e,p_\nu) F(Z',E_e) \sqrt{E_e^2-m_e^2c^4} E_e (E-E_e)^2,$$

which compares to relations in Enomoto (2005), Dye (2012), Fiorentini et al (2007), etc...

References