diff --git a/run.sh b/run.sh
index 615f08f23..e973ec5bd 100755
--- a/run.sh
+++ b/run.sh
@@ -5,8 +5,9 @@ cd test
 rm -R output-plasma
 rm -R output-plasma-refine2
 
-cp Restart/lowP/restart_output-plasma.LTE.h5 restart_output-plasma.sol.h5
-cp Restart/highP/restart_output-plasma-refine2.sol.h5 restart_output-plasma-refine2.sol.h5
+# cp Restart/lowP/restart_output-plasma.LTE.h5 restart_output-plasma.sol.h5
+# cp Restart/highP/restart_output-plasma-refine2.LTE.h5 restart_output-plasma-refine2.sol.h5
+cp Restart/highP/restart_output-plasma-refine2.NLTE.h5 restart_output-plasma-refine2.sol.h5
 
 
 
diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp
index 77b1c269f..566cdc946 100644
--- a/src/M2ulPhyS.cpp
+++ b/src/M2ulPhyS.cpp
@@ -156,9 +156,12 @@ void M2ulPhyS::initMixtureAndTransportModels() {
           gpu::instantiateDeviceChemistry<<<1, 1>>>(d_mixture, config.chemistryInput, chemistry_);
 #else
           chemistry_ = new Chemistry(mixture, config.chemistryInput);
+
 #endif
           break;
       }
+
+
       break;
     case WorkingFluid::LTE_FLUID:
       int table_dim;
@@ -3367,7 +3370,6 @@ void M2ulPhyS::parseReactionInputs() {
 
     if (model == "arrhenius") {
       config.reactionModels[r - 1] = ARRHENIUS;
-
       double A, b, E;
       tpsP->getRequiredInput((basepath + "/arrhenius/A").c_str(), A);
       tpsP->getRequiredInput((basepath + "/arrhenius/b").c_str(), b);
@@ -3386,6 +3388,9 @@ void M2ulPhyS::parseReactionInputs() {
       config.reactionModels[r - 1] = TABULATED_RXN;
       std::string inputPath(basepath + "/tabulated");
       readTable(inputPath, config.chemistryInput.reactionInputs[r - 1].tableInput);
+
+    } else if (model == "radiative_decay") {
+      config.reactionModels[r - 1] = RADIATIVE_DECAY;
     } else {
       grvy_printf(GRVY_ERROR, "\nUnknown reaction_model -> %s", model.c_str());
       exit(ERROR);
@@ -3457,6 +3462,7 @@ void M2ulPhyS::parseReactionInputs() {
 
       // energy conservation.
       // TODO(kevin): this will need an adjustion when radiation comes into play.
+      // if (config.reactionModels[r] != RADIATIVE_DECAY) {
       double reactEnergy = 0.0, prodEnergy = 0.0;
       for (int sp = 0; sp < numSpecies_; sp++) {
         // int inputSp = (*mixtureToInputMap_)[sp];
@@ -3474,9 +3480,12 @@ void M2ulPhyS::parseReactionInputs() {
         grvy_printf(GRVY_ERROR, "Difference = %.8E\n", delta);
         exit(-1);
       }
+      // }
+
     }
   }
 
+
   // Pack up chemistry input for instantiation.
   {
     config.chemistryInput.model = config.chemistryModel_;
@@ -3486,9 +3495,11 @@ void M2ulPhyS::parseReactionInputs() {
     } else {
       config.chemistryInput.electronIndex = -1;
     }
+    config.chemistryInput.speciesMapping = &config.speciesMapping;
+    config.chemistryInput.speciesNames = &config.speciesNames;
 
     size_t rxn_param_idx = 0;
-
+    bool RadiativeDecayReactionsIncluded = false;
     config.chemistryInput.numReactions = config.numReactions;
     for (int r = 0; r < config.numReactions; r++) {
       config.chemistryInput.reactionEnergies[r] = config.reactionEnergies[r];
@@ -3504,11 +3515,16 @@ void M2ulPhyS::parseReactionInputs() {
             config.equilibriumConstantParams[p + r * gpudata::MAXCHEMPARAMS];
       }
 
-      if (config.reactionModels[r] != TABULATED_RXN) {
+      if (config.reactionModels[r] == ARRHENIUS ||  config.reactionModels[r] == HOFFERTLIEN) {
         assert(rxn_param_idx < config.rxnModelParamsHost.size());
         config.chemistryInput.reactionInputs[r].modelParams = config.rxnModelParamsHost[rxn_param_idx].Read();
         rxn_param_idx += 1;
       }
+
+      if (config.reactionModels[r] == RADIATIVE_DECAY) RadiativeDecayReactionsIncluded = true;
+    }
+    if (RadiativeDecayReactionsIncluded) {
+      tpsP->getRequiredInput("reactions/characteristic_length", config.chemistryInput.char_length);
     }
   }
 }
@@ -4287,7 +4303,7 @@ void M2ulPhyS::updateVisualizationVariables() {
       Th = prim[1 + _nvel];
       Te = (in_mix->IsTwoTemperature()) ? prim[_num_equation - 1] : Th;
       double kfwd[gpudata::MAXREACTIONS], kC[gpudata::MAXREACTIONS];
-      in_chem->computeForwardRateCoeffs(Th, Te, kfwd);
+      in_chem->computeForwardRateCoeffs(nsp, Th, Te, kfwd);
       in_chem->computeEquilibriumConstants(Th, Te, kC);
       // get reaction rates
       double progressRates[gpudata::MAXREACTIONS];
diff --git a/src/chemistry.cpp b/src/chemistry.cpp
index 631139232..bf0e40268 100644
--- a/src/chemistry.cpp
+++ b/src/chemistry.cpp
@@ -37,7 +37,9 @@ using namespace std;
 
 Chemistry::Chemistry(GasMixture *mixture, RunConfiguration &config) : Chemistry(mixture, config.chemistryInput) {}
 
-MFEM_HOST_DEVICE Chemistry::Chemistry(GasMixture *mixture, const ChemistryInput &inputs) : mixture_(mixture) {
+MFEM_HOST_DEVICE Chemistry::Chemistry(GasMixture *mixture, const ChemistryInput &inputs) 
+: mixture_(mixture),
+  inputs_(inputs) {
   numEquations_ = mixture->GetNumEquations();
   numSpecies_ = mixture->GetNumSpecies();
   numActiveSpecies_ = mixture->GetNumActiveSpecies();
@@ -86,6 +88,12 @@ MFEM_HOST_DEVICE Chemistry::Chemistry(GasMixture *mixture, const ChemistryInput
       case TABULATED_RXN: {
         reactions_[r] = new Tabulated(inputs.reactionInputs[r].tableInput);
       } break;
+      case RADIATIVE_DECAY: {
+        reactions_[r] = new RadiativeDecay(inputs.char_length, inputs.speciesMapping, inputs.speciesNames, 
+                                            numSpecies_,
+                                            &reactantStoich_[0 + r * numSpecies_],
+                                            &productStoich_[0 + r * numSpecies_]);    
+      } break;
       default:
         printf("Unknown reactionModel.");
         assert(false);
@@ -106,13 +114,13 @@ MFEM_HOST_DEVICE Chemistry::~Chemistry() {
   }
 }
 
-void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, Vector &kfwd) {
+void Chemistry::computeForwardRateCoeffs(const mfem::Vector &ns, const double &T_h, const double &T_e, Vector &kfwd) {
   kfwd.SetSize(numReactions_);
 
   const double Thlim = max(T_h, min_temperature_);
   const double Telim = max(T_e, min_temperature_);
 
-  computeForwardRateCoeffs(Thlim, Telim, &kfwd[0]);
+  computeForwardRateCoeffs(&ns[0], Thlim, Telim, &kfwd[0]);
   // kfwd = 0.0;
   //
   // for (int r = 0; r < numReactions_; r++) {
@@ -123,8 +131,7 @@ void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, V
   return;
 }
 
-MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double &T_h, const double &T_e, double *kfwd) {
-  // kfwd.SetSize(numReactions_);
+MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double *ns, const double &T_h, const double &T_e, double *kfwd) {
   for (int r = 0; r < numReactions_; r++) kfwd[r] = 0.0;
 
   const double Thlim = max(T_h, min_temperature_);
@@ -132,7 +139,11 @@ MFEM_HOST_DEVICE void Chemistry::computeForwardRateCoeffs(const double &T_h, con
 
   for (int r = 0; r < numReactions_; r++) {
     bool isElectronInvolved = isElectronInvolvedAt(r);
-    kfwd[r] = reactions_[r]->computeRateCoefficient(Thlim, Telim, isElectronInvolved);
+    if (inputs_.reactionModels[r] == RADIATIVE_DECAY) {
+      kfwd[r] = reactions_[r]->computeRateCoefficient(ns, Thlim, Telim);
+    } else {
+      kfwd[r] = reactions_[r]->computeRateCoefficient(Thlim, Telim, isElectronInvolved);
+    }
   }
 
   return;
@@ -240,7 +251,7 @@ MFEM_HOST_DEVICE void Chemistry::computeCreationRate(const double *progressRate,
   for (int sp = 0; sp < numSpecies_; sp++) {
     for (int r = 0; r < numReactions_; r++) {
       creationRate[sp] +=
-          progressRate[r] * (productStoich_[sp + r * numSpecies_] - reactantStoich_[sp + r * numSpecies_]);
+        progressRate[r] * (productStoich_[sp + r * numSpecies_] - reactantStoich_[sp + r * numSpecies_]);
     }
     creationRate[sp] *= mixture_->GetGasParams(sp, GasParams::SPECIES_MW);
   }
diff --git a/src/chemistry.hpp b/src/chemistry.hpp
index d56394195..3448b5eb5 100644
--- a/src/chemistry.hpp
+++ b/src/chemistry.hpp
@@ -77,7 +77,7 @@ class Chemistry {
   double equilibriumConstantParams_[gpudata::MAXCHEMPARAMS * gpudata::MAXREACTIONS];
 
   ChemistryModel model_;
-
+  ChemistryInput inputs_;
   // /*
   //   From input file, reaction/reaction%d/equation will specify this mapping.
   // */
@@ -98,9 +98,9 @@ class Chemistry {
   // return Vector of reaction rate coefficients, with the size of numReaction_.
   // WARNING(marc) I have removed "virtual" qualifier here assuming these functions will not
   // change for child classes. Correct if wrong
-  void computeForwardRateCoeffs(const double &T_h, const double &T_e, Vector &kfwd);
-  MFEM_HOST_DEVICE void computeForwardRateCoeffs(const double &T_h, const double &T_e, double *kfwd);
-
+  void computeForwardRateCoeffs(const Vector &ns, const double &T_h, const double &T_e, Vector &kfwd);
+  MFEM_HOST_DEVICE void computeForwardRateCoeffs(const double *ns, const double &T_h, const double &T_e, 
+                                                 double *kfwd);
   void computeEquilibriumConstants(const double &T_h, const double &T_e, Vector &kC);
   MFEM_HOST_DEVICE void computeEquilibriumConstants(const double &T_h, const double &T_e, double *kC);
 
diff --git a/src/dataStructures.hpp b/src/dataStructures.hpp
index 252422db9..a98e4c419 100644
--- a/src/dataStructures.hpp
+++ b/src/dataStructures.hpp
@@ -74,7 +74,7 @@ enum TransportModel { ARGON_MINIMAL, ARGON_MIXTURE, CONSTANT, LTE_TRANSPORT, MIX
 
 enum ChemistryModel { /* CANTERA, */ NUM_CHEMISTRYMODEL };
 
-enum ReactionModel { ARRHENIUS, HOFFERTLIEN, TABULATED_RXN, NUM_REACTIONMODEL };
+enum ReactionModel { ARRHENIUS, HOFFERTLIEN, TABULATED_RXN, RADIATIVE_DECAY, NUM_REACTIONMODEL };
 
 enum RadiationModel { NONE_RAD, NET_EMISSION, P1_MODEL, NUM_RADIATIONMODEL };
 
@@ -629,7 +629,8 @@ struct ReactionInput {
 
 struct ChemistryInput {
   ChemistryModel model;
-
+  std::map<std::string, int> *speciesMapping; // Needed for rediative decay reactions
+  std::vector<std::string> *speciesNames;
   int electronIndex;
   int numReactions;
   double reactionEnergies[gpudata::MAXREACTIONS];
@@ -645,7 +646,7 @@ struct ChemistryInput {
   ReactionModel reactionModels[gpudata::MAXREACTIONS];
   double equilibriumConstantParams[gpudata::MAXCHEMPARAMS * gpudata::MAXREACTIONS];
   ReactionInput reactionInputs[gpudata::MAXREACTIONS];
-
+  double char_length;
   double minimumTemperature;
 };
 
diff --git a/src/equation_of_state.cpp b/src/equation_of_state.cpp
index cce96aba6..6232aa8bb 100644
--- a/src/equation_of_state.cpp
+++ b/src/equation_of_state.cpp
@@ -2100,46 +2100,4 @@ void PerfectMixture::GetSpeciesFromLTE(double *conserv, double *primit, TableInt
 
   conserv[iTh] = totalEnergy;
   primit[iTh] = T;
-
-
-    // double nh = 0.0;
-           
-    // for (int sp = 0; sp < nPop; sp++) nh = nh + n_sp[sp]; 
-    // nh = nh + n_e + n_e + nB;
-    // double p_Mycalc = nh*  UNIVERSALGASCONSTANT *T;
-
-    // double electronPressure = 0.0;
-    // double p_calc =  PerfectMixture::ComputePressure(conserv, &electronPressure);
-
-    // double rho_calc = 0.0;
-    // for (int sp = 0; sp < nPop; sp++) rho_calc = rho_calc + n_sp[sp]* GetGasParams(sp, GasParams::SPECIES_MW); 
-    // rho_calc = rho_calc + n_sp[iIon1]* GetGasParams(iIon1, GasParams::SPECIES_MW); 
-    // rho_calc = rho_calc + n_sp[iElectron]* GetGasParams(iElectron, GasParams::SPECIES_MW); 
-    // rho_calc = rho_calc + n_sp[iBackground]* GetGasParams(iBackground, GasParams::SPECIES_MW); 
-
-    // std::cout << p_0 << " " << p_calc << " " << p_Mycalc << " " <<  n_e << " " << T << " " << rho_calc << std::endl;
-
-    // std::cout << " " << n_sp[0]*GetGasParams(0, GasParams::SPECIES_MW)/GetGasParams(0, GasParams::SPECIES_DEGENERACY) << " " 
-    //           << " " << n_sp[1]*GetGasParams(1, GasParams::SPECIES_MW)/GetGasParams(1, GasParams::SPECIES_DEGENERACY) << " " 
-    //           << " " << n_sp[2]*GetGasParams(2, GasParams::SPECIES_MW)/GetGasParams(2, GasParams::SPECIES_DEGENERACY) << " " 
-    //           << " " << n_sp[iIon1]*GetGasParams(iIon1, GasParams::SPECIES_MW)/GetGasParams(iIon1, GasParams::SPECIES_DEGENERACY) << " " 
-    //           << " " << n_sp[iElectron]*GetGasParams(iElectron, GasParams::SPECIES_MW)/GetGasParams(iElectron, GasParams::SPECIES_DEGENERACY) << " " 
-    //           << " " << n_sp[iBackground]*GetGasParams(iBackground, GasParams::SPECIES_MW)/GetGasParams(iBackground, GasParams::SPECIES_DEGENERACY) << " "
-    //           << std::endl;
-
-    // std::cout << " " << GetGasParams(0, GasParams::SPECIES_MW)*1000.0 << " " 
-    //           << " " << GetGasParams(1, GasParams::SPECIES_MW)*1000.0 << " " 
-    //           << " " << GetGasParams(2, GasParams::SPECIES_MW)*1000.0 << " " 
-    //           << " " << GetGasParams(iIon1, GasParams::SPECIES_MW)*1000.0 << " " 
-    //           << " " << GetGasParams(iElectron, GasParams::SPECIES_MW)*1000.0 << " " 
-    //           << " " << GetGasParams(iBackground, GasParams::SPECIES_MW)*1000.0 << " "
-    //           << std::endl;
-
-
-    // exit(0);
-
-
-
-
-
 }
diff --git a/src/equation_of_state.hpp b/src/equation_of_state.hpp
index b25ca93c3..779849289 100644
--- a/src/equation_of_state.hpp
+++ b/src/equation_of_state.hpp
@@ -57,10 +57,11 @@ static const double AVOGADRONUMBER = 6.0221409e+23;    // mol^(-1)
 static const double BOLTZMANNCONSTANT = UNIVERSALGASCONSTANT / AVOGADRONUMBER;
 static const double PLANCKCONSTANT = 6.62607015e-34;  // m^2 kg / s
 static const double VACUUMPERMITTIVITY = 8.8541878128e-12;
-static const double ELECTRONCHARGE = 1.60218e-19;
+static const double ELECTRONCHARGE = 1.602176634e-19;
 static const double MOLARELECTRONCHARGE = ELECTRONCHARGE * AVOGADRONUMBER;
 static const double ELECTRONMASS = 9.1093837015e-31;  // kg
 static const double qeOverkB = ELECTRONCHARGE / BOLTZMANNCONSTANT;
+static const double SPEEDOFLIGHT = 299792458.0; // m/s
 
 static const double IonizationEnergy_Argon = 15.7596119;  // eV
 
@@ -630,6 +631,7 @@ MFEM_HOST_DEVICE inline double DryAir::ComputeTemperature(const double *state) {
 ////// Perfect Mixture GasMixture                     ////////////////////
 //////////////////////////////////////////////////////////////////////////
 
+
 class PerfectMixture : public GasMixture {
  private:
   // Vector specificHeatRatios_;
diff --git a/src/reaction.cpp b/src/reaction.cpp
index 54c8baa28..62b2d1a60 100644
--- a/src/reaction.cpp
+++ b/src/reaction.cpp
@@ -75,3 +75,198 @@ MFEM_HOST_DEVICE double Tabulated::computeRateCoefficient(const double &T_h, con
   double temp = (isElectronInvolved) ? T_e : T_h;
   return table_->eval(temp);
 }
+
+
+
+
+
+MFEM_HOST_DEVICE RadiativeDecay::RadiativeDecay(const double &_R, 
+                                                std::map<std::string, int> *_speciesMapping, std::vector<std::string> *_speciesNames,
+                                                const int _numSpecies, const double *_reactantStoich, const double *_productStoich) 
+  : Reaction(), R(_R) {
+
+  rank0_ = Mpi::Root();
+  L = 2.0 * R;
+
+
+  int numOfReactants = 0, numOfProducts = 0;
+  for (int sp = 0; sp < _numSpecies; sp++) {
+    numOfReactants = numOfReactants + _reactantStoich[sp];   
+    numOfProducts = numOfProducts + _productStoich[sp];   
+  }
+  assert(numOfReactants ==1);
+  assert(numOfProducts ==1);
+
+  for (int sp = 0; sp < _numSpecies; sp++) {
+    if (_reactantStoich[sp] > 0) upper_sp_name = (*_speciesNames)[sp];                                              
+    if (_productStoich[sp] > 0) lower_sp_name = (*_speciesNames)[sp];
+  }
+
+  iAr_u = (*_speciesMapping)[upper_sp_name];
+  iAr_l = (*_speciesMapping)[lower_sp_name];
+
+
+  bool flag = ((upper_sp_name == "Ar_r" && lower_sp_name == "Ar")    || 
+               (upper_sp_name == "Ar_p" && lower_sp_name == "Ar_r")  || 
+               (upper_sp_name == "Ar_p" && lower_sp_name == "Ar_m"));
+  assert(flag);
+
+  if (upper_sp_name == "Ar_r" || upper_sp_name == "Ar") {
+    NumOfInteral_lvl_u = E_lvl_r.size();
+    E_lvl_u = &E_lvl_r;
+    g_lvl_u = &g_lvl_r;
+    n_sp_lvl_u.resize(E_lvl_r.size());
+    E_lvl_l = &E_lvl_g;
+    g_lvl_l = &g_lvl_g;
+    n_sp_lvl_l.resize(E_lvl_g.size());
+    Aji = &Aji_r_g;
+  } else if (upper_sp_name == "Ar_p" && lower_sp_name == "Ar_m") {
+    NumOfInteral_lvl_u = E_lvl_4p.size();
+    E_lvl_u = &E_lvl_4p;
+    g_lvl_u = &g_lvl_4p;
+    n_sp_lvl_u.resize(E_lvl_4p.size());
+    E_lvl_l = &E_lvl_m;
+    g_lvl_l = &g_lvl_m;
+    n_sp_lvl_l.resize(E_lvl_m.size());
+    Aji = &Aji_4p_m;
+  } else if (upper_sp_name == "Ar_p" && lower_sp_name == "Ar_r") {
+    NumOfInteral_lvl_u = E_lvl_4p.size();
+    E_lvl_u = &E_lvl_4p;
+    g_lvl_u = &g_lvl_4p;
+    n_sp_lvl_u.resize(E_lvl_4p.size());
+    E_lvl_l = &E_lvl_r;
+    g_lvl_l = &g_lvl_r;
+    n_sp_lvl_l.resize(E_lvl_r.size());
+    Aji = &Aji_4p_r;
+  } else {
+    if (rank0_) std::cout << " Species " << upper_sp_name.c_str() << " not recognized for this reactive. " << std::endl;
+    MPI_Barrier(MPI_COMM_WORLD);    
+    assert(false);
+  }
+
+  effAcoef_lvl.resize(NumOfInteral_lvl_u);
+  std::fill(effAcoef_lvl.begin(), effAcoef_lvl.end(), 0.0);
+
+
+}
+
+MFEM_HOST_DEVICE RadiativeDecay::~RadiativeDecay() {
+
+}
+
+
+
+
+MFEM_HOST_DEVICE double RadiativeDecay::computeRateCoefficient(const double *nsp,
+                                                               const double &T_h, const double &T_e) {
+
+
+  double n_sp_u = nsp[iAr_u]; // Find the index corresponding to the correct species
+  double n_sp_l = nsp[iAr_l];
+
+  GetNumDensityOfInteralLevels(n_sp_lvl_u.size(),n_sp_u, T_e, &n_sp_lvl_u[0]);                                                       
+  GetNumDensityOfInteralLevels(n_sp_lvl_l.size(),n_sp_l, T_e, &n_sp_lvl_l[0]);                                                       
+
+  GetEinsteinACoefficient(T_h, T_e, &effAcoef_lvl[0]);                                                       
+
+  // The effective Einstein A Coefficient of a lumped level is number-density 
+  // weighted over the internal levels.
+  double effAcoef = 0.0;  
+  for (int i_lvl = 0; i_lvl < NumOfInteral_lvl_u; i_lvl++) {
+    effAcoef += n_sp_lvl_u[i_lvl] * effAcoef_lvl[i_lvl] ;  
+  }
+  effAcoef = effAcoef / (n_sp_u + small);
+
+  effAcoef = max(min(effAcoef,1.0),0.0);  
+  return effAcoef;
+}
+
+MFEM_HOST_DEVICE void RadiativeDecay::GetNumDensityOfInteralLevels(const int NumOfInteral_lvl, const double &n_sp, 
+                                                       const double &T_e, double *n_sp_internal) {
+
+  // n_sp in [mol/m^3] is the number concentration of a lumped state/level.  
+  // Here, we evaluate the number concentration of the internal levels assuming that they 
+  // are Boltzmann distributed.                                               
+  // E_lvl_u in [J/mol]
+  // T_e in [K] is the Electron / Excitation Temperature 
+  // Should we use T_h instead?
+
+  double Q_n = 0.0;  // partition function of internal levels
+  for (int i_lvl = 0; i_lvl < NumOfInteral_lvl; i_lvl++) {
+    Q_n += (*g_lvl_u)[i_lvl] * exp(-(*E_lvl_u)[i_lvl] / UNIVERSALGASCONSTANT / (T_e + small));  // E0 in J/mol, so e/kT = E0/NA/RT
+  }
+
+  for (int i_lvl = 0; i_lvl < NumOfInteral_lvl; i_lvl++) {  
+    n_sp_internal[i_lvl] = (*g_lvl_u)[i_lvl] * exp(-(*E_lvl_u)[i_lvl] / UNIVERSALGASCONSTANT / (T_e + small)) * n_sp / Q_n;  // E0 in J/mol, so e/kT = E0/NA/RT
+  }
+}
+
+
+MFEM_HOST_DEVICE void RadiativeDecay::GetEinsteinACoefficient(const double &T_h, const double &T_e, 
+                                                              double *effAcoef) {
+
+  for (int i_lvl = 0; i_lvl < NumOfInteral_lvl_u; i_lvl++) {
+    for (int itrans = 0; itrans < int((*Aji)[i_lvl].size()); itrans++) 
+    { 
+      double Acoef = (*Aji)[i_lvl][itrans];
+      double eta = escapeFactCalc(n_sp_lvl_l[itrans], 
+                  (*E_lvl_u)[i_lvl], (*E_lvl_l)[itrans], 
+                  (*g_lvl_u)[i_lvl], (*g_lvl_l)[itrans], Acoef, T_h);  
+      effAcoef[i_lvl] = effAcoef[i_lvl] + eta * Acoef;
+    } 
+
+
+  }
+}
+
+
+MFEM_HOST_DEVICE double RadiativeDecay::escapeFactCalc(const double &n_i, const double &E_j, const double &E_i, 
+                                                       const double &g_j, const double &g_i, const double &A_ji, 
+                                                       const double &T_g) { 
+  //  Calculations for escape factor
+  // i -> lower level
+  // j -> upper level
+  // Ar(j) -> Ar(i) + hv 
+  // R =  n_j * A_ji * eta        
+
+  //  n_i in [mol/m^3] is the number density of the upper level 
+  //  E_i in [J/mol] is the energy of the upper level 
+  //  g_i in [] is the degeneracy of the upper level 
+  //  A_ji in [1/s] is the Einstein A coeficient of transition j->i
+
+  //  T_g in [K] is the Temperature 
+
+  double eta = 1.0; 
+
+  // wavelength of transition                                                      
+  double lambda_0 = PLANCKCONSTANT * SPEEDOFLIGHT / 
+                    ( (E_j-E_i) / AVOGADRONUMBER );
+
+  // absorption coefficient at line center, for Doppler absorption 
+  double k0 = pow(lambda_0,3) * (n_i * AVOGADRONUMBER) * g_j * A_ji * sqrt(M_Ar) / 
+              (8.0 * PI * g_i * sqrt(2.0 * BOLTZMANNCONSTANT * PI * T_g));
+
+  double q0 = R;
+  double Lq = L/(2*q0);              
+
+  // compute escape factor   
+  // Mewe (1967)
+  // eta = ( 2.0 - exp(-1e-3 * k0 * R) ) / (1.0 + k0*R) ;
+  if ((k0*(L/2) > 1) && (k0*q0 > 1)) {
+    // Chai & Kwon Doppler lineshape
+    eta = ( 2.0 / (sqrt(PI * log(k0*L/2.0) ) * k0*L) / (2.0 * pow(Lq,2) + 2.0)
+            + 1.0 / ( sqrt( PI * log(k0*q0) ) * k0*2.0*q0) *
+            (Lq / (pow(Lq,2) + 1) + atan(Lq))); 
+    // Iordanova/Holstein    
+    // eta = 1.6 / (k0*R * pow((PI * log(k0*R)),2) ) ;
+    // Golubovskii et al. Lorentz lineshape 
+    // eta = (1.0 / ( sqrt(PI*k0*L/2) ) * (2.0/3.0 - 2.0 * pow(Lq,1.5) / (3.0 * pow((pow(Lq,2) + 1),0.75)) )
+            //  + 1.0/(2.0*sqrt(PI*k0*q0))* 4.0*Lq*sqrt(Lq) / (3.0 * pow((pow(Lq,2) + 1),0.75))); 
+  } else {
+    eta = 1.0;    
+  }    
+
+  eta = min(eta,1.0);  
+
+  return eta;
+}
\ No newline at end of file
diff --git a/src/reaction.hpp b/src/reaction.hpp
index a8507c2f2..cdffcbcc0 100644
--- a/src/reaction.hpp
+++ b/src/reaction.hpp
@@ -62,6 +62,12 @@ class Reaction {
     printf("computeRateCoefficient not implemented");
     return 0;
   }
+
+  MFEM_HOST_DEVICE virtual double computeRateCoefficient(const double *nsp, 
+                                                         const double &T_h, const double &T_e) {
+    printf("computeRateCoefficient not implemented");
+    return 0;    
+  }
 };
 
 class Arrhenius : public Reaction {
@@ -113,4 +119,193 @@ class Tabulated : public Reaction {
                                                          const bool isElectronInvolved = false);
 };
 
+
+
+ 
+class RadiativeDecay : public Reaction {
+ private:
+
+  const double big = 1.0e+50;
+  const double small = 1.0e-50;
+
+  // Length parameters used for the evaluation of the escape factors
+  const double R; // Radius of a cylinder
+  double L; // Length of a cylinder
+
+  std::string upper_sp_name;
+  std::string lower_sp_name;
+
+  int rank0_;
+
+  int iAr_u;
+  int iAr_l;
+
+  // Parameters
+  const double Mr_Ar = 0.039948;      // [kg/mol]
+  const double M_Ar = 6.63352088e-26; // Mr_Ar/spc.N_A # [kg] mass of argon atom (6.63352088e-26 kg)
+
+  // Variables
+  int NumOfInteral_lvl_u;
+  std::vector<double> *E_lvl_u;
+  std::vector<double> *g_lvl_u;
+  std::vector<double> n_sp_lvl_u;
+  std::vector<std::vector<double>> *Aji;
+
+  std::vector<double> effAcoef_lvl;
+
+  std::vector<double> *E_lvl_l;
+  std::vector<double> *g_lvl_l;
+  std::vector<double> n_sp_lvl_l;
+
+  // Data of internal levels of lumped species
+  // ground state
+  std::vector<double> E_lvl_g = {0.0}; 
+  std::vector<double> g_lvl_g = {1.0};
+
+  // 4s metastables
+  std::vector<double> E_lvl_m = {1114246.8116913952, 1131113.0237639823}; // [J/mol] // {11.54835442, 11.72316039}; // [eV]
+  std::vector<double> g_lvl_m = {5.0, 1.0};
+
+  // 4s resonacne
+  std::vector<double> E_lvl_r = {1121506.2040552883, 1141235.3742507447}; // [J/mol] // {11.62359272, 11.82807116}; // [eV]
+  std::vector<double> g_lvl_r = {3.0, 3.0};
+  std::vector<std::vector<double>>  Aji_r_g = {{132000000.0}, {532000000.0}};
+
+  // 4p states
+  std::vector<double> E_lvl_4p = {1245337.6579411437, 1280653.4893638478, 1261614.7730293325, 
+                                  1263463.1280640187, 1269085.454762629, 1270883.3460389085, 
+                                  1281579.837318737, 1283469.8354227678, 1285942.7139612488, 1300611.3568123293}; // [J/mol] 
+                                  // {12.90701530, 13.27303810, 13.07571571, 13.09487256, 13.15314387, 
+                                  // 13.17177770, 13.28263902, 13.30222747, 13.32785705, 13.47988682}; // [eV]
+  std::vector<double> g_lvl_4p  = {3.0, 1.0, 7.0, 5.0, 3.0, 5.0, 3.0, 5.0, 3.0, 1.0};
+  std::vector<std::vector<double>> Aji_4p_m = 
+    {
+    {18900000.0,  980000.0},
+    {33000000.0,  0.0},
+    {9300000.0,   0.0},
+    {5200000.0,   2430000.0},
+    {24500000.0,  0.0},
+    {0.0,         0.0},
+    {630000.0,    18600000.0},
+    {3800000.0,   0.0},
+    {6400000.0,   11700000.0},
+    {0.0,         0.0}    
+    };
+
+  std::vector<std::vector<double>> Aji_4p_r = 
+    {
+    {5400000.0,   190000.0},
+    {0.0,         0.0},
+    {21500000.0,  1470000.0},
+    {25000000.0,  1060000.0},
+    {4900000.0,   5000000.0},
+    {40000000.0,  8643.18384420115},
+    {22000.0,     13900000.0},
+    {8500000.0,   22300000.0},
+    {1830000.0,   15300000.0},
+    {236000.0,    45000000.0}    
+    };
+
+
+// Excited levels based on NIST database
+// Racah          Paschen   index  
+// Ar               1s1       n/a     Ground state
+
+// Ar(4s[3/2]2)     1s5       0       1st metastable
+// Ar(4s[3/2]1)     1s4       0       1st resonacne
+
+// Ar(4s'[1/2]0)    1s3       1       2nd metastable
+// Ar(4s'[1/2]1)    1s2       1       2nd resonacne
+
+// Ar(4p[1/2]1)     2p10	    0       4p states
+// Ar(4p[5/2]3)     2p9	      1          
+// Ar(4p[5/2]2)     2p8	      2
+// Ar(4p[3/2]1)     2p7	      3
+// Ar(4p[3/2]2)     2p6	      4
+// Ar(4p[1/2]0)     2p5	      5
+// Ar(4p'[3/2]1)    2p4	      6
+// Ar(4p'[3/2]2)    2p3	      7
+// Ar(4p'[1/2]1)    2p2       8
+// Ar(4p'[1/2]0)    2p1       9
+
+
+// Transitions based on NIST database
+// wavelength   upper level             lower level       A coef
+// 0            Ar(4s[3/2]2)      ->        Ar          0.0
+// 106.666      Ar(4s[3/2]1)      ->        Ar          132000000.0
+// 0            Ar(4s'[1/2]0)     ->        Ar          0.0
+// 104.822      Ar(4s'[1/2]1)     ->        Ar          532000000.0
+
+// 912.547      Ar(4p[1/2]1)      ->    Ar(4s[3/2]2)    18900000.0
+// 966.043      Ar(4p[1/2]1)      ->    Ar(4s[3/2]1)    5400000.0
+// 1047.292     Ar(4p[1/2]1)      ->    Ar(4s'[1/2]0)   980000.0
+// 1149.125     Ar(4p[1/2]1)      ->    Ar(4s'[1/2]1)   190000.0
+
+// 811.754      Ar(4p[5/2]3)      ->    Ar(4s[3/2]2)    33000000.0
+// 0            Ar(4p[5/2]3)      ->    Ar(4s[3/2]1)    0.0
+// 0            Ar(4p[5/2]3)      ->    Ar(4s'[1/2]0)   0.0
+// 0            Ar(4p[5/2]3)      ->    Ar(4s'[1/2]1)   0.0
+
+// 801.699      Ar(4p[5/2]2)      ->    Ar(4s[3/2]2)    9300000.0
+// 842.696      Ar(4p[5/2]2)      ->    Ar(4s[3/2]1)    21500000.0
+// 0            Ar(4p[5/2]2)      ->    Ar(4s'[1/2]0)   0.0
+// 978.719      Ar(4p[5/2]2)      ->    Ar(4s'[1/2]1)   1470000.0
+
+// 772.589      Ar(4p[3/2]1)      ->    Ar(4s[3/2]2)    5200000.0
+// 810.592      Ar(4p[3/2]1)      ->    Ar(4s[3/2]1)    25000000.0
+// 867.032      Ar(4p[3/2]1)      ->    Ar(4s'[1/2]0)   2430000.0
+// 935.678      Ar(4p[3/2]1)      ->    Ar(4s'[1/2]1)   1060000.0
+
+// 763.720      Ar(4p[3/2]2)      ->    Ar(4s[3/2]2)    24500000.0
+// 800.836      Ar(4p[3/2]2)      ->    Ar(4s[3/2]1)    4900000.0
+// 0            Ar(4p[3/2]2)      ->    Ar(4s'[1/2]0)   0.0
+// 922.703      Ar(4p[3/2]2)      ->    Ar(4s'[1/2]1)   5000000.0
+
+// 0            Ar(4p[1/2]0)      ->    Ar(4s[3/2]2)    0.0
+// 751.672      Ar(4p[1/2]0)      ->    Ar(4s[3/2]1)    40000000.0
+// 0            Ar(4p[1/2]0)      ->    Ar(4s'[1/2]0)   0.0
+// 858.042      Ar(4p[1/2]0)      ->    Ar(4s'[1/2]1)   8643.18384420115
+
+// 714.901      Ar(4p'[3/2]1)     ->    Ar(4s[3/2]2)    630000.0
+// 747.322      Ar(4p'[3/2]1)     ->    Ar(4s[3/2]1)    22000.0
+// 795.036      Ar(4p'[3/2]1)     ->    Ar(4s'[1/2]0)   18600000.0
+// 852.378      Ar(4p'[3/2]1)     ->    Ar(4s'[1/2]1)   13900000.0
+
+// 706.917      Ar(4p'[3/2]2)     ->    Ar(4s[3/2]2)    3800000.0
+// 738.601      Ar(4p'[3/2]2)     ->    Ar(4s[3/2]1)    8500000.0
+// 0            Ar(4p'[3/2]2)     ->    Ar(4s'[1/2]0)   0.0
+// 841.052      Ar(4p'[3/2]2)     ->    Ar(4s'[1/2]1)   22300000.0
+
+// 696.735      Ar(4p'[1/2]1)     ->    Ar(4s[3/2]2)    6400000.0
+// 727.494      Ar(4p'[1/2]1)     ->    Ar(4s[3/2]1)    1830000.0
+// 772.633      Ar(4p'[1/2]1)     ->    Ar(4s'[1/2]0)   11700000.0
+// 826.679      Ar(4p'[1/2]1)     ->    Ar(4s'[1/2]1)   15300000.0
+
+// 0            Ar(4p'[1/2]0)     ->    Ar(4s[3/2]2)    0.0
+// 667.912      Ar(4p'[1/2]0)     ->    Ar(4s[3/2]1)    236000.0
+// 0            Ar(4p'[1/2]0)     ->    Ar(4s'[1/2]0)   0.0
+// 750.59       Ar(4p'[1/2]0)     ->    Ar(4s'[1/2]1)   45000000.0
+
+ public:
+  MFEM_HOST_DEVICE RadiativeDecay(const double &_R, 
+                                  std::map<std::string, int> *_speciesMapping, std::vector<std::string> *_speciesNames,
+                                  const int _numSpecies ,const double *_reactantStoich,const double *_productStoich);
+
+  MFEM_HOST_DEVICE virtual ~RadiativeDecay();
+
+  MFEM_HOST_DEVICE virtual double computeRateCoefficient(const double *nsp, const double &T_h, const double &T_e); 
+
+  MFEM_HOST_DEVICE void GetEinsteinACoefficient(const double &T_h, const double &T_e, 
+                                                        double *effAcoef);
+
+  MFEM_HOST_DEVICE void GetNumDensityOfInteralLevels(const int NumOfInteral_lvl, const double &n_sp, 
+                                                     const double &T_e, double *n_sp_internal);
+
+  MFEM_HOST_DEVICE double escapeFactCalc(const double &n_i, const double &E_j, const double &E_i, 
+                                                 const double &g_j, const double &g_i, const double &A_ji, 
+                                                 const double &T_g);
+
+};
+
+
 #endif  // REACTION_HPP_
diff --git a/src/source_term.cpp b/src/source_term.cpp
index 69aa38bd8..8af4533c4 100644
--- a/src/source_term.cpp
+++ b/src/source_term.cpp
@@ -94,7 +94,6 @@ void SourceTerm::updateTerms(mfem::Vector &in) {
   TransportProperties *_transport = transport_; // Why do we do that? The pointers are already set in the class.
   Chemistry *_chemistry = chemistry_;
   Radiation *_radiation = radiation_;
-  const bool _enableRadiation = enableRadiation_;
 
   const int nnodes = vfes->GetNDofs();
   const int _dim = dim;
@@ -164,7 +163,7 @@ void SourceTerm::updateTerms(mfem::Vector &in) {
     double progressRates[gpudata::MAXREACTIONS], creationRates[gpudata::MAXSPECIES];
     if (_numSpecies > 1 && _numReactions > 0) {
       double kfwd[gpudata::MAXREACTIONS], kC[gpudata::MAXREACTIONS];
-      _chemistry->computeForwardRateCoeffs(Th, Te, kfwd);
+      _chemistry->computeForwardRateCoeffs(ns, Th, Te, kfwd);
       _chemistry->computeEquilibriumConstants(Th, Te, kC);
 
       // get reaction rates
@@ -203,17 +202,18 @@ void SourceTerm::updateTerms(mfem::Vector &in) {
 
     // TODO(kevin): energy sink for radiative reaction.
 
-    if (_enableRadiation) {
-      switch (radiation_->inputs.model) {
-        case NET_EMISSION:
-          srcTerm[1 + _nvel] += _radiation->computeEnergySink(Th);    
-          break;
-        case P1_MODEL:
-          srcTerm[1 + _nvel] += (*energySinkRad_)[n];
-          break;
-      }
-    }
 
+    // Radiation source term
+    switch (radiation_->inputs.model) {
+      case NET_EMISSION:
+        srcTerm[1 + _nvel] += _radiation->computeEnergySink(Th);    
+        break;
+      case P1_MODEL:
+        srcTerm[1 + _nvel] += (*energySinkRad_)[n];
+        break;
+      case NONE_RAD:
+        break;          
+    }
 
     if (_mixture->IsTwoTemperature()) {
       // energy sink from electron-impact reactions.
diff --git a/test/inputs/plasma.ini b/test/inputs/plasma.ini
index 27474a908..9a3a18d0c 100644
--- a/test/inputs/plasma.ini
+++ b/test/inputs/plasma.ini
@@ -31,7 +31,7 @@ order = 3
 integrationRule = 0
 basisType = 0
 refinement_levels = 0
-maxIters = 127600004 #551500010 #26500010  # Not used by cycle-avg-joule-coupled solver
+maxIters = 127600008 #551500010 #26500010  # Not used by cycle-avg-joule-coupled solver
 outputFreq = 2 #100000
 useRoe = 0
 enableSummationByParts = 0
@@ -226,7 +226,8 @@ mass = 5.48579908782496e-7 # [kg/mol]
 
 [species]
 numSpecies = 6
-background_index = 6
+background_index = 1
+#background_index = 2 # Which one should it be?
 # electron_index = 2
 
 [species/species1]
@@ -265,6 +266,15 @@ initialMassFraction = 0.0
 perfect_mixture/constant_molar_cv = 1.5
 perfect_mixture/constant_molar_cp = 2.5
 
+[species/species4]
+name = 'Ar_m'
+composition = '{Ar : 1}'
+formation_energy = 1116860.96186 # [J/mol]  =  11.5754 eV (weighted value such that n_m = the number densities of the two metastables) # (E_f = E_ev * spc.e * spc.N_A where E_f in [J/mol] and E_ev in [eV])   Exitation energy of metastables
+level_degeneracy = 6
+initialMassFraction = 0.0
+perfect_mixture/constant_molar_cv = 1.5
+perfect_mixture/constant_molar_cp = 2.5
+
 [species/species5]
 name = 'E'
 composition = '{E : 1}'
@@ -311,8 +321,8 @@ characteristic_length = 0.028
 [reactions/reaction1]
 equation = 'Ar + E <=> Ar.+1 + 2 E'
 reaction_energy = 1520571.3883 # [J/mol]
-reactant_stoichiometry = '0 0 0 0 1 1' 
-product_stoichiometry = '0 0 0 1 2 0'  
+reactant_stoichiometry = '1 1 0 0 0 0' 
+product_stoichiometry = '0 2 1 0 0 0'  
 model = arrhenius
 arrhenius/A = 74072.331348
 arrhenius/b = 1.511
@@ -321,32 +331,32 @@ arrhenius/E = 1176329.772504 # [J/mol]
 [reactions/reaction2]
 equation = 'Ar.+1 + 2 E <=> Ar + E'
 reaction_energy = -1520571.3883 # [J/mol] # Should it be negative?
-reactant_stoichiometry = '0 0 0 1 2 0' 
-product_stoichiometry = '0 0 0 0 1 1' 
+reactant_stoichiometry = '0 2 1 0 0 0' 
+product_stoichiometry = '1 1 0 0 0 0' 
 model = arrhenius
-arrhenius/A = 34126.4747526  
+arrhenius/A = 5.66683445516e-20
 arrhenius/b = 0.368
 arrhenius/E = -377725.908714 # [J/mol]
 
 [reactions/reaction3]
 equation = 'Ar_r <=> Ar'
 reaction_energy = -1130867.391486 # [J/mol] # Should it be negative?
-reactant_stoichiometry = '0 1 0 0 0 0' 
-product_stoichiometry = '0 0 0 0 0 1' 
+reactant_stoichiometry = '0 0 0 0 1 0' 
+product_stoichiometry = '1 0 0 0 0 0' 
 model = radiative_decay
 
 [reactions/reaction4]
 equation = 'Ar_p <=> Ar_r'
 reaction_energy = -139082.4944036866 # [J/mol]
-reactant_stoichiometry = '0 0 1 0 0 0' 
-product_stoichiometry = '0 1 0 0 0 0' 
+reactant_stoichiometry = '0 0 0 0 0 1' 
+product_stoichiometry = '0 0 0 0 1 0' 
 model = radiative_decay
 
 [reactions/reaction5]
 equation = 'Ar_p <=> Ar_m'
 reaction_energy = -153088.92402968672 # [J/mol]
-reactant_stoichiometry = '0 0 1 0 0 0' 
-product_stoichiometry = '1 0 0 0 0 0' 
+reactant_stoichiometry = '0 0 0 0 0 1' 
+product_stoichiometry = '0 0 0 1 0 0' 
 model = radiative_decay
 
 #---------------------------------------