PMNS_Sterile.cxx
//......................................................................
///
/// Build and solve the full Hamiltonian for eigenvectors and eigenvalues.
///
/// This uses GSL to solve the eigensystem.
///
/// Here we divide the mass squared matrix Hms by the 2E
/// to obtain the vacuum Hamiltonian in eV. Then, the matter
/// potential is added to the electron and sterile components.
///
/// The matter potential is implemented as in arXiv:hep-ph/0606054v3 (eq. 3.15).
///
Build and solve the full Hamiltonian for eigenvectors and eigenvalues.

This uses GSL to solve the eigensystem.

Here we divide the mass squared matrix Hms by the 2E
to obtain the vacuum Hamiltonian in eV. Then, the matter
potential is added to the electron and sterile components.

The matter potential is implemented as in arXiv:hep-ph/0606054v3 (eq.
Build and solve the full Hamiltonian for eigenvectors and eigenvalues . This uses GSL to solve the eigensystem . Here we divide the mass squared matrix Hms by the 2E to obtain the vacuum Hamiltonian in eV . Then the matter potential is added to the electron and sterile components . The matter potential is implemented as in arXiv : hep - ph / 0606054v3 ( eq .
void PMNS_Sterile::SolveHam()
{
  BuildHms();
  if(fGotES) return;
  if(TryCache()) return;
  double rho = fPath.density;
  double zoa = fPath.zoa;
  double lv = 2 * kGeV2eV*fEnergy;                 
  double kr2GNe = kK2*M_SQRT2*kGf * rho*zoa;       
  double kr2GNn = kK2*M_SQRT2*kGf * rho*(1-zoa)/2; 
  for(size_t i=0;i<size_t(fNumNus);i++){
    complexD buf = fHms[i][i]/lv;
    *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_rect(real(buf),0);
    for(size_t j=i+1;j<size_t(fNumNus);j++){
      buf = fHms[i][j]/lv;
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
    }
    if(i>2){
      if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) ,  kr2GNn);
      else          *gsl_matrix_complex_ptr(fGSL.H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,i,i) , -kr2GNn);;
    }
  }
  if(!fIsNuBar) *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) ,  kr2GNe);
  else          *gsl_matrix_complex_ptr(fGSL.H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(fGSL.H_GSL,0,0) , -kr2GNe);
  gsl_eigen_hermv(fGSL.H_GSL, fGSL.fEvalGSL, fGSL.fEvecGSL, fGSL.W_GSL);
  for(int i=0;i<fNumNus;i++){
    fEval[i] = gsl_vector_get(fGSL.fEvalGSL,i);
    for(int j=0;j<fNumNus;j++){
      gsl_complex buf = gsl_matrix_complex_get(fGSL.fEvecGSL,i,j);
      fEvec[i][j] = complexD( GSL_REAL(buf), GSL_IMAG(buf) );
    }
  }
  fGotES = true;
  FillCache();
}
void PMNS_Sterile :: SolveHam ( ) { BuildHms ( ) ; if ( fGotES ) return ; if ( TryCache ( ) ) return ; double rho = fPath . density ; double zoa = fPath . zoa ; double lv = 2 * kGeV2eV * fEnergy ; double kr2GNe = kK2 * M_SQRT2 * kGf * rho * zoa ; double kr2GNn = kK2 * M_SQRT2 * kGf * rho * ( 1 - zoa ) / 2 ; for ( size_t i = 0 ; i < size_t ( fNumNus ) ; i ++ ) { complexD buf = fHms [ i ] [ i ] / lv ; * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_rect ( real ( buf ) , 0 ) ; for ( size_t j = i + 1 ; j < size_t ( fNumNus ) ; j ++ ) { buf = fHms [ i ] [ j ] / lv ; if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , j , i ) = gsl_complex_rect ( real ( buf ) , - imag ( buf ) ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , j , i ) = gsl_complex_rect ( real ( buf ) , imag ( buf ) ) ; } if ( i > 2 ) { if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , i , i ) , kr2GNn ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , i , i ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , i , i ) , - kr2GNn ) ; ; } } if ( ! fIsNuBar ) * gsl_matrix_complex_ptr ( fGSL . H_GSL , 0 , 0 ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , 0 , 0 ) , kr2GNe ) ; else * gsl_matrix_complex_ptr ( fGSL . H_GSL , 0 , 0 ) = gsl_complex_add_real ( gsl_matrix_complex_get ( fGSL . H_GSL , 0 , 0 ) , - kr2GNe ) ; gsl_eigen_hermv ( fGSL . H_GSL , fGSL . fEvalGSL , fGSL . fEvecGSL , fGSL . W_GSL ) ; for ( int i = 0 ; i < fNumNus ; i ++ ) { fEval [ i ] = gsl_vector_get ( fGSL . fEvalGSL , i ) ; for ( int j = 0 ; j < fNumNus ; j ++ ) { gsl_complex buf = gsl_matrix_complex_get ( fGSL . fEvecGSL , i , j ) ; fEvec [ i ] [ j ] = complexD ( GSL_REAL ( buf ) , GSL_IMAG ( buf ) ) ; } } fGotES = true ; FillCache ( ) ; }
......................................................................
......................................................................
// Build Hamiltonian // Check if anything changed // Try caching if activated   // 2E in eV  // Electron matter potential in eV // Neutron matter potential in eV // Finish building Hamiltonian in matter with dimension of eV // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example.  // Add nue CC coherent forward scattering from sterile neutrinos.  // Solve Hamiltonian for eigensystem // Fill fEval and fEvec vectors from GSL objects // Fill cache if activated  