PMNS_Sterile.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////
2 // $Id: PMNS_Sterile.cxx,v 1.7 2014/04/17 21:34:38 jcoelho Exp $
3 //
4 // Implementation of oscillations of neutrinos in matter in a
5 // n-neutrino framework.
6 //
7 //......................................................................
8 //
9 // Throughout I have taken:
10 // - L to be the neutrino flight distance in km
11 // - E to be the neutrino energy in GeV
12 // - Dmij to be the differences between the mass-squares in eV^2
13 // - Ne to be the electron number density in mole/cm^3
14 // - theta_ij and delta_ij to be in radians
15 //
16 // joao.coelho@tufts.edu
17 ////////////////////////////////////////////////////////////////////////
18 
19 #include "OscLib/PMNS_Sterile.h"
20 
21 #include <iostream>
22 #include <cassert>
23 #include <stdlib.h>
24 
25 #include <gsl/gsl_complex.h>
26 #include <gsl/gsl_complex_math.h>
27 #include <gsl/gsl_matrix.h>
28 #include <gsl/gsl_eigen.h>
29 #include <gsl/gsl_math.h>
30 #include <gsl/gsl_const_num.h>
31 #include <gsl/gsl_const_mksa.h>
32 #include <gsl/gsl_const_cgsm.h>
33 
34 namespace osc
35 {
37  {
38  Priv(int numNu)
39  {
40  // Allocate memory for the GSL objects
41  fEval = gsl_vector_alloc(numNu);
42  fEvec = gsl_matrix_complex_alloc(numNu, numNu);
43  H_GSL = gsl_matrix_complex_alloc(numNu, numNu);
44  W_GSL = gsl_eigen_hermv_alloc(numNu);
45  }
46 
47  Priv(const Priv& p)
48  {
49  const int numNu = p.fEval->size;
50  fEval = gsl_vector_alloc(numNu);
51  fEvec = gsl_matrix_complex_alloc(numNu, numNu);
52  H_GSL = gsl_matrix_complex_alloc(numNu, numNu);
53  W_GSL = gsl_eigen_hermv_alloc(numNu);
54 
55  gsl_vector_memcpy(fEval, p.fEval);
56  gsl_matrix_complex_memcpy(fEvec, p.fEvec);
57  gsl_matrix_complex_memcpy(H_GSL, p.H_GSL);
58  // No way to copy. Hope leaving it empty is OK
59  }
60 
62  {
63  // Free memory from GSL objects
64  gsl_matrix_complex_free(H_GSL);
65  gsl_eigen_hermv_free(W_GSL);
66  gsl_matrix_complex_free(fEvec);
67  gsl_vector_free(fEval);
68  }
69 
70  gsl_vector* fEval;
71  gsl_matrix_complex* fEvec;
72  gsl_matrix_complex* H_GSL;
73  gsl_eigen_hermv_workspace* W_GSL;
74  };
75 }
76 
77 
78 // Some useful complex numbers
79 static std::complex<double> zero(0.0,0.0);
80 static std::complex<double> one (1.0,0.0);
81 
82 // Unit conversion constants
83 static const double kKm2eV = 5.06773103202e+09; ///< km to eV^-1
84 static const double kK2 = 4.62711492217e-09; ///< mole/GeV^2/cm^3 to eV
85 static const double kGeV2eV = 1.0e+09; ///< GeV to eV
86 
87 //G_F in units of GeV^-2
88 static const double kGf = 1.166371E-5;
89 
90 using namespace std;
91 
92 using namespace osc;
93 
94 //......................................................................
95 
97  : d(new PMNS_Sterile::Priv(NumNus))
98 {
99 
100  fNumNus = NumNus;
101  this->SetStdPars();
102  this->ResetToFlavour(1);
103  fCachedNe = 0.0;
104  fCachedE = 1.0;
105  fCachedAnti = 1;
106 }
107 
109 {
110  delete d;
111 }
112 
114 {
115 
116  fDm = vector<double>(fNumNus, 0);
117  fTheta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
118  fDelta = vector< vector<double> >(fNumNus, vector<double>(fNumNus,0));
119 
120  fNuState = vector<complex>(fNumNus, zero);
121  fHms = vector< vector<complex> >(fNumNus, vector<complex>(fNumNus,zero));
122 
123 }
124 
126 {
127 
128  this->InitializeVectors();
129 
130  if(fNumNus>2){
131  this->SetAngle(1,2,0.6);
132  this->SetAngle(1,3,0.16);
133  this->SetAngle(2,3,0.7);
134  this->SetDm(2,7.6e-5);
135  this->SetDm(3,2.4e-3);
136  }
137  else if(fNumNus==2){
138  this->SetAngle(1,2,0.7);
139  this->SetDm(2,2.4e-3);
140  }
141 
142 }
143 
144 //......................................................................
145 ///
146 /// Set mixing angles in radians. Requires i < j.
147 ///
148 void PMNS_Sterile::SetAngle(int i, int j, double th)
149 {
150 
151  if(i>j){
152  cout << "First argument should be smaller than second argument" << endl;
153  cout << "Setting reverse order (Theta" << j << i << "). " << endl;
154  int temp = i;
155  i = j;
156  j = temp;
157  }
158  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
159  cout << "Theta" << i << j << " not valid for " << fNumNus;
160  cout << " neutrinos. Doing nothing." << endl;
161  return;
162  }
163 
164  fTheta[i-1][j-1] = th;
165 
166  fBuiltHms = false;
167 
168 }
169 
170 //......................................................................
171 ///
172 /// Set CP violating phases in radians. Requires i+1 < j.
173 ///
174 void PMNS_Sterile::SetDelta(int i, int j, double delta)
175 {
176 
177  if(i>j){
178  cout << "First argument should be smaller than second argument" << endl;
179  cout << "Setting reverse order (Delta" << j << i << "). " << endl;
180  int temp = i;
181  i = j;
182  j = temp;
183  }
184  if(i<1 || i>fNumNus-1 || j<2 || j>fNumNus){
185  cout << "Delta" << i << j << " not valid for " << fNumNus;
186  cout << " neutrinos. Doing nothing." << endl;
187  return;
188  }
189  if(i+1==j){
190  cout << "Rotation " << i << j << " is real. Doing nothing." << endl;
191  return;
192  }
193 
194  fDelta[i-1][j-1] = delta;
195 
196  fBuiltHms = false;
197 
198 }
199 
200 //......................................................................
201 ///
202 /// Set the mass-splittings. These are m_j^2-m_1^2
203 ///
204 void PMNS_Sterile::SetDm(int j, double dm)
205 {
206 
207  if(j<2 || j>fNumNus){
208  cout << "Dm" << j << "1 not valid for " << fNumNus;
209  cout << " neutrinos. Doing nothing." << endl;
210  return;
211  }
212 
213  fDm[j-1] = dm;
214 
215  fBuiltHms = false;
216 
217 }
218 
219 //......................................................................
220 ///
221 /// Rotate the Hamiltonian by the angle theta_ij and phase delta_ij.
222 /// The rotations assume all off-diagonal elements with i > j are zero.
223 /// This is correct if the order of rotations is chosen appropriately
224 /// and it speeds up computation by skipping null terms
225 ///
226 void PMNS_Sterile::RotateH(int i,int j){
227 
228  // Do nothing if angle is zero
229  if(fTheta[i][j]==0) return;
230 
231  double fSinBuffer = sin(fTheta[i][j]);
232  double fCosBuffer = cos(fTheta[i][j]);
233 
234  double fHmsBufferD;
235  complex fHmsBufferC;
236 
237  // With Delta
238  if(i+1<j){
239  complex fExpBuffer = complex(cos(fDelta[i][j]), -sin(fDelta[i][j]));
240 
241  // General case
242  if(i>0){
243  // Top columns
244  for(int k=0; k<i; k++){
245  fHmsBufferC = fHms[k][i];
246 
247  fHms[k][i] *= fCosBuffer;
248  fHms[k][i] += fHms[k][j] * fSinBuffer * conj(fExpBuffer);
249 
250  fHms[k][j] *= fCosBuffer;
251  fHms[k][j] -= fHmsBufferC * fSinBuffer * fExpBuffer;
252  }
253 
254  // Middle row and column
255  for(int k=i+1; k<j; k++){
256  fHmsBufferC = fHms[k][j];
257 
258  fHms[k][j] *= fCosBuffer;
259  fHms[k][j] -= conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
260 
261  fHms[i][k] *= fCosBuffer;
262  fHms[i][k] += fSinBuffer * fExpBuffer * conj(fHmsBufferC);
263  }
264 
265  // Nodes ij
266  fHmsBufferC = fHms[i][i];
267  fHmsBufferD = real(fHms[j][j]);
268 
269  fHms[i][i] *= fCosBuffer * fCosBuffer;
270  fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
271  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
272 
273  fHms[j][j] *= fCosBuffer * fCosBuffer;
274  fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
275  fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j] * conj(fExpBuffer));
276 
277  fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j] * conj(fExpBuffer)) * fSinBuffer * fExpBuffer;
278  fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC) * fExpBuffer;
279 
280  }
281  // First rotation on j (No top columns)
282  else{
283  // Middle rows and columns
284  for(int k=i+1; k<j; k++){
285  fHms[k][j] = -conj(fHms[i][k]) * fSinBuffer * fExpBuffer;
286 
287  fHms[i][k] *= fCosBuffer;
288  }
289 
290  // Nodes ij
291  fHmsBufferD = real(fHms[i][i]);
292 
293  fHms[i][j] = fSinBuffer * fCosBuffer * (fHms[j][j] - fHmsBufferD) * fExpBuffer;
294 
295  fHms[i][i] *= fCosBuffer * fCosBuffer;
296  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
297 
298  fHms[j][j] *= fCosBuffer * fCosBuffer;
299  fHms[j][j] += fSinBuffer * fHmsBufferD * fSinBuffer;
300  }
301 
302  }
303  // Without Delta (No middle rows or columns: j = i+1)
304  else{
305  // General case
306  if(i>0){
307  // Top columns
308  for(int k=0; k<i; k++){
309  fHmsBufferC = fHms[k][i];
310 
311  fHms[k][i] *= fCosBuffer;
312  fHms[k][i] += fHms[k][j] * fSinBuffer;
313 
314  fHms[k][j] *= fCosBuffer;
315  fHms[k][j] -= fHmsBufferC * fSinBuffer;
316  }
317 
318  // Nodes ij
319  fHmsBufferC = fHms[i][i];
320  fHmsBufferD = real(fHms[j][j]);
321 
322  fHms[i][i] *= fCosBuffer * fCosBuffer;
323  fHms[i][i] += 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
324  fHms[i][i] += fSinBuffer * fHms[j][j] * fSinBuffer;
325 
326  fHms[j][j] *= fCosBuffer * fCosBuffer;
327  fHms[j][j] += fSinBuffer * fHmsBufferC * fSinBuffer;
328  fHms[j][j] -= 2 * fSinBuffer * fCosBuffer * real(fHms[i][j]);
329 
330  fHms[i][j] -= 2 * fSinBuffer * real(fHms[i][j]) * fSinBuffer;
331  fHms[i][j] += fSinBuffer * fCosBuffer * (fHmsBufferD - fHmsBufferC);
332 
333  }
334  // First rotation (theta12)
335  else{
336 
337  fHms[i][j] = fSinBuffer * fCosBuffer * fHms[j][j];
338 
339  fHms[i][i] = fSinBuffer * fHms[j][j] * fSinBuffer;
340 
341  fHms[j][j] *= fCosBuffer * fCosBuffer;
342 
343  }
344  }
345 
346 }
347 
348 //......................................................................
349 ///
350 /// Build Hms = H*2E, where H is the Hamiltonian in vacuum on flavour basis
351 /// and E is the neutrino energy in eV. Hms is effectively the matrix of
352 /// masses squared.
353 ///
354 /// This is a hermitian matrix, so only the
355 /// upper triangular part needs to be filled
356 ///
357 /// The construction of the Hamiltonian avoids computing terms that
358 /// are simply zero. This has a big impact in the computation time.
359 /// This construction is described in DocDB-XXXX (to be posted)
360 ///
362 {
363 
364  // Check if anything changed
365  if(fBuiltHms) return;
366 
367  for(int j=0; j<fNumNus; j++){
368  // Set mass splitting
369  fHms[j][j] = fDm[j];
370  // Reset off-diagonal elements
371  for(int i=0; i<j; i++){
372  fHms[i][j] = 0;
373  }
374  // Rotate j neutrinos
375  for(int i=0; i<j; i++){
376  this->RotateH(i,j);
377  }
378  }
379 
380  // Tag as built
381  fBuiltHms = true;
382 
383 }
384 
385 //......................................................................
386 ///
387 /// Solve the full Hamiltonian for eigenvectors and eigenvalues
388 /// This uses GSL to solve the eigensystem.
389 ///
390 /// The matter effect assumes # of neutrons = # of electrons.
391 ///
392 void PMNS_Sterile::SolveHam(double E, double Ne, int anti)
393 {
394 
395  // Check if anything has changed before recalculating
396  if(Ne!=fCachedNe || E!=fCachedE || anti!=fCachedAnti || !fBuiltHms ){
397  fCachedNe = Ne;
398  fCachedE = E;
399  fCachedAnti = anti;
400  this->BuildHms();
401  }
402  else return;
403 
404  double lv = 2 * kGeV2eV*E; // 2E in eV
405  double kr2GNe = kK2*M_SQRT2*kGf*Ne; // Matter potential in eV
406 
407  // Finish building Hamiltonian in matter with dimension of eV
408 
409  for(size_t i=0;i<size_t(fNumNus);i++){
410  complex buf = fHms[i][i]/lv;
411  *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_rect(real(buf),0);
412  for(size_t j=i+1;j<size_t(fNumNus);j++){
413  buf = fHms[i][j]/lv;
414  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf),-imag(buf));
415  else *gsl_matrix_complex_ptr(d->H_GSL, j, i) = gsl_complex_rect(real(buf), imag(buf));
416  }
417  if(i>2){
418  // Subtract NC coherent forward scattering from sterile neutrinos. See arXiv:hep-ph/0606054v3, eq. 3.15, for example.
419  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,i,i) , kr2GNe/2);
420  else *gsl_matrix_complex_ptr(d->H_GSL, i, i) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,i,i) , -kr2GNe/2);;
421  }
422  }
423  // Add nue CC coherent forward scattering from sterile neutrinos.
424  if(anti>0) *gsl_matrix_complex_ptr(d->H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,0,0) , kr2GNe);
425  else *gsl_matrix_complex_ptr(d->H_GSL, 0, 0) = gsl_complex_add_real(gsl_matrix_complex_get(d->H_GSL,0,0) , -kr2GNe);
426 
427  // Solve Hamiltonian for eigensystem
428  gsl_eigen_hermv(d->H_GSL, d->fEval, d->fEvec, d->W_GSL);
429 
430 }
431 
432 ///.....................................................................
433 ///
434 /// Propagate the current neutrino state over a distance L in km
435 /// with an energy E in GeV through constant matter of density
436 /// Ne in mole/cm^3.
437 /// @param anti - +1 = neutrino case, -1 = anti-neutrino case
438 ///
439 void PMNS_Sterile::PropMatter(double L, double E, double Ne, int anti)
440 {
441 
442  // Solve Hamiltonian
443  this->SolveHam(E, Ne, anti);
444 
445  // Store coefficients of propagation eigenstates
446  vector<complex> nuComp(fNumNus, zero);
447  for(int i=0;i<fNumNus;i++){
448  nuComp[i] = 0;
449  for(int j=0;j<fNumNus;j++){
450  gsl_complex buf = gsl_matrix_complex_get(d->fEvec,j,i);
451  complex evecji = complex( GSL_REAL(buf), GSL_IMAG(buf) );
452  nuComp[i] += fNuState[j] * conj( evecji );
453  }
454  }
455 
456  // Propagate neutrino state
457  for(int i=0;i<fNumNus;i++){
458  fNuState[i] = 0;
459  for(int j=0;j<fNumNus;j++){
460  gsl_complex buf = gsl_matrix_complex_get(d->fEvec,i,j);
461  complex evecij = complex( GSL_REAL(buf), GSL_IMAG(buf) );
462  complex iEval(0.0,gsl_vector_get(d->fEval,j));
463  fNuState[i] += exp(-iEval * kKm2eV*L) * nuComp[j] * evecij;
464  }
465  }
466 
467 }
468 
469 //......................................................................
470 ///
471 /// Do several layers in a row. L and Ne must have the same length
472 ///
473 void PMNS_Sterile::PropMatter(const std::list<double>& L,
474  double E,
475  const std::list<double>& Ne,
476  int anti)
477 {
478  if (L.size()!=Ne.size()) abort();
479  std::list<double>::const_iterator Li (L.begin());
480  std::list<double>::const_iterator Lend(L.end());
481  std::list<double>::const_iterator Ni (Ne.begin());
482  for (; Li!=Lend; ++Li, ++Ni) {
483  //cout << *Li << " km and " << *Ni << " mol/cm^3" << endl;
484  this->PropMatter(*Li, E, *Ni, anti);
485  }
486 }
487 
488 
489 //......................................................................
490 ///
491 /// Reset the neutrino state back to a pure flavour where
492 /// it starts
493 ///
495 {
496  int i;
497  for (i=0; i<fNumNus; ++i){
498  if (i==flv) fNuState[i] = one;
499  else fNuState[i] = zero;
500  }
501 }
502 
503 //......................................................................
504 ///
505 /// Compute oscillation probability of flavour flv
506 ///
507 /// 0 = nue, 1 = numu, 2 = nutau
508 ///
509 double PMNS_Sterile::P(int flv) const
510 {
511  assert(flv>=0 && flv<fNumNus);
512  return norm(fNuState[flv]);
513 }
514 
515 ////////////////////////////////////////////////////////////////////////
std::vector< double > fDm
m^2_i - m^2_1 in vacuum
Definition: PMNS_Sterile.h:84
virtual void SetDelta(int i, int j, double delta)
virtual void SetAngle(int i, int j, double th)
static const double kK2
mole/GeV^2/cm^3 to eV
Definition: PMNSOpt.h:42
double delta
Definition: runWimpSim.h:98
virtual double P(int flv) const
virtual void PropMatter(double L, double E, double Ne, int anti=1)
const char * p
Definition: xmltok.h:285
virtual void SetDm(int i, double dm)
gsl_matrix_complex * fEvec
double fCachedNe
Cached electron density.
Definition: PMNS_Sterile.h:91
PMNS_Sterile(int NumNus)
virtual void BuildHms()
gsl_matrix_complex * H_GSL
static const double kGeV2eV
GeV to eV.
Definition: PMNSOpt.h:43
virtual void RotateH(int i, int j)
Rotate the Hamiltonian by theta_ij and delta_ij.
static constexpr double L
gsl_eigen_hermv_workspace * W_GSL
Float_t E
Definition: plot.C:20
std::vector< complex > fNuState
The neutrino current state.
Definition: PMNS_Sterile.h:88
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
std::complex< double > complex
Definition: PMNS_Sterile.h:65
virtual void SetStdPars()
Set standard 3-flavor parameters.
const double j
Definition: BetheBloch.cxx:29
virtual void ResetToFlavour(int flv=1)
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
T sin(T number)
Definition: d0nt_math.hpp:132
virtual ~PMNS_Sterile()
int fCachedAnti
Cached nu/nubar selector.
Definition: PMNS_Sterile.h:93
Float_t norm
bool fBuiltHms
Tag to avoid rebuilding Hms.
Definition: PMNS_Sterile.h:94
std::vector< std::vector< double > > fTheta
theta[i][j] mixing angle
Definition: PMNS_Sterile.h:85
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
virtual void SolveHam(double E, double Ne, int anti)
double fCachedE
Cached neutrino energy.
Definition: PMNS_Sterile.h:92
std::vector< std::vector< complex > > fHms
matrix H*2E in eV^2
Definition: PMNS_Sterile.h:89
Float_t e
Definition: plot.C:35
static std::complex< double > one(1.0, 0.0)
std::vector< std::vector< double > > fDelta
delta[i][j] CP violating phase
Definition: PMNS_Sterile.h:86
virtual void InitializeVectors()
static const double kKm2eV
km to eV^-1
Definition: PMNSOpt.h:41
static const double kGf
Definition: PMNSOpt.h:46
static std::complex< double > zero(0.0, 0.0)