Application of integral equations to neutrino mass searches in beta decay
1 Application of integral equations to neutrino mass searches in beta decay
Thomas M. Semkow and Xin Li
Wadsworth Center, New York State Department of Health, Albany, NY 12201, USA and School of Public Health, University at Albany, State University of New York, Rensselaer, NY 12144, USA Email: [email protected]
Abstract
A new mathematical method for elucidating neutrino mass from beta decay is studied. It is based upon the solutions of transformed Fredholm and Volterra integral equations. In principle, theoretical beta-particle spectra can consist of several neutrino-mass eigenvalues. Integration of the theoretical beta spectrum with a normalized instrumental response function results in the Fredholm integral equation of the first kind. This equation is transformed in such a way that the solution of it is a superposition of the Heaviside step-functions, one for each neutrino mass eigenvalue. A series expansion leading to matrix linear equations is then derived to solve the transformed Fredholm equation. Another approach is derived when the theoretical beta spectrum is obtained by a separate deconvolution of the observed spectrum. It is then proven that the transformed Fredholm equation reduces to the Abel integral equation. The Abel equation has a general integral solution, which is proven in this work by using a specific function for the beta spectrum. As an example, a numerical solution of the Abel integral equation is also provided, which has a fractional sensitivity of about 10 ― for subtle neutrino eigenvalue searches, and can distinguish from experimental beta-spectrum discrepancies, such as shape and energy nonlinearities. Keywords: Fredholm equation, Abel equation, Heaviside function, heavy neutrino, shape correction, energy nonlinearity
1. Introduction
We start by outlining the physics of neutrino-mass searches. Since the neutrino was conceptually proposed by Pauli in 1930 and given its name by Fermi, followed by his development of the beta decay theory in 1933-1934 [1-3], studies of this elusive neutral particle have never ceased. Whether or not neutrino bears a rest mass involves the extension of the Standard Model, the composition of Dark Matter in the Universe, and the asymmetry between matter and antimatter [4]. An indirect evidence for neutrino mass has been provided by neutrino-oscillation experiments, which give us the differences between the masses of neutrinos with different flavors [5]. Three approaches are generally used to determine the absolute masses of neutrinos: investigation on the fluctuations of cosmic microwave background and large-scale structure, measurement of neutrino-less double-beta decay rate, and direct determinations of neutrino mass [6]. Currently, the first approach constraints the sum of three neutrino masses to less than 0.23 eV [7]. The second approach relies on the hypothesis that the neutrino is a Majorana particle, where a lower bound of the 0 νββ decay half-life of 5.3×10 years has been set [8], whereas an effective neutrino mass of 0.1-0.9 eV was claimed [9]. The third approach has been focusing on a precise measurement of neutrino kinetic energy, as well as an endpoint region of a beta particle spectrum. The former requires a very strong neutrino source: the analysis of the neutrino energy spectrum as a function of arrival time from SN1987A supernova data yielded an upper limit of neutrino mass at 5.7 eV [10]. In case of the latter, the Mainz and Troitsk experiments for H beta decay yielded the upper limits of electron antineutrino mass at 2.3 eV [11] and 2.05 eV [12], respectively, while the KATRIN experiment is expected to have a sensitivity down to 0.2 eV [13]. The
Re beta-decay bolometric experiments resulted in upper limits for neutrino mass of 26 eV [14] and 15 eV [15]. The
Ho (electron capture) neutrino experiments yielded higher upper limits of 225 eV [16] and 490 eV [17] and efforts are in progress to lower these limits in the ECHO experiment [6]. In the science of direct search for massive neutrino in beta decay away from the endpoint energy, there were claims of the discovery of 17-keV neutrino in H [18], C [19], and S [20,21]. The 17-keV component appeared to be present as a discontinuity in beta spectra at a fractional contribution of approximately 0.01. These effects were later explained as electrons scattering from collimators [22,23]. The remainder of the paper concentrates on the neutrino mass searches in β – decay (hereinafter referred to as beta decay). The electron antineutrino, ν(cid:3364) (cid:2915) (hereinafter referred to as neutrino) is emitted in the weak-interaction process of neutron decay inside the nucleus: n ⟶ p (cid:2878) (cid:3397) β (cid:2879) (cid:3397) ν(cid:3364) (cid:2915) . The emitted β -particle (electron) kinetic-energy spectrum is a continuous function, owing to the energy sharing between the beta particle and the neutrino. Possible neutrino mass (cid:1839) in the beta spectrum is included in the factor: (cid:3493)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) (cid:3398) (cid:1839) (cid:2870) , where (cid:1843) is the nuclear-recoil corrected (cid:1843) -value of beta decay and (cid:1846) is the emitted electron kinetic energy. The maximum beta energy (called the endpoint energy) is equal to (cid:1843) or (cid:1843) (cid:3398) (cid:1839) , when the neutrino mass is zero or positive, respectively. In principle, there can be several mixed neutrino eigenvalues. The details of the theoretical beta spectra are described in Section 2. The emitted beta spectrum is convoluted with the instrumental response to yield the observed beta spectrum. The existing methods of identification of the neutrino mass in beta decay, whether close to or away from the endpoint energy, employ two approaches: i) convolution of the theoretical beta spectrum with the normalized instrumental response function and comparison of the convolution with the experimental spectrum, and ii) deconvolution of the observed spectrum and comparison with the theoretical spectrum. Statistical measures, such as (cid:2031) (cid:2870) minimization or Bayesian likelihood, are used for the above-mentioned comparisons. Statistical measures can be non-specific, i.e ., there may be several factors other than the neutrino mass, which affect them. Therefore, we aimed in this work at developing mathematical method based on transformational properties, rather than statistical measures, to elucidate neutrino mass in beta decay. In Section 3, we derive a transformed Fredholm equation of the first kind, which includes both theoretical beta spectrum and instrumental response function, and show that the solution to it is a superposition of Heaviside step-functions with an abscissa of (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) , one for each possible neutrino mass eigenvalue. The neutrino mass can then be identified by the abscissa value at the raise of the step function, whereas its eigenvalue contribution is shown to be proportional to the ordinate value past the step. In Section 4, we outline a possible solution to the transformed Fredholm equation using series expansion and linear equations in a matrix notation. If method ii) of deconvoluting the beta spectrum first, before comparison with the theory, is employed, then the transformed Fredholm equation from Section 3 is shown in Section 5 to reduce to the Abel integral equation, also known as the Abel transformation. The Abel integral equation has a known integral solution. We prove in Section 5 that, if the function undergoing Abel transformation has a form proportional to the beta spectrum, then the solution to the Abel equation is a Heaviside step-function. In Section 6, we provide a numerical example of a solution of the Abel equation, when the beta spectrum contains several neutrino-mass eigenvalues. We also describe in detail how the Abel solution behaves in the presence of small experimental discrepancies, which can influence neutrino mass detection, such as β shape correction or small nonlinearities in the measured β -energy scale. Finally, in Section 7 we summarize the advantages and limitations of the proposed new mathematical method of elucidating neutrino mass in beta decay, and propose future directions of this work.
2. Beta-decay theory
There exist established formalisms for calculating allowed and superallowed beta spectra, which differ in some approaches [6,24,25]. For the purpose of this work, we define the beta-energy spectrum as follows: (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) (cid:3404) (cid:1827)(cid:1830)(cid:4666)(cid:1846)(cid:4667)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667)(cid:3493)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) (cid:3398) (cid:1839) (cid:2870) , (cid:4666)1(cid:4667) where (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) (cid:3404) (cid:1856)(cid:1840) (cid:1856)(cid:1846)⁄ represents a distribution of beta particles with energies between and (cid:1846) and (cid:1846) (cid:3397) (cid:1856)(cid:1846) , factor (cid:1827) contains the quantum-transition matrix element and an overall normalization, whereas (cid:1843) , (cid:1846) , and (cid:1839) have been defined in Section 1. Factor (cid:1830)(cid:4666)(cid:1846)(cid:4667) is defined as follows: (cid:1830)(cid:4666)(cid:1846)(cid:4667) (cid:3404) (cid:1868)(cid:1831)(cid:1832) (cid:2924) (cid:1832) (cid:2928) (cid:1829) (cid:2928) (cid:1845) , (cid:4666)2(cid:4667) where (cid:1868) is relativistic momentum, (cid:1831) is total (relativistic) energy, (cid:1832) (cid:2924) is a non-relativistic Fermi function, (cid:1832) (cid:2928) is a relativistic correction to the Fermi function [26], and (cid:1829) (cid:2928) is a radiative correction [24]. (cid:1845) is a shape correction [24], which corrects for small smooth discrepancies between the theoretical and measured spectra at low energies, is assumed to be a parabola with an argument of (cid:1843) (cid:3398) (cid:1846) . Also, a screening correction [27,28] is applied to (cid:1830)(cid:4666)(cid:1846)(cid:4667) . Additional terms involving the daughter’s final excitation states in low beta-energy emitters [6], as well as a possible Lorentz-invariance violation correction [29], have not been included. The examples of beta spectra from H decay (cid:4666)(cid:1843) (cid:3404) 18590.29 eV(cid:4667) [30] calculated according to Eqs. (1,2) are plotted in Fig. 1. The physical constants and nuclear parameters for the calculations were taken from Ref. [30]. It is seen that if neutrino mass is non-zero, the endpoint beta energy lies below (cid:1843) , while another neutrino mass eigenvalue is evidenced by the discontinuity of the first derivative in the spectrum. If neutrino mass is close to zero and any other eigenvalue has a very small contribution, these effects are difficult to discern from the beta spectrum.
Figure 1. Beta particle spectra from H decay. Spectrum (a): zero neutrino mass; spectrum (b): equal superposition of two neutrino eigenvalues with the masses of 3000 and 5000 eV.
3. Transformed Fredholm integral equation
The observed beta spectrum usually differs from the emitted (theoretical) beta spectrum, given by Eqs. (1,2), owing to an instrumental response. Let (cid:1846) represent the emitted kinetic energy while (cid:1846) (cid:4593) , the observed kinetic energy. The emitted and observed beta spectra are then (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) and (cid:1840)(cid:4662) (cid:4666)(cid:1846) (cid:4593) (cid:4667) , respectively. They are coupled with each other by the response function (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667) through the integral: (cid:1840)(cid:4662) (cid:4666)(cid:1846) (cid:4593) (cid:4667) (cid:3404) (cid:3505) (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667)(cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667)(cid:1856)(cid:1846). (cid:4666)3(cid:4667) (cid:3018)(cid:2879)(cid:3014)(cid:2868)
Equation (3) is a Fredholm integral equation of the first kind [31], where (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667) is a kernel. The variable (cid:1846) (cid:4593) can extend past (cid:1843) (cid:3398) (cid:1839) due to the abovementioned instrumental response. The response function is usually measured separately. Since the emission and detection of beta particles are stochastic processes, the normalized (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667) and (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) can be considered statistical probability density functions (pdf). The integral in Eq. (3) is then a mixture of distributions [32]. If, additionally, (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667) is a function of (cid:1846) (cid:4593) (cid:3398) (cid:1846) , e.g ., a Gaussian, then Eq. (3) is statistically a convolution [33]. The process of solving for (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) is referred to as deconvolution and is discussed in Section 4. The terms: convolution and deconvolution will be used here even if Eq. (3) is not a strictly mathematical convolution. If there existed (cid:1863) neutrino-mass eigenstates with the masses (cid:1839) (cid:3038) and normalization factors (cid:1827) (cid:3038) , then combining Eqs. (1,3) and summing over the eigenstates would result in the observed spectrum: (cid:1840)(cid:4662) (cid:4666)(cid:1846) (cid:4593) (cid:4667) (cid:3404) (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667)(cid:1830)(cid:4666)(cid:1846)(cid:4667)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667)(cid:3495)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) (cid:3398) (cid:1839) (cid:3038)(cid:2870) (cid:1856)(cid:1846). (cid:4666)4(cid:4667) (cid:3018)(cid:2879)(cid:3014) (cid:3286) (cid:2868) Subsequently, we introduce a mass variable (cid:1865) and utilize the Dirac (cid:2012) -function [34]. Then Eq. (4) transforms to: (cid:1840)(cid:4662) (cid:4666)(cid:1846) (cid:4593) (cid:4667) (cid:3404) (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1856)(cid:1865) (cid:2012)(cid:4666)(cid:1865) (cid:3398) (cid:1839) (cid:3038) (cid:4667) (cid:3018)(cid:2868) (cid:3505) (cid:1856)(cid:1846) (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667)(cid:1830)(cid:4666)(cid:1846)(cid:4667)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667)(cid:3493)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) (cid:3398) (cid:1865) (cid:2870) . (cid:4666)5(cid:4667) (cid:3018)(cid:2879)(cid:3040)(cid:2868)
New variables are defined: (cid:1876) (cid:3404) (cid:4666)(cid:1843) (cid:3398) (cid:1846) (cid:4593) (cid:4667) (cid:2870) and (cid:1878) (cid:3404) (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) . By substitution, Eq. (5) results in (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/2 (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1856)(cid:1865) (cid:2012)(cid:4666)(cid:1865) (cid:3398) (cid:1839) (cid:3038) (cid:4667) (cid:3018)(cid:2868) (cid:3505) (cid:1856)(cid:1878) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1878)(cid:4667)(cid:3493)(cid:1878) (cid:3398) (cid:1865) (cid:2870) . (cid:4666)6(cid:4667) (cid:3018) (cid:3118) (cid:3040) (cid:3118)
By using a property of the (cid:2012) -function [34], (cid:2012)(cid:4666)(cid:1865) (cid:3398) (cid:1839) (cid:3038) (cid:4667) (cid:3404) 2(cid:1865) (cid:2012)(cid:4666)(cid:1865) (cid:2870) (cid:3398) (cid:1839) (cid:3038)(cid:2870) (cid:4667) , (cid:4666)7(cid:4667) and substituting of the variables, (cid:1877) (cid:3404) (cid:1865) (cid:2870) , (cid:1853) (cid:3404) (cid:1843) (cid:2870) , and (cid:1853) (cid:3038) (cid:3404) (cid:1839) (cid:3038)(cid:2870) , one obtains from Eq. (6): (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/2 (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1856)(cid:1877) (cid:2012)(cid:4666)(cid:1877) (cid:3398) (cid:1853) (cid:3038) (cid:4667) (cid:3028)(cid:2868) (cid:3505) (cid:1856)(cid:1878) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1878)(cid:4667)(cid:3493)(cid:1878) (cid:3398) (cid:1877). (cid:4666)8(cid:4667) (cid:3028)(cid:3052) Changing the order of integration in Eq. (8), yields (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/2 (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1856)(cid:1878) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1878)(cid:4667) (cid:3028)(cid:2868) (cid:3505) (cid:1856)(cid:1877) (cid:2012)(cid:4666)(cid:1877) (cid:3398) (cid:1853) (cid:3038) (cid:4667)(cid:3493)(cid:1878) (cid:3398) (cid:1877). (cid:4666)9(cid:4667) (cid:3053)(cid:2868)
The inner integral in Eq. (9) can be integrated by parts. We also use the fact that the integral of the (cid:2012) -function on positive argument is the Heaviside step function: (cid:3505) (cid:2012)(cid:4666)(cid:1877) (cid:3398) (cid:1853) (cid:3038) (cid:4667) (cid:1856)(cid:1877) (cid:3404) (cid:1834) (cid:3028) (cid:3286) (cid:4666)(cid:1877)(cid:4667) (cid:3404) (cid:3420)0, (cid:1877) (cid:3409) (cid:1853) (cid:3038) (cid:3038) . (cid:4666)10(cid:4667)
The result is given below: (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/4 (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:3505) (cid:1856)(cid:1878) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1878)(cid:4667) (cid:3028)(cid:2868) (cid:3505) (cid:1856)(cid:1877) (cid:1834) (cid:3028) (cid:3286) (cid:4666)(cid:1877)(cid:4667) (cid:3493)(cid:1878) (cid:3398) (cid:1877)⁄ . (cid:4666)11(cid:4667) (cid:3053)(cid:2868)
The solution to Eq. (11) is the function (cid:1858)(cid:4666)(cid:1877)(cid:4667) (cid:3404) (cid:3533) (cid:1827) (cid:3038)(cid:3038) (cid:1834) (cid:3028) (cid:3286) (cid:4666)(cid:1877)(cid:4667), (cid:4666)12(cid:4667) which is a superposition of the step functions for all possible neutrino mass eigenvalues, scaled by their appropriate contributions. By inserting (12) into (11) one obtains (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/4 (cid:3505) (cid:1856)(cid:1878) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1878)(cid:4667) (cid:3028)(cid:2868) (cid:3505) (cid:1856)(cid:1877) (cid:1858)(cid:4666)(cid:1877)(cid:4667) (cid:3493)(cid:1878) (cid:3398) (cid:1877)⁄ . (cid:4666)13(cid:4667) (cid:3053)(cid:2868)
By reversing the order of integration in Eq. (13) back to that of Eq. (8) and by dividing both sides by (cid:1830)(cid:4666)(cid:1876)(cid:4667) , i.e ., the theoretical factor (cid:1830) given by Eq. (2) evaluated at the observed energy (cid:1846) (cid:4593) , one obtains (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667)(cid:1830)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:3505) (cid:1837)(cid:4666)(cid:1876), (cid:1877)(cid:4667)(cid:1858)(cid:4666)(cid:1877)(cid:4667)(cid:1856)(cid:1877), with (cid:3028)(cid:2868) (cid:4666)14a(cid:4667) (cid:1837)(cid:4666)(cid:1876), (cid:1877)(cid:4667) (cid:3404) 1/4 (cid:3505) (cid:1844)(cid:4666)(cid:1876), (cid:1878)(cid:4667)(cid:3493)(cid:1878) (cid:3398) (cid:1877) (cid:1830)(cid:4666)(cid:1878)(cid:4667)(cid:1830)(cid:4666)(cid:1876)(cid:4667) (cid:1856)(cid:1878). (cid:4666)14b(cid:4667) (cid:3028)(cid:3052) Equation (14a) is the transformed Fredholm equation. The kernel is given by Eq. (14b). It links the observed beta spectrum with the emitted (theoretical) spectrum and with the instrumental response. In the solution function given by Eq. (12), one does not need to assume the presence of neutrino mass or its eigenvalues. They may be revealed by the abscissa and ordinate of any step functions present in the solution, to within the accuracy of the solution, which is dependent on the detailed shape of the observed spectrum. Therefore, the proposed approach interrogates the shape of the observed spectrum using the transformation derived, in contrast to the methods based on statistical measures.
4. Solution of the Fredholm integral equation
There exist algorithms for solving the Fredholm equations given by Eqs. (3) and (14a) [35,36]. If Eq. (3) is a mathematical convolution, then it can also be solved by the Laplace transform method, taking the advantage of the theorem that the Laplace transform of the convolution is a product of Laplace transforms of the convoluting functions [37]. However, Eq. (14a) is not a convolution and we derive below a matrix method to solve it. We take into consideration the fact that the observed (cid:1846) (cid:4593) variable is discrete and incremented at equally spaced intervals (channels) owing to the instrumental digitizing of the spectrum (modern instrumentation can have (cid:1866) (cid:3404) 2 (cid:2869)(cid:2873) or more channels). Consequently, variable (cid:1876) is also digitized. Consider a general integral equation of the first kind, of which Eqs. (3,14a) are the examples: (cid:2038)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:3505) (cid:1837)(cid:4666)(cid:1876), (cid:1877)(cid:4667)(cid:2032)(cid:4666)(cid:1877)(cid:4667) (cid:1856)(cid:1877) . (cid:4666)15(cid:4667) If (cid:1876) is discrete, one has: (cid:2038)(cid:4666)(cid:1876) (cid:3036) (cid:4667) (cid:3404) (cid:3505) (cid:1837)(cid:4666)(cid:1876) (cid:3036) , (cid:1877)(cid:4667)(cid:2032)(cid:4666)(cid:1877)(cid:4667) (cid:1856)(cid:1877) , (cid:1861) (cid:3404) 1, … , (cid:1866). (cid:4666)16(cid:4667) Let us assume a trial solution in a form of a series: (cid:2032)(cid:4666)(cid:1877)(cid:4667) (cid:3404) (cid:3533) (cid:1837)(cid:3435)(cid:1876) (cid:3037) , (cid:1877)(cid:3439) (cid:3041)(cid:3037)(cid:2880)(cid:2869) (cid:1855) (cid:3037) . (cid:4666)17(cid:4667)
By inserting the trial solution into Eq. (16), one obtains a set of linear equations: (cid:2038)(cid:4666)(cid:1876) (cid:3036) (cid:4667) (cid:3404) (cid:3533) (cid:1837) (cid:3036)(cid:3037)(cid:3041)(cid:3037)(cid:2880)(cid:2869) (cid:1855) (cid:3037) , (cid:1861) (cid:3404) 1, … , (cid:1866), (cid:4666)18a(cid:4667) with (cid:1837) (cid:3036)(cid:3037) (cid:3404) (cid:1516) (cid:1837)(cid:4666)(cid:1876) (cid:3036) , (cid:1877)(cid:4667)(cid:1837)(cid:3435)(cid:1876) (cid:3037) , (cid:1877)(cid:3439) (cid:1856)(cid:1877). (cid:4666)18b(cid:4667)
In matrix notation, Eq. (18a) reads: (cid:2740) (cid:3404) (cid:1785)(cid:1777) , (cid:4666)19(cid:4667) where (cid:2740) and (cid:1777) are vectors and (cid:1785) is a symmetric square matrix. By solving the linear equations for (cid:1777) , the solution (17), continuous in (cid:1877) variable ( i.e ., in (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) ) can, in principle, be obtained.
5. Abel integral equation and its solution
In Section 3, the general treatment for elucidating of neutrino mass in beta decay was presented, leading to the transformed Fredholm equation (14a). This equation can be considerably simplified if the observed beta spectrum is first deconvoluted using Eq. (3). The limit to the transformed Fredholm equation (14a) can then be obtained by using the fact that, if the spectrum was deconvoluted, then for the purpose of (14a), (cid:1844)(cid:4666)(cid:1846) (cid:4593) , (cid:1846)(cid:4667) (cid:3404) (cid:2012)(cid:4666)(cid:1846) (cid:4593) (cid:3398) (cid:1846)(cid:4667) (cid:3404) (cid:2012)(cid:4670)(cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:3398) (cid:4666)(cid:1843) (cid:3398) (cid:1846) (cid:4593) (cid:4667)(cid:4671) (cid:3404) (cid:2012)(cid:3435)√(cid:1878) (cid:3398) √(cid:1876)(cid:3439) (cid:3404) 2√(cid:1878)(cid:2012)(cid:4666)(cid:1878) (cid:3398) (cid:1876)(cid:4667), where we used a version of Eq. (7). Then, starting from Eq. (13), (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667) (cid:3404) 1/4 (cid:3505) (cid:1856)(cid:1878) (cid:2012)(cid:4666)(cid:1878) (cid:3398) (cid:1876)(cid:4667)2√(cid:1878)(cid:1830)(cid:4666)(cid:1878)(cid:4667) (cid:3505) (cid:1856)(cid:1877) (cid:1858)(cid:4666)(cid:1877)(cid:4667) (cid:3493)(cid:1878) (cid:3398) (cid:1877) ⁄ (cid:3053)(cid:2868)(cid:3028)(cid:2868) (cid:3404) 1/2 (cid:1830)(cid:4666)(cid:1876)(cid:4667)√(cid:1876) (cid:3505) (cid:1858)(cid:4666)(cid:1877)(cid:4667) (cid:3493)(cid:1876) (cid:3398) (cid:1877)⁄ (cid:3051)(cid:2868) (cid:1856)(cid:1877). (cid:4666)20(cid:4667)
By defining a new function: (cid:1859)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:1840)(cid:4662) (cid:4666)(cid:1876)(cid:4667)(cid:1830)(cid:4666)(cid:1876)(cid:4667)√(cid:1876) , (cid:4666)21(cid:4667) we obtain from Eq. (20) (cid:1859)(cid:4666)(cid:1876)(cid:4667) (cid:3404) 12 (cid:3505) (cid:1858)(cid:4666)(cid:1877)(cid:4667)(cid:3493)(cid:1876) (cid:3398) (cid:1877) (cid:1856)(cid:1877). (cid:4666)22(cid:4667) (cid:3051)(cid:2868)
Equation (22) is a special case of the Abel transformation, also called the Abel integral equation [31,38], which is a special case of the Volterra integral equation of the first kind. The general form of the Abel equation is (cid:1859)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:3505)(cid:4666)(cid:1876) (cid:3398) (cid:1877)(cid:4667) (cid:2879)(cid:3080) (cid:1858)(cid:4666)(cid:1877)(cid:4667) (cid:1856)(cid:1877), 0 (cid:3407) (cid:2009) (cid:3407) 1, (cid:4666)23(cid:4667) (cid:3051)(cid:2868) with a solution of (cid:1858)(cid:4666)(cid:1876)(cid:4667) (cid:3404) sin (cid:2024)(cid:2009)(cid:2024) (cid:1856)(cid:1856)(cid:1876) (cid:3505)(cid:4666)(cid:1876) (cid:3398) (cid:1877)(cid:4667) (cid:3080)(cid:2879)(cid:2869) (cid:1859)(cid:4666)(cid:1877)(cid:4667) (cid:1856)(cid:1877). (cid:4666)24(cid:4667) (cid:3051)(cid:2868)
Therefore, the solution of Eq. (22) for (cid:2009) (cid:3404) 1/2 is (cid:1858)(cid:4666)(cid:1876)(cid:4667) (cid:3404) 2(cid:2024) (cid:1856)(cid:1856)(cid:1876) (cid:3505) (cid:1859)(cid:4666)(cid:1877)(cid:4667)(cid:3493)(cid:1876) (cid:3398) (cid:1877) (cid:1856)(cid:1877). (cid:4666)25(cid:4667) (cid:3051)(cid:2868) By using Eq. (1) for a single neutrino-mass eigenvalue (cid:1839) , and setting (cid:1877) (cid:3404) (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) as well as (cid:1854) (cid:3404) (cid:1839) (cid:2870) , we obtain from Eq. (21) (cid:1859)(cid:4666)(cid:1877)(cid:4667) (cid:3404) (cid:1827)(cid:3493)(cid:1877) (cid:3398) (cid:1854). (cid:4666)26(cid:4667)
By inserting Eq. (26) to Eq. (25), one obtains (cid:1858)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:3422) 0, (cid:1877) (cid:3409) (cid:1854) 2(cid:1827)(cid:2024) (cid:1856)(cid:1856)(cid:1876) (cid:3505) (cid:3493)(cid:1877) (cid:3398) (cid:1854)(cid:3493)(cid:1876) (cid:3398) (cid:1877) (cid:1856)(cid:1877) , (cid:1876) (cid:3408) (cid:1877) (cid:3408) (cid:1854). (cid:3051)(cid:3029) (cid:3505) (cid:3493)(cid:1877) (cid:3398) (cid:1854)(cid:3493)(cid:1876) (cid:3398) (cid:1877) (cid:1856)(cid:1877) (cid:2912)(cid:2935) (cid:2926)(cid:2911)(cid:2928)(cid:2930)(cid:2929) (cid:4657)(cid:1755)(cid:1755)(cid:1755)(cid:1755)(cid:4654) 12 (cid:4666)(cid:1876) (cid:3398) (cid:1854)(cid:4667) (cid:3505) (cid:1856)(cid:1877)(cid:3493)(cid:4666)(cid:1877) (cid:3398) (cid:1854)(cid:4667)(cid:4666)(cid:1876) (cid:3398) (cid:1877)(cid:4667) (cid:3051)(cid:3029) (cid:3404) (cid:2024)2 (cid:3051)(cid:3029) (cid:4666)(cid:1876) (cid:3398) (cid:1854)(cid:4667), since the last integral is equal to (cid:2024) [39]. Whence, (cid:1858)(cid:4666)(cid:1876)(cid:4667) (cid:3404) (cid:3421) 0, (cid:1876) (cid:3409) (cid:1854) (cid:1827) (cid:1856)(cid:1856)(cid:1876) (cid:4666)(cid:1876) (cid:3398) (cid:1854)(cid:4667) (cid:3404) (cid:1827), (cid:1876) (cid:3408) (cid:1854). (cid:3425) (cid:3404) (cid:1827)(cid:1834) (cid:3029) (cid:4666)(cid:1876)(cid:4667).
Therefore, the solution in Eq. (25) is a Heaviside step function for the beta spectrum. If there are (cid:1863) neutrino mass eigenstates (cid:1839) (cid:3038) , (cid:1853) (cid:3038) (cid:3404) (cid:1839) (cid:3038)(cid:2870) , with relative contributions (cid:1827) (cid:3038) , then the function (cid:1859) from Eq. (26) can be written as: (cid:1859)(cid:4666)(cid:1877)(cid:4667) (cid:3404) (cid:3533) (cid:1827) (cid:3038) (cid:3493)(cid:1877) (cid:3398) (cid:1853) (cid:3038) , (cid:4666)27(cid:4667) (cid:3038) and a generalization of the solution to the Abel equation is given by Eq. (12).
6. Numerical examples
In this Section, several numerical examples are provided for the solution to the Abel equation. In a possible applications of this method to the experimental spectrum, test function (cid:1859) is defined as the ratio of the deconvoluted β-energy spectrum (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) , obtained from Eq. (3), to calculated (cid:1830)(cid:4666)(cid:1846)(cid:4667) . The function (cid:1859) is given by Eq. (21). Since we do not know a priori how many eigenvalues there are and what their corresponding neutrino masses are, one needs to perform a numerical solution on test function (cid:1859) using Eq. (25). However, for the purpose of this demonstration, we assume one or more components in a simulated spectrum (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) , calculated according to Eq. (1), without deconvolution. For numerical integration and differentiation of the test function using Eq. (25), we apply the established high-performance algorithms [35] with an additional adaptive optimization of the differential step. In the first example, we demonstrate how the step-function solution (cid:1858) depends on the bin size of abscissa (cid:1876) , for mass-less neutrino. Several step functions for zero neutrino-mass H decay (cid:4666)(cid:1843) (cid:3404)18590.29 eV(cid:4667) [30] are plotted in Fig. 2 as functions of the (cid:1876) (cid:3404) (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) argument for different bin sizes. The leftmost points at (cid:1876) = 0 could not be plotted on a logarithmic scale, so they were placed on the abscissa at the smallest bin studied, (cid:1876) = 0.01 eV . Therefore, the position of the step for a given bin size can be taken as an upper limit of the square of neutrino mass. Figure 2. Numerical solutions of the Abel equation for zero neutrino mass H decay (except curve (a)), plotted for various abscissa bin sizes in eV indicated in the legend. Zero-abscissa values are arbitrarily placed at 0.01 eV . The relative contributions are scaled to 1.0, 0.9, 0.8, and 0.7. Curve (a) depicts the effect of energy nonlinearity for 0.01 eV bin. In the second example, we selected beta decay of S ( (cid:1843) = 167.319 keV [30]) and assumed 3 neutrino mass eigenvalues at 0, 10, and 15 keV as well as their contributions (normalization factors) as given in Table 1. Function (cid:1859) was calculated using Eq. (27), and subjected to the numerical solution using Eq. (25). The solution (cid:1858) is plotted as a function of argument (cid:4666)(cid:1843) (cid:3398) (cid:1846)(cid:4667) (cid:2870) in Fig. 3 (solid curve), emphasizing the vicinity of 1 on the ordinate axis. It is seen that the three step functions are easily identified. The raise of the step functions on the abscissa of the numerical solution correctly identify squares of the neutrino masses at 0, 100, and 225 keV, respectively. The identified contributions agree with the assumed contributions. While the method is in principle non-statistical, the relative uncertainties of the identified contributions are caused by rounding errors in the numerical solution. The relative uncertainty is negligible for the principal eigenvalue and between 1-2% for the minor eigenvalues. Parameter Eigenvalue 1 2 3 Assumed neutrino mass (keV) Assumed contribution Identified contribution Relative uncertainty 0 1 1.000E+00 2.2E-11 10 0.001 9.997E-04 1.2E-02 15 0.001 1.000E-03 1.7E-02
Figure 3. Numerical solution of the Abel equation for searching of neutrino mass in the beta decay of S. Solid curve: three neutrino components with masses of 0, 10, and 15 keV. Dashed curve: two neutrino components with masses 0 and 20 keV with a shape correction function. Table 1. Assumed and identified parameters for 3-neutrino mass eigenvalue test of the numerical solution of the Abel equation for searching of neutrino mass in beta decay of S. When comparing experimental and theoretical β spectra, small smooth discrepancies were observed, which increase with decreasing kinetic energy [24]. They can be corrected for with a shape-correction factor (cid:1845) , which is assumed to be a parabola with an argument (cid:1843) (cid:3398) (cid:1846) . An example of such a shape correction function is given in Fig. 4 (solid curve) for S decay. This correction is equal to zero at a maximum kinetic energy, and it is assumed 1.01 at zero energy. In the following, we analyzed the effect of this correction on the Abel solution. We assumed that (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) had a deviation from a theoretical shape according to the shape correction function in Fig. 4, included in (cid:1830)(cid:4666)(cid:1846)(cid:4667) in Eq. (2), and did not correct for it when dividing by (cid:1830)(cid:4666)(cid:1846)(cid:4667) in Eq. (21). We also assumed two neutrino mass eigenstates at 0 (contribution of 1) and 20 keV (contribution of 0.001). Subsequent application of the Abel solution resulted in a dashed curve in Fig. 3. It is seen that not applying shape correction results in the flat portions of the step function sloping upward, however, both the neutrino mass and its contribution were unaffected. Figure 4. Shape correction and energy nonlinearity plotted as functions of energy for S decay. Besides the shape correction, which affects the intensity of the β spectrum (in the ordinate direction), we study the effect of energy-scale nonlinearities (in the abscissa direction) on the Abel solution and neutrino-mass detection for S decay. Such minute nonlinearities may result from improper energy calibration or from nonlinear response of a digital signal processor. As a next example, we assumed the energy nonlinearity with functional dependence as depicted in Fig. 4 (dashed curve). We calculated (cid:1840)(cid:4662) (cid:4666)(cid:1846)(cid:4667) with the energy nonlinearities included and did not correct for them in (cid:1830)(cid:4666)(cid:1846)(cid:4667) for Eq. (21). We study the effects close to the maximum kinetic energy (endpoint energy, at zero neutrino mass) and at lower kinetic energies (for a possible presence of heavy neutrino) separately. For any nonlinearity that is positive at the end-point energy and is not corrected for, the measured energy would exceed (cid:1843) and thus would be identified as an experimental problem. Therefore, we assumed nonlinearity equal to –0.1 keV at the maximum kinetic energy (see Fig. 4). Such nonlinearity would result in the measured energy being shifted down by 0.1 keV (~ 0.06%) below (cid:1843) . This might simulate the presence of neutrino mass. The Abel solution for this case is depicted in Fig. 5. Without nonlinearities present (dashed curve), the upper limit on the square of neutrino mass would be placed at 0.001 keV . The above-mentioned nonlinearity appear to simulate a square of the neutrino mass at 0.01 keV (solid curve), however, it can be rejected because the step is not sharp but is a rather slowly raising function instead. The effect of the energy nonlinearity from Fig. 4 on the wider range of the Abel solution is depicted in Fig. 6 (solid curve), where we have also added a 100-keV neutrino component at 0.001 contribution. It is seen that the curve resembles an exact step, as it does in the case without nonlinearities depicted by the dashed curve. Therefore the heavy neutrino is identified at 100 keV and only its mass is shifted by 0.1 keV nonlinearity (see Fig. 4). Figure 5. Step function solutions for zero neutrino mass S decay, with (solid curve) and without (dashed curve) energy nonlinearities. We have also assessed the effect of energy nonlinearities on neutrino mass in H decay. We assumed a functional nonlinearity similar to the one in Fig. 4, except the energy shift at the endpoint was assumed to be –1 eV (~ 0.005%) below (cid:1843) . The Abel solution resulted in curve (a) in Fig. 2. It appears to simulate a square of the neutrino mass at 1 eV , however, it can be rejected following the arguments for S in Fig. 5. Figure 6. Step function solutions for S decay consisting of two neutrino components having mases of 0 and 100 keV, with (solid curve) and without (dashed curve) energy nonlinearities.
7. Summary and conclusions
We have developed a math-theoretical model for searches of neutrino mass in beta-decay energy spectra. The transformed Fredholm equation was derived (Eq. (14a) with the kernel given by Eq. (14b)), linking the observed beta spectrum with the theoretical one as well as with the instrumental response. The solution to this equation is a superposition of the Heaviside step-functions representing the possible neutrino mass eigenvalues, as well as their contributions (Eq. (12)). In this way, the shape of beta spectrum is tested for discrepancies, which could be associated with the neutrino mass. The step function is a very characteristic signature of such discrepancies. A possible solution to the Fredholm equation was derived based on series expansion (Eq. (17)) and solution of the matrix equation (Eq. (19)), taking advantage of the fact that the observed spectrum is discrete. Numerical verification of the transformed Fredholm equation deserves further investigation.
A limiting version of the transformed Fredholm equation was also derived. It is based upon regular deconvoluting of the observed spectrum from the instrumental response using Eq. (3). The factors dependent on neutrino mass (Eq. (27)) are then subject to Abel transformation (Abel integral equation, Eq. (22)), the solution of which (Eq. (25)) is proven in to be a superposition of the step functions (Eq. (12)). The Abel transformation was numerically tested on theoretically calculated beta spectra. It was shown that this method is sensitive down to a relative contribution of about 10 ― for heavy neutrinos and can set upper limits on the neutrino mass of 0.1 eV for H and 0.03 keV for S. We have also numerically tested two possible experimental deviations in beta spectra: shape correction and energy nonlinearities. It was shown that a smooth, moderate shape correction up to 1% at low energy changes the slope of the step function but does not affect the numerical recognition of either zero-mass or heavy neutrino. Likewise, heavy neutrino detection is not affected by energy nonlinearity. Nevertheless, the energy nonlinearity can simulate an artificial neutrino mass close to the endpoint of β spectrum. However, if the Abel method is applied to such a case, the step function is not sharp and, therefore, the hypothesis of neutrino mass can be rejected. If not neutrino mass, the proposed method can indicate some discrepancies between measured and calculated β spectra. The math-theoretical methods derived here are inherently non-statistical, however they have possible limitations arising from numerical errors, theoretical considerations, and statistical fluctuations present in the observed spectra. They are discussed below in that order. It was shown that at the 10 ― level, the Abel-transformation method resulted in relative fluctuations below 2%, originating from rounding errors in the calculations. Such rounding errors could be possibly lowered, resulting in a higher sensitivity, if higher-precision arithmetic were used in the numerical calculations. The theory of beta decay with all possible corrections, from which theoretical beta spectra are derived, can be given only to a finite accuracy, which is a limiting factor for neutrino mass searches in beta decay, regardless of the method used [24]. Finally, the statistical fluctuations in the observed beta spectra limit the sensitivity of neutrino-mass searches in beta decay by any method, including this one. Owing to the fact that the observed spectrum is discrete, there are Poisson fluctuations in each spectrum channel. The variation coefficient of the Poisson fluctuations is given by , where (cid:1875) is the number of observed counts in the channel [33]. As seen in Fig. 1, the fluctuations are expected to be lower at the low beta energy region and higher close to the endpoint. The fluctuations present in the observed spectrum will translate into the fluctuations of any solution derived from the spectrum. It is, therefore, imperative to have the observed statistics in the data as high as possible. One way to alleviate problems with low statistics could be fitting of a polynomial to such data and applying integral equations to the polynomial fit instead of the original data points. Acknowledgments
Thanks are due to C. J. Bradt for his valuable comments.
References [1]
E. Fermi, Tentativo di una teoria dell’emissione dei raggi « beta » , Ric. Scientifica (2) (1933) 491-495. [2] E. Fermi, Versuch einer theorie der β-strahlen, Z. Physik (1934) 161-171. [3] F.L. Wilson, Fermi’s theory of beta decay, Am. J. Phys. (1968) 1150-1160. [4] C. Giunti, C.W. Kim, Fundamentals of Neutrino Physics and Astrophysics (Oxford University Press, Oxford, 2007). [5]
O.G. Miranda, J.W.F. Valle, Neutrino oscillations and the seesaw origin of neutrino mass, Nucl. Phys. B (2016) 436-455. [6]
G. Drexlin, V. Hannen, S. Mertens, C. Weinheimer, Current direct neutrino mass experiments, Adv. High En. Phys. 293986 (2013) 1-39. [7]
Planck Collaboration, Planck 2015 results XIII. Cosmological parameters, Astron. Astrophys.
A13 (2016) 1-63. [8]
GERDA Collaboration, Background-free search for neutrinoless double-β decay of Ge with GERDA, Nature (2017) 47-52. [9]
H.V. Klapdor-Kleingrothaus, I.V. Krivosheina, A. Dietz, O. Chkvorets, Search for neutrinoless double beta decay with enriched Ge in Gran Sasso 1990–2003, Phys. Lett. B (2004) 198-212. [10]
T.J. Loredo, D.Q. Lamb, Bayesian analysis of neutrinos observed from supernova SN 1987A, Phys. Rev. D
C. Kraus, B. Bornschein, L. Bornschein et al ., Final results from phase II of the Mainz neutrino mass search in tritium (cid:2010) decay, Eur. Phys. J. C (2005) 447-468. [12] V.N. Aseev, A.I. Belesev, A.I. Berlev et al ., Upper limit on the electron antineutrino mass from the Troitsk experiment, Phys. Rev. D S. Mertens, Direct neutrino mass experiments, J. Phys.: Conf. Ser.
F. Gatti, Microcalorimeter measurements, Nucl. Phys. B (Proc. Suppl.) (2001) 293-296. [15] M. Sisti, C. Arnaboldi, C. Brofferio et al ., New limits from the Milano neutrino mass experiment with thermal microcalorimeters, Nucl. Instr. Meth. Phys. Res. A (2004) 125-131. [16]
P.T. Springer, C.L. Bennett, P.A. Baisden, Measurement of the neutrino mass using the inner bremsstrahlung emitted in the electron-capture decay of
Ho, Phys. Rev. A (1987) 679-689. [17] S. Yasumi, H. Maezawa, K. Shima et al ., The mass of the electron neutrino from electron capture in
Ho, Phys. Lett. B (1994) 229-233. [18]
A. Hime, J.J. Simpson, Evidence of the 17-keV neutrino in the β spectrum of H, Phys. Rev. D (1989) 1837-1850. [19] E.B. Norman, B. Sur, K.T. Lesko, M.M. Hindi, R.-M. Larimer, T.R. Ho, J.T. Witort, P.N. Luke, W.L. Hansen, E.E. Haller, Evidence for the emission of a massive neutrino in the nuclear beta decay, J. Phys. G: Nucl. Part. Phys. (1991) S291-S299. [20] J.J. Simpson, A. Hime, Evidence of the 17-keV neutrino in the β spectrum of S, Phys. Rev. D (1989) 1825-1836. [21] A. Hime, N.A. Jelley, New evidence for the 17 keV neutrino, Phys. Lett. B (1991) 441-449. [22]
A. Hime, Do scattering effects resolve the 17-keV conundrum? Phys. Lett. B (1993) 165-173. [23] F.E. Wietfeldt, E.B. Norman, The 17-keV neutrino, Phys. Rep. (1996) 149-197. [24]
M. Bahran, W.R. Chen, G.R. Kalbfleisch, Fermi theory of nuclear decay and heavy neutrino searches, Phys. Rev. D (1993) R759-R763. [25] X. Mougeot, Systematic comparison of beta spectra calculations using improved analytical screening correction with experimental shape factors, Appl. Radiat. Isot. (2016) 177-182. [26]
H. Behrens, J. Jӓnecke, Numerical Tables for Beta-Decay and Electron Capture, in: K.-H. Hellwege (Ed.), Landolt-Bӧrnstein Numerical Data and Functional Relationships in Science and Technology (Springer-Verlag, Berlin, 1969). [27]
M.E. Rose, A note on the possible effect of screening in the theory of beta-disintegration, Phys. Rev. (1936) 727-729. [28] L. Durand, Electron screening corrections to beta-decay spectra, Phys. Rev. (1964) B310-B313. [29]
J.M. Carmona, J.L. Cortés, Testing Lorentz invariance violations in the tritium beta-decay anomaly, Phys. Lett. B (2000) 75-80. [30]
W.J. Huang, G. Audi, M. Wang, F.G. Kondev, S. Naimi, Xing Xu, The Ame2016 atomic mass evaluation, (I). Evaluation of input data; and adjustment procedures, Chin. Phys. C (3) 030002 (2017) 1-344. [31] D. Porter, D.S.G. Stirling, Integral Equations (Cambridge University Press, Cambridge, 1993) Sects. 1.2, 1.3.4, 9.2. [32]
N.L. Johnson, S. Kotz, N. Balakrishnan, Continuous Univariate Distributions, Vol. 1 (Wiley-Interscience Publication, New York, 1994) Sect. 12.3. [33]
N.L. Johnson, A.W. Kemp, S. Kotz, Univariate Discrete Distributions (Wiley-Interscience, Hoboken, 2005) Sects. 1.2.11, 4.4. [34]
P.A.M. Dirac, The Principles of Quantum Mechanics (Oxford University Press, Oxford, 1988) Sect. III.15. [35]
Numerical Algorithms Group Ltd, Oxford, UK. [36]
Wolfram Research, Champaign, IL, USA. [37]
H. Margenau, G.M. Murphy, The Mathematics of Physics and Chemistry (Krieger, Huntington, 1976) Sect. 8.5. [38]