Angular Distributions of Emitted Electrons in the Two-Neutrino ββ Decay

: The two-neutrino double-beta decay (2 νββ -decay) process is attracting more and more attention of the physics community due to its potential to explain nuclear structure aspects of involved atomic nuclei and to constrain new (beyond the Standard model) physics scenarios. Topics of interest are energy and angular distributions of the emitted electrons, which might allow the deduction of valuable information about fundamental properties and interactions of neutrinos once a new generation of the double-beta decay experiments will be realized. These tasks require an improved theoretical description of the 2 νββ -decay differential decay rates, which is presented. The dependence of the denominators in nuclear matrix elements on lepton energies is taken into account via the Taylor expansion. Both the Fermi and Gamow-Teller matrix elements are considered. For nuclei of experimental interest, relevant phase-space factors are calculated by using exact Dirac wave functions with ﬁnite nuclear size and electron screening. The uncertainty of the angular correlation factor on nuclear structure parameters is discussed. It is emphasized that the effective axial-vector coupling constant g eff A can be determined more reliably by accurately measuring the angular correlation factor.


Introduction
There is an increased interest in the two-neutrino double-beta decay (2νββ-decay) mode with an emission of two electrons and two antineutrinos which was suggested by Maria Goeppert-Mayer in 1935 [1]. The 2νββ-decay is the second order process of weak interaction fully consistent with the Standard model (SM) of particle physics, where two neutrons are simultaneously transformed into two protons inside an atomic nucleus, and two pairs of electrons and antineutrinos are emitted. It is the rarest process measured so far in nature with a half-life above 10 18 years. The first direct observation of this process dates back to 1987 [2]. Currently, the 2νββ-decay has been observed in eleven even-even nuclei for transition to ground state and in two cases also for transition to the first excited 0 + 1 state of the final nucleus [3]. Several next-generation double-beta experiments with a variety of suitable isotopes are in different stages of construction or preparation and planning. Their main goal is the observation of the neutrinoless double-beta decay (0νββ-decay), a double-beta decay mode violating the total lepton number, in which only two electrons are emitted. It would prove that neutrinos are Majorana particles (i.e., neutrinos are their own antiparticles), as most theories beyond SM requires. Besides the Majorana nature of the neutrino, valuable information on the origin of neutrino masses, violation of CP in the lepton sector and possible existence of sterile neutrinos is expected, as well.
An important by-product of these experimental searches will be a precise measurement of the 2νββ-decay with high statistics. The study of this process is important on its own and represents a research program with many important questions waiting for the answers. The 2νββ-decay provides insights into nuclear structure [4] of nuclear systems under consideration which can then be used in a reliable calculation of the nuclear matrix elements (NMEs) for 0νββ-decay, necessary to extract the particle physics parameters responsible for the lepton number violation [5] or in fixing the effective value of the axial-vector coupling constant g eff A [6][7][8][9]. With the increasing statistics of the 2νββ-decay, the possibility of exploring new physics beyond SM with this process becomes more and more reliable. The 2νββ-decay can be, and it is used for the search of the right-handed neutrino interactions without having to assume their nature [10], the neutrino self-interactions [11], the sterile neutrinos with masses up to Q-value of the process [12], the bosonic neutrino component [13], the violation of the Lorentz invariance [14][15][16][17], the 0νββ-decays with Majoron(s) emission [18,19] and the quadruple-β decay [20].
A necessary condition to gain rich information about fundamental properties and interactions of neutrinos from the precisely measured differential characteristics of the 2νββ-decay is that all atomic and nuclear structure physics aspects of this process are well understood. Recently, the theoretical description of the shape of single and summed electron energy distributions of the 2νββ-decay has been significantly improved by taking into account the dependence on lepton energies from the energy denominators of nuclear matrix elements [6], which was commonly neglected. However, their importance was already manifested by the study of 2νββ-decay energy and angular distributions within the framework of the Single State Dominance (SSD) versus Higher State Dominance (HSD) hypothesis in [21,22]. In this paper, the theoretical study of Ref. [6] is extended to analyze the angular correlations of outgoing electrons by considering both the Fermi and Gamow-Teller components of the 2νββ-decay nuclear matrix elements. For nuclei of experimental interest, all phase-space factors are calculated by using exact Dirac wave functions with finite nuclear size and electron screening.

The 2νββ Angular Distribution
The inverse half-life of the 2νββ-decay transition to the 0 + ground state of the final nucleus is defined as where Γ 2ν is the decay rate. If the contribution of higher partial waves of electron and higher order terms of the nucleon current are neglected as they are strongly suppressed [23], the angular distribution of emitted electrons takes the form [10] dΓ 2ν d(cos θ) where is the angular correlation factor and θ is the angle between the emitted electrons. The decay rates are given by with Here, G β = G F cos θ C (G F is the Fermi constant and θ C is the Cabbibo angle [24]) and m e is the mass of electron. E i , E f and E e i (E e i = p 2 e i + m 2 e , i = 1, 2) are the energies of initial and final nuclei and electrons, respectively.
are radial components of electron wave functions evaluated at nuclear surface with radius R, which will be introduced in the following sections. A 2ν and B 2ν consist of products of the Fermi and Gamow-Teller nuclear matrix elements, which depend on lepton energies: where with Here, |0 + i , |0 + f are the 0 + ground states of the initial and final even-even nuclei, respectively, and |1 + n (|0 + n ) are all possible states of the intermediate nucleus with angular momentum and parity J π = 1 + (J π = 0 + ) and energy E n (1 + ) (E n (0 + )). The lepton energies enter in the factors The maximal value of |ε K | and |ε L | is the half of Q value of the process (ε K,L ∈ (−Q/2, Q/2)). For 2νββ decay with energetically forbidden transition to intermediate nucleus ( The Fermi (Gamow-Teller) operator governing the matrix elements M F (n) (M GT (n)) given in Equation (10) is a generator of an isospin SU(2) (spin-isospin SU(4)) multiplet symmetry. In the case both isospin and spin-isospin symmetries would be exact in nuclei the 2νββ-decay would be forbidden as ground states of initial and final nuclei would belong to different multiplets. Usually, it is assumed that the contribution from Fermi matrix elements to the 2νββ-decay amplitude is negligibly small as the isospin is known to be a reasonable approximation in nuclei. The main contribution is due to Gamow-Teller matrix elements. This fact is partially confirmed by the nuclear structure studies, which are, however, model dependent. There is a question whether it is possible to prove experimentally the dominance of the Gamow-Teller over Fermi contribution in the 2νββ-decay.

The Standard Approximation and the HSD Hypothesis
Commonly, the calculation of the 2νββ-decay distributions and decay rate is simplified by neglecting ε K,L in energy denominators of NME's by which a separation of phase-space factors and NMEs is achieved. We get where Fermi and Gamow-Teller matrix elements are given by As a consequence of Equation (12) the angular correlation factor K 2ν (see Equations (4) and (5)) is not affected by the nuclear structure as the dependence on A 2ν and B 2ν is eliminated. We note that this approximation scheme is justified when the contribution from higher lying states (0 + n , 1 + n ) of the intermediate nucleus to the 2νββ-decay rate dominates over the low lying states. Such theoretical assumption is denoted as the higher state dominance (HSD) hypothesis. The experimental observation shows that Fermi strength distribution from the initial to the intermediate nucleus is concentrated to the Isobaric Analogue State region with excitation energy above 10 MeV, that is, the isospin symmetry prevents any significant fragmentation of the Fermi transition. Contrary, the Gamow-Teller strength distribution is fragmented over many daughter states covering a region of Gamow-Teller resonance and a region of low-lying states of the intermediate nucleus. Thus, the HSD assumption is justified for Fermi matrix elements M K,L F and it is an open question whether transitions through low-lying or higher-lying (region of the GT resonance) 1 + states of the intermediate nucleus give the main contribution to M K,L GT , or whether there is a mutual cancellation of these contributions.

The Angular Distribution within the SSD Hypothesis
A long time ago, it was proposed by Abad et al. [25] that the 2νββ-decay is governed by the Gamow-Teller transitions connecting the 0 + ground states of initial (A,Z) and final (A,Z+2) even-even nuclei with the first 1 + 1 state of the intermediate nucleus (A,Z+1). This assumption is known as the SSD hypothesis. Then, we have and M K,L F = 0. In this case the normalized to the full width single and summed electron differential decay rates and the angular correlation factor K 2ν are free of M GT (n = 1), but they are affected by the lepton energies incorporated in ε K,L and the energy difference [21,22]. The experimental study of energy distributions performed for 2νββ-decay of 82 Se [26] and 100 Mo [19] showed a clear preference for the SSD scenario over the HSD scenario. However, a better interpretation of experimental data requires an improved theoretical description the 2νββ-decay process.

The Dependence of Angular Distribution on Lepton Energies via Taylor Expansion
The improved formulas for the 2νββ differential decay rates can be obtained by taking advantage of the Taylor expansion over the parameters ε K,L containing the lepton energies in the NMEs denominators. This procedure allows the factorization of NMEs and phasespace factors. In [6] it was manifested that additional terms due to Taylor expansion in the decay rate are significant. It indicates that the effect of nuclear structure on the angular distribution of emitted electrons in the 2νββ-decay might be important as well.
Our assumptions are as follows: (i) The Fermi matrix element M K,L F is governed by the transitions through isobaric analogue state with energy above 10 MeV and the effect of lepton energies in the energy denominators can be neglected, that is, M K,L F M F ; (ii) The dependence of Gamow-Teller matrix element M K,L GT on lepton energies cannot be neglected and will be taken into account by Taylor expansion of energy denominators over the ratio By limiting our consideration to the fourth power in ε K,L , we get Here, the leading Γ 2ν 0 (Λ 2ν 0 ), next to leading Γ 2ν 2 (Λ 2ν 2 ) and next-to-next to leading Γ 2ν 22 and Γ 2ν 4 (Λ 2ν 22 and Λ 2ν 4 ) orders in Taylor expansion are given by and The phase-space factors are defined as with N = {0, 2, 22, 4}. In the phase-space expressions, The products of nuclear matrix elements are given by and where nuclear matrix elements take the forms By introducing ratios of NMEs, for the angular distribution we get where and The factor ξ 2ν 31 is a free parameter which can be determined from the measured single electron, summed electron energy distributions or the angular correlation factor K 2ν . The corresponding analysis can be performed by exploiting the SSD relation for the factor ξ 2ν 51 : Here, E 1 (1 + ) is the energy of the lowest 1 + state of the intermediate nucleus. This approximation reflects the higher power of energy denominators in matrix elements M 2ν and M 2ν GT−5 , which suppresses strongly the contributions from higher lying states of the intermediate nucleus.
The subject of interest might be also dependence of the angular distribution on energies of the emitted electrons,

Electron Wave Function
An important ingredient for the evaluation of the phase space factors, both in single β and double β decay, is the distortion of the outgoing particle wave function in the electrostatic potential of the daughter nucleus, V(r). If the potential is spherically symmetric V(r) ≡ V(r), the time-independent Dirac equation is satisfied by [27] Ψ(E e , r) = ∑ κµ g κ (E e , r)χ where g κ (E e , r) and f κ (E e , r) are radial functions, indexed by the relativistic quantum number κ, eigenvalue of the operator K = β(σ · L + 1) which takes negative or positive integer values (|κ| = k, k = 1, 2, 3, . . . ). Here, 1 is the 4 × 4 unit matrix, L is the orbital angular momentum operator, σ stands for the Pauli matrices in four dimensions, and β is equal to the γ 0 Dirac matrix. E e and r are the energy and the position vector of the electron, respectively. χ µ κ are the usual spin-angular momentum wave functions with and χ m the traditional spinor. The total angular momentum of the electron is j κ = |κ| − 1/2 and the orbital angular momentum can be defined by The radial wave functions g κ (E e , r) and f κ (E e , r) satisfy the Dirac radial equation and must be normalized so that they have for large values of pr the asymptotic behavior where is the Sommerfeld parameter and δ κ is the phase shift. Here Z f is the charge of the final nucleus which generates the potential V(r).
The wave function of the outgoing electron, can be expanded in spherical waves as where the spherical waves index is displayed in the spectroscopic notation. For the case of 2νββ decay, 0 + → 0 + transitions, we are interested to study wherep = p/p, p is the electron momentum vector and σ stands for the Pauli matrices in two dimensions.

The Approximation Scheme A
We adopt the relativistic electron wave function in a uniform charge distribution in the final nucleus, generating the potential where R is the radius of the daughter nucleus, R = r 0 A 1/3 with r 0 = 1.2 fm. By keeping the lowest power of the expansion of r, the radial wave functions for the s 1/2 wave are given by [28] The normalization constant can be expressed in a good approximation as where and For this scheme, the functions F ss (E e ) and E ss (E e ), entering in the decay rate, can be safely approximated with

The Approximation Scheme B
We consider the analytical solution of the Dirac equation for the Coulomb potential of the pointlike daughter nucleus with charge Z f , V(r) = −αZ f /r. The radial wave functions are expressed as [29] with Here, Γ(z) is the Gamma function and 1 F 1 (a, b, z) is the confluent hypergeometric function.

The Approximation Scheme C
In this scheme, we consider the numerical solutions of the radial Dirac equation for a uniform charged distribution of the final nucleus. Using the numerical wave functions as solutions of the radial Dirac equation for the potential described by Equation (37), we consider the finite size effect of the final nucleus. Moreover, the screening effect of atomic electrons is considered by the Thomas-Fermi approximation. The universal screening function used φ(r) is the solution of the Thomas-Fermi equation with the boundary conditions φ(0) = 1 and φ(∞) = 0. In Equation (45) f and a 0 is the Bohr radius. For solving the Thomas-Fermi equation we implied the numerical Majorana method described in [30]. The Majorana method was also used in the 2νββ phase-space factors calculation in Refs. [31][32][33].
To describe the physical context of 2νβ − β − decay, the boundary condition on infinity should describe the final nucleus as a positive ion with electrical charge +2. Thus, the input potential for the radial Dirac equation is The radial wave functions are evaluated with the subroutine package RADIAL [34,35]. Following the subroutine procedure, the truncation errors can be avoided entirely, and the radial wave functions can be obtained with the desired accuracy. Thus, the numerical solutions can be considered as exact for the potential described by Equation (46).
The approximation schemes presented above differ by the treatment of the Coulomb potential of the final nucleus and the precision of the electronic wave functions. Although in approximation scheme A, the finite nuclear size correction is taken into account, the radial wave functions for the s 1/2 wave are approximate by keeping the lowest power of the expansion of r. In scheme B, the radial wave functions are exact, but the final nucleus is considered a point-like nucleus. The most precise treatment presented in the paper is the approximation scheme C. Besides the fact that the radial wave functions are exact, we also consider the finite nuclear size correction. Moreover, in this scheme, the screening effect of atomic electrons is considered. At this point, it is important to mention that further improvements were made in the treatment of the Coulomb potential for double beta decay in Refs. [32,33], by taking into account the diffuse nuclear surface correction. This correction is a small one compared to the finite nuclear size correction, and it can be safely neglected for the purpose of our paper.
In Figure 1, we present the radial wave functions of an electron in s 1/2 spherical wave state from the double-β emitter 136 Xe, as functions of kinetic energy E e − m e and evaluated on the surface of the daughter nucleus. We observe that the approximation scheme A, corresponding to leading finite-size Coulomb, is in good agreement with the other two approaches for g −1 (E e , R) but fails badly for f +1 (E e , R), especially at low energies of the electron. We note that approximation schemes B and C are in good agreement in the representation intervals as a final remark.

Results and Discussion
The integration over all leptons energies in the decay rate is essential when predicting the 2νββ half-life and the angular correlation coefficient of the emitted electrons. We performed the integration over the energies of the electrons numerically, with a 15-points Gauss-Kronrod quadrature. The integration over neutrino energy is performed analytically in the Appendix A.
In Table 1, we present the phase-space factors entering in the decay rate, that is, Equation (18), calculated within approximation schemes A, B, and C, for 2νββ decay of 100 Mo, 136 Xe and 150 Nd. As was already pointed out in [6], the G N phase-space factors calculated with exact relativistic electron wave functions (scheme C) are smaller than those provided by approximation scheme A. We can see that the results for G N obtained with scheme B are between the results obtained with schemes A and C. In the angular phasespace factors H N , we observe a good agreement between schemes A and B, and again smaller results from scheme C. The different behavior of G N and H N , for the same approximations schemes for the electronic wave functions, indicates that the angular correlation is sensitive to the treatment of the Coulomb interaction for the emitted electrons.
In Table 2, we update the values of the phase-space factors with the approximation scheme C for 11 nuclei of experimental interest. We took the Q-value for each nucleus from the experiment with the smallest uncertainty when available or from tables of recommended values [36], as it is the case of 2νββ-decay of 124 Sn. In what follows, we used just those values for the maximum sum of the electrons' kinetic energies.
Considering ξ 2ν 31 =−0.2, 0, 0.2, 0.4 and 0.6, we present in Table 3 the values of the angular correlation coefficient K 2ν , defined in Equation (26), for multiple nuclei of interest. For each nucleus, we display the ξ 2ν 51 fixed by the SSD relation and the energy difference between the lowest 1 + state of the intermediate nucleus and the ground state of the initial nucleus (see Equation (27)). Columns 5 and 6 in Table 3 represent the angular correlation coefficient results using the electronic wave function in approximation schemes A and C, respectively. Besides the fact that we obtain smaller results with scheme C compared with scheme A, we can also observe that the dependence of K 2ν on ξ 2ν 31 is not linear. To emphasize the dependence on the ratio of the nuclear matrix elements we depict, in Figure 2, the angular correlation coefficient K 2ν as a function of ξ 2ν 31 . The representation is done from ξ 2ν 31 = −0.8 to ξ 2ν 31 = 0.8 for all nuclei. We can see that K 2ν is a quadratic function of ξ 2ν 31 for all nuclei, at least in the representation interval. If we assume the SSD in the 2νββ process, then the unknown parameter ξ 2ν 31 is fixed to the ξ 2ν 31 SSD value. We present with filled blue circles, in Figure 2, the values of the angular correlation coefficient evaluated at ξ 2ν 31 SSD . The integration over all leptons energies is crucial in calculating the global observables. Still, the integrand itself provides valuable information about the energy and angular distributions of the emitted electrons. The analytical approach for the integration over the neutrino energy (see Appendix A) ensures that we can evaluate the obtained distributions at any energetical point, not only in the points fixed by the numerical method of integration.
From the analytical expressions of the partial double distributions in Equation (A8), we can see that they exhibit a completely different behavior as functions of the electron energies, E e 1 and E e 2 . This fact is highlighted in Figure 3, where we present, for 2νββ-decay of 100 Mo, the contour plots of the partial double electron distribution normalized to the corresponding partial decay rate. The normalization procedure ensures that the normalized partial distributions do not depend on the ratios of nuclear matrix elements ξ 2ν 31 and ξ 2ν 51 .  136 Xe and 150 Nd. The results are obtained using different approximations for the radial wave functions g −1 (E e ) and f 1 (E e ) of the electron: (A) The standard approximation of Doi et al. [28]. (B) The analytical solution of the Dirac equation for a point-like nucleus [29]. (C) The exact solution of the Dirac equation for a uniform charge distribution of the nucleus screened by the atomic electronic cloud. The results are presented in yr −1 . For each nucleus and each phase space factor we present the percent deviation between approximation schemes A and B, δ AB = 100(X A − X B )/X B , and the percent deviation between approximation schemes A and C, δ AB = 100(X A − X C )/X C , with X = G 2ν N or X = H 2ν N .  The filled blue circles indicate the values of K 2ν for the ratio of the nuclear matrix elements fixed by the SSD assumption, ξ 2ν 31 SSD . In the case of 82 Se and 116 Cd, the right filled circle correspond to 82 Se and the left one to 116 Cd. We used the approximation scheme C for the relativistic wave function of the emitted electrons. Table 3. The values of the angular correlation coefficient for different values of ξ 2ν 31 . We assume the approximation scheme A and C for the relativistic wave function of the emitted electrons. For each nuclei, we display the energy difference E(1 + ) − E i in MeV, necessary to calculate ξ 2ν 51 via Equation (27).  The angular correlation distribution κ as a function of the energies of the electrons and the ratios of the nuclear matrix elements can be found in Equation (A10). In Figure 4, we show, for 100 Mo, the contour plots of the angular distributions for different values of ξ 2ν 31 as functions of the electron energies. The left, middle and right panels are obtained for ξ 2ν 31 = −1, 0 and 1, respectively. The middle panel corresponds to the standard approximation angular distribution.

Nucleus
We compared the improved formalism presented in this paper with the exact SSD formalism of the 2νββ-decay [21,22]. If we assume the SSD, the ratios ξ 2ν 31 and ξ 2ν 51 are fixed. We present the single electron differential decay rates normalized to the full width, in Figure 5, for 82 Se, 100 Mo and 150 Nd. With a solid black line, we present the normalized single electron spectra within the exact SSD formalism [21,22], and with the blue dashed line, the spectra presented in Equation (A11). The value of the ratio ξ 2ν 31 fixed in SSD assumption can be found in the plots' legend. In the lower panels, the residuals between exact and improved formalisms are plotted. We can conclude that the spectra are in a good agreement. A small disagreement appears only for low electron energy, which is not accessible by the double-beta experiments due to large background events.   To compare the angular correlation distribution between the exact SSD formalism and the improved one, we display in Figure 6, the distribution reproduced within the exact SSD formalism. It should be noted that there are no significant deviations of this distribution from the one with positive ξ 2ν 31 from Figure 4. Considering the comparison of the single electron spectra and angular distributions, we can conclude that the improved formalism presented in this paper, with ξ 2ν 31 fixed by SSD assumption, is in excellent agreement with the exact SSD formalism. Mo. The distribution is obtained using the exact SSD formalism presented in [21,22] and the screened exact finite-size Coulomb wave functions for s 1/2 electron state.

A
It was already pointed out in [6] that the calculation of M 2ν GT−3 can be more reliable than that of M 2ν GT−1 , because M 2ν GT−3 is saturated by contributions through the lightest states of the intermediate nucleus. The 2νββ-decay half-life can be expressed with M 2ν GT−3 , ξ 2ν 31 and ξ 2ν 51 as follows: that is, without explicit dependence on nuclear structure factor (M 2ν GT − (g V /g eff A ) 2 M 2ν F ). The nuclear structure parameter ξ 2ν 31 can be deduced from the energy distribution of the emitted electrons or from the angular correlation factor K 2ν as solution of a quadratic equation with coefficients A, B and C, which are functions of the measured angular correlation factor: The fact that for a measured angular correlation coefficient there are two possible values of ξ 2ν 31 , it can also be seen graphically from Figure 2. To determine the physical solution, a cross-check with the determination of ξ 2ν 31 from the energy distribution is necessary.
GT−3 is reliably calculated, and ξ 2ν 31 is precisely obtained from the angular and energetic measurements. In that case, we can determine g eff A from the measured 2νββdecay half-life expression via Equation (47). A disagreement between ξ 2ν 31 deduced from the energy and angular distributions might indicate a new physics scenario observed in the data analysis of the 2νββ-decay.

Summary and Conclusions
The 2νββ-decay has been a subject of theoretical and experimental research for more than 85 years and remains an important topic in modern nuclear and particle physics. In the presented paper, the theoretical description of the angular distribution of outgoing electrons is achieved by considering the effect of lepton energies in energy denominators of the nuclear matrix elements via the Taylor expansion. It is claimed that by a precise measurement of angular correlation factor K 2ν of the emitted electrons, the nuclear structure parameter ξ 31 can be fixed, which might allow determining the effective axial-vector coupling constant g eff A from the measured 2νββ-decay half-life, once the nuclear matrix element M GT−3 is calculated reliably. A non-zero ξ 31 is a signature of the dominance of the Gamow-Teller over Fermi transitions in the 2νββ-decay. A more accurate treatment of the nuclear physics aspects of the 2νββ-decay will allow more reliability to address possible new physics scenarios associated with neutrino properties and interactions.
In the case of 100 Mo and 150 Nd, the NEMO-3 experiment already recorded a large number of 2νββ-events describing all kinds of differential characteristics of emitted electrons. It would be difficult to reanalyze the recorded data within the improved formalism of this paper. However, many running and planed double beta decay experiments, for example, GERDA, CUORE, KamLAND-Zen, EXO, CUPID, LEGEND, and so forth, will achieve sufficiently large statistics to fix nuclear structure parameter ξ 31 from the energy distribution of emitted electrons. The potential of the SuperNEMO experiment, in which the first module Demonstrator is in a construction phase, will be able to measure this quantity even from the angular distribution. The angular information might also be obtained with high-pressure gaseous time-projection technology and pixelated detectors; however, the present methods have significant limitations. It goes without saying that there is continuous progress in double beta decay technologies, the ultimate goal of which is to register emitted electrons with high energy and angular resolutions. Acknowledgments: The figures for this article have been created using the SciDraw scientific figure preparation system [47].

Conflicts of Interest:
The authors declare no conflict of interest.