diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index c8d34cf1..bc9726d0 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -17,7 +17,29 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     forAll(Y, i)
     {
-        if (Y[i].name() != inertSpecie)
+        if (Y[i].name() == electronSpecie)
+        {
+            fvScalarMatrix neEqn
+                (
+                 fvm::ddt(ne)
+                 + mvConvection->fvmDiv(ve, ne)
+                 - mvConvection->fvmDiv(fvc::interpolate(De/Te*fvc::grad(Te)) & mesh.Sf(), ne)
+                 - fvm::laplacian(De, ne)
+                 ==
+                 fvOptions(ne)
+                );
+
+            neEqn.relax();
+
+            fvOptions.constrain(neEqn);
+
+            neEqn.solve(mesh.solver("ne"));
+
+            fvOptions.correct(ne);
+
+            ne.max(0.0);
+        }
+        else if (Y[i].name() != inertSpecie)
         {
             volScalarField& Yi = Y[i];
 
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 50a4bb1d..e8466faf 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -55,10 +55,15 @@ basicMultiComponentMixture& composition = thermo.composition();
 
 //- Elementary charge (default in [C])
 const dimensionedScalar eCharge = constant::electromagnetic::e;
+//- Avogadro number (default in [1/mol])
+const dimensionedScalar NA = constant::physicoChemical::NA;
+//- Universal gas constant (default in [J/mol/K])
+const dimensionedScalar R = constant::physicoChemical::R;
 
 PtrList<volScalarField>& Y = composition.Y();
 
 word inertSpecie(thermo.lookup("inertSpecie"));
+word electronSpecie("E-");
 
 volScalarField rho
 (
@@ -164,19 +169,7 @@ volScalarField ne
 );
 
 Info<< "Creating field gas number density\n" << endl;
-volScalarField ng
-(
-    IOobject
-    (
-        "ng",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh,
-    dimensionedScalar("ng", ne.dimensions(), 1e-25)
-);
+volScalarField ng("ng", p / R / T * NA);
 
 Info<< "Creating field reduced electric field\n" << endl;
 volScalarField En ("En", mag(E) / (ne+ng));
diff --git a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
index 5efbb114..0d766c36 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
+++ b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
@@ -67,6 +67,7 @@ int main(int argc, char *argv[])
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
         #include "PhiEqn.H"
+        En = mag(E) / (ne+ng);
         #include "rhoEqn.H"
 
         while (pimple.loop())