PMNS_DensityMatrix.cxx

Go to the documentation of this file.

 
//
// Base class for methods using density matrices
//
// This  class inherits from the PMNS_Base class
//
// jcoelho\@apc.in2p3.fr
 
#include <cassert>
 
#include "MatrixDecomp/zheevh3.h"
 
#include "PMNS_DensityMatrix.h"
#include "exceptions.h"
 
using namespace OscProb;
 
using namespace std;
 
//.............................................................................
PMNS_DensityMatrix::PMNS_DensityMatrix()
    : PMNS_Fast(), fRho(3, vectorC(3, 0)), fMBuffer(3, vectorC(3, 0))
{
}
 
//.............................................................................
PMNS_DensityMatrix::~PMNS_DensityMatrix() {}
 
//.............................................................................
void PMNS_DensityMatrix::RotateState(bool to_basis, matrixC U)
{
  // buffer = rho . U
  for (int i = 0; i < fNumNus; i++) {
    for (int j = 0; j < fNumNus; j++) {
      fMBuffer[i][j] = 0;
      for (int k = 0; k < fNumNus; k++) {
        if (to_basis)
          fMBuffer[i][j] += fRho[i][k] * U[k][j];
        else
          fMBuffer[i][j] += fRho[i][k] * conj(U[j][k]);
      }
    }
  }
 
  // rho = U^\dagger . buffer = U^\dagger . rho . U
  // Final matrix is Hermitian, so copy upper to lower triangle
  for (int i = 0; i < fNumNus; i++) {
    for (int j = i; j < fNumNus; j++) {
      fRho[i][j] = 0;
      for (int k = 0; k < fNumNus; k++) {
        if (to_basis)
          fRho[i][j] += conj(U[k][i]) * fMBuffer[k][j];
        else
          fRho[i][j] += U[i][k] * fMBuffer[k][j];
      }
      if (j > i) fRho[j][i] = conj(fRho[i][j]);
    }
  }
}
 
//.............................................................................
void PMNS_DensityMatrix::ResetToFlavour(int flv)
{
  assert(flv >= 0 && flv < fNumNus);
 
  PMNS_Base::ResetToFlavour(flv);
 
  for (int i = 0; i < fNumNus; ++i) {
    for (int j = 0; j < fNumNus; ++j) {
      if (i == flv && i == j)
        fRho[i][j] = one;
      else
        fRho[i][j] = zero;
    }
  }
}
 
//.............................................................................
double PMNS_DensityMatrix::P(int flv)
{
  assert(flv >= 0 && flv < fNumNus);
  return abs(fRho[flv][flv]);
}
 
//.............................................................................
void PMNS_DensityMatrix::SetPureState(vectorC nu_in)
{
  assert(nu_in.size() == fNumNus);
 
  for (int i = 0; i < fNumNus; i++) {
    for (int j = 0; j < fNumNus; j++) {
      fRho[i][j] = conj(nu_in[i]) * nu_in[j];
    }
  }
}
 
//.............................................................................
void PMNS_DensityMatrix::SetInitialRho(matrixC rho_in)
{
  assert(rho_in.size() == fNumNus);
  assert(rho_in[0].size() == fNumNus);
 
  double eps = 1e-12;
 
  double   trace = 0;
  complexD rho_array[3][3];
  for (int i = 0; i < fNumNus; i++) {
    trace += abs(rho_in[i][i]);
    for (int j = i; j < fNumNus; j++) {
      // Ensure rho_in is hermitian
      assert(abs(rho_in[i][j] - conj(rho_in[j][i])) < eps);
      rho_array[i][j] = rho_in[i][j];
    }
  }
 
  // Ensure rho_in has trace 1
  assert(abs(trace - 1) < eps);
 
  double   fEvalGLoBES[3];
  complexD fEvecGLoBES[3][3];
 
  // Find the eigenvalues of rho_in
  zheevh3(rho_array, fEvecGLoBES, fEvalGLoBES);
 
  // Ensure rho_in is positive definite
  for (int i = 0; i < fNumNus; i++) assert(fEvalGLoBES[i] >= -eps);
 
  fRho = rho_in;
}
 
//.............................................................................
matrixD PMNS_DensityMatrix::ProbMatrix(int nflvi, int nflvf)
{
  THROW_ON_INVALID_ARG(nflvi >= 0, nflvi);
  THROW_ON_INVALID_ARG(nflvi <= fNumNus, nflvi, fNumNus);
  THROW_ON_INVALID_ARG(nflvf >= 0, nflvf);
  THROW_ON_INVALID_ARG(nflvf <= fNumNus, nflvf, fNumNus);
 
  // Output probabilities
  matrixD probs(nflvi, vectorD(nflvf));
 
  // List of states
  vector<matrixC> allstates(nflvi, matrixC(fNumNus, vectorC(fNumNus)));
 
  // Reset all initial states
  for (int i = 0; i < nflvi; i++) {
    ResetToFlavour(i);
    allstates[i] = fRho;
  }
 
  // Propagate all states in parallel
  for (int i = 0; i < int(fNuPaths.size()); i++) {
    for (int flvi = 0; flvi < nflvi; flvi++) {
      fRho = allstates[flvi];
      PropagatePath(fNuPaths[i]);
      allstates[flvi] = fRho;
    }
  }
 
  // Get all probabilities
  for (int flvi = 0; flvi < nflvi; flvi++) {
    for (int flvj = 0; flvj < nflvf; flvj++) {
      probs[flvi][flvj] = abs(allstates[flvi][flvj][flvj]);
    }
  }
 
  return probs;
}