/////////////////////////////////////////////////////////////////////////////// /// class OscProb::PMNS_LIV /// /// Implementation of neutrino oscillations in matter in a /// three-neutrino framework with LIV as modelled by the SME. /// The SME coefficients are included up to the 8th order, /// following the approach described in /// https://doi.org/10.1103/PhysRevD.85.096005. /// /// This developement is part of the QGRANT project with /// ID: 101068013, /// founded by the HORIZON-MSCA-2021-PF-01-01 programme. /// /// \author Nafis R. K. Chowdhury - nrkhanchowdhury\@km3net.de /// \author Joao Coelho - jcoelho\@apc.in2p3.fr /// \author Alba Domi - alba.domi\@fau.de /////////////////////////////////////////////////////////////////////////////// #include #include "PMNS_LIV.h" using namespace OscProb; using namespace std; //............................................................................. /// Constructor. /// This class is restricted to 3 neutrino flavours. /// By default, all LIV coefficients are set to zero. /// PMNS_LIV::PMNS_LIV() : PMNS_Fast() { SetStdPath(); for (int dim = 3; dim < 8; dim += 2) { for (int flvi = 0; flvi < 3; flvi++) { for (int flvj = flvi; flvj < 3; flvj++) { SetaT(flvi, flvj, dim, 0, 0); SetcT(flvi, flvj, dim + 1, 0, 0); } } } } //............................................................................. /// /// Destructor. /// PMNS_LIV::~PMNS_LIV() {} //............................................................................. /// Set any given aT parameter of a chosen dimension. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - First flavour index /// @param flvj - Second flavour index /// @param dim - Dimension of the coefficient: can be 3,5,7. /// @param val - Absolute value of the parameter /// @param phase - Complex phase of the parameter in radians /// void PMNS_LIV::SetaT(int flvi, int flvj, int dim, double val, double phase) { if (flvi > flvj) { cerr << "WARNING: First argument should be smaller or equal to second " "argument" << endl << "Setting reverse order (aT_" << flvj << flvi << "). " << endl; int temp = flvi; flvi = flvj; flvj = temp; } if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) { cerr << "WARNING: aT_" << flvi << flvj << " not valid for " << fNumNus << " neutrinos. Doing nothing." << endl; return; } if (dim != 3 && dim != 5 && dim != 7) { cerr << "WARNING: Invalid aT coefficient dimension dim=" << dim << " not in [3,5,7].\n"; return; } int pos = (dim - 3) / 2; complexD h = val; if (flvi != flvj) h *= complexD(cos(phase), sin(phase)); bool isSame = (faT[flvi][flvj][pos] == h); if (!isSame) ClearCache(); fGotES *= isSame; faT[flvi][flvj][pos] = h; } //............................................................................. /// Set any given cT parameter of a chosen dimension. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - First flavour index /// @param flvj - Second flavour index /// @param dim - Dimension of the coefficient: can be 4,6,8. /// @param val - Absolute value of the parameter /// @param phase - Complex phase of the parameter in radians /// void PMNS_LIV::SetcT(int flvi, int flvj, int dim, double val, double phase) { if (flvi > flvj) { cerr << "WARNING: First argument should be smaller or equal to second " "argument" << endl << "Setting reverse order (cT_" << flvj << flvi << "). " << endl; int temp = flvi; flvi = flvj; flvj = temp; } if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) { cerr << "WARNING: cT_" << flvi << flvj << " not valid for " << fNumNus << " neutrinos. Doing nothing." << endl; return; } if (dim != 4 && dim != 6 && dim != 8) { cerr << "WARNING: Invalid cT coefficient dimension dim=" << dim << " not in [4,6,8].\n"; return; } int pos = dim / 2 - 2; complexD h = val; if (flvi != flvj) h *= complexD(cos(phase), sin(phase)); bool isSame = (fcT[flvi][flvj][pos] == h); if (!isSame) ClearCache(); fGotES *= isSame; fcT[flvi][flvj][pos] = h; } //............................................................................. /// /// Get any given aT parameter of a chosen dimension. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - The first flavour index /// @param flvj - The second flavour index /// @param dim - Dimension of the coefficient: can be 3,5,7. /// /// @return - The aT parameter value /// complexD PMNS_LIV::GetaT(int flvi, int flvj, int dim) { if (flvi > flvj) { cerr << "WARNING: First argument should be smaller or equal to second " "argument" << endl << "Setting reverse order (aT_" << flvj << flvi << "). " << endl; int temp = flvi; flvi = flvj; flvj = temp; } if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) { cerr << "WARNING: aT_" << flvi << flvj << " not valid for " << fNumNus << " neutrinos. Returning 0." << endl; return zero; } if (dim != 3 && dim != 5 && dim != 7) { cerr << "WARNING: Invalid aT coefficient dimension dim=" << dim << " not in [3,5,7].\n"; return zero; } int pos = (dim - 3) / 2; return faT[flvi][flvj][pos]; } //............................................................................. /// /// Get any given cT parameter of a chosen dimension. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - The first flavour index /// @param flvj - The second flavour index /// @param dim - Dimension of the coefficient: can be 4,6,8. /// /// @return - The cT parameter value /// complex PMNS_LIV::GetcT(int flvi, int flvj, int dim) { if (flvi > flvj) { cerr << "WARNING: First argument should be smaller or equal to second " "argument" << endl << "Setting reverse order (cT_" << flvj << flvi << "). " << endl; int temp = flvi; flvi = flvj; flvj = temp; } if (flvi < 0 || flvi > 2 || flvj < flvi || flvj > 2) { cerr << "WARNING: cT_" << flvi << flvj << " not valid for " << fNumNus << " neutrinos. Returning 0." << endl; return zero; } if (dim != 4 && dim != 6 && dim != 8) { cerr << "WARNING: Invalid cT coefficient dimension dim=" << dim << " not in [4,6,8].\n"; return zero; } int pos = dim / 2 - 2; return fcT[flvi][flvj][pos]; } //............................................................................. /// /// Build the full LIV Hamiltonian in matter /// void PMNS_LIV::UpdateHam() { double lv = 2 * kGeV2eV * fEnergy; // 2*E in eV // Set the vacuum Hamiltonian for (int i = 0; i < fNumNus; i++) { for (int j = i; j < fNumNus; j++) { fHam[i][j] = fHms[i][j] / lv; } } // Add matter potential double kr2GNe = kK2 * M_SQRT2 * kGf * fPath.density * fPath.zoa; if (!fIsNuBar) fHam[0][0] += kr2GNe; else fHam[0][0] -= kr2GNe; // Add LIV terms for all dimensions double energy_pow = kGeV2eV; for (int dim = 3; dim < 9; dim++) { // Increase energy power if (dim > 3) energy_pow *= fEnergy; // Set parameter index int coeff_idx = (dim - 3) / 2; for (int i = 0; i < fNumNus; i++) { for (int j = i; j < fNumNus; j++) { complexD liv_term = energy_pow; if (dim % 2 == 1) { // aT is CPT-odd if (fIsNuBar) liv_term *= -faT[i][j][coeff_idx]; else liv_term *= faT[i][j][coeff_idx]; } else { // cT is CPT-even // Add 4/3 coefficient for backward compatibility if (dim == 4) liv_term *= -4 / 3. * fcT[i][j][coeff_idx]; else liv_term *= -fcT[i][j][coeff_idx]; } fHam[i][j] += liv_term; } } } // Conjugate Hamiltonian for antineutrinos if (fIsNuBar) { for (int i = 0; i < fNumNus; i++) { for (int j = i + 1; j < fNumNus; j++) { fHam[i][j] = conj(fHam[i][j]); } } } } /////////////////////////////////////////////////////////////////////////////// // // Obsolete functions for backward compatibility... // /////////////////////////////////////////////////////////////////////////////// //............................................................................. /// /// Set all LIV parameters at once. /// /// @param aT_ee - The absolute value of the parameter aT_ee /// @param aT_mumu - The absolute value of the parameter aT_mumu /// @param aT_tautau - The absolute value of the parameter aT_tautau /// @param aT_emu - The absolute value of the parameter aT_emu /// @param aT_etau - The absolute value of the parameter aT_etau /// @param aT_mutau - The absolute value of the parameter aT_mutau /// @param delta_aT_emu - The phase of the complex parameter aT_emu in radians /// @param delta_aT_etau - The phase of the complex parameter aT_etau in /// radians /// @param delta_aT_mutau - The phase of the complex parameter aT_mutau in /// radians /// @param cT_ee - The absolute value of the parameter cT_ee /// @param cT_mumu - The absolute value of the parameter cT_mumu /// @param cT_tautau - The absolute value of the parameter cT_tautau /// @param cT_emu - The absolute value of the parameter cT_emu /// @param cT_etau - The absolute value of the parameter cT_etau /// @param cT_mutau - The absolute value of the parameter cT_mutau /// @param delta_cT_emu - The phase of the complex parameter cT_emu in radians /// @param delta_cT_etau - The phase of the complex parameter cT_etau in /// radians /// @param delta_cT_mutau - The phase of the complex parameter cT_mutau in /// radians /// void PMNS_LIV::SetLIV(double aT_ee, double aT_mumu, double aT_tautau, double aT_emu, double aT_etau, double aT_mutau, double cT_ee, double cT_mumu, double cT_tautau, double cT_emu, double cT_etau, double cT_mutau, double delta_aT_emu, double delta_aT_etau, double delta_aT_mutau, double delta_cT_emu, double delta_cT_etau, double delta_cT_mutau) { SetaT(0, 0, aT_ee, 0); SetaT(1, 1, aT_mumu, 0); SetaT(2, 2, aT_tautau, 0); SetaT(0, 1, aT_emu, delta_aT_emu); SetaT(0, 2, aT_etau, delta_aT_etau); SetaT(1, 2, aT_mutau, delta_aT_mutau); SetcT(0, 0, cT_ee, 0); SetcT(1, 1, cT_mumu, 0); SetcT(2, 2, cT_tautau, 0); SetcT(0, 1, cT_emu, delta_cT_emu); SetcT(0, 2, cT_etau, delta_cT_etau); SetcT(1, 2, cT_mutau, delta_cT_mutau); } //............................................................................. /// /// Set any given aT parameter with dimension 3. /// /// This will check if value is changing to keep track of whether /// the eigensystem needs to be recomputed. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - The first flavour index /// @param flvj - The second flavour index /// @param val - The absolute value of the parameter /// @param phase - The complex phase of the parameter in radians /// void PMNS_LIV::SetaT(int flvi, int flvj, double val, double phase) { SetaT(flvi, flvj, 3, val, phase); } //............................................................................. /// /// Set any given cT parameter with dimension 4. /// /// This will check if value is changing to keep track of whether /// the eigensystem needs to be recomputed. /// /// Flavours are:\n /// - 0 = nue, 1 = numu, 2 = nutau /// /// Requires that flvi < flvj. Will notify you if input is wrong. /// If flvi > flvj, will assume reverse order and swap flvi and flvj. /// /// @param flvi - The first flavour index /// @param flvj - The second flavour index /// @param val - The absolute value of the parameter /// @param phase - The complex phase of the parameter in radians /// void PMNS_LIV::SetcT(int flvi, int flvj, double val, double phase) { SetcT(flvi, flvj, 4, val, phase); } //............................................................................. /// /// Set aT_ee parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetaT_ee(double a) { SetaT(0, 0, a, 0); } //............................................................................. /// /// Set aT_mumu parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetaT_mumu(double a) { SetaT(1, 1, a, 0); } //............................................................................. /// /// Set aT_tautau parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetaT_tautau(double a) { SetaT(2, 2, a, 0); } //............................................................................. /// /// Set aT_emu parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetaT_emu(double a, double phi) { SetaT(0, 1, a, phi); } //............................................................................. /// /// Set aT_etau parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetaT_etau(double a, double phi) { SetaT(0, 2, a, phi); } //............................................................................. /// /// Set aT_mutau parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetaT_mutau(double a, double phi) { SetaT(1, 2, a, phi); } //............................................................................. /// /// Set cT_ee parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetcT_ee(double a) { SetcT(0, 0, a, 0); } //............................................................................. /// /// Set cT_mumu parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetcT_mumu(double a) { SetcT(1, 1, a, 0); } //............................................................................. /// /// Set cT_tautau parameter /// /// @param a - The absolute value of the parameter /// void PMNS_LIV::SetcT_tautau(double a) { SetcT(2, 2, a, 0); } //............................................................................. /// /// Set cT_emu parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetcT_emu(double a, double phi) { SetcT(0, 1, a, phi); } //............................................................................. /// /// Set cT_etau parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetcT_etau(double a, double phi) { SetcT(0, 2, a, phi); } //............................................................................. /// /// Set cT_mutau parameter /// /// @param a - The absolute value of the parameter /// @param phi - The complex phase of the parameter in radians /// void PMNS_LIV::SetcT_mutau(double a, double phi) { SetcT(1, 2, a, phi); } ///////////////////////////////////////////////////////////////////////////////