This Mathematica notebook presents some of the results and discussions of the paper,
 R. Saravanan, A. Chakraborty, Chemical Physics Letters 731 (2019): 136567 (Click).

The survival probability for the excited state and ground state,

ClearAll["Global`*"];
ϵe=((ce cg)/(4 de dg))Sqrt[(dg/de)];ϵe=((ce cg)/(4 de dg))Sqrt[(de/dg)];Qe[t_]=1+Exp[-ϵe*xe+ϵe^2(σe^2+de*t)]/2(Sign[xe/(2(σe^2+de*t))-ϵe]-Erf[(xe-ϵe*2*(σe^2+de*t))/Sqrt[4(σe^2+de*t)]])-Erfc[xe/Sqrt[4(σe^2+de*t)]]+Exp[ϵe*xe+ϵe^2(σe^2+de*t)]/2(Sign[xe/(2(σe^2+de*t))+ϵe]-Erf[(xe+ϵe*2*(σe^2+de*t))/Sqrt[4(σe^2+de*t)]])-Exp[-ϵe^2(σe^2+de*t)]*HeavisideTheta[-ϵe*2*(σe^2+de*t)-xe]+Exp[-ϵe^2(σe^2+de*t)]*HeavisideTheta[ϵe*2*(σe^2+de*t)-xe];
Qg[t_]={1+Exp[-ϵg*xg+ϵg^2(σg^2+dg*t)]/2(Sign[xg/(2(σg^2+dg*t))-ϵg]-Erf[(xg-ϵg*2*(σg^2+dg*t))/Sqrt[4(σg^2+dg*t)]])-Erfc[xg/Sqrt[4(σg^2+dg*t)]]+Exp[ϵg*xg+ϵg^2(σg^2+dg*t)]/2(Sign[xg/(2(σg^2+dg*t))+ϵg]-Erf[(xg+ϵg*2*(σg^2+dg*t))/Sqrt[4(σg^2+dg*t)]])-Exp[-ϵg^2(σg^2+dg*t)]*HeavisideTheta[-ϵg*2*(σg^2+dg*t)-xg]+Exp[-ϵg^2(σg^2+dg*t)]*HeavisideTheta[ϵg*2*(σg^2+dg*t)-xg];

The above analytical results are benchmarked with numerical solutions and some of the following are theobservations of the two-state open systems;

observation 1: Pg(t) is zero and Pe(t) is non zero analytically when only excited states are populated. THis says that the model does not account ground state being populated when excited state is decaying. This observation matches and both the Pg(t) and Pe(t) matches quantitatively with numerical.

Observation 2: introducing Pg(x,0) doesn’t affect the Pe(t). Analytically this is seen. numerically also this fact is verified.

Observation 3: making either of the sink coefficients is suffice enough to stop the decay. Analytically this deduced for both cases (B1=0 and B1=A1). Numerically the fact is true when B1=and for  B1=A1, 

Supporting statement 4: When the states are equally populated the quantitative/qualitative expression are same when the parameters of |e> and |g> are same/different.numerically it has been verified.

Protect[aa,ab]
Remove["Global`*"]
ClearAll["Global'*"]

ClearAll["Global'*"];
Remove["Global'*"];
Unprotect["Global'*"];
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
xe=3;
xg=4;
σe=1/Sqrt[40];
σg=1/Sqrt[35];
de=0.5;
dg=0.6;
ke=3;
kg=2;
dThick=40;
dTime=100;
eqSimplee=D[Pe[x,t],t]==de*D[Pe[x,t],{x,2}]-ke*S[x]*Pg[x,t];
eqSimpleg=D[Pg[x,t],t]==dg*D[Pg[x,t],{x,2}]-kg*S[x]*Pe[x,t];
βe=1/(4σe^2);
βg=1/(4σg^2);
eqInie=Pe[x,0]==√(βe/π)*Exp[-βe*(x+xe)^2];
eqInig=Pg[x,0]==0;(*√(βg/π)*Exp[-βg*(x+xg)^2];*)
eqBlefte=Pe[-dThick,t]==0;
eqBleftg=Pg[-dThick,t]==0; 
eqBrighte=Pe[dThick,t]==0;
eqBrightg=Pg[dThick,t]==0;
h=-.4;
y=-1.5;
ys=0;
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqSimpleg,eqInie,eqInig,eqBrightg,eqBlefte,eqBleftg,eqBrighte},{Pe,Pg},{x,-dThick,dThick},{t,0,dTime}]
p1=Plot[Evaluate[Pe[x,12]/.res],{x,-dThick,dThick},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["probab distribution",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dashed}]
ϵe=Sqrt[((ke*kg)/(4*de*dg))Sqrt[dg/de]];
Panal[x_,t_]=(1/(2Sqrt[π(σe^2+de*t)]))Exp[-(((x+xe)^2)/(4(σe^2+de*t)))]-((ϵe)/4)Exp[ϵe^2(σe^2+de*t)-ϵe(Abs[x]+xe)](-Sign[((Abs[x]+xe)/(2(σe^2+de*t)))-ϵe]+Erf[((Abs[x]+xe-2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+((ϵe)/4)Exp[ϵe^2(σe^2+de*t)+ϵe(Abs[x]+xe)](-Sign[((Abs[x]+xe)/(2(σe^2+de*t)))+ϵe]+Erf[((Abs[x]+xe+2ϵe(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))]);+(ke/(2de))*Sign[x](-((Exp[ϵe^2(σe^2+de*t)])/4)Sum[((Exp[(2j-1)ϵe(Abs[x]+xe)])/(ϵe))(-Sign[((Abs[x]+xe)/(2(σe^2+de*t)))+(2j-1)ϵe]+Erf[((Abs[x]+xe+(2j-1)ϵe*2(σe^2+de*t))/(2Sqrt[(σe^2+de*t)]))])+(2j-1)Exp[(2j-1)*ϵe*2(σe^2+de*t)]*S[Abs[x]+xe+(2j-1)*ϵe*2(σe^2+de*t)],{j,0,1}]);
p2=Plot[Panal[x,12],{x,-dThick,dThick},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["probab distribution",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dashed}]
Show[p1,p2]
(*aa1=Plot[NIntegrate[Evaluate[Pg[x,t]/.res],{x,-dThick,dThick}],{t,0,dTime},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dashed}]
Plot[Evaluate[Pe[x,5]/.res],{x,-dThick,dThick},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["probab distribution",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dashed}]
ab2=Plot[NIntegrate[Evaluate[Pe[x,t]/.res],{x,-dThick,dThick}],{t,0,100},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
ϵe=Sqrt[((ke*kg)/(4*de*dg))Sqrt[dg/de]];
P1e[t_]=1+Exp[-ϵe*xe+ϵe^2(σe^2+de*t)]/2(Sign[xe/(2(σe^2+de*t))-ϵe]-Erf[(xe-ϵe*2*(σe^2+de*t))/Sqrt[4(σe^2+de*t)]])-Erfc[xe/Sqrt[4(σe^2+de*t)]]+Exp[ϵe*xe+ϵe^2(σe^2+de*t)]/2(Sign[xe/(2(σe^2+de*t))+ϵe]-Erf[(xe+ϵe*2*(σe^2+de*t))/Sqrt[4(σe^2+de*t)]])-Exp[-ϵe^2(σe^2+de*t)]*HeavisideTheta[-ϵe*2*(σe^2+de*t)-xe]+Exp[-ϵe^2(σe^2+de*t)]*HeavisideTheta[ϵe*2*(σe^2+de*t)-xe];
ϵg=Sqrt[((ke*kg)/(4*de*dg))Sqrt[de/dg]];
P2g[t_]=0;(*1+Exp[-ϵg*xg+ϵg^2(σg^2+dg*t)]/2(Sign[xg/(2(σg^2+dg*t))-ϵg]-Erf[(xg-ϵg*2*(σg^2+dg*t))/Sqrt[4(σg^2+dg*t)]])-Erfc[xg/Sqrt[4(σg^2+dg*t)]]+Exp[ϵg*xg+ϵg^2(σg^2+dg*t)]/2(Sign[xg/(2(σg^2+dg*t))+ϵg]-Erf[(xg+ϵg*2*(σg^2+dg*t))/Sqrt[4(σg^2+dg*t)]])-Exp[-ϵg^2(σg^2+dg*t)]*HeavisideTheta[-ϵg*2*(σg^2+dg*t)-xg]+Exp[-ϵg^2(σg^2+dg*t)]*HeavisideTheta[ϵg*2*(σg^2+dg*t)-xg];*)p1=Plot[P1e[t],{t,0,100},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[P2g[t],{t,0,100},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dotted}]
Show[ab1,p1]
Show[aa1,p2]
Show[ab1,ab2]*)

ϵe

Protect[ab1]

Show[ab1,p1]
Show[aa1,p2]

Smoluchowski Dynamics for a free particle

ClearAll["Global'*"]
Remove["Global`*"]
Unprotect["Global`*"]

ClearAll["Global'*"];
S[z_]=Sqrt[1000000/π]Exp[-1000000* (z)^2];
k0=5;
x0=3;
σ=1/Sqrt[40];
dDSim=0.5;
dThick=40;
dTime=100;
eqSimplee=D[P1[x,t],t]==dDSim*D[P1[x,t],{x,2}]-k0*S[x]*P2[x,t];
eqSimpleg=D[P2[x,t],t]==dDSim*D[P2[x,t],{x,2}]-k0*S[x]*P1[x,t];
β=1/(4σ^2);
eqInie=P1[x,0]==√(β/π)*Exp[-β*(x+x0)^2];
ϵ=0.001;
eqInig=P2[x,0]==0;
eqBlefte=P1[-dThick,t]==0;
eqBleftg=P2[-dThick,t]==0;
eqBrighte=P1[dThick,t]==0;
eqBrightg=P2[dThick,t]==0;
h=-.4;
y=-1.5;
ys=0;
mol[n_,o_]:={"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","MaxPoints"->n,"MinPoints"->n,"DifferenceOrder"->n}};
res=NDSolve[{eqSimplee,eqSimpleg,eqInie,eqInig,eqBrightg,eqBlefte,eqBleftg,eqBrighte},{P1,P2},{x,-dThick,dThick},{t,0,dTime}]
aa=Plot[NIntegrate[Evaluate[P2[x,t]/.res],{x,-dThick,dThick}],{t,0,dTime},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dashed}]
ab=Plot[NIntegrate[Evaluate[P1[x,t]/.res],{x,-dThick,dThick}],{t,0,100},PlotLegends->Placed[{Style["Simulation",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]

Protect[aa,ab];
ClearAll["Global'*"];
Remove["Global`*"];
a=ϵ0;
P1[t_]=1+Exp[-a*x0+a^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))-a]-Erf[(x0-a*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Erfc[x0/Sqrt[4(σ^2+d*t)]]+Exp[a*x0+a^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))+a]-Erf[(x0+a*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Exp[-a^2(σ^2+d*t)]*HeavisideTheta[-a*2*(σ^2+d*t)-x0]+Exp[-a^2(σ^2+d*t)]*HeavisideTheta[a*2*(σ^2+d*t)-x0];
P2[t_]=0;
k0=5;
x0=3;
σ=1/Sqrt[40];
ϵ0=k0/(2d);
a=ϵ0;
d=.5;
h=-.4;
y=-1.5;
ys=0;
p1=Plot[P1[t],{t,0,100},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[P2[t],{t,0,100},PlotLegends->Placed[{Style["Analytical",13]},{Scaled[{0,0.5}],{y-4.5ys,h-1}}],PlotRange->{-0.07,1.03},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, g])(t)",14]},PlotStyle->{Dotted}]
p3=Show[aa,p2]
Export["ground.eps",p3,ImageResolution->1500]
p4=Show[ab,p1]
Export["excited.eps",p4,ImageResolution->1500]

ClearAll["Global`*"]

ClearAll["Global"];
ϵ0=k0/(2d);
Pe[t_,k0_]=1+Exp[-ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))-ϵ0]-Erf[(x0-ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Erfc[x0/Sqrt[4(σ^2+d*t)]]+Exp[ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))+ϵ0]-Erf[(x0+ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[-ϵ0*2*(σ^2+d*t)-x0]+Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[ϵ0*2*(σ^2+d*t)-x0];
P2[t_]=0;
σ=0.5;
d=1;
x0=2;
h=1;
y=-1.3;
ys=0;
p1=Plot[Pe[t,0],{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=0 nm/ps",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[Pe[t,.1],{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=0.1 nm/ps",13]},{Scaled[{0,0.5}],{y+0.1,h-1}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Thick}]
p3=Plot[Pe[t,1],{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=1 nm/ps",13]},{Scaled[{0,0.5}],{y-4.5ys,h-0}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
p4=Plot[Pe[t,10],{t,0,100},PlotLegends->Placed[{Style["(Subscript[k, 0])=10 nm/ps",13]},{Scaled[{0,0.5}],{y
+.1,h+1}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Gray}]
p5=Show[p1,p2,p3,p4]
Export["k0effect.eps",p5,ImageResolution->1500,"CompressionLevel"->0]

ClearAll["Global"]
Remove["Global`*"]
Unprotect["Global`*"]

ClearAll["Global"];
ϵ0=k0/(2d);
Pe[t_,d_]=1+Exp[-ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))-ϵ0]-Erf[(x0-ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Erfc[x0/Sqrt[4(σ^2+d*t)]]+Exp[ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))+ϵ0]-Erf[(x0+ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[-ϵ0*2*(σ^2+d*t)-x0]+Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[ϵ0*2*(σ^2+d*t)-x0];
P2[t_]=0;
σ=0;
k0=10;
x0=5;
h=.2;
y=-1.4;
ys=0;p1=Plot[Pe[t,0.5],{t,0,100},PlotLegends->Placed[{Style["D=0.5 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y+0.2,h-2}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[Pe[t,1],{t,0,100},PlotLegends->Placed[{Style["D=1 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y,h-1}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Thick}]
p3=Plot[Pe[t,2],{t,0,100},PlotLegends->Placed[{Style["D=2 (nm^2)/ps",13]},{Scaled[{0,0.5}],{y-4.5ys,h-0}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
p5=Show[p1,p2,p3]
Export["deffecttwo.eps",p5,ImageResolution->1500,"CompressionLevel"->0]
SystemOpen[%]

Show[,p9]

ClearAll["Global"]
Remove["Global`*"]
Unprotect["Global`*"]

ClearAll["Global"];
ϵ0=k0/(2d);
Pe[t_,σ_]=1+Exp[-ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))-ϵ0]-Erf[(x0-ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Erfc[x0/Sqrt[4(σ^2+d*t)]]+Exp[ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))+ϵ0]-Erf[(x0+ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[-ϵ0*2*(σ^2+d*t)-x0]+Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[ϵ0*2*(σ^2+d*t)-x0];
P2[t_]=0;
k0=10;
x0=5;
d=2;
h=.2;
y=-1.9;
ys=0;p1=Plot[Pe[t,1],{t,0,100},PlotLegends->Placed[{Style["σ=1 nm",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[Pe[t,2],{t,0,100},PlotLegends->Placed[{Style["σ=2 nm",13]},{Scaled[{0,0.5}],{y,h-1}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Thick}]
p3=Plot[Pe[t,5],{t,0,100},PlotLegends->Placed[{Style["σ=5 nm",13]},{Scaled[{0,0.5}],{y-4.5ys,h-0}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
p5=Show[p1,p2,p3]
Export["sigeffecttwo.eps",p5,ImageResolution->1500,"CompressionLevel"->0]

ClearAll["Global"]
Remove["Global`*"]
Unprotect["Global`*"]

ClearAll["Global"];
ϵ0=k0/(2d);
Pe[t_,x0_]=1+Exp[-ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))-ϵ0]-Erf[(x0-ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Erfc[x0/Sqrt[4(σ^2+d*t)]]+Exp[ϵ0*x0+ϵ0^2(σ^2+d*t)]/2(Sign[x0/(2(σ^2+d*t))+ϵ0]-Erf[(x0+ϵ0*2*(σ^2+d*t))/Sqrt[4(σ^2+d*t)]])-Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[-ϵ0*2*(σ^2+d*t)-x0]+Exp[-ϵ0^2(σ^2+d*t)]*HeavisideTheta[ϵ0*2*(σ^2+d*t)-x0];
P2[t_]=0;
k0=10;
d=2;
σ=0;
h=.1;
y=-1.8;
ys=0;p1=Plot[Pe[t,1],{t,0,100},PlotLegends->Placed[{Style["(Subscript[x, 0])=1 nm",13]},{Scaled[{0,0.5}],{y-4.5ys,h-2}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dotted}]
p2=Plot[Pe[t,2],{t,0,100},PlotLegends->Placed[{Style["(Subscript[x, 0])=2 nm",13]},{Scaled[{0,0.5}],{y,h-1}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Thick}]
p3=Plot[Pe[t,5],{t,0,100},PlotLegends->Placed[{Style["(Subscript[x, 0])=5 nm",13]},{Scaled[{0,0.5}],{y-4.5ys,h-0}}],PlotRange->{0,1.1},PlotTheme->{"Scientific","Monochrome"},PlotLabel->{Style["Survival probability vs time",14]},FrameLabel->{Style["t (ps)",14],Style["(Subscript[Q, e])(t)",14]},PlotStyle->{Dashed}]
p5=Show[p1,p2,p3]
Export["x0effecttwo.eps",p5,ImageResolution->1500,"CompressionLevel"->0]

Show[,p15]

x0effecttwo.eps"