A detailed discussion on the closed two-state model for reversible reactions

Abstract: A model for reversible reactions in condensed phases is solved. The model assumes population transfer between molecular potential surfaces (reactant and product). The population is represented by a probability distribution of the molecular configuration at a particular instant t. The probability distribution is governed by the time-dependent Smoluchowski equations. The population distribution between both states explicitly makes it a coupled PDE and has been solved in time-domain for the first time. Using the derived analytical results, a theory for electronic relaxation and reversible chemical reactions is developed.

Verifications of the analytical solution  by checking the boundary conditions and the original equation;

The verification on the obtained analytical solutions when both the state are initially populatedSOlubility and boundary conditions test when both the states are initially populated.

ClearAll["Global`*"];
Pe[x_,t_]=(Exp[I k (x)]+((dg kg (ce-cg))/(-k cg de -ce dg kg+2I k kg de dg ))Exp[I k (Abs[x])])(Exp[I k (xe)]Exp[-(σe^2+de  t) k^2]);Pg[x_,t_]=(Exp[I kg (x)]+((de k (-ce+cg))/(-k cg de -ce dg kg+2I k kg de dg ))Exp[I kg (Abs[x])])(Exp[I kg (xg)]Exp[-(σg^2+dg  t) kg^2]);(*(1/2)(Exp[I kg (xe)]Exp[-(σe^2+dg  t) kg^2]+Exp[I kg (xg)]Exp[-(σg^2+dg  t) kg^2]);*)
k=Sqrt[S/de]
kg=Sqrt[S/dg]
Pep[x_,t_]=Pe[x,t]/.Abs[x]->x;
Pem[x_,t_]=Pe[x,t]/.Abs[x]->-x;
Pgp[x_,t_]=Pg[x,t]/.Abs[x]->x;
Pgm[x_,t_]=Pg[x,t]/.Abs[x]->-x;
de=dg;
xe=xg;
σe=σg;
Print["Solubility check"]
Simplify[de D[Pep[x,t],{x,2}]-D[Pep[x,t],t]]
Simplify[dg D[Pgm[x,t],{x,2}]-D[Pgm[x,t],t]]
Print["Boundary check"]
Simplify[(D[Pgp[x,t],x]-D[Pgm[x,t],x])-cg/dg*Pg[x,t]+ce*Pe[x,t]/dg/.x->0]
Simplify[(D[Pep[x,t],x]-D[Pem[x,t],x])-ce/de*Pe[x,t]+cg/de*Pg[x,t]/.x->0]
Print["Initial check"]
Pe[x,0]
Pg[x,0]

SOlubility and boundary conditions test when only the excited state is initially populated;

ClearAll["Global`*"];
e1[x_]:=a Exp[I k x]+b Exp[-I k x];
e2[x_]:=c*Exp[I k x];
g1[x_]:=d*Exp[I kg x]+e Exp[-I kg x];
g2[x_]:=f Exp[I kg x];
a=1;
d=0;
Solve[{e1[0]==e2[0],g1[0]==g2[0],de(e2'[0]-e1'[0])==ce e1[0]-cg g1[0],dg(g2'[0]-g1'[0])==cg g1[0]-ce e1[0]},{b,c,e,f}]

ClearAll["Global`*"];
Pe[x_,t_]=(Exp[I k (x)]+((dg kg (ce))/(-k cg de -ce dg kg+2I k kg de dg ))Exp[I k (Abs[x])])(Exp[I k (xe)]Exp[-(σe^2+de  t) k^2]);Pg[x_,t_]=(0Exp[I kg (x)]+((de k (-ce))/(-k cg de -ce dg kg+2I k kg de dg ))Exp[I kg (Abs[x])])(Exp[I kg (xe)]Exp[-(σe^2+dg  t) kg^2]);(*(1/2)(Exp[I kg (xe)]Exp[-(σe^2+dg  t) kg^2]+Exp[I kg (xg)]Exp[-(σg^2+dg  t) kg^2]);*)
k=Sqrt[S/de]
kg=Sqrt[S/dg]
de=dg;
Pep[x_,t_]=Pe[x,t]/.Abs[x]->x;
Pem[x_,t_]=Pe[x,t]/.Abs[x]->-x;
Pgp[x_,t_]=Pg[x,t]/.Abs[x]->x;
Pgm[x_,t_]=Pg[x,t]/.Abs[x]->-x;
Print["Solubility check"]
Simplify[de D[Pep[x,t],{x,2}]-D[Pep[x,t],t]]
Simplify[dg D[Pgm[x,t],{x,2}]-D[Pgm[x,t],t]]
Print["Boundary check"]
Simplify[dg(D[Pgp[x,t],x]-D[Pgm[x,t],x])-cg*Pg[x,t]+ce*Pe[x,t]/.x->0]
Simplify[de(D[Pep[x,t],x]-D[Pem[x,t],x])-ce*Pe[x,t]+cg*Pg[x,t]/.x->0]
Print["Initial check"]
Pe[x,0]
Pg[x,0]

WHen both the states are initially populated, the model mimics a reversible chemical reaction

Plot when both the states are initially populated; every considerations are numerically verified; Some of the results related to the closed ness of the system
({{Considerations,                                                                        Observations,                                                                                                           Inferences}})
({{ce=cg, Conditions for the dynamic equilibrium irrespective of the other values, Illustrates dynamic equilibrium of reversible reactions}, {de= dg, ce>cg, System obeys closed ness, till when (∑_j)(Subscript[K, ij])=0 condition satisfies for the coupling matrix, Can represent a reversible chemical reaction}, {de=dg, ce<cg, system is always a closed system, till when (∑_j)(Subscript[K, ij])=0 condition satisfies for the coupling matrix, ,, ,,}})
Effect of de and dg (de,dg)=(1,2)
({{ce≠ cg&&(de≠dg
|| xe≠ xg
|| σe≠ σg), No closedness if more than one parameters change, violates the particle conservation. However, more explanation can be given only by
Experimenting more with the parameters }, {ce> cg&&(de>dg), Initially total population decay and than grow more than unity
for smaller time, excited state decay more than the ground state growth
for larger time, excited state decay less than the ground state growth, true}, {ce< cg&&(de<dg), Initially total population decay and than grow more than unity (quantitatively same as above)
for smaller time, ground state decay more than the excited state growth
for larger time,  ground state decay less than the excited state growth, true}})
({{ce> cg&&(de<dg), Initially total population grows and than decay less than unity(reduces as the same amount it increased)
for smaller time, ground state grow more and excited state decay less
for larger time, excited state decays more and ground state grow less, true}, {ce< cg&&(de>dg), Initially total population grows and than decay less than unity (quantitatively same as above)
for smaller time, ground state decay less and excited state grow more
for larger time,ground state decays more and excited state grow less, true}})
Effect of xe, and xg. (xe,xg)= (4,8)
({{ce> cg&&(|xe|>|xg|),  total population grows and than decay, always more than unity
For all time ground state grows more than the excited state. Rate is more in small time and less for large times, true}, {ce< cg&&(|xe|<|xg|), total population grows and than decay, always more than unity(quantitatively same as above)
For all time excited state grows more than the ground state decay. Rate is more in small time and less for large times, true}})
({{ce> cg&&((|xe|<|xg|), Initially total population decays and then growss. Always less than unity(reduces as the same amount it increased)
For all time excited state decay more than the ground state growth. Rate is more in small time and less for large times, true}, {ce< cg&&((|xe|>|xg|, Initially total population decays and then growss. Always less than unity(quantitatively same as above)
For all time ground state decay more than the excited state growth. Rate is more in small time and less for large times, true}})
Effect of σe, and σg. (σe, σg)=(4,0)
({{ce> cg&&(σe>σg),  total population decay,grow, decay, and then grow to unity. Always lesser than unity.
|g> grows, |e> decays., true}, {ce< cg&&(σe<σg),  total population decay,grow, decay, and then grow to unity. Always lesser than unity.
|g> decays, |e> grows., true}})
({{ce> cg&&(σe<σg), total population grow, decay, and then grow ,decay to unity. Always more than unity.
|g> grows, |e>decays., true}, {ce< cg&&(σe>σg), total population grow, decay, and then grow ,decay to unity. Always more than unity.
|g> decays, |e>grows., true}})

Let us generate graphs on the effect of the rates

ClearAll["Global`*"];
Qe[t_]=(1/2)+(((ce*dg-cg*dg))/((2*de*dg*2ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](Sign[(xe/(2(σe^2+de*t)))+ϵe]-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);Qg[t_]=(1/2)+(((cg*de-ce*de))/((2*de*dg*2ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](Sign[(xg/(2(σg^2+dg*t)))+ϵg]-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
de=2;
dg=1;
xe=4;
xg=4;
σe=0;
σg=0;
h=-1;
y=-1.3;
ce=1;
cg=2;
Plot[{Qe[t],Qg[t]},{t,0,100}];
ab1=Plot[{Qe[t]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)",13]},{Scaled[{0,0.5}],{y-.6,h-.7}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dashed,Thick}];
aa1=Plot[{Qg[t]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y-0.6,h}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted}];
ac1=Plot[{Qg[t]+Qe[t]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)+(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+0.05,h+.7}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Gray}];
Show[ac1,aa1,ab1]

The origin of the discontinuity

ClearAll["Global`*"];
Qe[t_]=(1/2)+(((ce*dg-cg*dg))/((2*de*dg*2ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);Qg[t_]=(1/2)+(((cg*de-ce*de))/((2*de*dg*2ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
Pe1[x_,t_]=(1/(4Sqrt[π(σe^2+de*t)]))*Exp[((-(x+xe)^2)/(4(σe^2+de*t)))]+(((ce*dg-cg*dg))/((8*de*dg)))*Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](-Sign[((xe+Abs[x])/(2(σe^2+de*t)))+ϵe]+Erf[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);
Pg1[x_,t_]=(1/(4Sqrt[π(σg^2+dg*t)]))*Exp[((-(x+xg)^2)/(4(σg^2+dg*t)))]+(((cg*de-ce*de))/((8*de*dg)))*Exp[ϵg^2(σg^2+dg*t)+ϵg(Abs[x]+xg)](-Sign[((xg+Abs[x])/(2(σg^2+dg*t)))+ϵg]+Erf[((Abs[x]+xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))]);
Pe[x_,t_]=(1/(4Sqrt[π(σe^2+de*t)]))*Exp[((-(x+xe)^2)/(4(σe^2+de*t)))]-(((ce*dg-cg*dg))/((8*de*dg)))*Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](Erfc[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);
Pg[x_,t_]=(1/(4Sqrt[π(σg^2+dg*t)]))*Exp[((-(x+xg)^2)/(4(σg^2+dg*t)))]-(((cg*de-ce*de))/((8*de*dg)))*Exp[ϵg^2(σg^2+dg*t)+ϵg(Abs[x]+xg)](Erfc[((Abs[x]+xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))]);
Pep[x_,t_]=Pe[x,t]/.Abs[x]->x;
Pem[x_,t_]=Pe[x,t]/.Abs[x]->-x;
Pgp[x_,t_]=Pg[x,t]/.Abs[x]->x;
Pgm[x_,t_]=Pg[x,t]/.Abs[x]->-x;
Simplify[de D[Pep[x,t],{x,2}]-D[Pep[x,t],t]]
Simplify[dg D[Pgm[x,t],{x,2}]-D[Pgm[x,t],t]]
Simplify[dg(D[Pgp[x,t],x]-D[Pgm[x,t],x])-cg*Pg[x,t]+ce*Pe[x,t]/.x->0]
Simplify[de(D[Pep[x,t],x]-D[Pem[x,t],x])-ce*Pe[x,t]+cg*Pg[x,t]/.x->0]
de=2;
dg=2;
xe=3;
xg=3;
σe=1/Sqrt[40];
σg=1/Sqrt[40];
h=-.4;
y=-1.5;
ce=2;
cg=1;
ab1=Plot[{Qe'[t],-ce*Pe[0,t]+cg*Pg[0,t]},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, e]) (t))/(𝕕t))",13],Style["-(Subscript[c, e])(Subscript[P, e])(0,t)+(Subscript[c, g])(Subscript[P, g])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, e])(0,t)",14]}]
aa1=Plot[{Qg'[t],-cg*Pg[0,t]+ce*Pe[0,t]},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, g]) (t))/(𝕕t))",13],Style["-(Subscript[c, g])(Subscript[P, g])(0,t)+(Subscript[c, e])(Subscript[P, e])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, g])(0,t)",14]}]
ac1=Plot[{Qg'[t]+Qe'[t],0},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["((𝕕 ((Subscript[Q, e]) (t) + (Subscript[Q, g]) (t)))/(𝕕t))",13],Style["0",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Transition flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, 1 -> 2])(x,t)",14]}]

Let us try to find the invariant in the closed two-state system

ClearAll["Global`*"];
Qe[t_]=(1/2)+(((ce*dg-cg*dg))/((2*de*dg*2ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);Qg[t_]=(1/2)+(((cg*de-ce*de))/((2*de*dg*2ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
Pe[x_,t_]=(1/(4Sqrt[π(σe^2+de*t)]))*Exp[((-(x+xe)^2)/(4(σe^2+de*t)))]+(((ce*dg-cg*dg))/((8*de*dg)))*Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](-Sign[((xe+Abs[x])/(2(σe^2+de*t)))+ϵe]+Erf[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);
Pg[x_,t_]=(1/(4Sqrt[π(σg^2+dg*t)]))*Exp[((-(x+xg)^2)/(4(σg^2+dg*t)))]+(((cg*de-ce*de))/((8*de*dg)))*Exp[ϵg^2(σg^2+dg*t)+ϵg(Abs[x]+xg)](-Sign[((xg+Abs[x])/(2(σg^2+dg*t)))+ϵg]+Erf[((Abs[x]+xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))]);
de=2;
dg=1;
xe=3;
xg=3;
σe=1/Sqrt[40];
σg=1/Sqrt[40];
h=-.4;
y=-1.5;
ce=2;
cg=1;
ab1=Plot[{Qe'[t],-ce*Pe[0,t]+cg*Pg[0,t]},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, e]) (t))/(𝕕t))",13],Style["-(Subscript[c, e])(Subscript[P, e])(0,t)+(Subscript[c, g])(Subscript[P, g])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, e])(0,t)",14]}];
aa1=Plot[{Qg'[t],-cg*Pg[0,t]+ce*Pe[0,t]},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, g]) (t))/(𝕕t))",13],Style["-(Subscript[c, g])(Subscript[P, g])(0,t)+(Subscript[c, e])(Subscript[P, e])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, g])(0,t)",14]}];
ac1=Plot[{Qg'[t]+Qe'[t],0},{t,0,100},PlotStyle->{Dashed,{Dotted,Red}},PlotLegends->{Placed[{Style["((𝕕 ((Subscript[Q, e]) (t) + (Subscript[Q, g]) (t)))/(𝕕t))",13],Style["0",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Transition flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, 1 -> 2])(x,t)",14]}]

Plot when only the excited state is populated; every considerations are numerically verified; Some of the results related to the closed ness of the system
({{Considerations,                                                                        Observations,                                                                                                           Inferences}})
({{ce=cg, System is closed. one-directional reaction because of the population imbalance, Illustrates dynamic equilibrium}, {de= dg, ce>cg, System obeys closed ness, till when (∑_j)(Subscript[K, ij])=0 condition satisfies for the coupling matrix, Can represent a reversible chemical reaction}, {de=dg, ce<cg, system is always a closed system, till when (∑_j)(Subscript[K, ij])=0 condition satisfies for the coupling matrix, ,, ,,}})
Effect of de and dg (de,dg)=(1,2)
({{ce≠ cg&&(de≠dg
|| xe≠ xg
|| σe≠ σg), No closedness if more than one parameters change, violates the particle conservation. However, more explanation can be given only by
Experimenting more with the parameters }, {ce> cg&&(de>dg), Initially total population decay and than grow more than unity
for smaller time, excited state decay more than the ground state growth
for larger time, excited state decay less than the ground state growth, true}, {ce< cg&&(de<dg), Initially total population decay and than grow more than unity (quantitatively same as above)
for smaller time, ground state decay more than the excited state growth
for larger time,  ground state decay less than the excited state growth, true}})
({{ce> cg&&(de<dg), Initially total population grows and than decay less than unity(reduces as the same amount it increased)
for smaller time, ground state grow more and excited state decay less
for larger time, excited state decays more and ground state grow less, true}, {ce< cg&&(de>dg), Initially total population grows and than decay less than unity (quantitatively same as above)
for smaller time, ground state decay less and excited state grow more
for larger time,ground state decays more and excited state grow less, true}})
Effect of x0
({{ce> cg&&(|xe|>|xg|),  total population grows and than decay, always more than unity
For all time ground state grows more than the excited state. Rate is more in small time and less for large times, true}, {ce< cg&&(|xe|<|xg|), total population grows and than decay, always more than unity(quantitatively same as above)
For all time excited state grows more than the ground state decay. Rate is more in small time and less for large times, true}})
({{ce> cg&&((|xe|<|xg|), Initially total population decays and then growss. Always less than unity(reduces as the same amount it increased)
For all time excited state decay more than the ground state growth. Rate is more in small time and less for large times, true}, {ce< cg&&((|xe|>|xg|, Initially total population decays and then growss. Always less than unity(quantitatively same as above)
For all time ground state decay more than the excited state growth. Rate is more in small time and less for large times, true}})
Effect of σ
({{ce> cg&&(σe>σg),  total population decay,grow, decay, and then grow to unity. Always lesser than unity.
|g> grows, |e> decays., true}, {ce< cg&&(σe<σg),  total population decay,grow, decay, and then grow to unity. Always lesser than unity.
|g> decays, |e> grows., true}})
({{ce> cg&&(σe<σg), total population grow, decay, and then grow ,decay to unity. Always more than unity.
|g> grows, |e>decays., true}, {ce< cg&&(σe>σg), total population grow, decay, and then grow ,decay to unity. Always more than unity.
|g> decays, |e>grows., true}})

Let us generate graphs on the effect of the diffusivity

ClearAll["Global`*"];
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
Qe[t_,de_,dg_,ce_,cg_]=1+(((ce*dg))/((2*de*dg*ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);
Qg[t_,de_,dg_,ce_,cg_]=(((-ce*de))/((2*de*dg*ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
xe=4;
xg=4;
σe=0;
σg=0;
h=-.4;
y=-1.5;
a1=Plot[{Qe[t,.5,.5,2,2]},{t,0,100},PlotRange->{-0.02,1.01},PlotLegends->Placed[{Style["D=0.5 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed,Black}];
a2=Plot[{Qe[t,1,1,2,2]},{t,0,100},PlotLegends->Placed[{Style["D=1 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Thick,Black}];
a3=Plot[{Qe[t,2,2,2,2]},{t,0,100},PlotLegends->Placed[{Style["D=2 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
Show[a1,a2,a3]
b1=Plot[{Qe[t,1,1,2,2]},{t,0,100},PlotRange->{-0.02,1.01},PlotLegends->Placed[{Style["(Subscript[D, g])-(Subscript[D, e])=0 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dashed,Black}];
ab2=Plot[{Qe[t,1,1.5,2,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[D, g])-(Subscript[D, e])=0.5  (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Thick,Black}];
ab3=Plot[{Qe[t,1,2,2,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[D, g])-(Subscript[D, e])=1  (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
ac1=Plot[{Qg[t,2,1,1,1]+Qe[t,2,1,1,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)+(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h-2}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Gray}];
Show[ac1,ab2,ab3,b1]
Clear[σe,σg,xe,xg];
xe=3;
xg=3;
σe=1/Sqrt[40];
σg=1/Sqrt[40];
aa1=Plot[{Qg[t,2,2,2,3]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
ab1=Plot[{Qe[t,2,2,2,3]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dashed,Black}];
ac1=Plot[{Qg[t,2,2,2,3]+Qe[t,2,2,2,3]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)+(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+0.5,h+1}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Gray,Thick}];
Show[ac1,aa1,ab1]

Effect of x0

ClearAll["Global`*"];
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
xe=x0;
xg=x0;
Qe[t_,x0_]=1+(((ce*dg))/((2*de*dg*ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);
Qg[t_,x0_]=(((-ce*de))/((2*de*dg*ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
σe=0;
σg=0;
ce=2;
cg=2;
de=1;
dg=1;
h=-.4;
y=-1.5;
ab1=Plot[{Qe[t,1]},{t,0,100},PlotRange->{-0.02,1.01},PlotLegends->Placed[{Style["(Subscript[x, 0])=1 nm",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dashed,Black}];
ab2=Plot[{Qe[t,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[x, 0])=2 nm",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Thick,Black}];
ab3=Plot[{Qe[t,5]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[x, 0])=5 nm",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->{-0.01,1.01},PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted,Black}];
Show[ab3,ab2,ab1]
aa1=Plot[{Qg[t,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted}];
ac1=Plot[{Qg[t,2]+Qe[t,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])[t]+Qe[t]",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->{-0.01,1.01},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["Q(t)",14]},PlotStyle->{Red}];
Show[ac1,aa1,ab1]

Effect of σ

ClearAll["Global`*"];
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
xe=4;
xg=4;
σe=σ;
σg=σ;
Qe[t_,σ_]=1+(((ce*dg))/((2*de*dg*ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);
Qg[t_,σ_]=(((-ce*de))/((2*de*dg*ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
ce=2;
cg=2;
de=1;
dg=1;
h=-.4;
y=-1.5;
ab1=Plot[{Qe[t,1]},{t,0,100},PlotRange->{-0.02,1.01},PlotLegends->Placed[{Style["σ=1 nm",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed,Black}];
ab2=Plot[{Qe[t,2]},{t,0,100},PlotLegends->Placed[{Style["σ=2 nm",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Thick,Black}];
ab3=Plot[{Qe[t,7]},{t,0,100},PlotLegends->Placed[{Style["σ=7 nm",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
Show[ab1,ab2,ab3]
aa1=Plot[{Qg[t,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted}];
ac1=Plot[{Qg[t,2]+Qe[t,2]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])[t]+Qe[t]",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Red}];
Show[ac1,aa1,ab1]

Let us generate graphs on the effect of the rates

ClearAll["Global`*"];
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
Qe[t_,de_,dg_,ce_,cg_]=1+(((ce*dg))/((2*de*dg*ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);
Qg[t_,de_,dg_,ce_,cg_]=(((-ce*de))/((2*de*dg*ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
xe=4;
xg=4;
σe=0;
σg=0;
h=-.4;
y=-1.5;
ab1=Plot[{Qe[t,.5,.5,.1,.1]},{t,0,100},PlotRange->{-0.02,1.01},PlotLegends->Placed[{Style["(Subscript[k, 0])= 0.1 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed,Black}];
ab2=Plot[{Qe[t,.5,0.5,1,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])= 1 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Thick,Black}];
ab3=Plot[{Qe[t,.5,0.5,7,7]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])= 7 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
Show[ab1,ab2,ab3]
a1=Plot[{Qe[t,.5,.5,1,1]},{t,0,100},PlotRange->{-0.02,1.05},PlotLegends->Placed[{Style["(Subscript[c, e])-(Subscript[c, g])= 0 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h-1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Population probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed,Black}];
a2=Plot[{Qe[t,.5,0.5,1.5,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[c, e])-(Subscript[c, g])= 0.5 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Thick,Black}];
a3=Plot[{Qe[t,.5,0.5,3,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[c, e])-(Subscript[c, g])= 2 nm/ps",13]},{Scaled[{0,0.5}],{y+.1,h+1}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted,Black}];
ac1=Plot[{Qg[t,.5,0.5,3,1]+Qe[t,.5,0.5,3,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, e])(t)+(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h+2}}],PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Gray}];
Show[a1,a2,a3,ac1]
aa1=Plot[{Qg[t,.5,0.5,1,1]},{t,0,100},PlotLegends->Placed[{Style["(Subscript[Q, g])(t)",13]},{Scaled[{0,0.5}],{y+.1,h+2}}],PlotRange->All,PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]},PlotStyle->{Dotted}];
Show[ac1,aa1,ab1]

The origin of the discontinuity

ClearAll["Global`*"];
Qe[t_]=1+(((ce*dg))/((2*de*dg*ϵe)))*(Exp[ϵe*xe+ϵe^2(σe^2+de*t)](1-Erf[((xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+2HeavisideTheta[-xe-2ϵe(σe^2+de*t)]*Exp[-ϵe^2(σe^2+de*t)]-Erfc[xe/Sqrt[4(σe^2+de*t)]]);Qg[t_]=(((-ce*de))/((2*de*dg*ϵg)))*(Exp[ϵg*xg+ϵg^2(σg^2+dg*t)](1-Erf[((xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))])+2HeavisideTheta[-xg-2ϵg(σg^2+dg*t)]*Exp[-ϵg^2(σg^2+dg*t)]-Erfc[xg/Sqrt[4(σg^2+dg*t)]]);
ϵe=((Sqrt[dg*de]*cg+ce*dg)/(2*de*dg));
ϵg=((Sqrt[dg*de]*ce+cg*de)/(2*de*dg));
Pe1[x_,t_]=(1/(2Sqrt[π(σe^2+de*t)]))*Exp[((-(x+xe)^2)/(4(σe^2+de*t)))]+(((ce*dg))/((4*de*dg)))*Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](-Sign[((xe+Abs[x])/(2(σe^2+de*t)))+ϵe]+Erf[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);
Pg1[x_,t_]=(((-ce*de))/((4*de*dg)))*Exp[ϵg^2(σg^2+dg*t)+ϵg(Abs[x]+xg)](-Sign[((xg+Abs[x])/(2(σg^2+dg*t)))+ϵg]+Erf[((Abs[x]+xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))]);
Pe[x_,t_]=(1/(2Sqrt[π(σe^2+de*t)]))*Exp[((-(x+xe)^2)/(4(σe^2+de*t)))]+(((ce*dg))/((4*de*dg)))*Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](-1+Erf[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);
Pg[x_,t_]=(((-ce*de))/((4*de*dg)))*Exp[ϵg^2(σg^2+dg*t)+ϵg(Abs[x]+xg)](-1+Erf[((Abs[x]+xg+2ϵg(σg^2+dg*t))/(2Sqrt[(σg^2+dg*t)]))]);
Pgplus[x_,t_]=Pg[x,t]/.{Abs[x]->x,Sign[((xg+Abs[x])/(2(σg^2+dg*t)))+ϵg]->1};
Pgminus[x_,t_]=Pg[x,t]/.{Abs[x]->-x,Sign[((xg+Abs[x])/(2(σg^2+dg*t)))+ϵg]->1};
Peplus[x_,t_]=Pe[x,t]/.{Abs[x]->x,Sign[((xe+Abs[x])/(2(σe^2+de*t)))+ϵe]->1};
Peminus[x_,t_]=Pe[x,t]/.{Abs[x]->-x,Sign[((xe+Abs[x])/(2(σe^2+de*t)))+ϵe]->1};
Print["Boundary check"]
σe=σg;
xe=xg;
$Assumptions=de>0&&dg >0;
b=Simplify[dg(D[Pgplus[x,t],x]-D[Pgminus[x,t],x])-cg*Pg[x,t]+ce*Pe[x,t]/.x->0]
c=Simplify[de(D[Peplus[x,t],x]-D[Peminus[x,t],x])-ce*Pe[x,t]+cg*Pg[x,t]/.x->0]
de=2;
dg=2;
xe=3;
xg=3;
σe=1/Sqrt[40];
σg=1/Sqrt[40];
ce=2;
cg=2;
h=-.4;
y=-1.5;
a1=Plot[{b,c},{t,0,100},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, e]) (t))/(𝕕t))",13],Style["-(Subscript[c, e])(Subscript[P, e])(0,t)+(Subscript[c, g])(Subscript[P, g])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotStyle->{{Dashed,Black},{Dotted,Red}},PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, e])(0,t)",14]}]
ab1=Plot[{Qe'[t],-ce*Pe[0,t]+cg*Pg[0,t]},{t,0,100},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, e]) (t))/(𝕕t))",13],Style["-(Subscript[c, e])(Subscript[P, e])(0,t)+(Subscript[c, g])(Subscript[P, g])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotStyle->{{Dashed,Black},{Dotted,Red}},PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, e])(0,t)",14]}]
aa1=Plot[{Qg'[t],-cg*Pg[0,t]+ce*Pe[0,t]},{t,0,100},PlotStyle->{{Dashed,Black},{Dotted,Red}},PlotLegends->{Placed[{Style["(((Subscript[𝕕Q, g]) (t))/(𝕕t))",13],Style["-(Subscript[c, g])(Subscript[P, g])(0,t)+(Subscript[c, e])(Subscript[P, e])(0,t)",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["population flux vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[j, g])(0,t)",14]}]
ac1=Plot[{Qg'[t]+Qe'[t],0},{t,0,100},PlotStyle->{Dashed,Dotted},PlotLegends->{Placed[{Style["((𝕕 ((Subscript[Q, e]) (t) + (Subscript[Q, g]) (t)))/(𝕕t))",13],Style["0",13]},{Scaled[{0,0.5}],{y+.1,h-1},Scaled[{0,0.5}],{y+.1,h}}]},PlotRange->All,PlotTheme->{"Scientific"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14]}]