Seesaw neutrinos with one right-handed singlet field and a second Higgs doublet

We study parameters of an extension of the Standard Model. The neutrino sector is enlarged by one right-handed singlet field, allowing for the seesaw mechanism type-I, and the Higgs sector contains one additional doublet, which contributes to light neutrino masses through one-loop radiative corrections. Employing an approximation for the effective light neutrino mass matrix we express the masses of the light neutrinos analytically, allowing us to parametrize the Yukawa couplings to neutrinos by the experimental measurements on the neutrino sector and only two free parameters. We focus on a CP-conserving Higgs potential for which we present the allowed ranges of the input parameters and a statistical overview over the possible values of the Yukawa couplings.


Introduction
The precise interpretation of the neutral lepton fields in the particle physics Lagrangian is not settled yet, owing to the very small mass of the known neutrinos and the weakness of their interaction with other particles [1]. The observed neutrino oscillations support the notion that neutrinos have non-vanishing masses, calling for a modification of the Standard Model (SM). The size of the neutrino mass is not the only puzzle to solve. Absence of an electrical charge allows neutrinos to be their own antiparticles. The nature of the neutrinos -whether they are Dirac or Majorana particles -might be determined by future experiments [2]. For the experimental constraints see [3][4][5].
The Standard Model considers neutrinos as massless. Adding heavy right-handed neutral singlets and additional Higgs doublets, the authors of ref. [6] combined the seesaw mechanism (type-I) with the radiative mass generation. The spontaneous symmetry breaking of the SM gauge group leads to a Dirac mass term for neutrinos. The assumption that neutrinos are Majorana particles allows an additional term in the Lagrangian, namely, the Majorana mass term for the heavy singlets.
The model parameters allow small masses of the light neutrinos that are compatible with the experimental observations. We use this model in the formulation of Grimus and Lavoura [7,8], restricting the number of additional Higgs doublets to one. The case of three additional heavy neutrino fields was studied e.g. in refs. [9,10]. We assume only one heavy neutrino field and consider only 1-loop corrections to the neutrino mass matrix. Ibarra and Simonetto [11] analysed this scenario in the decoupling limit by renormalization group methods and predicted qualitatively our quantitative results. Our preliminary results were presented at several conferences [12][13][14][15][16]. This paper provides a more complete description of the performed numerical analysis. We reduce the number of free model parameters by linking the model predictions with experimental neutrino observables.
Our extended model has several subsets of parameters. The neutrino sector is characterized by the mass of the heavy neutrino and the strength of the coupling to the neutral Higgs fields. The masses of the three light neutrinos are the result of our model parameters. They are subject to experimental constraints, namely the experimental neutrino mass differences, ∆m 2 21 and ∆m 2 31 , as well as the experimental neutrino oscillation angles θ 12 , θ 13 , and θ 23 [17]. We follow the ideas from ref. [18] on neutrino oscillation angle estimation from the neutrino mixing matrix. More details are given in appendix B. It should be noted that experimental data is usually interpreted in the "3 × 3" neutrino mixing model [1,17], i.e. three flavoured neutrinos are considered as mixed states of three neutrino mass eigenstates. We do not attempt to reinterpret the experimental results in the context of an extended neutrino model.
We parameterize the Higgs sector following the analysis of Haber and O'Neil [19]. The Yukawa couplings are parameterized similarly to Grimus and Lavoura [7,8], which coincide with [19] in the Higgs sector. For the numerical analysis we take the mass of the SM-like Higgs boson as m h = 125.18 GeV [1] and allow the masses of two other neutral Higgs bosons to vary in the range from m h to 3000 GeV.
Using cosmological arguments the PLANCK collaboration finds [20] that the sum of all light neutrino masses is limited by m ν < 0.12 eV. The earlier upper bound estimate was notably larger: m ν < 0.23 eV [21]. If the new bound is correct, the overall scale of the neutrino masses must be smaller, and the mass of the lightest neutrino could be much smaller than the masses of the other neutrinos, especially for the inverted hierarchy. As a matter of fact, the lightest neutrino has no mass in the model of ref. [6] with only one heavy neutrino. We call this setup the Grimus-Neufeld model. However, this Grimus-Neufeld model is compatible with the results of [20,22] and is fully consistent with the current experimental neutrino data.
The outline of the paper is the following. Section 2 reviews the seesaw mechanism and the formalism of the two-Higgs-doublet model as it is used in our analysis. Section 3 shows the analytic determination of the neutrino masses that can be used to replace free model parameters by the measured neutrino mass differences and mixing angles. Section 4 describes our main results, namely, the analysis of the free model parameters and restrictions for the Higgs sector. Our findings are summarized in section 5. For completeness, appendix A describes the features of the weight vectors b i that relate the scalar Higgs fields to their mass eigenfields, appendix B gives the details of the oscillation angle calculation, and appendix C summarizes the restrictions that we apply to the parameters of the 2HDM potential.

Description of the model
We discuss an extension of the Standard Model with enlarged Higgs and neutrino sectors. Our main interest is the neutrino sector. Since we need the Higgs sector for the radiative neutrino masses, we give a short overview of the properties of the Higgs sector that we use in our calculations.

The Higgs sector
The authors of ref. [23] discuss the basis independent formulation of the general two-Higgsdoublet model (2HDM). Using their definition of the Higgs basis, we can write the two complex doublets of our model in a unique way where the vacuum expectation value (VEV) v 246 GeV and the Goldstone bosons G 0 and G + appear only in the first Higgs doublet φ 1 . The tree-level relations between the basis independent parameters defining the Higgs potential and the parameters describing the physical states are linear and can be easily inverted. This feature allows us to use the VEV, the masses of the physical Higgs bosons, m H 0 1 , m H 0 2 , m H 0 3 , and m H + , and their mixing angles ϑ 12 and ϑ 13 as input parameters.
The mass eigenstate for the charged Higgs boson corresponds directly to the field H + with the mass m H + , but the mass eigenstates for the neutral Higgs bosons with the masses m H 0 1 , m H 0 2 , and m H 0 3 , respectively, are linear superpositions of the neutral fields H 0 1r , H 0 2r , and H 0 2i . Following the formulation of Grimus and Lavoura [7,8] these linear superpositions are conveniently expressed by where φ 0 are the neutral parts of the Higgs doublets without the VEV: φ where n H is the number of Higgs doublets, i.e. n H = 2 in the 2HDM. We discuss those vectors in the general case in appendix A. There we also show how to obtain the following parametric values for the vectors b k : Table 1. Basis-independent conditions for a CP-conserving 2HDM scalar potential and vacuum [19]. ϑ ij are the mixing angles of the neutral Higgses and β − α is the invariant angle constructed from the angle α which mixes the CP-even Higgs bosons and the angle β which relates the values of the VEV's; ε ≡ sgn(β − α) is a pseudo-invariant quantity; m H and m A denote the masses for the CP-even and CP-odd Higgses. Relations between neutral Higgs fields and angular factors are explained in more detail in appendix C and in ref. [23]. Our case I corresponds to the case I of [19], whereas our case II corresponds to the case IIa of [19].
where c 1j = cos ϑ 1j and s 1j = sin ϑ 1j (j = 2, 3) are determined by the angles ϑ 1j that describe the mixing of the neutral Higgs fields. Restricting ourselves to the CP conserving case we use the analysis of ref. [19], where the authors discuss the CP-invariant Higgs potential in the 2HDM framework under various basis-independent conditions. The possible overall phase, that can be written in front of the second Higgs doublet and that acts like a mixing angle ϑ 23 between H 0 2r and H 0 2i , is used to define the CP-property of the mass eigenstates, corresponding to their coupling to gauge bosons. The choice is H 0 to be CP-even and A 0 to be CP-odd.
This assignment does not order the masses of the neutral Higgs bosons: both m H < m A and m H > m A are possible, giving us two conditions (case I and case II), which are listed in table 1. We still assume the fixed SM-like Higgs mass m H 0 1 ≡ m h to be smaller than the other two: m h < m H,A . Authors of [23] argue that one can assume − π 2 ϑ 12 , ϑ 13 < π 2 without the loss of generality. We perform the numerical analysis of the neutrino mass spectrum considering the named two cases, but using only the single mixing angle (β − α).

The Yukawa couplings
Using the vector-and-matrix notation, the Yukawa Lagrangian for the leptons is expressed [7,8] as The quantities R and ν R are the vectors of the right-handed charged leptons and the right-handed projection of the neutrino singlets, respectively. L and ν L form the lepton doublet under the weak interactions and combine with the Higgs doublets φ k to form SU (2) weak -invariant terms. They are also vectors in the generation space of dimension n L = 3. The Yukawa coupling matrices Γ k have the dimension n L × n L , while ∆ k have the dimension n R × n L , where n R is the number of the singlet neutrino fields, n R = 1 in our case.
Taking the bilinear terms of eq. (2.4), which means taking only the VEV from the Higgs doublets, we get the Dirac mass terms for charged leptons and neutrinos, assuming the charged leptons to be in their mass eigenstates: and These matrices have to be diagonalized using the singular-value decomposition (SVD) like in the SM to get the correct definition for the mass eigenstates that will describe the physical particles. Having done this transformation to the mass eigenstates, which we write down as the fields appearing in eq. (2.4), the respective transformation matrices reappear in two unique combinations, V CKM and V PMNS , in the interactions with the charged gauge bosons W ∓ or the charged scalar bosons H + and G + , giving the charged current Lagrangian where g is the SU (2) gauge coupling constant and ζ stands for the neutrino mass eigenstates. We give this part of the Lagrangian only as a reference, to show what neutrino experiments measure, as this PMNS matrix V PMNS is the basis for the interpretation of experimental data in the "3 × 3" neutrino mixing model [1].

Neutrinos at tree level
The singlet neutrinos, added to the SM, are neutral with respect to all gauge groups of the SM. This offers the possibility that they are Majorana particles, allowing to write a Majorana mass term for them. Since the Lagrangian has to be a scalar with respect to Lorentz transformations, we have to combine a spinor with itself in a Lorentz invariant way. The Dirac spinors can only be combined using the charge conjugation matrix C, which also appears in the definition of the Lorentz covariant conjugation 1 where Ψ is a Dirac spinor. The Majorana condition can now be written aŝ where η Ψ is the Majorana phase. Assuming ν R to be n R Majorana fermions we can write down the Majorana mass term as where the order of M R and C is irrelevant, as these matrices act on different indices of the spinor ν R : C is a 4 × 4 matrix, connecting the spinor indices of ν R , whereas M R is a symmetric n R × n R matrix, acting on the "generation" index of ν R . Since in our case ν R = 1, the Majorana mass matrix of the heavy singlet M R is just a number. The mass terms for the neutrinos, including the Dirac mass terms originating from the Yukawa terms in eq. (2.4), can be written as and can be written in a compact form by introducing the (n L + n R ) × (n L + n R ) symmetric neutrino mass matrix The Majorana mass matrix of the light neutrinos is vanishing at tree level, M L = 0. The neutrino mass matrix M ν can be diagonalized [6][7][8] using the properties of the singular-value decomposition of a symmetric matrix, or Takagi factorization [25] U M ν U =m = diag (m 1 , m 2 , m 3 , m 4 ) , (2.13) where m i are real and non-negative. Following the conventions of [17] we adopt the massordering m 1 ≤ m 2 < m 3 m 4 for the normal hierarchy and m 3 ≤ m 1 < m 2 m 4 for the inverted hierarchy of the neutrino mass spectrum. In order to implement the seesaw mechanism [26,27]  and the mass of one light neutrino that is generated by the seesaw mechanism. We will refer to it as the "seesaw neutrino" ζ tree s with the mass m tree s . (The neutrino states in the mass basis are denoted as ζ to distinguish them from the flavour eigenstates denoted as ν.) The remaining two neutrino states are massless at tree-level. Since the radiative corrections [6] generate only one mass, one of these two states will stay massless. We call this state ζ o with the mass m o = 0. The seesaw neutrino ζ s has the mass m s . The remaining third light neutrino ζ r has the mass m r . As argued in ref. [6], the loop generated (i.e. radiative) mass m r can be of the same order as the seesaw generated mass m tree s . Hence we do not impose an ordering between these two states (m s and m r ). Combining these two possibilities of the ordering with the normal or inverted hierarchy we  Table 2. Index arrangements between the naming and numbering of the light neutrino states. The overbarred scenarios describe the case, when the loop-generated mass m r becomes bigger than the loop-corrected seesaw mass m s . The mass m o is always 0 in our model.
can have four arrangements of indices between the names o, r, and s, and the numbers 1, 2, and 3, as displayed in Table 2. Since the formulation of the theoretical basis does not care about the numbering, we stay with the names and refer to Table 2 only when implementing the physical values.
It is useful to decompose the (n L + n R ) × (n L + n R ) unitary matrix U from eq. (2.13) into two submatrices [6][7][8] where the submatrix U L is of size n L × (n L + n R ) and the submatrix U R is n R × (n L + n R ). These submatrices obey certain unitarity relations: (2.15) Combining with eq. (2.13), we can obtain the following relations: With these submatrices of U , the left-and right-handed neutrinos can be written as linear superpositions of the n L + n R physical Majorana neutrino fields ζ α (to the remainder of this section, we omit the superscript "tree"): where P L and P R are the projectors of chirality.
Switching to the physical Majorana mass states ζ, we have to express the field couplings using the matrices U L and U R . Neutrino interaction with the Z boson is given by where c w is the cosine of the Weinberg angle. The Yukawa couplings for the neutral scalars take the form where we treat the Goldstone boson G 0 as H 0 4 . The Yukawa coupling ∆ b k is the result of rewriting the Yukawa Lagrangian eq. (2.4) using the physical Higgs fields defined in eq. (2.2): (2.20) The tree level quantities are used to calculate 1-loop corrections.

Loop corrections to the neutrino masses
We are interested in radiatively generated neutrino masses at one-loop level [7]. The light neutrino Majorana mass term δM L has the largest influence from the corrections to the neutrino mass matrix, since this submatrix is zero at tree level, M L | tree = 0. The contributions to the masses from charge-changing currents are subdominant [7,8,28].
Once the one-loop corrections are taken into account, the neutral fermion mass matrix is given by [7] The one-loop corrections to δM L originate via the self-energy functions Σ S(X) L (0) (where X = Z, G 0 , H 0 k , k = 1, 2, 3) that arise from the self-energy Feynman diagrams. The contributions Σ S L (p 2 ) are evaluated at zero external momentum squared (p 2 = 0). The neutrino couplings to the Z, Higgs H 0 k and Goldstone G 0 bosons are determined by eqs. (2.18) and (2.19). Each diagram contains a divergent piece but the sum of the three contributions yields a finite result. The expression for these one-loop corrections is given by (see e.g. [7]) where the sum index k runs over all neutral physical Higgses H 0 k . The 1-loop corrections are defined in terms of tree level quantities.

Parameters of the model
As the Grimus-Neufeld model is a minimal extension of the Standard Model, the only additions to the Lagrangian of the Standard Model are the heavy singlet Majorana mass term, eq. (2.10), the Yukawa couplings to the heavy singlet fermion, ∆ j , the Yukawa couplings of the second Higgs doublet to the charged leptons, Γ 2 , both given in eq. (2.4), and the Higgs potential of the two Higgs doublets, that replaces the Higgs potential of the Standard Model. That gives us as the primary parameters of our model. p i,SM denotes the SM parameters like the masses of the charged leptons or the Fermi coupling constant G F . p i,2HDM stands for the parameterization of the 2HDM potential and can be either the potential parameters m 2 ij and λ j or, following the idea of [29,30], the masses and physical couplings of the Higgs fields. It means also that we assume the charged fermion fields to be in their mass eigenstates, making Γ 1 =

Reducing parameters by neutrino measurements
The main goal of this section is to show, how we can replace the 6 complex parameters in ∆ 1 and ∆ 2 by the measured mass differences of the light neutrinos, the entries of the PMNS matrix and two additional real parameters. Of course, this works only because not all of the 6 complex parameters in ∆ j are physically independent.
Using the approximation to the contributions of the 1-loop corrections to the neutrino mass matrix, eq. (2.22), we can relate the calculated neutrino masses to the measured neutrino mass differences.
Following [7] we treat only the effective 3 × 3 light neutrino mass matrix M ν , which is a rank 1 matrix at tree level and equals Similarly to the treatment in [33], we can write the diagonalization of the tree-level neutrino mass matrix as The equation for V s , and corresponds to the Dirac mass term of the effective 2 × 2 seesaw between ζ tree s and ζ tree 4 . Using the notation of M D , eq. (2.6), we can express the Yukawa coupling ∆ 1 as We would like to write ∆ 2 in terms of the vectors V i as well. We can assume 2 that the massless neutrino state ζ o does not couple to the second Higgs doublet, either: This condition ensures that the lightest neutrino only couples to the electroweak sector. Then we can express ∆ 2 in terms of the parameters d and d and the vectors V i as where we choose the phase of V r in such a way, that the coefficient d becomes real and positive. The coefficient d may be a complex number. Our goal is to express these coefficients d and d in terms of the other model parameters.
The neutrino mass matrix, corrected for 1-loop contributions written in eq. (2.22), gives an effective 3 × 3-matrix which has to be diagonalized like eq. (2.13). This diagonalization gives a vanishing neutrino mass m o = 0 and two positive masses m s and m r , which can provide the two measured neutrino mass squared differences. Note, that m s can differ from the tree-level value m tree s obtained from the diagonalization of eq. (3.1). The tree-level diagonalization matrix V partially diagonalizes the effective light neutrino mass matrix M ν , eq. (3.10), and we see explicitly, that it is rank 2: Since there are two Yukawa couplings that couple the fermionic singlet νR to the three generations of neutral leptons, they can be viewed as two 3-vectors in generation space. But two 3-vectors always have a single 3-vector that is orthogonal to both of them. This orthogonal state corresponds to the massless neutrino state ζo = ζ tree o and is therefore the justification of our assumption. It is our choice to consider specific intermediate neutrino states: (1) the state ζ tree s is aligned to one Yukawa coupling, eq. (3.7), and (2) one state, ζo, is orthogonal to the other two states. The former is affected by the mixing due to R3, eq. (3.21), and the later is not. where and f 3 is defined to contain the tree-level contribution, too: The values of a, b, c,f 3 , and f i (i = 1, 2, 3) are complex in the general case, as can be seen from the complex entries in the vectors b k , eq. (2.3). If the Higgs potential is CPconserving, the entries in the vectors b k , table 1, become either real or purely imaginary, hence giving real functions f 1 and f 3 . For getting the masses and the mass eigenstates, we use the Takagi Factorization [25] for eq. (3.11) with the unitary matrix R 3 R 3 only mixes the massive states ζ r and ζ s , hence we can parameterize it as where the parameters β and γ describe effectively only a 2 × 2 unitary matrix. The phases α i have to be determined together with the possible Majorana phases of the light neutrinos. β and γ can be determined from the linear relation and The masses are most easily obtained as the eigenvalues of the squared matrix The masses then are given by The phases α r and α s have to be extracted from the relation linear in M 2×2 , eq. (3.20), and as they drop out in the squared relations. One additional relation for the phases can be obtained from the determinant which can serve as a numerical consistency condition for the extraction of the phases from eqs. (3.27) and (3.28).
With the rotation matrix R 3 , eq. (3.21), we have now the transformation matrix between the flavour eigenstates ν L and the light neutrino mass eigenstates ζ which allows us to identify our vectors V i with columns of the PMNS matrix, eq. (2.7).
Since we chose to identify ζ o with the massless neutrino, ζ r with the neutrino, that gets its mass only with radiative corrections, and ζ s with the neutrino that already has a mass from the seesaw mechanism, we have to take the corresponding columns from the PMNS matrix to determine our vectors, that we want to use for the definition of the Yukawa couplings: where the numbers for the columns have to be taken according to table 2. Relating the measured mass squared differences ∆m 2 21 and ∆m 2 31 to the masses of the three light neutrinos m i we can express the parameters of the Yukawa couplings of the neutrinos to the second Higgs doublet by measured quantities.
Inserting the definitions of the matrix elements a, b, and c (eqs. (3.12)-(3.14)) into the relation eq. (3.29), we can derive: Taking the modulus we get the functional expression for d 2 where we treat m 2 D as a free parameter, since in general m s = m tree s . To get an expression for the modulus of d we take the trace of [eq. (3.20)] · [eq. (3.20)] † , which gives m 2 r + m 2 s on the r.h.s. and (|a| 2 + |b| 2 ) + (|b| 2 + |c| 2 ) on the l.h.s. By reversing the sides, we write a fourth order polynomial in |d |: The general expressions for the coefficients a i are simpler in our CP conserving case with the b-vectors having the form of eq. (2.3). In this case, the values of f 1 and f 3 , given in eqs. (3.15) and (3.19), are real numbers, leading to The value of |d | is then given as a real positive solution to the fourth order equation Therefore |d | has the dependence In order to find a real and positive solution, the value of the phase φ = arg(d ) can be restricted.
Using d and |d | we have analytically replaced two of our input parameters with the neutrino masses. The replacement of the Yukawa couplings to the fermionic singlet is done by eqs. (3.7) and (3.9), using the determined parameters d and |d |, and the vectors V i , eqs. (3.31), (3.32), and (3.33).

The choice of the initial parameters for the numerical analysis
The primary choice are the parameters appearing in the Lagrangian, which define the model. More convenient is a choice, where some of the parameters can be directly related to measured quantities. At the tree-level we achieve this simplification, by using part of the guidelines in [31,32]. Taking for the 2HDM potential the masses of the Higgs particles and their mixing angles and ignoring the 2HDM parameters that do not enter our calculations we get the list: This is the general and basic parameter list for scans of the parameter space of the model. When using our analytical result for the neutrino masses, we can reduce this parameter list by replacing the Yukawa couplings ∆ j by their values, eqs. (3.7) and (3.9). This means, we have to take the neutrino masses as input, also replacing M R with m 4 , as the seesaw mechanism, eq. (2.13), gives the relation even if we do not expect to measure the mass of the heavy neutrino. Since for simplicity we assumed a CP-conserving Higgs potential, we can also simplify the mixing angles of the neutral Higgs bosons, either s 12 or s 13 , as given in table 1, to a single angle s β−α . That leaves us with the same parameter list as the dependencies of |d |, eq. (3.44). From these the only parameter, that does not have an immediate physical meaning is m 2 D . We know from the tree-level seesaw relation eq. (3.6) that but that does not tell us the value of m tree s . Assuming that our model has a sensible loop expansion, we can make the educated guess, that m tree s should be of the same order as the physical mass m s , which we identify with one of the light neutrino masses. For simplicity we parameterize the change from m tree s to m s as a multiplicative parameter that we call λ D , as it enters at the place of m 2 D . Since in our model the lightest neutrino stays massless, the measured neutrino squared mass differences [17], Using the assignments of Table 2 we connect the value of the parameter m D to the masses obtained in eqs. (3.50) and (3.51): We assume that the one loop corrections do not invalidate the tree level assumptions for the seesaw. This allows us to restrict the scaling parameter to the range 1 2 ≤ λ D ≤ 2. Having made these adjustments to the parameters, we arrive at three separate sets of parameters for our analysis. (1) There are several input parameters that (a) are not affected by our calculations, like the Standard Model parameters p i,SM , (b) the parameters of the 2HDM, that do not enter in the calculation of the neutrino masses, like the Higgs potential parameters λ i that do not enter the tree-level Higgs masses, and (c) the Yukawa coupling of the second Higgs doublet to the charged fermions. We summarize the first set of parameters with the namẽ (3.56) Then (2) there are the parameters that are always used as input for the calculation of the neutrino masses, where we use the same symbol s ϑ for both angles s ϑ 1j , table 1, as we have only one nonvanishing mixing angle due to our simplification of taking only a CP-conserving Higgs sector. Comparing to [23] we have c ϑ = s β−α . It is easy to generalize our calculation to a CP non-conserving Higgs potential. We do not expect additional difficulties. In principle, just the intermediate parameter f 1 , defined in eq. (3.15), will become complex. The biggest difficulty would be to present our results in the extended parameter space, while the conclusions of our study would not change.
And (3) there are parameters, that can be both input and output of our calculations. For example, the neutrino parameters are an input, if we use the procedure of this section. But they become an output, if we stay with the Lagrangian parameters {∆ j } as input, as is the starting point of Section 2.

Numerical analysis
Usually, the model parameters are the quantities that are defined in the Lagrangian, and the predictions are the measureable quantities, like masses and cross sections. Therefore, the Yukawa couplings and the parameters of the Higgs potential should be treated as our input parameters. Since our interest in the Higgs sector is limited, we take the masses and the mixing angle of the neutral Higgs bosons as input parameters.
Considerations are different in the neutrino sector. On one hand, using the approximations of Grimus and Lavoura [7] we treat the Yukawa couplings ∆ j , eq. (2.4), together with the Majorana mass M R , the Higgs masses and the Higgs mixing angle as input parameters and "predict" the neutrino masses and mixings. Analysing the model in this way, one can fit the input parameters to obtain the physically measured neutrino mass differences and the neutrino mixing matrix. For this approach one has to construct a minimization function that allows to find the global minimum, which should give the model parameters that correspond to the physically measured values.
On the other hand we can use our analytic results for the neutrino masses to directly determine the Yukawa couplings ∆ j , eqs. (3.7) and (3.9), from the measured neutrino parameters and other input parameters via evaluation of the orthonormal vectors V i . We determine d, d , and R 3 from eqs.  [17] of the oscillation angles θ 12 , θ 13 , θ 23 , and the Dirac phase δ CP as input for the PMNS matrix. Using thus obtained values of the Yukawa couplings ∆ j , we can again go back and "predict" the neutrino masses and mixings, compare the result with the measured masses, and hopefully save a lot of time by having to sample over a much smaller parameter space: we have to vary only one phase φ and the scaling parameter λ D , eq. (3.48), instead of 6 complex entries in ∆ j , eq. (2.4). This is the procedure we adopt in first subsection, 4.1, to check the consistency of our approach.
The second subsection discusses the allowed parameter space by showing various distributions of parameters and interpreting the restrictions that can be seen in the plots. In the third subsection we argue that our analytical approach has advantages over the "blind" systematic scanning of the allowed parameters that go beyond the simple saving of computer time.
All the numerical analysis was performed using data points of the Higgs sector that were subjected to additional theoretical and experimental constraints similar to [34,35], as described in appendix C: the CP-conserving 2HDM potential should be stable, guarantee tree-level unitarity of the S matrix, be bounded from below, have a global minimum, and fulfill the experimental restrictions of the Peskin-Takeuchi S, T , and U parameters. Additionally, the SM Higgs boson has the mass m h = 125.18 GeV [1] and the masses m H and m A of the other two neutral Higgs bosons vary in the range from m h to 3000 GeV. The mixing angle between h 0 and H 0 varies in the range from −π/2 to π/2, where we assume that h 0 to corresponds to the SM Higgs boson. Progress in the experimental particle physics program at the LHC limits the Higgs sector parameter space. Haller et al [36] summarized the restrictions on the 2HDM parameters. We checked that most of the Higgs potential points that pass our theoretical restrictions will also fulfill the more restrictive and specialized constraints of the "typed" 2HDM models.
The behavior of some numerical solutions in our model is illustrated using the benchmark point B1 of [37]. Originally, this point was used for the 2HDM type-II studies [38] (there it was named H-1) and was recently excluded [36]. However, the values would still be valid for the 2HDM type-I model [36,39]. Since we make no distinction for the type of the 2HDM model, we use this point having updated the mass of the lightest Higgs boson h 0 . More details are provided below.

Numerical consistency of the model
As the first test of our approach we calculate the Yukawa couplings, eqs. (3.7) and (3.9), using the input parameters eq. (3.57). For that we have to calculate also the matrix R 3 , eq. (3.21), in order to use the correct orthonormal basis ( V o , V r , V s ) that defines the Yukawa couplings. These numerical values of the Yukawa couplings we treat as input in the sense of eq. (3.45) and calculate the masses and the mixing matrix between the neutrino mass eigenstates and the interaction states: as expected, the mass differences agree between input and output. For the neutrino mixing angles it makes a difference, whether we adopt our complicated procedure, described in sec. 3, or we just take the measured PMNS matrix as the orthonormal basis and ignore the difficulty of calculating R 3 . We get the neutrino mixing angles back in the first case, whereas in the second case the range of φ , that allows solutions to eq. (3.43), is reduced to few points where the angles of the obtained mixing matrix lie in the 3σ bands of the experimentally allowed values. In the second case it can even happen, that we cannot find any values of φ that allow suitable angles of the calculated PMNS matrix.

Distributions of the model parameters
As a first overview we show the distribution of masses of the heavier scalar and pseudoscalar Higgs bosons and the cosine of the mixing angle β − α between the two CP-even Higgs bosons h 0 and H 0 in figure 1. The allowed Higgs potential points are calculated with the procedure described in appendix C. The density of points in (m H , m A ) plane is equalized in the non-logarithmic scale to have a more uniform representation of different m H and m A combinations. Figure 1 clearly shows the restriction on the mixing angle β − α when the masses become large, indicating the onset of the decoupling regime: |β − α| → π/2, when m H , m A 700 GeV. This agrees with the experimental constraints from the LHC measurements [36] suggesting cos(β − α) 0.4. Nearly all considered points satisfy this limit, as indicated by the relative frequency distribution of cos(β − α).  Their order is not related to their magnitude or the numeric nature (whether the value is real or complex). The left panel of figure 3 shows the solutions for a large variation in λ D . When λ D = 0.5 or λ D = 1, we find one physical solution |d | > 0, one non-physical solution |d | < 0, and two complex |d | solutions for every phase value φ . But we see also at certain phase values φ = π ± π/2, that the index numbers of the solutions switch: at φ = π/2 the 1 st and 2 nd solutions become complex and the earlier complex 3 rd and 4 th solutions become real. When λ D = 1.5 or λ D = 2, we find two pairs of real solutions, but only in a very limited range of the phase φ . These two physically allowed values of |d |, the 3 rd and 4 th solutions, will give different Yukawa couplings, for the same point in the Higgs sector. The plot also contains colored dots that are later used to illustrate the value of the Higgs potential parameters in the Higgs basis 0.264299 0.082535 3.67689 -1.75837 -1.72924 0.0842752 0.0699585 Table 3. Benchmark point B1 of [37]. The value of the lightest Higgs boson h 0 is updated to the newest PDG value [1]. The vacuum expectation value v 2 = G −1 F / √ 2, needed to calculate the potential parameters in the generic or in the Higgs basis, is defined in the same way as in [37], but the value for G F = 1.1663787(6) × 10 −5 GeV −2 is taken from [1]. The bilinear parameters m 2 jk or Y i are given in GeV 2 .      We continue our discussion with the statistical description of the whole parameter space of our model. In figure 6 we show the distribution of the size of individual components of the Yukawa couplings for λ D = 1 fixed and m 4 varying between 10 2 to 10 12 on a logarithmic scale. The components of the first Yukawa coupling ∆ 1 show only little variation and a linear dependence in the logarithmic plot on m 4 . But the second Yukawa coupling ∆ 2 exhibits a much larger variation and also a structure at low values of m 4 . This comes from the definition of the Yukawa couplings, eq. (3.7) and (3.9). Whereas ∆ 2 carries the whole variation of d and d , ∆ 1 only sees the dependence of R 3 , eq. (3.21), and of m D , eq. (3.48).   |∆ jk |, seen as the dips in figure 7.
The striped regions in figure 7 depict the values of |∆ jk | for the inverted hierarchy. One can notice the similarity of the vertical thickness of (a) the striped areas for |∆ j2 | and |∆ j3 | with the colored area of |∆ j1 |, and (b) the striped area of |∆ j1 | with the colored areas of |∆ j2 | and |∆ j3 |. The behavior reflects the exchange of the related neutrino states: in the inverted hierarchy the two heavier states are more similar whereas in the normal hierarchy the two lighter states are closer related. In some way |∆ 11 | represents the decoupled state for the normal hierarchy and by that the vector of the PMNS matrix that stands for the massless neutrino. |∆ 12 | and |∆ 13 | give the states mixed from the seesaw mechanism and radiative mass generation. In the inverted hierarchy it is |∆ 13 | that represents the massless neutrino and |∆ 11 | and |∆ 12 | that give the mixed states. This behavior is not so pronounced in ∆ 2k , since this Yukawa coupling is the superposition of the PMNS vectors with the complex numbers d and d , giving a much larger spread of values, as could already be seen in figure 6. Figure 8 shows the wave-like behaviour of the median of |∆ 21 | that comes from the interplay between the scale of the Higgs boson masses and the scale of the Majorana mass term. A hint for this interesting behavior is already seen in the bump of |∆ 21 | for low values of m 4 in figure 6. We obtain very similar plots for |∆ 22 | and |∆ 23 |, as both elements of the second Yukawa coupling have a similar dependence on the parameters of the model.

Numerical advantage of the analytic approach
In our analysis we can find observables which satisfy the experimental bounds by scanning over only one parameter, for example the phase φ . But the usual way for the calculation of observables in such a model (or more sophisticated models) is fitting the parameters by using some global minimization algorithm. Due to the small number of parameters and observables in our study there is a good possibility to compare these two different methods of calculation. In order to find the numerical values for the parameters we construct a minimization function χ 2 where n is the number of observables to be fitted. In this case we fit the neutrino masses and the oscillation parameters θ 12 , θ 13 , θ 23 , and δ CP . H is the Heaviside step function,Ō i denotes the central value of each observable O i , δ ± O i are the upper and lower experimental errors of the observable, and O v i is the calculated value of the observable. The data is fitted by minimizing χ 2 with respect to the Yukawa couplings ∆ 1 and ∆ 2 in eq. (2.4), which means that there are twelve real parameters to be fitted. The central values (i.e. the experimental best fits)Ō i and the 1σ errors δ ± O i are taken from [17].
For the numerical minimization of the χ 2 function we have used the "differential evolution" algorithm which is expensive with respect to computer resources but quite effective. The calculations were carried out for the 1.000 points in the Higgs sector that were used also for Fig. 2, by running forty to sixty separate minimizations on each data point. 3 The parameter sets were saved if χ 2 < 10 −15 , meaning that the calculated value of the observables coincide to a high accuracy with the respective experimental central value. The relative frequency of the value of |∆ 23 | is obtained by counting how often the value appears in the selected bin, divided by the total number of points. We assume the normal hierarchy and take M R ∼ m 4 = 10 10 GeV fixed. The yellow area represents data points obtained with the minimization algorithm, while the black and red histograms represent the distributions of |∆ 23 | calculated with our analytical procedure, described in sec. 3. The black histogram uses λ D = 1 and the red histogram λ D = 1.5. The green histogram shows a part of the red histogram that is determined by taking only the third solution of the fourth order equation (3.43).
In figure 9 we compare the distribution of the Yukawa coupling ∆ 23 calculated using different methods. The histogram in the figure shows the statistical distribution of more than 40.000 points. The yellow area represents data obtained by the minimization algorithm while the black and red histograms represent the distributions of ∆ 23 calculated with the analytical procedure, described in sec. 3, for λ D = 1 or λ D = 1.5, respectively. We see that the black histogram almost coincides with the fitted data, but the red histogram is moved to larger values of |∆ 23 |. As can be seen in Fig. 3, λ D = 1.5 restricts the phase φ to a rather small interval, but gives larger values of |d | at the same time. Since |∆ 23 | depends strongly on |d |, the larger value of |d | explains the shift in the distributions. But the minimization algorithm just looks for any solution and therefore finds the points more probably in larger areas of the parameter space. Since λ D = 1 has a larger phase space than λ D = 1.5, the distribution coming from the minimization algorithm should be more similar to the histogram with λ D = 1, which is what we see in figure 9. For λ D = 1.5 approximately half of the values of ∆ 23 come from the third solution of the fourth order equation (3.43). These values are shown by the green histogram in figure 9. Approximately another half of the values come from the fourth solution and only a very small fraction of values is given by other solutions.
From this example we can guess that using the minimization algorithm in the general case, i.e. by varying twelve free parameters, we could miss some regions of the parameter space, if we do not repeat the minimization often enough. But the main difference between the fitting procedure and calculations using our analytical method is the usage of computational resources. The calculations described in this example took about 430 times longer -25 -using the minimization method than using our analytical method.

Summary
The seesaw mechanism is one of the most successful extensions of the SM which explains neutrino masses. In the usual setup, one adds a heavy singlet fermion for each light neutrino. The Grimus-Neufeld model adds only a single Majorana fermion to the fermion content of the SM, producing only a single seesaw mass for the SM-like neutrinos. Finite corrections to the neutrino mass matrix arise from one-loop diagrams mediated by the heavy neutrino. In the Grimus-Neufeld model these loop corrections produce a radiative mass for one SM-like neutrino. In order to allow this radiative mass the Higgs sector of the Grimus-Neufeld model is constructed from two Higgs doublets, giving two Yukawa couplings to the heavy neutrino. These Yukawa couplings have to be linearly independent, thus characterizing the Higgs sector of the model as a general type. For simplicity we assume a CP-invariant Higgs potential. For the numerical calculations, we take the masses of the neutral Higgs bosons and their mixing angle as input parameters.
We parameterize the Yukawa couplings to the heavy neutrino and calculate the neutrino masses and oscillation parameters following the approximations of Grimus and Lavoura [7]. Since we obtain analytical solutions for the neutrino masses, and the Grimus-Neufeld model has the lightest neutrino massless at one loop level [40], we can use the two measured mass differences as input to determine the Yukawa couplings. With this approach we also retain the neutrino mixing matrix as an unchanged input for our calculation. This change in the parameterization is a new feature compared to previous treatments of seesaw models and has the major advantage, that it reduces the undetermined parameters of the model.
After the distribution of tree-level heavy Higgs masses and the physical tree-level mixing angle between h 0 and H 0 in figure 1 we show the distribution of the parameters d and |d |, that parameterize the second Yukawa coupling, in dependence of the heavy Majorana mass m 4 in figure 2. Taking the benchmark point B1 from [37] as a reference, we show in figure 3 the behavior of the solutions of the fourth order equation that we need to solve to obtain |d | and in figure 4 we show the corresponding Yukawa couplings.
We present a statistical analysis of the modulus of the Yukawa couplings in dependence on m 4 in figure 6 and in dependence on λ D in figure 7. As a final plot 8 in the presentation of the parameter space we also show the wave of the median of |∆ 21 | depending on the mass of the heavy scalar Higgs and discuss its origin, finishing the overview over the parameter space of the Grimus-Neufeld model. The last subsection illustrates with figure 9 the numerical advantage of finding an analytical solution.
In summary, we parameterized and discussed the Grimus-Neufeld model in terms of mostly physically measured low energy scale quantities. The only two "non-physical" parameters that are used in our model are (1) the phase of the Yukawa coupling d of the tree-level "seesaw" neutrino mass state to the second Higgs doublet, denoted as φ , and (2) the proportionality between the tree-level "seesaw" neutrino mass and its mass after the 1-loop radiative correction, denoted as λ D . The other parameters are directly measureable quantities. Having only two not directly measureable parameters increases the testability of our model: a few measurements that restrict the neutrino Yukawa couplings can confirm or rule out our model.
Our study of the Grimus-Neufeld model does not end here. This paper discussed only the Higgs and neutrino sectors. We aim to study the full model with all particle sectors included and get additional restrictions on the Grimus-Neufeld model parameter space from the estimated predictions of rare processes.
Physical Higgs fields φ 0 b k =G 0 must be orthogonal to the Goldstone field G 0 which follows from (A.2). This leads to the condition To study the unit vectors b, introduced in eq. (A.1) (which are the same as eq. (2.2) in the text) and corresponding to the Higgs fields other than the Goldstone boson G 0 , lets define them in the following form: From the orthogonality relations (A.2 -A.4) and due to the fixed value of b G 0 (A.5) it is possible to write the orthogonality equations for the vector components in the following manner: Introducing three sign-parameters s 32im , s 11 , and s 22 (they can take values ±1), we can write where c ij ≡ cos(ϑ ij ) and s ij ≡ sin(ϑ ij ).

B Parameterization of the mixing matrix
Neutrino oscillation angles are introduced using the neutrino mass diagonalization matrix U (2.13) and factorizing it to contain the ordinary Pontecorvo-Maki-Nakagawa-Sakata (PMNS) neutrino mixing matrix [1]. We introduce the formalism by discussing the 3 × 3 neutrino mixing case, where the relationships are simpler; then we expand it to the 4 × 4 case.
The simplest case (3×3) considers only the light neutrinos, assuming they are Majorana particles. This case is discussed in ref. [41] in a slightly different notation of the matrix elements. Factorization of the rotation matrix with the PMNS matrix included explicitly in the case 3 + 3 is discussed in ref. [18]. Here we give formulas for the 3 + 1 case.
The neutrino masses and their mixing angles are predicted from a given neutrino mass matrix (the "top-down" method, as discussed in [41]). Exact analytical expressions for the mixing angles, Dirac and Majorana phases, and formulas for the non-physical phases can be given for the 3-and 4-dimensional cases. Only numerical solutions are possible in the case of 2 or 3 additional neutrinos (i.e. 5-or 6-dimensional [18] cases).

The 3-dimensional case
First we parameterize the neutrino diagonalisation matrix by including explicitly the PMNS mixing matrix for the 3 × 3 mixing [41]. The neutrino mass matrix can be diagonalised by a unitary transformation U , obtained by the singular value decomposition method, see eq. (2.13). Lets denote the complex matrix elements in the following way: This matrix could be factorized into three terms where V PMNS is the standard PMNS mixing matrix [1] for Dirac neutrinos: There are 9 parameters: 3 mixing angles (θ 12 , θ 13 , θ 23 ); 1 Dirac phase δ 13 = δ CP ; 2 Majorana phases α 21 and α 31 ; and the matrixÛ (3) φ containing 3 non-physical and unmeasurable phases φ i (i = 1, 2, 3).
Comparing eqs. (B.1) and (B.2) we can find the relations between the elements of the rotation matrix in a general form and its parameters in the factorized form: These relations are obtained by comparing eq. (B.2) with the corresponding matrix elements from eq. (B.1) forming the upper-triangular matrix: x 1 , x 2 , x 3 , y 2 , y 3 , and z 3 . Other (numerically identical) solutions are possible, using the diagonal elements and the matrix elements from the lower-triangular matrix (y 1 , z 1 , and z 2 ). It should be noted that the Dirac phase can be evaluated using the Jarlskog invariant, for example expressed in the "standard" parameterization [1] J CP = Im (y 3 x * 3 x 2 y * 2 ) = 1 8 cos θ 13 sin 2θ 12 sin 2θ 23 sin 2θ 13 sin δ 13 . (B.11) However, using this equation we need to be careful because J CP has the same value for sin(δ 13 ) and sin(π − δ 13 ), which gives a degeneracy of the δ 13 values.

4-dimensional case
If there is one additional Majorana neutrino, decomposition of the neutrino mass diagonalization matrix into factors including the PMNS neutrino mixing matrix is more complicated. Lets define the 2-dimensional rotation matrices in the 4-dimensional complex space, similarly to ref. [18], α = diag 1, e iα 21 /2 , e iα 31 /2 , 1 . Note that a shorter notation can be used to define the elements of the rotation matrices: where δ b a equals 1, when a = b, or 0, otherwise. This notation is not restricted to the 4-dimensional case.
For the model with n R = 1 the diagonalization matrix U (2.13) is calculated numerically. Defining its elements as x 1 x 2 x 3 x 4 y 1 y 2 y 3 y 4 z 1 z 2 z 3 z 4 t 1 t 2 t 3 t 4      (B. 16) and comparing to eq. (B.14) we find the relations: As U (4×4) is unitary, there are relations between the elements. The expressions for the angles do not contain all entries of the rotation matrix U (4×4) , defined in eq. (B.16). The relations used are obtained comparing eq. (B.14) with the matrix elements from eq. (B.16) forming the upper-triangular matrix: x 1 , x 2 , x 3 , x 4 , y 2 , y 3 , y 4 , z 3 , z 4 , and t 4 . Other (numerically identical) solutions are possible using the diagonal elements x 1 , y 2 , z 3 , and t 4 , and the matrix elements z 1 , z 2 , t 1 , t 2 , and t 3 .

C The two-Higgs-doublet model
The most general 2HDM scalar potential of two doublets φ 1 and φ 2 is where the parameters m 2 11 , m 2 22 , and λ 1−4 are real numbers, whereas the remaining parameters λ 5 , λ 6 , λ 7 and m 2 12 in general can be complex. Since our main purpose of the paper is the analysis of the neutrino sector we restrict our analysis of the Higgs sector a CP-conserving Higgs potential with a softly broken Z 2 symmetry, where λ 5 and m 2 12 are real, but λ 6 = λ 7 = 0.
To ensure a stable vacuum, the scalar potential has to be bound from below (BFB), i.e. there should be no direction in the field space along which the potential tends to minus infinity. Necessary and sufficient conditions for the most general 2HDM scalar potential to be BFB were first derived in ref. [43] and later in ref. [44]. The procedure of ref. [44] can only be handled numerically. For 2HDM scalar potentials, that are more constrained by symmetries, like in potentials where one has λ 6 = λ 7 = 0, the necessary and sufficient BFB conditions can be derived 4 as simple analytical expressions: (C.6) We also apply the condition from ref. [44], which guarantees that the vacuum state has a lower value than all the other possible stability points of the potential: