Creators: J. de Graaf and S. Samin
Created: 24.04.2019
Version: Mathematica 11.2.0.0

In this Mathematica notebook we solve the doubly linearized equation system governing self-thermo(di)electrophoresis for the intermediate regime of ionic strength. We also consider the two limiting cases, obtaining analytic expressions that we report on in the main text.

Introducing Useful Functions

In the following we will make use of spherical polar coordinates. It is therefore useful to introduce the gradient, divergence, and Laplacian for this system.

GradSPH[fun_]:=D[fun,r]*{1,0}+(1/r)*D[fun,θ]*{0,1}
DivSPH[fvec_]:=(1/(r^2))*D[(r^2)*fvec[[1]],r]+(1/(r*Sin[θ]))*D[Sin[θ]*fvec[[2]],θ]
LapSPH[fun_]:=DivSPH[GradSPH[fun]];

(* Check the forms of the derivatives. *)

Simplify[GradSPH[f[r,θ]]]
Simplify[DivSPH[{R[r,θ],T[r,θ]}]]
Simplify[LapSPH[f[r,θ]]]

The Equilibrium Equations

Here we solve for the reduced equilibrium electrostatic potential, which we will need in constructing and solving equations for the non-equilibrium fields. Throughout we will assume a homogeneous boundary condition for the electrostatic potential.

(* The equation for the equilibrium electrostatic potential. *)

phidiff=(LapSPH[pheq[r]]==κ^2*pheq[r]);
phieq=Expand[Simplify[DSolve[phidiff,pheq[r],r],Assumptions->{κ>0}][[1,1,2]]]

(* The part that diverges can be eliminated. *)

phieqn=phieq/.{C[2]->0}

(* The electrostatic potential for a conducting surface. *)

subcond=Solve[(phieqn/.{r->a})==ϕ0,C[1]][[1]];
ϕceq[r_]:=phieqn/.subcond;
ϕceq[r]

(* The electrostatic potential for an insulating surface. *)

dphieqn=Simplify[GradSPH[phieqn]][[1]];
subiso=Solve[(dphieqn/.{r->a})==-σ0*el/(ϵ0*kB*T0),C[1]][[1]];
ϕseq[r_]:=phieqn/.subiso;
ϕseq[r]

(* Note that these two only differ in terms of the prefactor. Hence and for convenience we will assume the generic form. *)

ϕeq[r_,θ_]:=ϕs*(a/r)*Exp[-κ*r];
ϕeq[r,θ]

The Non-Equilibrium Equations: The Temperature

The temperature field is easily solved for, because it is a simple solution to the Laplace equation.

(* Assume the following form for the temperature field. *)

T[r_,θ_]:=tr[r]*Cos[θ];

(* Specify the differential equation. *)

tdiff=Simplify[LapSPH[T[r,θ]]==0,Assumptions->{r>0,Cos[θ]>0}]

(* Solve for this equation. *)

tdsol=DSolve[tdiff,tr[r],r][[1,1,2]]

(* So that the full solution is given by. *)

tds=tdsol*Cos[θ]

(* Now we must study the limit to infinity of this equation. *)

Limit[tds,r->Infinity]

(* Clearly we must drop the C[1] prefactor term. *)

tdsn=tds/.{C[1]->0}

Limit[tdsn,r->Infinity]

(* The integration constant can be set by using the boundary condition. *)

tsubs=Solve[Simplify[(tdsn/.{r->a})==Cos[θ]],C[2]][[1]]

(* Solution to the reduced heat equation. *)

tdsn/.tsubs

(* Note that we have to play a trick to get the symbols recognised properly by Mathematica. Also we introduce a minus sign to account for the hot cap being on the bottom (negative values of z). *)

t[p_,q_]:=-(tdsn/.tsubs/.{r->p,θ->q});
t[r,θ]

The Non-Equilibrium Equations: The Ionic Strength

Next we take the steps necessary to solve for the non-equilibrium fields. We start with the local ionic strength which is decoupled from the other equations, other than through the boundary conditions.

(* Assume the following form for the local non-equilibrium ionic strength *)

X[r_,θ_]:=x[r]*Cos[θ];

(* Then the associated differential equation becomes. *)

diffXlhs=Simplify[LapSPH[X[r,θ]]]

diffXrhs=Simplify[β*Dot[GradSPH[t[r,θ]],GradSPH[ϕeq[r,θ]]]]

(* With this the piece of interest. *)

xdiff=FullSimplify[diffXlhs==diffXrhs,Assumptions->{r>0,a>0,Cos[θ]>0}]

(* Now solve this differential equation. *)

xdsol=DSolve[xdiff,x[r],r][[1,1,2]]

(* So that the full solution is given by. *)

xds=Simplify[Expand[xdsol*Cos[θ]]]

(* Now we must study the limit to infinity of this equation. *)

Limit[xds,r->Infinity,Assumptions->{a>0,κ>0}]

(* We conclude that we must drop the C[1] prefactored term. *)

xdsn=xds/.{C[1]->0}

Limit[xdsn,r->Infinity,Assumptions->{a>0,κ>0}]

(* We cannot as of yet determine the integration constant so that we arrive at. *)

Xs[p_,q_]:=(xdsn/.{C[2]->Xint,r->p,θ->q});
Xs[r,θ]

The Non-Equilibrium Equations: The Ion Excess

As we argue in the main text, the local ion excess for may be isolated from the non-equilibrium electrostatic potential, which makes it an ideal target to solve for next.

(* Assume the following form for the local non-equilibrium ionic strength. *)

dX[r_,θ_]:=dx[r]*Cos[θ];

(* Then the associated differential equation becomes. *)

diffdXlhs=Simplify[Expand[LapSPH[dX[r,θ]]-(κ^2)*dX[r,θ]]]

diffdXrhs=Simplify[(1+γ+ϵst)*Dot[GradSPH[t[r,θ]],GradSPH[ϕeq[r,θ]]]+(κ^2)*(1+ϵst)*t[r,θ]*ϕeq[r,θ]]

(* With this the piece of interest. *)

dxdiff=FullSimplify[diffdXlhs==diffdXrhs,Assumptions->{r>0,a>0,Cos[θ]>0}]

(* Now solve this differential equation. *)

dxdsol=DSolve[dxdiff,dx[r],r][[1,1,2]];
sdxdsol=Simplify[dxdsol,Assumptions->{a>0,κ>0,r>0}]

(* So that the full solution is given by. *)

dxds=Simplify[Expand[sdxdsol*Cos[θ]]]

(* Now we must study the limit to infinity of this equation. *)

Limit[dxds,r->Infinity,Assumptions->{a>0,κ>0}]

(* So that we must choose C[1] such that this term drops. *)

dxdsn=FullSimplify[dxds/.{C[1]->I*C[2]},Assumptions->{a>0,κ>0,r>0}]

Limit[dxdsn,r->Infinity,Assumptions->{a>0,κ>0}]

(* Note that the above implies that we are left with complex values in our expression. This means that C[2] must be some complex number as well. This is surprising from a physics perspective. Hence we fix this in assigning the as-of-yet-unknown integration constant. We are always free to make this choice. *)

dXs[p_,q_]:=(dxdsn/.{C[2]->I*dXint,r->p,θ->q});
dXs[r,θ]

The Non-Equilibrium Equations: The Poisson Equation

Now that we have a solution for the local ion excess, we can solve for the non-equilibrium part of the Poisson Equation, which depends on this quantity.

(* Assume the following form for the local non-equilibrium ionic strength. *)

P[r_,θ_]:=p[r]*Cos[θ];

(* Then the associated differential equation becomes. *)

diffPlhs=Simplify[LapSPH[P[r,θ]]]

diffPrhs=Simplify[-ϵst*DivSPH[t[r,θ]*GradSPH[ϕeq[r,θ]]]-(κ^2)*dXs[r,θ]]

(* With this the piece of interest. *)

pdiff=FullSimplify[diffPlhs==diffPrhs,Assumptions->{r>0,a>0,κ>0,Cos[θ]>0}]

(* Now solve this differential equation. *)

pdsol=DSolve[pdiff,p[r],r][[1,1,2]];
spdsol=Simplify[pdsol,Assumptions->{a>0,κ>0,r>0}]

(* So that the full solution is given by. *)

pds=Simplify[Expand[spdsol*Cos[θ]]]

(* Now we must study the limit to infinity of this equation. *)

Limit[pds,r->Infinity,Assumptions->{a>0,κ>0}]

(* Clearly we must drop the C[1] prefactor. *)

pdsn=FullSimplify[pds/.{C[1]->0},Assumptions->{a>0,κ>0,r>0}]

Limit[pdsn,r->Infinity,Assumptions->{a>0,κ>0}]

(* We must still determine the integration constant. *)

Ps[p_,q_]:=(pdsn/.{C[2]->Pint,r->p,θ->q});
Ps[r,θ]

The Non-Equilibrium Equations: Fixing the Integration Constants

We now have to use the flux boundary conditions on the particle for X and dX to fix the integration constants. The integration constant for the potential must be specified separately for conducting and insulating surfaces. Here we use a two-vector notation, where the first component corresponds to the radial contribution and the second to the polar.

(* The flux of positive ions is given by. *)

fluxp=Simplify[(β+γ)*(ϕeq[r,θ]-1)*GradSPH[t[r,θ]]-GradSPH[Xs[r,θ]]-GradSPH[dXs[r,θ]]+t[r,θ]*GradSPH[ϕeq[r,θ]]-GradSPH[Ps[r,θ]]];
fluxp//MatrixForm

(* The flux of negative ions is given by. *)

fluxm=Simplify[(β-γ)*(ϕeq[r,θ]+1)*GradSPH[t[r,θ]]-GradSPH[Xs[r,θ]]+GradSPH[dXs[r,θ]]-t[r,θ]*GradSPH[ϕeq[r,θ]]+GradSPH[Ps[r,θ]]];
fluxm//MatrixForm

(* Now evaluate these at the surface dotted with the normal to the surface *)

fpa=Dot[fluxp,{1,0}]/.{r->a}
fma=Dot[fluxm,{1,0}]/.{r->a}

(* Intriguingly this expression does not depend on dXint. One might have expected such a dependence as the fluxes only depend on X and dX and hence their integration constants. However this is an expression of the interconnectedness of dX and ϕ. We proceed by solving for zero flux through the boundary. *)

fsub=Solve[{fpa==0,fma==0},{Pint,Xint}][[1]];
fsub//MatrixForm

(* Now we substitute this solution into the original expressions. *)

Print["The Local Non-Equilibrium Ionic Strength:"]
Xsol=Simplify[Expand[Xs[r,θ]/.fsub]]
Print["The Local Non-Equilibrium Ion Excess:"]
dXsol=Simplify[Expand[dXs[r,θ]/.fsub]]
Print["The Non-Equilibrium Potential:"]
Psol=Simplify[Expand[Ps[r,θ]/.fsub]]

(* Next we must determine the remaining integration constant dXint from the electrostatic boundary conditions. We do this for a conducting surface first by equating the potential to zero at the surface. Note that the equilibrium electrostatic potential already accounted for the entirety of the electrostatic boundary condition. *)

Print["The Integration Constant for a Conducting Surface:"]
solcond=Solve[Simplify[Psol//.{r->a}]==0,dXint][[1]]
Print["The Local Non-Equilibrium Ionic Strength:"]
Xctot=Simplify[Expand[Xsol]]
Print["The Local Non-Equilibrium Ion Excess:"]
dXctot=Simplify[Expand[dXsol/.solcond]]
Print["The Non-Equilibrium Potential:"]
Pctot=Simplify[Expand[Psol/.solcond]]

(* For the insulating surface we must set the gradient of the potential equal to zero to ensure that no surface charge is added. *)

Print["The Integration Constant for an Insulating Surface:"]
soliso=Solve[Simplify[(Dot[GradSPH[Psol],{1,0}])//.{r->a}]==0,dXint][[1]]
Print["The Local Non-Equilibrium Ionic Strength:"]
Xstot=Xctot
Print["The Local Non-Equilibrium Ion Excess:"]
dXstot=Simplify[Expand[dXsol/.soliso]]
Print["The Non-Equilibrium Potential:"]
Pstot=Simplify[Expand[Psol/.soliso]]

The Swim Speed using Teubner’s Formalism

We now have solutions for the non-equilibrium potential and ion excess, as well as for the temperature and equilibrium potential, we can determine the speed of the particle.

(* We first determine the force acting on the fluid. Here we need to distinguish between conducting and onsulating surfaces. For convenience we also strip off the dielectrophoretic contribution which is equal in both cases barring the obvious difference in ϕs. *)

peq=ϕeq[r,θ];
temp=t[r,θ];
Fc=-2*kB*T0*n0*(dXctot*GradSPH[peq]-peq*GradSPH[Pctot]);
Fs=-2*kB*T0*n0*(dXstot*GradSPH[peq]-peq*GradSPH[Pstot]);
Fϵ=-2*kB*T0*n0*((1/(2*(κ^2)))*ϵst*Dot[GradSPH[peq],GradSPH[peq]]*GradSPH[temp]);

(* The integration kernel for Teubner's formalism is given by. Here we also introduce a prefactor that we will use as a multiplicative constant towards the end of the calculation. We also incorporate the azimuthal integration into this prefactor. *)

Kpref=(1/(6*Pi*eta*a))*2*Pi;
Tker={((3*a)/(2*r)-((a^3)/(2*(r^3)))-1)*Cos[θ],-((3*a)/(4*r)+((a^3)/(4*(r^3)))-1)*Sin[θ]};
Tker//MatrixForm

(* So that the integrand is given by the following expressions. Here we also take into account the Jacobian of the transformation to spherical polar coordinates. Lastly the variables are exchanged again to ensure all terms are recognized by Mathematica's integration routine. *)

Intc=(Dot[Tker,Fc]*(r^2)*Sin[θ])//.{r->ri,θ->ti};
Ints=(Dot[Tker,Fs]*(r^2)*Sin[θ])//.{r->ri,θ->ti};
Intϵ=(Dot[Tker,Fϵ]*(r^2)*Sin[θ])//.{r->ri,θ->ti};

(* Unfortunately the above presents a set of extremely challenging intgrands. We therefore break these up into terms and integrate the terms separately. *)

exIntc=Expand[Intc];
Length[exIntc]

exInts=Expand[Ints];
Length[exInts]

exIntϵ=Expand[Intϵ];
Length[exIntϵ]

(* The routine to integrate a single term *)

IntTerm[IIin_,Lin_]:=Module[{II=IIin,L=Lin,p0,p0a,p1,xi,yi,res},
p0=Simplify[L[[II]]];
p0a=Simplify[p0//.{ri->a}];
p1=Simplify[p0/p0a];
yi=Simplify[Integrate[p0a,{ti,0,Pi}]];
xi=Integrate[p1,{ri,a,Infinity},Assumptions->{a>0,κ>0}];
res=Simplify[xi*yi];
res];

The Speed resulting from the Dielectrophoretic Force Contribution

Let us consider the speed generated by the force term that scales with (SuperStar[ϵ]) first. Note that through the dependence of the other fields on this parameter, there will be other contributions, but this term is relatively short and will give us a feeling for the way the integration routines perform. There are dielectric contributions to the other integrands as well. Thus, even knowing these terms here, does not give us sufficient information to analyse the Hückel and Smoluchowski limits.

(*** Sum the full series for the dielectrophoretic term ***)

(* Uncomment the following lines to perform the calculation *)
(*exSϵ=Sum[IntTerm[k,exIntϵ],{k,1,Length[exIntϵ]}];*)
(*SϵR=Kpref*Simplify[Expand[exSϵ]]*)

(* The calculations are long so we present the result below. *)

Sϵ=(1/(630 a eta (κ^2)))(E^(-2 a κ)) kB n0 T0 ϵst ((ϕs)^2) (-12-24 a κ+4 (a^2) (κ^2)+5 (a^3) (κ^3)-22 (a^4) (κ^4)-6 (a^5) (κ^5)+12 (a^6) (κ^6)+8 (a^5) (E^(2 a κ)) (κ^5) (-7+3 (a^2) (κ^2)) ExpIntegralEi[-2 a κ]);

(* Dielectrophoretic term for conducting surface *)

Sϵc=Simplify[Sϵ//.{ϕs->ϕ0*Exp[a*κ],n0->κ^2*ϵ0*kB*T0/(2*el^2)}]

(* Dielectrophoretic term for conducting surface *)

Sϵs=Simplify[Sϵ//.{ϕs->a*el*σ0*Exp[a*κ]/(kB*T0*ϵ0*(1+κ*a)),n0->κ^2*ϵ0*kB*T0/(2*el^2)}]

The Speed resulting from the Other Forces for a Conducting Surface

Now that we have covered the dielectrophoretic force contribution, we turn our attention to the more complicated force expressions for the conductor. Here we add the dielectrophoretic contribution.

(*** Sum the full series for the conducting surface ***)

(* Uncomment the following lines to perform the calculation *)
(*exSc= Sum[IntTerm[k,exIntc],{k,1,Length[exIntc]}];*)
(*ScR=Kpref*Simplify[Expand[exSc]]*)

(* The calculations are long so we present the result below. *)

Sc=Simplify[(-(1/(4320 eta (1+a κ)))a (E^(-2 a κ)) kB n0 T0 ϕs (720 (E^(a κ)) β-480 a (E^(a κ)) β κ-1320 (a^2) (E^(a κ)) β (κ^2)+120 (a^4) (E^(a κ)) β (κ^4)+84 ϕs-276 γ ϕs-228 a κ ϕs+372 a γ κ ϕs-1204 (a^2) (κ^2) ϕs+676 (a^2) γ (κ^2) ϕs-900 (a^3) (κ^3) ϕs-60 (a^3) γ (κ^3) ϕs-37 (a^4) (κ^4) ϕs+133 (a^4) γ (κ^4) ϕs-19 (a^5) (κ^5) ϕs+211 (a^5) γ (κ^5) ϕs+5 (a^6) (κ^6) ϕs-5 (a^6) γ (κ^6) ϕs-5 (a^7) (κ^7) ϕs+5 (a^7) γ (κ^7) ϕs+5 (a^6) (E^(2 a κ)) (-1+γ) (κ^6) (-12-12 a κ+(a^2) (κ^2)+(a^3) (κ^3)) ϕs ((ExpIntegralEi[-a κ])^2)+32 (a^2) (E^(2 a κ)) (κ^2) ExpIntegralEi[-2 a κ] (360 (E^(a κ)) β+(165-30 a κ-150 (a^2) (κ^2)-49 (a^3) (κ^3)-4 (a^4) (κ^4)-60 ϵst (-5-2 a κ+3 (a^2) (κ^2))+γ (-45+90 a κ-90 (a^2) (κ^2)+61 (a^3) (κ^3)+16 (a^4) (κ^4))) ϕs+45 (a^4) (E^(a κ)) (-1+γ) (κ^4) ϕs ExpIntegralEi[-a κ])+11520 (a^2) (E^(3 a κ)) β (κ^2) Gamma[0,2 a κ]+5280 (a^2) (E^(2 a κ)) (κ^2) ϕs Gamma[0,2 a κ]-1440 (a^2) (E^(2 a κ)) γ (κ^2) ϕs Gamma[0,2 a κ]+9600 (a^2) (E^(2 a κ)) ϵst (κ^2) ϕs Gamma[0,2 a κ]+960 (a^3) (E^(2 a κ)) (κ^3) ϕs Gamma[0,2 a κ]+2880 (a^3) (E^(2 a κ)) γ (κ^3) ϕs Gamma[0,2 a κ]+3840 (a^3) (E^(2 a κ)) ϵst (κ^3) ϕs Gamma[0,2 a κ]-2880 (a^4) (E^(2 a κ)) (κ^4) ϕs Gamma[0,2 a κ]-2880 (a^4) (E^(2 a κ)) γ (κ^4) ϕs Gamma[0,2 a κ]-5760 (a^4) (E^(2 a κ)) ϵst (κ^4) ϕs Gamma[0,2 a κ]-1440 (a^5) (E^(2 a κ)) (κ^5) ϕs Gamma[0,2 a κ]+1440 (a^5) (E^(2 a κ)) γ (κ^5) ϕs Gamma[0,2 a κ]+10 (a^2) (E^(a κ)) (κ^2) ExpIntegralEi[-a κ] ((1+a κ) (12 (E^(a κ)) β (-12+(a^2) (κ^2))+(-1+γ) (84-12 a κ+6 (a^2) (κ^2)-10 (a^3) (κ^3)-(a^4) (κ^4)+(a^5) (κ^5)) ϕs)+144 (a^4) (E^(2 a κ)) (-1+γ) (κ^4) ϕs Gamma[0,2 a κ])))//.{ϕs->ϕ0*Exp[a κ],n0->(κ^2)*ϵ0*kB*T0/(2*el^2)}]+Sϵc;

(* Large Debye length limit. Note that the functional dependencies are slightly different than for the splitting approach. This makes sense because of the limitations of this method. *)

Normal[Series[Sc,{κ,0,2},Assumptions->{a>0,κ>0}]]/.{(κ^2)->(2*el^2)*n0/(ϵ0*kB*T0)}

(* Small Debye length limit. *)

U1=Limit[Sc[[1]],κ->Infinity,Assumptions->{a>0,el>0,eta>0,kB>0,T0>0,ϵ0>0,ϵst∈Reals,ϕ0∈Reals}];
U2=Limit[Sc[[2]],κ->Infinity,Assumptions->{a>0,el>0,eta>0,kB>0,T0>0,ϵ0>0,ϵst∈Reals,ϕ0∈Reals,β∈Reals,γ∈Reals}];
Usmall=Simplify[U1+U2]

The Speed resulting from the Other Forces for an Insulating Surface

Last we integrate out the contributions for the insulating surface.

(*** Sum the full series for the insulating surface ***)

(* Uncomment the following lines to perform the calculation *)
(*exSs=Sum[IntTerm[k,exInts],{k,1,Length[exInts]}];*)
(*SsR=Kpref*Simplify[Expand[exSs]]*)

(* The calculations are long so we present the result below. *)

Ss=Simplify[(-(1/(4320 eta (2+2 a κ+(a^2) (κ^2))))a (E^(-2 a κ)) kB n0 T0 ϕs (1440 (E^(a κ)) β-960 a (E^(a κ)) β κ-1920 (a^2) (E^(a κ)) β (κ^2)-1200 (a^3) (E^(a κ)) β (κ^3)+120 (a^4) (E^(a κ)) β (κ^4)+120 (a^5) (E^(a κ)) β (κ^5)+168 ϕs-552 γ ϕs-456 a κ ϕs+744 a γ κ ϕs-2324 (a^2) (κ^2) ϕs+1076 (a^2) γ (κ^2) ϕs-2112 (a^3) (κ^3) ϕs+528 (a^3) γ (κ^3) ϕs-966 (a^4) (κ^4) ϕs+294 (a^4) γ (κ^4) ϕs-46 (a^5) (κ^5) ϕs+334 (a^5) γ (κ^5) ϕs-19 (a^6) (κ^6) ϕs+211 (a^6) γ (κ^6) ϕs-5 (a^8) (κ^8) ϕs+5 (a^8) γ (κ^8) ϕs+32 (a^2) (E^(2 a κ)) (κ^2) (720 (E^(a κ)) β+(330-60 a κ-135 (a^2) (κ^2)-203 (a^3) (κ^3)-8 (a^4) (κ^4)-4 (a^5) (κ^5)-60 ϵst (-13-7 a κ+(a^2) (κ^2)+3 (a^3) (κ^3))+γ (-90+180 a κ-45 (a^2) (κ^2)-13 (a^3) (κ^3)+32 (a^4) (κ^4)+16 (a^5) (κ^5))) ϕs) ExpIntegralEi[-2 a κ]+10 (a^2) (E^(a κ)) (κ^2) (2+2 a κ+(a^2) (κ^2)) (12 (E^(a κ)) β (-12+(a^2) (κ^2))+(-1+γ) (84-12 a κ+6 (a^2) (κ^2)-10 (a^3) (κ^3)-(a^4) (κ^4)+(a^5) (κ^5)) ϕs) ExpIntegralEi[-a κ]+5 (a^6) (E^(2 a κ)) (-1+γ) (κ^6) (-24-24 a κ-10 (a^2) (κ^2)+2 (a^3) (κ^3)+(a^4) (κ^4)) ϕs ((ExpIntegralEi[-a κ])^2)+23040 (a^2) (E^(3 a κ)) β (κ^2) Gamma[0,2 a κ]+10560 (a^2) (E^(2 a κ)) (κ^2) ϕs Gamma[0,2 a κ]-2880 (a^2) (E^(2 a κ)) γ (κ^2) ϕs Gamma[0,2 a κ]+24960 (a^2) (E^(2 a κ)) ϵst (κ^2) ϕs Gamma[0,2 a κ]+1920 (a^3) (E^(2 a κ)) (κ^3) ϕs Gamma[0,2 a κ]+5760 (a^3) (E^(2 a κ)) γ (κ^3) ϕs Gamma[0,2 a κ]+13440 (a^3) (E^(2 a κ)) ϵst (κ^3) ϕs Gamma[0,2 a κ]-480 (a^4) (E^(2 a κ)) (κ^4) ϕs Gamma[0,2 a κ]-1440 (a^4) (E^(2 a κ)) γ (κ^4) ϕs Gamma[0,2 a κ]-1920 (a^4) (E^(2 a κ)) ϵst (κ^4) ϕs Gamma[0,2 a κ]-4320 (a^5) (E^(2 a κ)) (κ^5) ϕs Gamma[0,2 a κ]-1440 (a^5) (E^(2 a κ)) γ (κ^5) ϕs Gamma[0,2 a κ]-5760 (a^5) (E^(2 a κ)) ϵst (κ^5) ϕs Gamma[0,2 a κ]))//.{ϕs->a*el*σ0*Exp[a*κ]/(kB*T0*ϵ0*(1+κ*a)),n0->κ^2*ϵ0*kB*T0/(2*el^2)}]+Sϵs;

(* Large Debye Length Limit *)

Normal[Series[Ss,{κ,0,2},Assumptions->{a>0,κ>0}]]/.{(κ^2)->(2*el^2)*n0/(ϵ0*kB*T0)}

(* Small Debye Length Limit, which is vanishing as expected *)

U1=Limit[Ss[[1]],κ->Infinity,Assumptions->{a>0,el>0,eta>0,kB>0,T0>0,ϵ0>0,ϵst∈Reals,σ0∈Reals}];
U2=Limit[Ss[[2]],κ->Infinity,Assumptions->{a>0,el>0,eta>0,kB>0,T0>0,ϵ0>0,ϵst∈Reals,σ0∈Reals,β∈Reals,γ∈Reals}];
Usmall=Simplify[U1+U2]*((2*el^2)*n0*λ0^2/(ϵ0*kB*T0))