Neutrinoless double beta decay

software

on github: pybbsens

(ignore 2vbb background shape)

calculation steps

draw spectrum

• 2vbb energy spectrum shape
• detector resolution convolution with spectrum

define sensitivity

• signal-to-noise = 1
• draw graph: R(T(2v)/T(0v)) - resolution
• draw graph: T(0v) - resolution

calculate mbb

• phase space factor
• nuclear matrix elements: IBM, QRPA, ISM ...
• draw graph: mbb - resolution

normal hierarchy or inverted hierarchy

• different hierarchy and mbb - resolution
• neutrino mixing matrix

data acquisition

• det-mass * time
• redefine sensitivity: poisson dist. theory
• add background from other processes

Two neutrino double beta decay spectrum:

$\dfrac{dN}{dK} \sim K(T_0-K)^5(1 + 2K + \dfrac{4K^2}{3} + \dfrac{K^3}{3} + \dfrac{K^4}{30})$

%matplotlib inline

import numpy as np  
from matplotlib import pyplot as plt

# Q-value of 76Ge double beta decay
e0 = 2.039

# total energy of two electrons
eres = 0.0001  
et = np.arange(0., 2.5, eres)  
# 2vbb
s2v = et * (e0 - et)**5 * (1. + 2.*et + 4./3.*et**2 + 1./3.*et**3 + 1./30.*et**4)  
s2v[np.where(et > e0)] = 0.0  
s2v = s2v / sum(s2v)  # normalize spectrum  
# 0vbb
s0v = np.zeros(et.shape)  
s0v[np.where(et==e0)] = 1.0

# draw
plt.plot(et, s2v*1e5, label='2vbb') # times 10^5 for illutration  
plt.plot(et, s0v, label='0vbb')  
plt.legend()  
plt.xlabel('Energy (MeV)')  
plt.ylabel('a.u.')  
plt.show()  

Convolute with a gaussian function, set FWHM = 2 keV, conresponding to a resolution of 0.1 % @2.039 MeV.

from math import *

# define resolution
FWHM = 0.001 * 2  
sigma = FWHM / 2.355  
# construct gaussian signal
gx = np.arange(-5*sigma, 5*sigma, eres)  
gaussian = np.exp(-(gx/sigma)**2 / 2.)

# convolution
s2v_conv = np.convolve(s2v, gaussian, mode="same")  
s0v_conv = np.convolve(s0v, gaussian, mode="same")

# assume half-life ratio 2vbb / 0vbb
ratio = 1e14

# draw 
region = np.where(abs(et - e0) <= 5*sigma )  
plt.plot(et[region], ratio * s2v_conv[region], label='2vbb')  
plt.plot(et[region], s0v_conv[region], label='0vbb')  
plt.plot(et[region], ratio * s2v_conv[region] + s0v_conv[region], label='total')  
plt.legend()  
plt.xlabel('Energy (MeV)')  
plt.ylabel('a.u.')  
plt.show()  

Define Sensitivity: $ R = \dfrac{N{0\nu\beta\beta}}{N{2\nu\beta\beta}} = 1 $ where $N$ corresponds to counts within FWHM centered at $Q_{\beta\beta}$

# define functions to simplify calculation

def convSpec(s, FWHM):  
    # define gaussian
    sigma = FWHM / (2.*sqrt(2.*log(2.)))
    # construct gaussian signal
    gx = np.arange(-5*sigma, 5*sigma, eres)
    gaussian = np.exp(-(gx/sigma)**2 / 2.)
    # convolution
    return np.convolve(s, gaussian, mode="same")

def countsROI(s, FWHM, Q):  
    return sum(s[np.where(abs(et-Q) <= FWHM/2.)])

def getHalflifeRatio(n0v, n2v):  
    return n0v / n2v

def getHalflife(ratio, T2v):  
    return ratio * T2v


# spectrum from last section
# Q-value of 76Ge double beta decay
e0 = 2.039  
T2v = 1.83e21  # from GERDA 2013  
# total energy of two electrons
eres = 0.0001  
et = np.arange(0., 2.5, eres)  
# 2vbb
s2v = et * (e0 - et)**5 * (1. + 2.*et + 4./3.*et**2 + 1./3.*et**3 + 1./30.*et**4)  
s2v[np.where(et > e0)] = 0.0  
s2v = s2v / sum(s2v)  # normalize spectrum  
# 0vbb
s0v = np.zeros(et.shape)  
s0v[np.where(et==e0)] = 1.0


# calculate ratio and half-life under different resolution

resolutions = np.arange(0.1, 10., 0.1) / 100.  
FWHMs = resolutions * e0  
ratios = []  
taus = []  
for FWHM in FWHMs:  
    s2v_conv = convSpec(s2v, FWHM)
    s0v_conv = convSpec(s0v, FWHM)
    n2v = countsROI(s2v_conv, FWHM, e0)
    n0v = countsROI(s0v_conv, FWHM, e0)
    r = getHalflifeRatio(n0v, n2v)
    t = getHalflife(r, T2v)
    ratios.append(r)
    taus.append(t)

# draw graphs
plt.plot(resolutions*100., np.array(ratios))  
plt.xlabel(r'FWHM resolution at $Q_{\beta\beta}$ (%)')  
plt.ylabel(r'Ratio: $T^{0\nu}_{1/2}/T^{2\nu}_{1/2}