We start by summarizing 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    studies on 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. According to the theory of neutrinos, mass eigenstates mix to result in the flavor eigenstates, the latter bear masses close to zero   . In addition, a massive sterile neutrino has been sought   . In this work, we will abbreviate neutrino “eigenstate” any measurable neutrino, which can have zero or positive mass, without referring to a specific generation mechanism.
One of the approaches to the neutrino mass searches has been studying an endpoint region of a beta-particle spectrum  . The Mainz and Troitsk experiments for 3H beta decay yielded the upper limits of electron antineutrino mass at 2.3 eV  and 2.05 eV  , respectively, while the KATRIN experiment is expected to have a sensitivity down to 0.2 eV   .
In beta ( ) decay, β-particle (electron) and electron antineutrino, (hereinafter referred to as neutrino) are emitted in the weak-interaction process of neutron decay inside the nucleus . The emitted β-particle kinetic-energy spectrum is a continuous function, owing to the energy sharing between the beta particle and the neutrino. Possible neutrino mass M in the beta spectrum is included in the factor , where Q is the nuclear-recoil corrected Q-value of beta decay and T is the emitted electron kinetic energy. The maximum beta energy (called the endpoint energy) is equal to Q or , when the neutrino mass is zero or positive, respectively. In principle, there can be several mixed neutrino eigenstates. 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: 1) convolution of the theoretical beta spectrum with the normalized instrumental response function and comparison of the convolution with the observed spectrum, and 2) deconvolution of the observed spectrum and comparison with the theoretical spectrum  -  . Statistical measures, such as 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 neutrino mass, which affect them.
In order to advance beyond the difficulties and uncertainties associated with the statistical measures used for neutrino mass detection in beta decay of the existing methods, we aimed in this work at developing a novel mathematical theory based on transformational properties, rather than statistical measures. 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 , one for each possible neutrino mass eigenstate. The neutrino mass can then be identified by the abscissa value at the raise of the step function, whereas its eigenstate contribution is shown to be proportional to the ordinate increment past the step. In Section 4, we outline a possible solution to the transformed Fredholm equation using series expansion and linear algebra.
If method 2) of deconvoluting the beta spectrum first is employed, before comparison with the theory, then the transformed Fredholm equation from Section 3 is shown to reduce to the Abel integral equation in Section 5, also known as the Abel transformation. The Abel integral equation, which is a special case of the Volterra integral equation of the first kind, has a known integral solution. We prove as a Lemma 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 numerical examples of solutions of the Abel equation, when the beta spectrum contains several neutrino-mass eigenstates. We also describe how the Abel solution behaves in the presence of small experimental discrepancies, which can influence neutrino mass detection, such as beta shape correction or small nonlinearities in the measured beta-energy scale. We summarize the advantages and limitations of the proposed new mathematical method of elucidating neutrino mass in beta decay in Section 7, as well as 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    . For the purpose of this work, we define the beta-energy spectrum as follows:
where represents a distribution of beta particles with energies between and T and , factor A contains the quantum-transition matrix element and an overall normalization, whereas Q, T and M have been defined in Section 1.
Factor in Equation (1) is defined as follows:
where p is relativistic momentum, E is total (relativistic) energy, is a non-relativistic Fermi function, is a relativistic correction to the Fermi function  and is a radiative correction  . is a shape correction  , which corrects for small smooth discrepancies between the theoretical and observed spectra at low energies, and is assumed to be a parabola with an argument . Also, a screening correction   is applied to . Additional term involving the progeny’s final excitation states in low beta-energy emitters  as well as a possible Lorentz-invariance violation correction  have not been included.
The examples of beta spectra from 3H decay , calculated according to Equations ((1), (2)) are plotted in Figure 1. The physical constants and nuclear parameters for the calculations were taken from  . It is seen that
Figure 1. Beta particle spectra from 3H decay. Black spectrum: zero neutrino mass; red spectrum: equal superposition of two neutrino eigenstates with the masses of 3000 and 5000 eV.
if neutrino mass is non-zero, the endpoint beta energy lies below Q, while another neutrino mass eigenstate is evidenced by the non-differentiable point in the spectrum. If neutrino mass is close to zero and any other eigenstate has a very small contribution, these effects are difficult to discern from the beta spectrum.
3. Transformed Fredholm Integral Equation
The observed beta spectrum usually differs from the emitted (theoretical) beta spectrum, given by Equations ((1), (2)), owing to an instrumental response. Let T represent the emitted kinetic energy, while the observed kinetic energy. The emitted and observed beta spectra are then and , respectively. They are coupled with each other by the response function through the integral:
Equation (3) is a Fredholm integral equation of the first kind  , where is a kernel. The variable can extend past due to the abovementioned instrumental response. The response function is measured separately. Since the emission and detection of beta particles are stochastic processes, the normalized and can be considered statistical probability density functions (pdf). The integral in Equation (3) is then a mixture of distributions  . If, additionally, is a function of e.g., a Gaussian, then Equation (3) is a mathematical convolution  . The process of solving Equation (3) for is referred to as deconvolution and is discussed in Section 4. The terms: convolution and deconvolution will be used here even if Equation (3) is not a strictly mathematical convolution.
If there existed k neutrino-mass eigenstates with the masses and normalization factors , then combining Equations ((1), (3)) and summing over the eigenstates would result in the observed spectrum:
Subsequently, we introduce a mass variable m and utilize the Dirac δ-function  . Then Equation (4) transforms to:
New variables are defined: and . By substitution, Equation (5) results in
By using a property of the δ-function  ,
and substituting of the variables, , , and , one obtains from Equation (6):
Changing the order of integration in Equation (8), yields
The inner integral in Equation (9) can be integrated by parts. We also use the fact that the integral of the δ-function on positive argument is the Heaviside step function:
The result is given below:
The solution to Equation (11) is the function
which is a superposition of the step functions for all possible neutrino mass eigenstates, scaled by their appropriate contributions. By inserting Equation (12) into Equation (11) one obtains
By reversing the order of integration in Equation (13) back to that of Equation (8) and by dividing both sides by i.e., the theoretical factor D given by Eq. (2) evaluated at the observed energy , one obtains
Equation (14a) is the transformed Fredholm equation. The kernel is given by Equation (14b). It links the observed beta spectrum with the emitted (theoretical) spectrum and with the instrumental response. In the solution function given by Equation (12), one does not need to assume the presence of neutrino mass or its eigenstates. 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 established algorithms for solving Fredholm equations given by Equations (3) and (14a)   . If Equation (3) was a mathematical convolution, then it could 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 convoluted functions  .
However, Equation (14a) is not a convolution and we derive below a matrix method to solve it. We take into consideration the fact that the observed variable is discrete and incremented at equally spaced intervals (channels) owing to the instrumental digitizing of the spectrum (modern instrumentation can have or more channels). Consequently, variable x is also digitized. Consider a general integral equation of the first kind, of which Equations ((3), (14a)) are the examples:
If x is discrete, one has:
Let us assume a trial solution in a form of a series:
By inserting the trial solution into Equation (16), one obtains a set of linear equations:
In matrix notation, Equation (18a) reads:
where and are vectors and is a symmetric square matrix. By solving the linear equations for , the solution (17), continuous in y variable (i.e., in ) 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 was first deconvoluted using Equation (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 Equation (14a),
where we used a version of Equation (7). Then, starting from Equation (13),
By defining a new function:
we obtain from Equation (21)
Equation (23) is a special case of the Abel transformation, also called the Abel integral equation   , which is a special case of the Volterra integral equation of the first kind. The general form of the Abel equation is
with a solution of
Therefore, the solution of Equation (23) for is
By using Equation (1) for a single neutrino-mass eigenstate M, and setting as well as , we obtain from Equation (22)
Lemma. If is given by Equation (27), then , where H is the Heaviside step function given by Equation (10).
Proof. By inserting Equation (27) to Equation (26), we obtain
since the last integral is equal to π  . Whence,
which proves the Lemma.
If there are k neutrino mass eigenstates , , with relative contributions , then the function g from Equation (27) can be written as:
and, according to the Lemma, generalization of the solution to the Abel equation is given by Equation (12).
6. Numerical Examples
In this Section, several numerical calculations are provided for the solution to the Abel equation. In a possible application of this method to the observed spectrum, function g in Equation (22) is defined as the ratio of the deconvoluted beta-energy spectrum , obtained from Equation (3), to calculated using Equation (2). Since we do not know a priori how many neutrino eigenstates there are and what are their corresponding masses, one needs to perform a numerical solution on function g using Equation (26). However, for the purpose of this calculation, we assumed one or more eigenstate components in a simulated spectrum , calculated according to Equation (1), without deconvolution. For numerical integration and differentiation of the test function using Equation (26), we applied the established high-performance algorithms  with an additional adaptive optimization of the differential step.
In the first calculation, we demonstrate how the step-function solution f depends on the bin size of abscissa x, for mass-less neutrino. Several numerical solutions for zero neutrino-mass 3H decay are plotted in Figure 2 as functions of the argument for different bin sizes. The leftmost points at could not be plotted on a logarithmic scale, so they were placed on the abscissa at the smallest bin studied, eV2. Therefore, the raise of the step for a given bin size can be taken as an upper limit of the square of neutrino mass.
In the second calculation, we selected beta decay of 35S (Q = 167.319 keV  )
Figure 2. Numerical solutions of the Abel equation for zero neutrino mass 3H decay (except curve (a)), plotted for various abscissa bin sizes in eV2 indicated in the legend. Zero-abscissa values are arbitrarily placed at 0.01 eV2. The relative contributions are scaled to 1.0, 0.9, 0.8, and 0.7 to enable better visualization. Curve (a) depicts the effect of energy nonlinearity for 0.01 eV2 bin.
Table 1. Assumed and identified parameters for 3-neutrino mass eigenstate test of the numerical solution of the Abel equation for searching of neutrino mass in beta decay of 35S.
and assumed 3 neutrino eigenstates with the masses of 0, 10, and 15 keV as well as their contributions (normalization factors) of 1, 0.001, and 0.001, respectively, as given in Table 1.
Function g was calculated using Equation (31), and subjected to the numerical solution using Equation (26). The solution f is plotted as a function of argument in Figure 3 (black 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 given in Table 1 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 eigenstate and between 1% - 2% for the minor eigenstates.
When comparing experimental and theoretical beta spectra, small smooth discrepancies were observed, which increase with decreasing kinetic energy  . They can be corrected for with a shape-correction factor S, which is assumed to be a parabola with an argument . An example of such a shape correction function is given in Figure 4 (black curve) for 35S decay. This correction factor is equal to one 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 had a deviation from a theoretical shape according to the shape correction function in Figure 4, included in in Equation (2), and did not correct for it when dividing by in Equation (22). We also
Figure 3. Numerical solution of the Abel equation for searching of neutrino mass in the beta decay of 35S. Black curve: three neutrino components with masses of 0, 10, and 15 keV. Red curve: two neutrino components with masses 0 and 20 keV with a shape correction function.
Figure 4. Shape correction and energy nonlinearity plotted as functions of energy for 35S decay.
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 red curve in Figure 3. It is seen that not correcting for shape results in the flat portions of the step function sloping upward, however, both the neutrino masses and their contributions were unaffected.
Besides the shape correction, which affects the intensity of the beta 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 35S 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 a functional dependence as depicted in Figure 4 (red curve). We calculated with the energy nonlinearities included and did not correct for them in for Equation (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 deconvoluted energy would exceed Q and thus would be identified as an experimental problem. Therefore, we assumed nonlinearity equal to −0.1 keV at the maximum kinetic energy (see Figure 4). Such nonlinearity would result in the observed energy being shifted down by 0.1 keV (~0.06%) below Q. This might simulate the presence of neutrino mass. The Abel solution for this case is depicted in Figure 5. Without nonlinearities present (black curve), the upper limit on the square of neutrino mass would be placed at 0.001 keV2. The above-mentioned nonlinearity appear to simulate a square of the neutrino mass at 0.01 keV2 (red curve), however, it can be rejected because the step is not sharp but is a rather slowly raising function instead.
Figure 5. Step function solutions for zero neutrino mass 35S decay, with (red curve) and without (black curve) energy nonlinearities.
Figure 6. Step function solutions for 35S decay consisting of two neutrino components having mases of 0 and 100 keV, with (red curve) and without (black curve) energy nonlinearities.
The effect of the energy nonlinearity from Figure 4 on the wider range of the Abel solution is depicted in Figure 6 (red 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 at 104 eV2, as it does in the case without nonlinearities depicted by the black curve. Therefore the heavy neutrino is identified at 100 keV and only its mass is shifted by the 0.1 keV nonlinearity (see Figure 4).
We have also assessed the effect of energy nonlinearities on neutrino mass in 3H decay. We assumed a functional nonlinearity similar to the one in Figure 4, except the energy shift at the endpoint was assumed to be ?1 eV (~0.005%) below Q. The Abel solution resulted in curve (a) in Figure 2. It appears to simulate a square of the neutrino mass at 1 eV2, however, it can be rejected following the arguments for 35S in Figure 5.
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 (Equation (14a) with the kernel given by Equation (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 eigenstates, as well as their contributions (Equation (12)). In this way, the shape of beta spectrum is tested for non-differentiable points, which could be associated with the neutrino mass. The step function is a very characteristic signature of such non-differentiabilities. A possible solution to the Fredholm equation was derived based on a series expansion (Equation (17)) and solution of the matrix equation (Equation (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 prior normal deconvolution of the observed spectrum from the instrumental response using Equation (3). The factors dependent on neutrino mass (Equation (31)) are then subjected to the Abel transformation (Abel integral equation, Equation (23)), the solution to which (Equation (26)) is proven in the Lemma to be a superposition of the step functions. The Abel transformation was numerically tested on theoretically calculated beta spectra. It was shown that this method is sensitive down to a fractional contribution of about 10−3 for heavy neutrinos and can set upper limits on the neutrino mass of 0.1 eV for 3H and 0.03 keV for 35S.
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. However, the energy nonlinearity can simulate an artificial neutrino mass close to the endpoint of beta spectrum. Nevertheless, 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 observed and calculated beta spectra.
The novelty of the model presented in this work is that the solutions with the step functions are very sensitive and characteristic signatures of neutrino mass, and they can in many cases distinguish between experimental discrepancies and true neutrino-mass signatures. Therefore, this model may have some advantages over the statistical methods, which may be less specific and less sensitive for neutrino mass detection in beta decay.
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−3 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  .
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 N is the number of observed counts in the channel  . As seen in Figure 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 original data points.
Thanks are due to C. J. Bradt for his valuable comments.
 Drexlin, G., Hannen, V., Mertens, S. and Weinheimer, C. (2013) Current Direct Neutrino Mass Experiments. Advances in High Energy Physics, 2013, Article ID: 293986.
 Kraus, C., Bornschein, B., Bornschein, L., et al. (2005) Final Results from Phase II of the Mainz Neutrino Mass Search in Tritium β Decay. The European Physical Journal C, 40, 447-468.
 Aseev, V.N., Belesev, A.I., Berlev, A.I., et al. (2011) Upper Limit on the Electron Antineutrino Mass from the Troitsk Experiment. Physical Review D, 84, 112003.
 Mougeot, X. (2016) Systematic Comparison of Beta Spectra Calculations Using Improved Analytical Screening Correction with Experimental Shape Factors. Applied Radiation and Isotopes, 109, 177-182.
 Behrens, H. and Janecke, J. (1969) Numerical Tables for Beta-Decay and Electron Capture. In: Hellwege, K.-H., Ed., Landolt-Bornstein Numerical Data and Functional Relationships in Science and Technology, Springer-Verlag, Berlin.
 Huang, W.J., Audi, G., Wang, M., Kondev, F.G., Naimi, S. and Xu, X. (2017) The Ame2016 Atomic Mass Evaluation, (I). Evaluation of Input Data; and Adjustment Procedures. Chinese Physics C, 41, Article ID: 030002.