Complex neutrino Yukawa

[FStaub]

Both is necessary.

[SDrau]

Thanks, it's nice to have a confirmation!

However, I now return to problems of having complex neutrino Yukawa couplings...I am trying to implement a model with two Higgs doublets and heavy Majorana neutrino. I have attached the model files.

I have edited the src/Makefile, recompiled SPheno and also set the 5th flag to 2 in the MODSEL block, yet I am having problems with appearing NaNs and negative masses squared. I did try to run the BLSM with the parameter point from this topic with QP enabled and it worked out fine.

The errors that I get if the Yukawa coupling is fully imaginary are:
with QP

Code: Select all

NaN appearing in CalculateMFv                            
				Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
				STOP Subroutine TerminateProgram 	

and without QP

Code: Select all

Calculating mass spectrum
							  1 .-iteration
				   ... reached precision:                  Infinity
							  2 .-iteration
				   ... reached precision:   5.1981170444690511E-006
				 Calculate loop corrected masses 
				 Warning from Subroutine OneLoopMFv                              
				 a mass squarred is negative:            1  -1.0623512217328707E-025
				Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_OVERFLOW_FLAG IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
				STOP Subroutine TerminateProgram

But I do get normal output if I use a purely real coupling with the same numerical values. I would also expect to get normal output with the purely imaginary Yukawa since that is what I get when I use FlexibleSUSY.

I have attached the model files for SARAH (Flexible312 and SPheno312 folders), SPheno (in SPheno312) and FlexibleSUSY( in Flexible312/Flexible312/), LesHouches.in.<model> files for FlexibleSUSY and SPheno, and a little file, in which I have specified the used Yukawa coupling as well as the masses of neutrinos (or errors) that I get out in all cases. SARAH model files differ only by the definition of VEVs, because FlexibleSUSY is able to handle the explicit 0 in one of the doublets and in SPheno I set one of the VEVs to 0 when solving the tadpole equations.

I hope that I have provided all the needed information and thank you in advance!

Sincerely,

S.Drauksas

[FStaub]

Hi,

the problem was the following: once CP violation is turned on, the Neutrino matrix is no longer symmetric. Therefore, one needs to diagonalise M^\dagger M to get the eigenvalues squared. However, some neutrinos are massless in your matrix what returns tiny random numbers O(10^-30) for the mass squared with both signs. When taking the sqrt, that leads to the NaN. I attach a patch for this where I set a cut-off of |10^-25| for m^2. Everything below is a numerical zero which is set to +10^-25.
I'll send a copy to Werner Porod to include that in SPheno.

A second place where tiny negative masses give a problem is this line in TreeMasses_MODEL generated by SARAH.

Code: Select all

If (MaxVal(Abs(mat2(i1,(i1+1):4))).gt.Abs(mat2(i1,i1))) SecondDiagonalisationNeeded = .True. 

You should replace that by

Code: Select all

If (MaxVal(Abs(mat2(i1,(i1+1):4))).gt.Max(Abs(mat2(i1,i1)),1.0E-15_dp)) SecondDiagonalisationNeeded = .True. 

[SDrau]

Dear Florian,

thank you for the fixes and looking into it! Though, I am still bumping into problems.
I did replace the MathematicsQP.f90 file, the lines in TreeMasses_MODEL and recompiled everything. I no longer get NaNs, but with purely imaginary Yukawa I get this error:

Code: Select all

 Calculating mass spectrum
              1 .-iteration
  Problem in RGE Running. Errorcode:       -1006
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
STOP Subroutine TerminateProgram

and in Messages.out:

Code: Select all

 sugra            1 .-iteration
 Warning from Subroutine CalculateMFv, ierr =       -1006
 Warning from Subroutine CalculateMFv, ierr =       -1006
   
 ErrorLevel, Iname:           0           0
 The error has occured in the following chain of subroutines:
 Match_and_Run                           
 OneLoopMasses                           
 TreeMassesSPheno312                     
 CalculateMFv                            
 RealEigenSystem_DP                      
 RealEigenSystem  

With purely real Yukawa coupling I get the same output as before.

Also, I have tried another coupling. Real coupling seems to work fine, but the imaginary coupling gave out Sqrt[10^-25] value as masses of the neutrinos.
I did use a considerably smaller coupling:

Code: Select all

Block IMYV1IN #
1   5.000000E-10         # Yv1(1)
2   4.000000E-08         # Yv1(2)
3   2.000000E-07         # Yv1(3)

With a larger coupling than the first one that also gave output when given as the real part:

Code: Select all

Block IMYV1IN #
1   5.000000E-3         # Yv1(1)
2   4.000000E-2         # Yv1(2)
3   2.000000E-1         # Yv1(3)

I got a NaN:

Code: Select all

 Calculating mass spectrum
              1 .-iteration
 NaN appearing in CalculateMFv                            
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO IEEE_UNDERFLOW_FLAG IEEE_DENORMAL
STOP Subroutine TerminateProgram

and in Messages.out:

Code: Select all

 Warning from Subroutine CalculateMFv, ierr =       -1006
 sugra            1 .-iteration
   
 ErrorLevel, Iname:           0           4
 The error has occured in the following chain of subroutines:
 Match_and_Run                           
 BoundaryBSM                             
 TreeMassesSPheno312                     
 CalculateMFv                            
 ComplexEigenSystem_DP                   
 RealEigenSystem 

[FStaub]

Hi,

it seems that the thresholds were too optimistic.
You can change in MatheamticsQP.f90 the check

Code: Select all

   Do i1=1,n1
    If (Abs(wr(2*i1-1)).gt.1.0E-25_dp) Then
      eigenvalues(i1) = Sqrt( wr(2*i1-1) )
    Else   
      eigenvalues(i1) = 1.0E-25_dp
    End If

and set other values. Maybe, one can even think of calculating the threshold dynamical:
In Quadruple precision you have about 30 digits which come out correctly. This, anything what is below "10^-30 x (biggest entry in mass matrix)^2" is a random number. So, for MR ~ 10^4, a proper threshold is more about 10^-22.

Similarly,

Code: Select all

If (MaxVal(Abs(mat2(i1,(i1+1):4))).gt.Max(Abs(mat2(i1,i1)),1.0E-15_dp)) SecondDiagonalisationNeeded = .True. 

also that should be adjusted. The value which I put in the check is the minimal value which is not considered to be zero in the off-diagonal entries in the diagonalised mass matrix (without squares).

Unfortunately, dealing with neutrino masses in the complex case is numerical quite delicate as you can see. I just wonder what FlexibleSUSY does use for routines for that?

Cheers,
Florian

[SDrau]

Dear Florian,

so in essence, the minimal mass of neutrinos in the complex case is limited by 10^(-30)*MR^2 instead of 10^(-30)*MR?

I took a little peek at the code in FlexibleSUSY. As far as I can tell, it does not try to deal with asymmetric matrices. Specifically, only the upper triangle of the Mv mass matrix is written down and then it is copied to the lower one, so FS is able to deal with complex entries with the same precision as in the real case. Though, FlexibleSUSY is using double precision and probably all of the numbers smaller than ~10^-16 are random as well.

Regards,

Simonas

[FStaub]

Hi,

thanks for sharing your observations.
I'm just surprised about the FS approach: hey should get complex eigenvalues since the matrix is no longer hermitian, or?

Cheers
Florian

[SDrau]

Dear Florian,

I did some more digging and in their code documentation they reference this paper: https://doi.org/10.1016/j.amc.2014.01.170 , which I could not access for free. I think that by using Takagi's factorization they are able to bypass actually having M^{Dagger}M, but they do get square roots of eigenvalues of M^{Dagger}M.

Regards,

Simonas

[AlexanderVoigt]

Hi Simonas and Florian,

you are correct, for symmetric matrices FlexibleSUSY implements Takagi's factorization, which matters in case some singular values are degenerate. In this approach M^{Dagger}M is not diagonalized.

As you say, for your Flexible312 model FlexibleSUSY symmetrizes the mass matrix, i.e. only the upper triangle is calculated and then copied to the lower one. As far as I can see, the mass matrix must be symmetric, otherwise you'd have a physical problem. So, if I've not overlooked something, the result of FlexibleSUSY should in principle be correct, except for the numerical instabilities you've discussed, because only double precision is used.

Best regards,
Alex