Numerical analysis - KL model - Heat current - Gapless- for fixed n different temperature

Quit[];

(*In this section we analysed the linear temperature dependence of quantum current in the gapless regime*)

dω=0.00025;ωtab=Table[i,{i,-N[1.8],N[1.8],dω}];
ωntab=Drop[ωtab,{Position[ωtab,0.][[1,1]]}];
zl[m_, r_,ω_,k_,v1_,γ_]:= -m (r^2) (ω^2)+k (v1^2)- I *(γ *ω);
zr[m_, r_,ω_,k_,v2_,γ_]:= -m (r^2) (ω^2)+k (v2^2)- I *(γ *ω);
z[m_,r_,ω_,k_,v1_,v2_]:= -m (r^2) (ω^2)+k ((v1^2) + (v2^2));
m=1.0;r=2.0;a=1.0; k=1.0;γ=0.2;η=0.2;
f[ω_,T_]:= (1/(Exp[((ℏ ω)/(kb T))]-1));kb=1.0;ℏ=1.0; n=32;

Θ=N[0.0];

v1 = (((2r Sin[Θ]+a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));
v2=  (((2r Sin[Θ]-a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));

v1

v2

(*dtl1=0.025;
tl1tab=Table[i,{i,0.0101,0.1001,dtl1}];
dtl2=0.1;
tl2tab=Table[i,{i,0.1001,2.1001,dtl2}];
tltab=Join[tl1tab,tl2tab];
dtr1=0.025;
tr1tab=Table[i,{i,0.01,0.10,dtr1}];
dtr2=0.1;
tr2tab=Table[i,{i,0.10,2.10,dtr2}];
trtab=Join[tr1tab,tr2tab];
Ttab=((tltab+trtab)/2);
tltab-trtab;*)

(*Ttab*)

(*dtl=0.1;tltab=Table[i,{i,2.9,3.0,dtl}];trtab=Table[i,{i,1.9,2.0,dtl}];
Ttab=((tltab+trtab)/2);*)

dtl=0.005;tltab=Table[i,{i,0.0051,0.2001,dtl}];trtab=Table[i,{i,0.005,0.2,dtl}];
Ttab=((tltab+trtab)/2);

Ttab

GList={};
transcofflist={};
currlist={};
currclasslist={};
currQlist={};
gmatslist={};
mats= Table[0,{n},{n}];
For[l=1,l<=Length[ωntab],l++,
For[i=1,i<=n-1,i++,
mats[[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];mats[[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];mats[[n,n]]=zr[m,r,ωntab[[l]],k,v2,η];mats[[i,i+1]]=-k * Abs[v1 v2];mats[[i+1,i]]=-k * Abs[v1 v2]; Z=mats;];GList=Join[GList,{{ωntab[[l]],Inverse[Z]}}];gmatslist=Join[gmatslist,{{ωntab[[l]],((Abs[GList[[l,2,1,n]]])^2)}}];
transcofflist=Join[transcofflist,{{ωntab[[l]],(((Abs[GList[[l,2,1,n]]])^2))*((ωntab[[l]])^2)*γ*η}}];]

For[p=1,p<=Length[Ttab],p++,For[l=1,l<=Length[ωntab],l++,currlist=Join[currlist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*((ℏ ((ωntab[[l]])^3))/π)*γ*η*(f[ωntab[[l]],tltab[[p]]]-f[ωntab[[l]],trtab[[p]]])}}];
currclasslist=Join[currclasslist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*(((ωntab[[l]])^2)/π)*γ*η*(tltab[[p]]-trtab[[p]])}}];
currQlist=Join[currQlist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*(((ℏ^2) ((ωntab[[l]])^4))/((((2*kb*Ttab[[p]]))^2)*π))*γ*η*kb*(tltab[[p]]-trtab[[p]])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])}}];]]

τ={};
For[j=1,j<=Length[Ttab],j++,τ=Join[τ,{{Ttab[[j]],Sum[currQlist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

ν={};
For[j=1,j<=Length[Ttab],j++,ν=Join[ν,{{Ttab[[j]],Sum[currlist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

η={};
For[j=1,j<=Length[Ttab],j++,η=Join[η,{{Ttab[[j]],Sum[currclasslist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

Export["CHC_tempDep_gapless.dat",η]

τ//N

ν

η

(*This plot shows the linear temperature dependence of the quantum current in the gapless regime*)

ListPlot[{τ,ν,η},PlotRange->{0,All}]

(*Linear fitting*)

τ1=Table[{τ[[i]][[1]],τ[[i]][[2]]},{i,1,15}]

τ1={{0.005050000000000001`,1.0333955711047692`*^-7},{0.01005`,2.0577737427840304`*^-7},{0.015050000000000001`,3.084796412100333`*^-7},{0.02005`,4.103784064571412`*^-7},{0.025050000000000003`,5.111831027430532`*^-7},{0.03005`,6.105786380622139`*^-7},{0.03505`,7.082086676732406`*^-7},{0.04005`,8.036627711722059`*^-7},{0.04505`,8.964828469590764`*^-7},{0.05005`,9.861942293242202`*^-7},{0.05505`,1.0723501870822204`*^-6},{0.06005`,1.1545732066012608`*^-6}}

Fit[τ1,{1,x},x]

fit1=Table[{x,2.135877682500529*^-8+0.000019239383785902408* x},{x,0,0.2,0.01}]

Export["QHC_tempDep_gaplessF.dat",fit1]

p1=Show[{ListPlot[{τ,η},PlotRange->{0,All}],Plot[4.3046956238498174*^-8+0.000018376455287362042*x,{x,0,0.22}]}]

p1=Show[ListLogPlot[{τ,η,fit1},PlotRange->{0,All}]]

(*ListLinePlot[gmatslist,PlotRange->All]*)

(*ListLinePlot[transcofflist,PlotRange->All]*)

(*ListLinePlot[transcofflist,PlotRange->{{-0.7,0.7},{-0.25,0.25}}]*)

m=1.0;r=1;a=0.025; k=1.0;γ=0.2;η=0.2;

m=1.0;r=2.0;a=1.0; k=1.0;γ=0.2;η=0.2;

Numerical analysis - KL model - Heat current - Gapless - for fixed temperature and different system size

(*In this section we analysed the system size dependence of quantum current in the gapless regime and showed that the transport is ballistic*)

dω=0.00025;ωtab=Table[i,{i,-N[1.8],N[1.8],dω}];
ωntab=Drop[ωtab,{Position[ωtab,0.][[1,1]]}];
zl[m_, r_,ω_,k_,v1_,γ_]:= -m (r^2) (ω^2)+k (v1^2)- I *(γ *ω);
zr[m_, r_,ω_,k_,v2_,γ_]:= -m (r^2) (ω^2)+k (v2^2)- I *(γ *ω);
z[m_,r_,ω_,k_,v1_,v2_]:= -m (r^2) (ω^2)+k ((v1^2) + (v2^2));m=1.0;r=2.0;a=1.0; k=1.0;γ=0.2;η=0.2;
f[ω_,T_]:= (1/(Exp[((ℏ ω)/(kb T))]-1));kb=1.0;ℏ=1.0;

ntab={2,4,8,16,24,32,40,48};

Θ=N[0.0];

v1 = (((2r Sin[Θ]+a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));
v2=  (((2r Sin[Θ]-a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));

v1

v2

(*dtl1=0.025;
tl1tab=Table[i,{i,0.0101,0.1001,dtl1}];
dtl2=0.1;
tl2tab=Table[i,{i,0.1001,2.1001,dtl2}];
tltab=Join[tl1tab,tl2tab];
dtr1=0.025;
tr1tab=Table[i,{i,0.01,0.10,dtr1}];
dtr2=0.1;
tr2tab=Table[i,{i,0.10,2.10,dtr2}];
trtab=Join[tr1tab,tr2tab];
Ttab=((tltab+trtab)/2);
tltab-trtab;*)

(*Ttab*)

(*dtl=0.1;tltab=Table[i,{i,2.9,3.0,dtl}];trtab=Table[i,{i,1.9,2.0,dtl}];
Ttab=((tltab+trtab)/2);*)

dtl=0.015;trtab=Table[i,{i,0.05,0.06,dtl}];tltab=Table[i,{i,0.06,0.07,dtl}];
Ttab=((tltab+trtab)/2);

Ttab

GListtab={};
transcofflisttab={};
currlisttab={};
currclasslisttab={};
currQlisttab={};
gmatslisttab={};
matstab= Table[Table[0,{ntab[[i]]},{ntab[[i]]}],{i,1,Length[ntab]}];
Ztab=Table[Table[0,{ntab[[i]]},{ntab[[i]]}],{i,1,Length[ntab]}];

(*For[l=1,l<=Length[ωntab],l++,For[j=1,j<=Length[ntab],j++,
For[i=1,i<=ntab[[j]]-1,i++,
matstab[[j]][[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];matstab[[j]][[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];matstab[[j]][[ntab[[j]],ntab[[j]]]]=zr[m,r,ωntab[[l]],k,v2,η];matstab[[j]][[i,i+1]]=-k * Abs[v1 v2];
matstab[[j]][[i+1,i]]=-k * Abs[v1 v2]];Ztab=Join[Ztab, {matstab[[j]]}];];]*)

For[l=1,l<=Length[ωntab],l++,For[j=1,j<=Length[ntab],j++,
For[i=1,i<=ntab[[j]]-1,i++,
matstab[[j]][[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];matstab[[j]][[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];matstab[[j]][[ntab[[j]],ntab[[j]]]]=zr[m,r,ωntab[[l]],k,v2,η];matstab[[j]][[i,i+1]]=-k * Abs[v1 v2];
matstab[[j]][[i+1,i]]=-k * Abs[v1 v2]];
Ztab[[j]]=matstab[[j]];
GListtab=Join[GListtab,{{ωntab[[l]],Inverse[Ztab[[j]]]}}];];]

combineByLabel[list_]:={#[[1,1]],#[[All,2]]}&/@GatherBy[list,First];

GListtabnew=combineByLabel[GListtab];

For[p=1,p<=Length[Ttab],p++,For[j=1,j<=Length[ntab],j++,For[l=1,l<=Length[ωntab],l++,currQlisttab=Join[currQlisttab,{{ntab[[j]],{Ttab[[p]],((Abs[GListtabnew[[l,2,j,1,ntab[[j]]]]])^2)*(((ℏ^2) ((ωntab[[l]])^4))/((((2*kb*Ttab[[p]]))^2)*π))*γ*η*kb*(tltab[[p]]-trtab[[p]])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])}}}];]]]

currQlisttabnew={#[[1,1]],#[[All,2]]}&/@GatherBy[currQlisttab,First];

tauForSize[sizeIndex_]:=Table[{Ttab[[j]],Sum[currQlisttabnew[[sizeIndex,2,i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]},{j,Length[Ttab]}]

tauAll=Table[{currQlisttabnew[[s,1]],(*system size n*)tauForSize[s]},{s,Length[currQlisttabnew]}];

tauAll

(*This plot shows the ballistic nature of the quantum current in the gapless regime*)

ListPlot[Table[{tauAll[[i]][[1]],Abs[tauAll[[i,2,1,2]]]},{i,1,Length[ntab]}],PlotRange->{0,All}]

Numerical analysis - KL model - Heat current - Gapped- for fixed n different temperature

m=1.0;r=1;a=0.025; k=1.0;γ=0.2;η=0.2;

m=1.0;r=2.0;a=1.0; k=1.0;γ=0.2;η=0.2;Θ=N[(π/6)];

(*In this section we analysed the linear temperature dependence of quantum current in the gapped regime*)

dω=0.00025;ωtab=Table[i,{i,-N[1.8],N[1.8],dω}];
ωntab=Drop[ωtab,{Position[ωtab,0.][[1,1]]}];
zl[m_, r_,ω_,k_,v1_,γ_]:= -m (r^2) (ω^2)+k (v1^2)- I *(γ *ω);
zr[m_, r_,ω_,k_,v2_,γ_]:= -m (r^2) (ω^2)+k (v2^2)- I *(γ *ω);
z[m_,r_,ω_,k_,v1_,v2_]:= -m (r^2) (ω^2)+k ((v1^2) + (v2^2));m=1.0;r=2.0;a=1.0; k=1.0;γ=0.2;η=0.2;
f[ω_,T_]:= (1/(Exp[((ℏ ω)/(kb T))]-1));kb=1.0;ℏ=1.0; n=32;

Θ=N[(π/6)];

v1 = (((2r Sin[Θ]+a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));
v2=  (((2r Sin[Θ]-a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));

v1

v2

ω0=(Sqrt[(k/(m (r^2)))])(v1-v2)

(*dtl1=0.025;
tl1tab=Table[i,{i,0.0101,0.1001,dtl1}];
dtl2=0.1;
tl2tab=Table[i,{i,0.1001,2.1001,dtl2}];
tltab=Join[tl1tab,tl2tab];
dtr1=0.025;
tr1tab=Table[i,{i,0.01,0.10,dtr1}];
dtr2=0.1;
tr2tab=Table[i,{i,0.10,2.10,dtr2}];
trtab=Join[tr1tab,tr2tab];
Ttab=((tltab+trtab)/2);
tltab-trtab;*)

(*Ttab*)

(*dtl=0.1;tltab=Table[i,{i,2.9,3.0,dtl}];trtab=Table[i,{i,1.9,2.0,dtl}];
Ttab=((tltab+trtab)/2);*)

dtl=0.005;tltab=Table[i,{i,0.0051,0.2001,dtl}];trtab=Table[i,{i,0.005,0.2,dtl}];
Ttab=((tltab+trtab)/2);

dtl=0.015;tltab=Table[i,{i,0.051,0.21,dtl}];trtab=Table[i,{i,0.05,0.2,dtl}];
Ttab=((tltab+trtab)/2);

Ttab

GList={};
transcofflist={};
currlist={};
currclasslist={};
currQlist={};
gmatslist={};
mats= Table[0,{n},{n}];
For[l=1,l<=Length[ωntab],l++,
For[i=1,i<=n-1,i++,
mats[[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];mats[[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];mats[[n,n]]=zr[m,r,ωntab[[l]],k,v2,η];mats[[i,i+1]]=-k * Abs[v1 v2];mats[[i+1,i]]=-k * Abs[v1 v2]; Z=mats;];GList=Join[GList,{{ωntab[[l]],Inverse[Z]}}];gmatslist=Join[gmatslist,{{ωntab[[l]],((Abs[GList[[l,2,1,n]]])^2)}}];
transcofflist=Join[transcofflist,{{ωntab[[l]],(((Abs[GList[[l,2,1,n]]])^2))*((ωntab[[l]])^2)*γ*η}}];]

For[p=1,p<=Length[Ttab],p++,For[l=1,l<=Length[ωntab],l++,currlist=Join[currlist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*((ℏ ((ωntab[[l]])^3))/π)*γ*η*(f[ωntab[[l]],tltab[[p]]]-f[ωntab[[l]],trtab[[p]]])}}];
currclasslist=Join[currclasslist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*(((ωntab[[l]])^2)/π)*γ*η*(tltab[[p]]-trtab[[p]])}}];
currQlist=Join[currQlist,{{Ttab[[p]],((Abs[GList[[l,2,1,n]]])^2)*(((ℏ^2) ((ωntab[[l]])^4))/((((2*kb*Ttab[[p]]))^2)*π))*γ*η*kb*(tltab[[p]]-trtab[[p]])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])}}];]]

τ={};
For[j=1,j<=Length[Ttab],j++,τ=Join[τ,{{Ttab[[j]],Sum[currQlist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

ν={};
For[j=1,j<=Length[Ttab],j++,ν=Join[ν,{{Ttab[[j]],Sum[currlist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

η={};
For[j=1,j<=Length[Ttab],j++,η=Join[η,{{Ttab[[j]],Sum[currclasslist[[i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]}}]]

τ

ν

η

Export["QHC_tempDep_gappedF.dat",fitd]

(*This plot shows the temperature dependence of the quantum current in the gapped regime*)

ListPlot[{τ,ν,η},PlotRange->{0,All}]

ω0=(Sqrt[(k/(m (r^2)))])(v1-v2)

ListLogPlot[{τ,ν,η},PlotRange->{0,All}]

9.495806577848959`*^-7

T1=0.10005;0.0000001*(((Exp[((-ℏ ω0)/(kb T1))])/(Sqrt[T1])))

LogPlot[0.0000005*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.005,0.1}]

ListLogPlot[τ,PlotRange->{0,All}]

p2=Show[LogPlot[0.0000005*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.001,0.1}],ListLogPlot[τ,PlotRange->{0,All}]]

p2=Show[LogPlot[0.000001*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.005,0.2}],ListLogPlot[τ,PlotRange->{0,All}]]

(*In this plot we compared the temperature dependence of quantum current with the analytic expression*)

p2=Show[Plot[0.000001*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.005,0.2}],ListPlot[τ,PlotRange->{0,All}]]

fitd=Table[{T,0.000001*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T]))},{T,0.001,0.2,0.001}];

(*ListLinePlot[gmatslist,PlotRange->All]*)

(*ListLinePlot[transcofflist,PlotRange->All]*)

(*ListLinePlot[transcofflist,PlotRange->{{-0.7,0.7},{-0.25,0.25}}]*)

Numerical analysis - KL model - Heat current - Gapped - for fixed temperature and different system size

(*In this section we analysed the system size dependence of quantum current in the gapped regime and showed that the transport is exponentially suppressed*)

dω=0.00025;ωtab=Table[i,{i,-N[1.8],N[1.8],dω}];
ωntab=Drop[ωtab,{Position[ωtab,0.][[1,1]]}];
zl[m_, r_,ω_,k_,v1_,γ_]:= -m (r^2) (ω^2)+k (v1^2)- I *(γ *ω);
zr[m_, r_,ω_,k_,v2_,γ_]:= -m (r^2) (ω^2)+k (v2^2)- I *(γ *ω);
z[m_,r_,ω_,k_,v1_,v2_]:= -m (r^2) (ω^2)+k ((v1^2) + (v2^2));
m=1.0;r=1;a=0.025; k=1.0;γ=0.2;η=0.2;
f[ω_,T_]:= (1/(Exp[((ℏ ω)/(kb T))]-1));kb=1.0;ℏ=1.0;

ntab={64};

Θ=N[(π/3)];

v1 = (((2r Sin[Θ]+a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));
v2=  (((2r Sin[Θ]-a)r Cos[Θ])/(Sqrt[(a^2)+ 4 (r^2) ((Cos[Θ])^2)]));

v1

v2

ω0=(Sqrt[(k/(m (r^2)))])(v1-v2)

(*dtl1=0.025;
tl1tab=Table[i,{i,0.0101,0.1001,dtl1}];
dtl2=0.1;
tl2tab=Table[i,{i,0.1001,2.1001,dtl2}];
tltab=Join[tl1tab,tl2tab];
dtr1=0.025;
tr1tab=Table[i,{i,0.01,0.10,dtr1}];
dtr2=0.1;
tr2tab=Table[i,{i,0.10,2.10,dtr2}];
trtab=Join[tr1tab,tr2tab];
Ttab=((tltab+trtab)/2);
tltab-trtab;*)

(*Ttab*)

(*dtl=0.1;tltab=Table[i,{i,2.9,3.0,dtl}];trtab=Table[i,{i,1.9,2.0,dtl}];
Ttab=((tltab+trtab)/2);*)

dtl=0.015;trtab=Table[i,{i,0.005,0.01,dtl}];tltab=Table[i,{i,0.01,0.02,dtl}];
Ttab=((tltab+trtab)/2);

Ttab=((tltab+trtab)/2);

Ttab

GListtab={};
transcofflisttab={};
currlisttab={};
currclasslisttab={};
currQlisttab={};
gmatslisttab={};
matstab= Table[Table[0,{ntab[[i]]},{ntab[[i]]}],{i,1,Length[ntab]}];
Ztab=Table[Table[0,{ntab[[i]]},{ntab[[i]]}],{i,1,Length[ntab]}];

(*For[l=1,l<=Length[ωntab],l++,For[j=1,j<=Length[ntab],j++,
For[i=1,i<=ntab[[j]]-1,i++,
matstab[[j]][[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];matstab[[j]][[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];matstab[[j]][[ntab[[j]],ntab[[j]]]]=zr[m,r,ωntab[[l]],k,v2,η];matstab[[j]][[i,i+1]]=-k * Abs[v1 v2];
matstab[[j]][[i+1,i]]=-k * Abs[v1 v2]];Ztab=Join[Ztab, {matstab[[j]]}];];]*)

For[l=1,l<=Length[ωntab],l++,For[j=1,j<=Length[ntab],j++,
For[i=1,i<=ntab[[j]]-1,i++,
matstab[[j]][[i,i]]=z[m,r,ωntab[[l]],k,v1,v2];matstab[[j]][[1,1]]=zl[m,r,ωntab[[l]],k,v1,γ];matstab[[j]][[ntab[[j]],ntab[[j]]]]=zr[m,r,ωntab[[l]],k,v2,η];matstab[[j]][[i,i+1]]=-k * Abs[v1 v2];
matstab[[j]][[i+1,i]]=-k * Abs[v1 v2]];
Ztab[[j]]=matstab[[j]];
GListtab=Join[GListtab,{{ωntab[[l]],Inverse[Ztab[[j]]]}}];];]

combineByLabel[list_]:={#[[1,1]],#[[All,2]]}&/@GatherBy[list,First];

GListtabnew=combineByLabel[GListtab];

For[p=1,p<=Length[Ttab],p++,For[j=1,j<=Length[ntab],j++,For[l=1,l<=Length[ωntab],l++,currQlisttab=Join[currQlisttab,{{ntab[[j]],{Ttab[[p]],((Abs[GListtabnew[[l,2,j,1,ntab[[j]]]]])^2)*(((ℏ^2) ((ωntab[[l]])^4))/((((2*kb*Ttab[[p]]))^2)*π))*γ*η*kb*(tltab[[p]]-trtab[[p]])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])*(Csch[((ℏ ωntab[[l]])/(2kb*Ttab[[p]]))])}}}];]]]

currQlisttabnew={#[[1,1]],#[[All,2]]}&/@GatherBy[currQlisttab,First];

tauForSize[sizeIndex_]:=Table[{Ttab[[j]],Sum[currQlisttabnew[[sizeIndex,2,i,2]]*dω,{i,(j-1)*Length[ωntab]+1,j*Length[ωntab]}]},{j,Length[Ttab]}]

tauAll=Table[{currQlisttabnew[[s,1]],(*system size n*)tauForSize[s]},{s,Length[currQlisttabnew]}];

tauAll

tauAll

tauAll={{2,{{0.0075`,0.000019320477766911355`}}},{4,{{0.0075`,0.00001869271711650402`}}},{8,{{0.0075`,0.000016734380797479692`}}},{16,{{0.0075`,0.000012366084928175385`}}},{24,{{0.0075`,8.762773218183325`*^-6}}},{32,{{0.0075`,6.1270182202976994`*^-6}}},{40,{{0.0075`,4.321689645470178`*^-6}}},{48,{{0.0075`,3.1440239882379405`*^-6}}},{64,{{0.0075`,1.9376000260691236`*^-6}}}}

(*This plot shows the insulating nature of the quantum current in the gapped regime*)

ListLogPlot[Table[{tauAll[[i]][[1]],Abs[tauAll[[i,2,1,2]]]},{i,1,Length[ntab]}],PlotRange->All]

LogPlot[0.00015*Exp[((-2(x-1)Abs[v1-v2])/(a(Sqrt[v1 v2])))],{x ,2,32}]

data=Table[{tauAll[[i]][[1]],tauAll[[i,2,1,2]]},{i,1,Length[ntab]}]

nlm=NonlinearModelFit[data,c* Exp[-b*x],{c,b},x]

Show[LogPlot[nlm[x],{x,1,50},PlotRange->All],ListLogPlot[Table[{tauAll[[i]][[1]],Abs[tauAll[[i,2,1,2]]]},{i,1,Length[ntab]}]]]

Show[{Plot[0.00095*Exp[((-2(x-1)Abs[v1-v2])/(a(Sqrt[v1 v2])))],{x ,0,25}],ListPlot[Table[{tauAll[[i]][[1]],tauAll[[i,2,1,2]]},{i,1,Length[ntab]}],PlotRange->{0,All}]}]

0.00012968*16

(π/3)//N

data//N

data={{2,Log[0.000019320477766911355]},{4,Log[0.00001869271711650402]},{8,Log[0.000016734380797479692]},{16,Log[0.000012366084928175385]},{24,Log[8.762773218183325*^-6]},{32,Log[6.1270182202976994*^-6]},{40,Log[4.321689645470178*^-6]},{48,Log[3.1440239882379405*^-6]}}

data={{8,Log[0.000016734380797479692]},{16,Log[0.000012366084928175385]},{24,Log[8.762773218183325*^-6]},{32,Log[6.1270182202976994*^-6]},{40,Log[4.321689645470178*^-6]},{48,Log[3.1440239882379405*^-6]},{64,Log[1.9376000260691236`*^-6]}};nlm=FindFit[data,c1+b1*x,{c1,b1},x]

Log[1.9376000260691236`*^-6]

Show[Plot[-10.7041-0.0396275*x,{x,1,64},PlotRange->All],ListPlot[data]]

data3={{2,0.000019320477766911355},{4,0.00001869271711650402},{8,0.000016734380797479692},{16,0.000012366084928175385},{24,8.762773218183325*^-6},{32,6.1270182202976994*^-6},{40,4.321689645470178*^-6},{48,3.1440239882379405*^-6}}

data2={{2,0.00011482340786421806},{4,0.00010724100043431044},{8,0.00010714780997200385},{16,0.00010714914358558576},{24,0.0001071474519802976},{32,0.00010715075401504113},{40,0.00010714637948420999},{48,0.00010714792747682814}}

Log[0.000107147]

Show[Plot[-10.6409-0.0423984*x,{x,1,50}],Plot[-9.14131,{x,1,50}],ListLogPlot[data2],ListPlot[data],PlotRange->All]

Show[LogPlot[0.000107147,{x,1,50},PlotRange->{{1,50},{2*(10^(-6)),2*(10^(-4))}}],LogPlot[Exp[-10.6409]*Exp[-0.0423984*x],{x,1,50}],ListLogPlot[data2],ListLogPlot[data3]]

tab1=Table[{x,Exp[-10.6409]*Exp[-0.0423984*x]},{x,0,50,1}]

tab1=Table[{x,Exp[-10.7041]*Exp[-0.0396275*x]},{x,0,64,1}]

Export["FitGappedL.dat",tab1]

dat0={{0.005050000000000001`,9.4751269997291`*^-8},{0.01005`,1.5091184451111626`*^-7},{0.015050000000000001`,2.2395884857978085`*^-7},{0.02005`,3.178435862060106`*^-7},{0.025050000000000003`,4.212276853802442`*^-7},{0.03005`,5.285607710072747`*^-7},{0.03505`,6.37829738331765`*^-7},{0.04005`,7.481913338807035`*^-7},{0.04505`,8.592029619536827`*^-7},{0.05005`,9.705937299920762`*^-7},{0.05505`,1.082181846670779`*^-6},{0.06005`,1.1938371035011948`*^-6},{0.06505`,1.3054610708233953`*^-6},{0.07005`,1.4169757156400702`*^-6},{0.07505`,1.5283164847912737`*^-6},{0.08005000000000001`,1.6394279112925707`*^-6},{0.08505000000000001`,1.750260709974525`*^-6},{0.09005000000000002`,1.8607697831223342`*^-6},{0.09505`,1.9709128002364458`*^-6},{0.10005`,2.0806491534736527`*^-6}}

Export["CHC_tempDep_gapped.dat",dat1]

dat1={{0.005050000000000001`,9.295608827842504`*^-6},{0.01005`,9.29560882784258`*^-6},{0.015050000000000001`,9.29560882784258`*^-6},{0.02005`,9.29560882784243`*^-6},{0.025050000000000003`,9.29560882784243`*^-6},{0.03005`,9.29560882784243`*^-6},{0.03505`,9.295608827842757`*^-6},{0.04005`,9.295608827842757`*^-6},{0.04505`,9.295608827842757`*^-6},{0.05005`,9.295608827842757`*^-6},{0.05505`,9.295608827842757`*^-6},{0.06005`,9.295608827842757`*^-6},{0.06505`,9.29560882784144`*^-6},{0.07005`,9.29560882784144`*^-6},{0.07505`,9.29560882784144`*^-6},{0.08005000000000001`,9.295608827842757`*^-6},{0.08505000000000001`,9.295608827842757`*^-6},{0.09005000000000002`,9.295608827842757`*^-6},{0.09505`,9.29560882784144`*^-6},{0.10005`,9.29560882784144`*^-6}}

ListLogPlot[{dat0,dat1},PlotRange->{0,All}]

ListLogPlot[{τ,dat0}]

Plot[0.000001*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.005,0.1}]

p2=Show[Plot[0.000001*((Exp[((-ℏ ω0)/(kb T))])/(Sqrt[T])),{T,0.005,0.1}],ListPlot[dat0]]