Quark Condensate Seesaw Mechanism for Neutrino Mass

We study a mechanism of generation of Majorana neutrino mass due to spontaneous breaking of chiral symmetry (SBChS) accompanied by the formation of a quark condensate. The effect of the condensate is transmitted to the neutrino sector via Lepton-Number Violating (LNV) lepton-quark dimension-$7$ operators known in the literature as an origin of the neutrino-mass-independent mechanism of neutrinoless double-beta ($0 \nu \beta \beta$) decay. The smallness of neutrino masses is due to a large ratio between the LNV scale and the scale of the SBChS. This is a new realization of the seesaw mechanism, which we dub the Quark Condensate SeeSaw (QCSS). We show that its dominance requires a softly-broken symmetry suppressing the Weinberg operator. In the QCSS, $u$ and $d$ quarks receive their masses at (multi-)loop level and the smallness of neutrino masses is intimately related to the smallness of the $u$ and $d$ current-quark masses. We examine the predictions of the QCSS for $0 \nu \beta \beta$-decay and reveal the importance of the nuclear mater effects in this case. We find that the genuine no-fine-tuning QCSS favors the normal neutrino-mass ordering and rather narrow ranges of the neutrino masses $m_{1, 2, 3}$. We propose a realization of the QCSS with $\mathbb{Z}_4$ symmetry which predicts these features of the neutrino spectrum.


I. INTRODUCTION
The smallness of neutrino masses, in comparison with the other Standard Model (SM) fermions, remains a mystery of particle physics theory. A common wisdom suggests that this smallness is related to some broken symmetry. One of the most natural candidates is U(1) L symmetry of the Lepton (L) Number, broken at some high-energy scale Λ. Then at the electroweak scale there appears the ∆L = 2 Weinberg operator which, after the electroweak symmetry breaking (EWSB) √ 2 H 0 = v = 246 GeV, leads to the Majorana neutrino mass For a generic case with f ∼ 1 and for m ν at a sub-eV scale one estimates Λ ∼ 10 [14][15] GeV, putting Lepton Number Violating (LNV) physics far beyond the experimental reach. This happens in the tree-level realizations of the Weinberg operator (1) in the celebrated seesaw Type I, II and III, where Λ is equal to the masses M of the corresponding seesaw messengers which, being very heavy, have no phenomenological value. In order to escape this situation and open up the possibility for a non-trivial phenomenology, various models have been proposed in the literature (for a recent review see, e.g., Ref. [1]) relaxing the above-mentioned limitation on the LNV scale Λ. Introducing new symmetries (softly-broken), one can forbid the operator (1) at the tree level, while allowing it at certain loop level l, so that in (2) there appears a loop suppression factor f ∼ (1/16π 2 ) l . With the appropriate l the LNV scale Λ can be reduced down to phenomenologically interesting values in the TeV ballpark (see, e.g., Refs. [2][3][4][5][6][7] and references therein). Another possibility is to resort to symmetries forbidding (1) at all, but allowing higher dimension-(5 + n) operators which after EWSB provide an extra suppression factor (v/Λ) n . As in the loop-based models, here, for sufficiently large n, the LNV scale Λ can be made as low as the current experimental limits. In some models both loop and higher-dimension suppressions can be combined.
In the present paper we consider another class of the SM gauge-invariant effective operators Here, all the possible SU(2) L contractions are assumed. It was observed in Ref. [8] that this operator contributes to the Majorana-neutrino mass matrix due to chiralsymmetry breaking (ChSB) via the light-quark conden-sate qq = −ω 3 = 0. The latter sets the ChSB scale, so that after the EWSB and ChSB one arrives at the contribution to the Majorana mass matrix of active neutrinos where g αβ = g u αβ + g d αβ and qq ≡ uu ≈ dd ≈ 2 u L u R ≈ 2 d R d L . This is a kind of seesaw formula relating the smallness of the Majorana masses of neutrino with the large ratio between the scale Λ of Lepton-Number Violation (LNV) and the scale of chiralsymmetry breaking ω = − qq 1/3 . We call the relation (4) Quark Condensate Seesaw (QCSS) formula. Taking from a renormalized lattice QCD within the MS scheme at a fixed scale µ = 2 GeV [9] and Λ ∼ a few TeV we get the neutrino mass in the sub-eV ballpark.
In the next section we study implications of the requirement of the dominance of the operator (3) for Ul-traViolet (UV) model building and certain phenomenological aspects of the QCSS. In Section III we extract limits on the couplings of nonstandard neutrino-quark contact interactions appearing in QCSS. Then, we analyze contributions of the operators (3) to neutrinoless double-beta (0νββ) decay and derive strong limitations on the QCSS mechanism from this LNV process.

II. DOMINANCE OF QCSS AND LIGHT-QUARK MASSES
Here we discuss the conditions for the dominance of the operator (3) in the Majorana-neutrino mass matrix. As usual, this can be guaranteed by imposing on the theory an appropriate symmetry group G which could be either continuous or discrete. General properties of this kind of symmetries were studied in Ref. [8]. This symmetry must forbid the Weinberg operator (1), but allow the operator O q 7 in Eq. (3). Therefore, the lepton bilinear LL must be a G non-singlet. Requiring that G remains a good symmetry after the EWSB and still forbids any contribution to the Majorana-neutrino mass term while allowing the quark-lepton coupling implies that we claim the SM Higgs H to be a G singlet. Thus, the condition of G invariance of the operator (3) requires that one of the quark bilinears (Q u R ) and (d R Q) or both be G non-singlets. The latter implies that the Yukawa couplings of the u and d quarks are not G-invariant and forbidden by this symmetry. Therefore, the light quarks do not receive their masses as a result of the EWSB. This is in line with the fact that the light quarks u and d are particular among other quarks by being much lighter than the others. However, the statement of vanishing masses m u,d = 0, or even one of them, contradicts the experimental observations. Therefore, small m u,d = 0 must be generated in some way to make the scenario we are discussing viable. In principle, for this scenario it is not necessary that both Yukawa couplings in Eq. (8) are forbidden. As seen from (3), it is sufficient that only one of them, say the u-quark Yukawa coupling, is forbidden as suggested in Ref. [8].
First, we note that in a generic effective theory the light-quark masses can be generated due to chiral symmetry breaking via the effective SM-invariant operators with q = u, d. Here, Λ q is some scale of the physics underlying these operators. All the three operators in Eq. (9) among their components have necessary for the light-quark mass m u,d generation. Thus, after the chiral symmetry breaking we get Note that the sign in the last equality is irrelevant and can be absorbed by a redefinition of κ.
The four-quark operators like (9) are well known in context of the Nambu-Jona-Lasinio model considered as chiral low-energy effective theory of QCD. Recall that in this approach the one-gluon exchange diagram with the amplitude converts to a point-like quark-lepton operator in a truncated theory where the gluon propagator D (g) (k 2 ) is replaced with g µν /Λ H . Here, Λ 2 H is some characteristic hadronic scale. After Fierz rearrangement in (12) one finds the operators shown in Eq. (9) with Λ 2 q ∼ α s /Λ 2 H . After the spontaneous breaking of chiral symmetry they render a contribution (11) to the masses of the u and d quarks, converting them to the so-called constituents quarks. However, chiral-symmetry breaking cannot be the only source of the quark masses. It is well known that quarks must acquire the so-called current-quark masses, not related to the chiral symmetry. These masses are needed for explicit breaking of chiral symmetry, otherwise, according to the Gell-Mann-Oakes-Renner relation (B8), pions as Goldstone bosons would be massless. Therefore, we need a source of the u and d current-quark masses.
In the present scenario, this source is provided by soft breaking of the G symmetry. The latter allows the couplings (8) at certain loop level, depending on the selected UV completion of the effective theory. In this case, the operator (8) acquires the corresponding loop suppression like the one mentioned in the Introduction. This loop contribution is finite (UV convergent), since the divergence-compensating counterterms of the pointlike form (8) are still forbidden by G. The loop-generated m u,d would be in this case proportional to a dimensionful soft-symmetry breaking parameter µ, which could be made as small as necessary to suppress m u,d in comparison with the other quarks. The smallness of µ would be technically natural due to the G-symmetry restoration in the limit µ → 0. On the other hand, a softly broken G symmetry also allows loop generation of the Weinberg operator (1) and, as a result, generation of small Majorana neutrino mass after EWSB. The loop suppression allows for the LNV scale Λ to lie in the TeV region, as we discussed in the Introduction. The idea of the loopgenerated hierarchy of the SM fermion masses where all the quark and lepton masses, except for the t quark, are generated at loop level underlies a popular trend in model building (for a recent realization see, e.g., Ref. [7,[10][11][12]). Here, we assume that there is a class of UV models where the Weinberg operators (1) are strongly loop-suppressed, so that the contribution of the operators (3) to the neutrino mass matrix (4) dominates.

III. LIMITS ON LNV LEPTON-QUARK INTERACTIONS
Let us examine phenomenological limits on strength of the effective LNV lepton-quark interactions from (3) and (7). The literature lacks the limits on this class of nonstandard interactions. We will comment on the existing limits at the end of this section.
Let us introduce the dimensionless parameters giving a measure of the relative strength of the fourfermion interactions (7) with respect to the Fermi con- Here, assuming the dominance of the QCSS in the Majorana mass of neutrinos, we can extract limits on ε + αβ from the neutrino-oscillation data since, according to Eq. (4), they are directly related to the elements of the neutrino mass matrix Using Eqs. (A5)-(A10), we relate these LNV leptonquark parameters to the neutrino-oscillation parameters.  15) vs. the lightestneutrino mass m0. The allowed gray regions between the curved lines are derived using current best-fit values of the neutrino-oscillation parameters [13] and CP phases δ ∈ [0, 2π) and α1, α2 ∈ [0, π). The solid and dashed lines refer to the normal (NO) and inverted (IO) orderings of the neutrino masses, respectively. The vertical gray bands represent the regions excluded by Planck measurements at 95% C.L. [14]. The horizontal gray bands in the first plot correspond to the regions excluded by KamLAND-Zen 90% C.L. limits on 0νββdecay [15]. The solid and dashed horizontal limiting lines refer to the scenarios with and without the nuclear-matter effect discussed in Section IV, respectively.
The current values of the latter we take from Ref. [13]. Then, varying the CP phases in the intervals δ ∈ [0, 2π) and α 1 , α 2 ∈ [0, π), we find the exclusion plots in the plane (m 0 , |ε αβ |) shown in Fig. 1 for the best-fit values of the neutrino oscillation parameters θ ij and ∆m 2 As is known, the upper cosmological limit on the sum of neutrino masses i m i < 0.12 eV at 95% C.L., set by the Planck measurements [14], imposes the limit on the mass m 0 of the lightest neutrino. Applying Eqs. (A3)-(A4), one finds for the normal (NO) and inverted (IO) neutrino-mass orderings, respectively. These limits are shown in Fig. 1 as vertical bands.
It is instructive to show allowed ranges (ε min αβ ; ε max αβ ) of the LNV lepton-quark interaction parameters ε αβ from (15). These ranges extracted from the exclusion plots in Fig. 1 The following comment is in order. Assuming that there is no significant cancellation between the two terms in (14), we can extend the above limits to the parameters ε u and ε d separately. However, these individual limits should be treated as order-of-magnitude estimates, since a priori there is no reason to expect the absence of the said cancellation.
To the best of our knowledge, the only analysis of phenomenological limits on the lepton-quark interaction strength ε from Eq. (13) existing in the literature is given in Refs. [8,16] where the SN 1987A and meson decays were studied. These limits are in the range ε < 10 −3 , which is an order of magnitude weaker than our limits in Eqs. (18)- (19).

IV. QUARK CONDENSATE SEESAW IN NEUTRINOLESS DOUBLE-BETA DECAY
The operators (3) were previously studied in the literature as a source of ∆L = 2 interactions able to induce the 0νββ decay with no explicit dependence on the Majorana neutrino mass [17][18][19]. After EWSB the interactions relevant to the 0νββ decay generated by these operators are where ε u,d ee are defined in Eq. (13). Let us examine the contribution of the effective ∆L = 2 interaction terms (20) and (21) to the 0νββ decay. The first term (20) combined with the SM weak charged-current (CC) interaction leads to the contribution shown in Fig. 2(b), which is independent of neutrino mass in the propagator due to chiralities in the vertices P L ( / q + m ν ) P R = / q. This is a manifestation of the fact that the ∆L = 2 is not provided by m ν , but solely by the upper vertex in Fig. 2(b). In our QCSS model, the term (21) also contributes to the 0νββ decay via the neutrino-mass mechanism shown in Fig. 2(a) due to the chiral symmetry breaking and formation of the quark condensate. As we discussed in the previous sections, the term (21) is the only source of the neutrino mass in the present model. However, there is a Figure 2. Contributions of the effective operators (3) to the 0νββ decay.
subtlety with the diagram in Fig. 2(a). It describes a process taking place in the nuclear environment, where the chiral quark condensate qq N is suppressed with respect to the one in the vacuum qq . As briefly discussed in Appendix B, this suppression is estimated to be a factortwo effect (B9). In what follows, we take this fact into account.
The inverse 0νββ-decay half-life in the QCSS model reads where g A ≈ 1.27 and m e are the axial-vector weak nucleon coupling constant and electron mass, respectively. The standard kinematical phase-space factor G 0νββ (Z, Q) is given in Table I. We introduced the ratio of the nuclear matrix elements (NME) which are defined in Appendix C. In the second row of Eq. (22) we used the prediction of the QCSS model Note that in the QCSS m ββ = m ν ee in contrast to the conventional neutrino mass models where these two quantities are equal in the diagonal charged-lepton basis. In the QCSS we have instead where m ν ee is given by (4) with the chiral condensate in the vacuum. We shall discuss the possible implications of this prediction in a while.
With the values of NME for 136 Xe given in Appendix C, quark condensate in matter (B9) with (5) and other known parameters we find From the currently most stringent upper bound on the 0νββ-decay half-life obtained for 136 Xe by the KamLAND-Zen experiment [15] T 0νββ 1/2 > 1.07 × 10 26 yr at 90% C.L., we find, using (22), an upper bound on the model parameters Let us consider several characteristic QCSS scenarios.
(i) In this scenario, the diagram in Fig. 2(a) dominates. This implies that a ε − ee b ε + ee and happens in the case of strong fine-tuning of the model parameters as follows from (14) and (26). Then, from (28) we find The region excluded by this limit is shown in the first plot of Fig. 1 as a horizontal band delimited by the solid line. The dashed line shows a limit without the nuclearmatter effect, which is apparently stronger by a factor of two.
It is also instructive to show the effect of nuclear matter in the conventional plane (m 0 , |m ν ee |), and this is done in Fig. 3 including constraints on the neutrino-oscillation parameters for both NO and IO. From (30) with (24) one finds |m ββ | < 116 meV.
Neglecting the nuclear-matter effects discussed above, one would have m ββ = m ν ee and, as a result, the same upper limit for |m ν ee | as for |m ββ | in (31); this limit is shown in Fig. 3 as the dashed horizontal line. However, the presence of nuclear matter is responsible for a factor-two effect (B9), and consequently the limit on |m ν ee | is twice as large as (31); the corresponding limit is displayed as the solid horizontal line. This is one of the predictions of the QCSS-model: the 0νββ-decay constraints in the plane (m 0 , |m ν ee |) are in fact by about a factor of two weaker than expected in other models of the neutrino mass. However, this effect is diluted by the comparable uncertainties in the theoretical values of nuclear matrix elements (for a recent review see, for instance Ref. [20]) and electron-shell overlap factors [21].
(ii) Let us examine the genuine QCSS scenario with no fine-tuning as in (29). In this case, an order-of-magnitude relation ε + ee ∼ ε − ee holds and the diagram in Fig. 2(b) is dominant because a |b| according to (26). Then, from (28) we find  Thus, in most of the model-parameter space, excluding the narrow fine-tuning domain (29), the contributions of both diagrams in Fig. 2 are roughly determined by the unique parameter ε + ee . In a while, we will discuss a possible realization of such scenario.
Comparing the limit (32) for ε + ee with the excluded regions in the first plot of Fig. 1, derived from the neutrinooscillation data, we conclude that the genuine QCSS predicts the NO and the range for the lightest-neutrino mass in a rather narrow range for the lightest neutrino mass 2.65 meV < m 0 = m 1 < 6.85 meV.
With these ranges for the neutrino masses, we find the corresponding range for the cosmological neutrino parameter We also determine the genuine QCSS range for the singlebeta-decay parameter which is beyond the reach of the current and near future tritium beta-decay experiments (for a recent review see, for instance, Refs. [22]). At last, the limit (32) translated into the 0νββ-decay parameter gives the upper bound |m ββ | < 6.8 × 10 −4 meV. (38) It is worth recalling that in the considered case of the QCSS this parameter characterizes the subdominant contribution to the 0νββ decay shown in Fig. 2(a), while the dominant one is given by the diagram Fig. 2(b).
(iii) Finally, let us consider a particular model in the QCSS framework, realizing the scenario (ii) discussed above. Instead of letting both quark bilinears (Q u R ) and (d R Q) to be G non-singlets, we allow this assignment only for one of them. We choose the (Q u R ) to be a non-singlet while leaving (d R Q) to be a singlet. Then, only one of the two operators in Eq. (3) is allowed, namely L L H Q u R . Therefore, only the ε u ee parameter survives in all the above equations, so that for the parameters in (13)- (15) we have ε + ee = ε − ee ≡ ε ee . Thus, in this particular model there is one parameter ε ee which determines the neutrino mass-matrix element m ν ee and the 0νββ-decay half-life according to Eqs. (4) and (22), respectively. Consequently, due to the limit (32), this model predicts the NO of the neutrino mass spectrum and the range of the neutrino masses shown in Eqs. (33)-(35).

V. PARTICULAR REALIZATION OF QCSS MODEL
There is one potential flaw in the model described above: the quark bilinear (d R Q), being a G singlet, allows the d-quark Yukawa term shown in (8). It leads to a tree-level d-quark mass after the EWSB, which makes the smallness of the d-quark mass rather weird. If we wish for both light quarks, u and d, to acquire their masses at the loop level and, on the other hand, only one lepton-quark operator O u 7 in Eq. (3) to survive, then we should take special care choosing an appropriate symmetry G. Under this symmetry, both quark bilinears should be nonsinglets while it must forbid the operators O W and O d Let us give an example of such a group G. One of the simplest cases resorts to the cyclic group G = Z 4 with the field assignments Here, we limit ourselves only to the first generation of the fermions. As seen from (39), this Z 4 model features the following properties. The Yukawa couplings of the light quarks (8) in Eq. (9) and the similar lepton-quark operator These operators lead to the corresponding contributions of the quark condensate to the u-and d-quark masses (11) and the electron mass. As we already argued in Section II, there must be additional contributions, at least to the u-and d-quark masses. Therefore, we postulated the Z 4 subgroup to be softly broken, so that these masses could receive small loop-level contributions.

VI. CONCLUSION AND DISCUSSIONS
We studied the Quark Condensate SeeSaw (QCSS) mechanism of generation of the Majorana mass matrix of neutrinos due to the spontaneous breaking of chiral symmetry and transmission of the effect of formation of the chiral condensate to the neutrino sector via the dimension-7 quark-lepton operators O u,d 7 in Eq (3). These operators can originate in the low-energy limit of certain class of UV models. We analyzed generic conditions in these models and found that they need to have a softly-broken symmetry, G, suppressing the Wilson operator down to some loop level in order for the QCSS mechanism to dominate over the ordinary tree-level Majorana neutrino mass generated by the electroweak symmetry breaking. We recognized that in this scenario the u and d quarks are naturally light in comparison with the other quarks and acquire their masses at loop level, as a result of the soft breaking of the model symmetry G. The order of the loop suppression depends on the particular UV model. We postpone the study of possible UV completions of the QCSS scenario for a subsequent publication. We also noted that the u-and d-quark masses always receive a contribution proportional to the chiral condensate via four-quark operators generated by nonperturbative QCD effects which convert these quarks to the constituent ones.
We calculated the predictions of the QCSS scenario for the LNV lepton-quark couplings (13)- (14) of the operators in Eq. (3). These couplings characterize both nonstandard neutrino-quark and charged lepton-quark interactions, which could be relevant for further studies of the phenomenological and astrophysical implications of the QCSS mechanism of the Majorana-neutrino-mass generation.
We analyzed predictions of the QCSS model for the neutrinoless double-beta (0νββ) decay. We calculated the corresponding nuclear matrix elements within the Quasiparticle Random Phase Approximation (QRPA) method with partial restoration of the isospin symmetry. Then, we studied particular QCSS scenarios with the dominance of either the neutrino-mass mechanism, Fig. 2(a), or the neutrino-mass-independent mechanism, Fig. 2(b). The first one implies a strong fine-tuning of the model parameters, while the second should be treated as a genuine QCSS model with no fine-tuning. In the first case we stressed the importance of taking into account the nuclear-matter effects, in order to properly compare the 0νββ-decay constraints on |m ν ee | with the corresponding limits derived from the neutrino-oscillation data. We showed that the genuine QCSS scenario, taking place in the largest part of the model-parameter space, disfavors the inverted ordering (IO) of the neutrino-mass spectrum and predicts the normal ordering (NO) with the lightestneutrino mass m 0 constrained to a rather narrow range (33). This is in accord with recent global analysis of the neutrino-oscillation parameters which favors NO over IO at more than 3σ [13]. In view of this fact we constructed an example of the QCSS model with softly-broken Z 4 symmetry, which predicts NO and neutrino masses in the ranges (33)-(35). We postpone a more detailed study of other possible realizations of the QCSS model as well as its UV completions.
We employ these relations for the analysis of the nonstandard neutrino-quark couplings (7) and (13), which is done in Section III. In this analysis, we use the neutrinooscillation data from Ref. [13].

Appendix B: Quark Condensate
Chiral symmetry is approximate invariance of the QCD Lagrangian under the global SU(3) L × SU(3) R transformations in the space of the lightest quark flavors q = u, d, s. Below the chiral scale 4π f π ∼ 1 GeV, this symmetry is spontaneously broken by the light-quark condensates 0|qq|0 ≡ qq = 0 in the QCD ground state (vacuum) |0 . The corresponding Goldstone bosons are the pions π. Their nonzero mass originates from the explicit breaking of the chiral symmetry by the light currentquark mass terms in the QCD Lagrangian/Hamiltonian: where we separated the isospin-singlet and isospin-triplet quark combinations. In what follows, we retain the singlet and consider only the lightest u and d quarks, denoting m q = 1 2 (m u + m d ) and qq = 1 2 (uu + dd). Here, we will examine the effect of nuclear environment on the formation of light quark condensate. Following Ref. [23], we use the Hellmann-Feynman theorem, allowing one to analyze the condensates in a modelindependent way to the first order in nucleon density. The Hellmann-Feynman theorem states where |ψ(λ) and E(λ) are, respectively, the normalized energy eigenstates and eigenvalues of the Hamiltonian H(λ) with explicit dependence on the parameter λ. Choosing λ = m q and H = d 3 x H m we get where both parts of this equation are multiplied by m q to ensure renormalization-group invariance of this relation [24]. Let us consider two cases |ψ(m q ) = |0 , |ρ N , where |0 is the QCD vacuum and |ρ N refers to the ground state of (uniform) nuclear matter at rest with nucleon density ρ N . Writing Eq. (B3) for these two cases, we subtract one from the other and obtain where E N is the energy density of nuclear matter. Provided the kinetic and potential energy of nucleons are known to be small, one has On the other hand, there is a current-algebra relation [25] where σ N is the pion-nucleon sigma term measuring the nucleon-mass m N shift from the chiral limit m u,d → 0. Then, using Eqs. (B4), (B5) and (B6), one finds a modelindependent relation [23] qq N qq characterizing the amount of chiral-symmetry restoration in dense medium. Here, we denoted qq N ≡ ρ N |qq|ρ N and qq ≡ 0|qq|0 . The Gell-Mann-Oakes-Renner relation [26] 2m q qq = −f 2 π m 2 π , has been used to derive Eq. (B7). In order to estimate the nuclear-matter effect on the quark condensate on the basis of Eq. (B7), we adopt the usual value for the nucleon density ρ = ρ p + ρ n = 0.17 fm −3 , the recent large value of σ = 64 MeV from a partial-wave analysis of the π-N scattering [27], f π = 92 MeV and the charged-pion mass m π = 140 MeV. Then, Eq. (B7) yields qq N ≈ 0.5 qq , demonstrating a substantial suppression of the quark condensate in the nuclear matter. The value qq N can be interpreted as the sum of scalar densities of the u (or d) quarks in vacuum and inside nucleons. The nucleon component of qq N is estimated in Ref. [16] to be ≈ (100 MeV) 3 and needs not be included separately. The sign of the nucleon component is opposite to the sign of the vacuum component; the latter is also numerically higher.
In the Quasiparticle Random-Phase Approximation (QRPA) method, M ν,ε are written via sums over the virtual intermediate states labeled by their angular momenta and parities J π and indices k i and k f [20,28] q m π , h ε GT (q 2 ) = 1 3 In the above formulas, m p and m π denote the proton and pion masses, respectively. For the nucleon form factors we use the following dipole parameterization: with q 2 = q · q (a small energy transfer in the nucleon vertex can be safely neglected), M V = 850 MeV and M V = 1086 MeV, and the normalization constants g V = 1, g A = 1.27 and g M = (µ p − µ n ) g V = 3.70. We used the values F P (0) = 4.41 from Ref. [29]. The induced pseudoscalar form factor is given by the PCAC relation g P (q 2 ) = 2m p q 2 + m 2 π g A (q 2 ). (C9) In Table I, we show the nuclear matrix elements calculated within the QRPA method with partial restoration of the isospin symmetry. More details on the calculation of these and other NMEs relevant for physics beyond the SM will be published elsewhere. It is worth noting that the main contribution to M ε comes from the Fermi nuclear matrix element M ε F . Table I. Nuclear matrix elements M ε , M ε F , M ε GT and M ε T calculated within the Quasiparticle Random-Phase Approximation (QRPA) method with partial restoration of the isospin symmetry, Argonne V18 nucleon-nucleon potential and the unquenched value of the axial-vector coupling constant gA = 1.27 [28]. The NMEs M ν , M ν F , M ν GT and M ν T were derived in Ref. [28]. The phase-space factors G 0νββ (Z, Q), taken from Ref. [30], are expressed in the units of 10 −14 yr −1 .