(*correction to the free Hamiltonian for a chain Hamiltonian for a chain with two ends*)
h1matchain[tl_,tr_,ϵ0_,t0_,ϵvec_,τvec_]:=Module[{ret,NT,i},
NT=Length[ϵvec];
ret = DiagonalMatrix[Join[{0},ϵvec-Table[ϵ0,{i,1,NT}],{0}]];
ret[[1,2]]=tl-t0;
ret[[2,1]]=tl-t0;
ret[[NT+1,NT+2]]=tr-t0;
ret[[NT+2,NT+1]]=tr-t0;
For[i=1,i<NT,++i,
ret[[i+1,i+2]]=τvec[[i]]-t0;
ret[[i+2,i+1]]=τvec[[i]]-t0;
];
ret
];
g0matcrchain[ϵ0_,t0_,ω_,NT_]:=Table[-I/Sqrt[4 t0^2 -(ω-ϵ0)^2]((ω-ϵ0)/(2 t0)-I Sqrt[1-(ω-ϵ0)^2/(4 t0^2)])^Abs[i-j],{i,1,NT+2},{j,1,NT+2}];
gfmatcrchain[tl_,tr_,ϵ0_,t0_,ϵvec_,τvec_,ω_]:=Inverse[IdentityMatrix[Length[ϵvec]+2]-g0matcrchain[ϵ0,t0,ω,Length[ϵvec]].h1matchain[tl,tr,ϵ0,t0,ϵvec,τvec]].g0matcrchain[ϵ0,t0,ω,Length[ϵvec]];
transmissionchain[tl_,tr_,ϵ0_,t0_,ϵvec_,τvec_,ω_]:=Which[ω>ϵ0-2t0 && ω<ϵ0+2t0,Abs[gfmatcrchain[tl,tr,ϵ0,t0,ϵvec,τvec,ω][[2,Length[ϵvec]+1]]]^2 Abs[tl]^2/t0^2 Sqrt[4t0^2-(ω-ϵ0)^2]Abs[tr]^2/t0^2 Sqrt[4t0^2-(ω-ϵ0)^2],True,0.0];

ListPlot[Monitor[Table[{w,transmissionchain[Sqrt[0.1]/Sqrt[2],Sqrt[0.1]/Sqrt[2],0,1.0*0.1,{0,0,0,0,0,0,0,0,0,0,0,0,0},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`,0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}/0.02,w]},{w,-1,+1,0.01}],w],PlotRange->All]

Sort[Eigenvalues[Take[h1matchain[tl,tr,0,0,{0,0,0,0,0,0,0,0,0,0,0,0,0},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`,0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}/0.02],{9,14},{9,14}]]]

Export["test.dat",Monitor[Table[{w,transmissionchain[Sqrt[100]/Sqrt[2],Sqrt[100]/Sqrt[2],0,1.0*100.,{0,0,0,0,0,0,0,0,0,0,0,0,0},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`,0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}/0.02,w]},{w,-1,+1,0.01}],w]]

Reverse[{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`,0.007269560988958051`,0.007238948047382406`}]

(*Tilmanns parameter:
(Subscript[β, L])=(Subscript[β, R])=β=1
Spektrale Kopplungsstärke Γ=0.02 (wideband)
ϵ=0
Von außen nach innen in der Kette
(Subscript[λ, L])=0.009999959154289148`
(Subscript[λ, 2L])=0.00781394527921633`
(Subscript[λ, 3L])=0.0074508497062710805`
(Subscript[λ, 4L])= 0.007326447541946075`
(Subscript[λ, 5L])=0.007269560988958051`
(Subscript[λ, 6L])=0.007238948047382406`
*)
(*
Anpassung wideband limes:
tl=tr=t und t0 sehr groß mit
2 t^2/t0=0.02 
*)
(* optimale Parameter für direkte Anpassung der SD
(Subscript[λ, L])=0.44798874545216505` Γ0
(Subscript[λ, 2L])=0.3720949697259449` Γ0
(Subscript[λ, 3L])=0.3666741914433289` Γ0
(Subscript[λ, 4L])=0.36586118013042696` Γ0
(Subscript[λ, 5L])=0.3656755715952492` Γ0
(Subscript[λ, 6L])=0.365617222323006` Γ0

*)

transmissionchain[1.,1.,0.,100.,{0.},{},1.]

(*1. Optimierung -> Arbeit Tilmann, Korrektur am 7.6.*)
Plot[{transmissionchain[1.,1.,0.,100.,{0.},{},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.},{0.009999959154289148`,0.009999959154289148`},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.009999959154289148`},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`},w]
},{w,-0.02,+0.02},PlotRange->All,PlotLabels->{"N=1","N=3","N=5","N=7","N=9","N=11","N=13"}]
(*2. Optimierung*)
Plot[{transmissionchain[1.,1.,0.,100.,{0.},{},w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.},{0.44798874545216505`,0.44798874545216505`}*0.02,w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3720949697259449`,0.44798874545216505`}*0.02,w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.36586118013042696`,0.36586118013042696`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.36586118013042696`,0.3656755715952492`,0.3656755715952492`,0.36586118013042696`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w],
transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.36586118013042696`,0.3656755715952492`,0.365617222323006`,0.365617222323006`,0.3656755715952492`,0.36586118013042696`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w]},{w,-0.02,+0.02},PlotRange->All,PlotLabels->{"N=1","N=3","N=5","N=7","N=9","N=11","N=13"}]

(*direkte Optimierung*)
diff[wtg_?NumericQ,tl_?NumericQ,tr_?NumericQ,ϵ0_?NumericQ,t0_?NumericQ,ϵvec_?(VectorQ[#,NumericQ]&),τvec_?(VectorQ[#,NumericQ]&)]:=NIntegrate[(HeavisideTheta[wtg^2-w^2]-transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w])^2,{w,-Infinity,-wtg,0,+wtg,+Infinity},MaxRecursion->20];

diff[0.75*0.02,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`, 0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.009999959154289148`}]

NMinimize[diff[0.75*0.02,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{λ6,λ5,λ4,λ3,λ2,λ1,λ1,λ2,λ3,λ4,λ5,λ6}],{λ1,λ2,λ3,λ4,λ5,λ6}]

(*compare same level transmissions: they both approximate the rectangular limit, the first scheme is better at the edges, the second better in the center.*)
Plot[{transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`, 0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.009999959154289148`},w],transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.36586118013042696`,0.3656755715952492`,0.365617222323006`,0.365617222323006`,0.3656755715952492`,0.36586118013042696`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w],transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.0099609897622874`,0.03163361777640106`,1.066290253346536`,0.4286598205548415`,1.1612881965805946`,3.0474653711824247`,3.0474653711824247`,1.1612881965805946`,0.4286598205548415`,1.066290253346536`,0.03163361777640106`,0.0099609897622874`},w]

},{w,-0.03,+0.03},PlotRange->All,PlotLabels->{"first","second","third"}]

Export["test.dat",Table[{w/0.02,transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`, 0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.009999959154289148`},w]},{w,-0.02,+0.02,0.0001}]]

Export["test.dat",Table[{w/0.02,transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.44798874545216505`,0.3720949697259449`,0.3666741914433289`,0.36586118013042696`,0.3656755715952492`,0.365617222323006`,0.365617222323006`,0.3656755715952492`,0.36586118013042696`,0.3666741914433289`,0.3720949697259449`,0.44798874545216505`}*0.02,w]},{w,-0.02,+0.02,0.0001}]]

Sort[Eigenvalues[Take[h1matchain[0.,0.,0.,0.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`, 0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.009999959154289148`}],{2,14},{2,14}]]]

Plot[transmissionchain[1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.0074508497062710805`, 0.007326447541946075`, 0.007269560988958051`,0.007238948047382406`,0.007238948047382406`,0.007269560988958051`,0.007326447541946075`,0.0074508497062710805`,0.009999959154289148`},w],{w,-0.02,+0.02},PlotRange->All]

(*for simplicity, we hardcode the some additional sampling points manually by throwing 40 points over the interval plotted above*)
entprod[β_,V_,tl_,tr_,ϵ0_,t0_,ϵvec_,τvec_]:=Module[{ret,cur},
cur=NIntegrate[transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w]/(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)-1/(Exp[β(w-(-V/2))]+1)),Join[{w,ϵ0-2t0},Table[x,{x,-0.02,+0.02,0.001}],{ϵ0+2 t0}],MaxRecursion->20];
β cur V
];
uncrel[β_,V_,tl_,tr_,ϵ0_,t0_,ϵvec_,τvec_]:=Module[{ret,cur,noise},
cur=NIntegrate[transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w]/(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)-1/(Exp[β(w-(-V/2))]+1)),Join[{w,ϵ0-2t0},Table[x,{x,-0.02,+0.02,0.001}],{ϵ0+2 t0}],MaxRecursion->20];
noise= NIntegrate[transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w]/(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)*(1-1/(Exp[β(w-(+V/2))]+1))+1/(Exp[β(w-(-V/2))]+1)*(1-1/(Exp[β(w-(-V/2))]+1))),Join[{w,ϵ0-2t0},Table[x,{x,-0.02,+0.02,0.001}],{ϵ0+2 t0}],MaxRecursion->20]+NIntegrate[transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w]*(1-transmissionchain[tl,tr,ϵ0,t0,ϵvec,τvec,w])/(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)-1/(Exp[β(w-(-V/2))]+1))^2,Join[{w,ϵ0-2t0},Table[x,{x,-0.02,+0.02,0.001}],{ϵ0+2 t0}],MaxRecursion->20];
β V noise/cur
];
uncrelid[β_,V_,wmin_,wmax_]:=Module[{ret,cur,noise},
cur=NIntegrate[1./(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)-1/(Exp[β(w-(-V/2))]+1)),{w,wmin,wmax},MaxRecursion->20];
noise= NIntegrate[1./(2 Pi)*(1/(Exp[β(w-(+V/2))]+1)*(1-1/(Exp[β(w-(+V/2))]+1))+1/(Exp[β(w-(-V/2))]+1)*(1-1/(Exp[β(w-(-V/2))]+1))),{w,wmin,wmax},MaxRecursion->20];
β V noise/cur
];

Export["test.dat",Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.0099609897622874`,0.03163361777640106`,1.066290253346536`,0.4286598205548415`,1.1612881965805946`,3.0474653711824247`,3.0474653711824247`,1.1612881965805946`,0.4286598205548415`,1.066290253346536`,0.03163361777640106`,0.0099609897622874`}]},{V,-30,+30,0.2}],V]]

Export["test.dat",Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051,0.007238948047382406,0.007238948047382406,0.007269560988958051,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]},{V,-30,+30,0.2}],V]]

Export["test.dat",Monitor[Table[{V,uncrelid[1.,V,-0.015,+0.015]},{V,-30,+30,0.1}],V]]

(*we do not exactly have the wideband limit with Γ(ω) = 2 τ^2/τ0 Sqrt[1-(ω-ϵ0)^2/(4 τ0^2)] but for large τ0 it works, which roughly corresponds to β Γ = 1.0 *)
tab1=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.},{}]},{V,-30.001,+30,1.}],V];
tab3=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.},{0.009999959154289148`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];
tab5=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];
tab7=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];
tab9=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];
tab11=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051,0.007269560988958051,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];
tab13=Monitor[Table[{V,uncrel[1.,V,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051,0.007238948047382406,0.007238948047382406,0.007269560988958051,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]},{V,-30.001,+30,1.}],V];

tabor1={{-30,1},{+30,1}};
tabor2={{-30,2},{+30,2}};

ListPlot[{tab1,tab3,tab5,tab7,tab9,tab11,tab13,tabor1,tabor2},Joined->True,PlotRange->{0,3},PlotLabels->{"N=1","N=3","N=5","N=7","N=9","N=11","N=13"}]

(*Tilmanns Ergebnisse vs. FindMinimum:
Die jeweiligen Minima der Unschärferelation mit den dazugehörigen Parametern (3,5,7,9,11,13) Quantenpunkte
{1.2458788244725605`,{λ->0.009999959154289148`,V->-7.221064554765229`}} vs. {1.2458788229026372`,{V->7.221067702630215`,λ->0.010000199061175588`}}
{0.7554305430129272`,{λ->-0.00781394527921633`,V->9.393004827858011`}} vs. {0.7558164712527939`,{V->9.391065913646672`,λ->0.007813725282389295`}}
{0.572847683176261`,{λ->0.0074508497062710805`,V->10.386440207427597`}} vs. {0.5726255406381631`,{V->10.38777281253604`,λ->0.007451002109568454`}}
{0.48965280097811276`,{λ->0.007326447541946075`,V->10.91473000281113`}} vs. {0.48959177828867484`,{V->10.915141270700905`,λ->0.007326444792376225`}}
{0.44564368710726454`,{λ->0.007269560988958051`,V->11.22239888748624`}} vs. {0.4455816929554629`,{V->11.222846950721534`,λ->0.007269559371554312`}}
{0.4197865144938278`,{λ->-0.007238948047382406`,V->11.41449491980886`}} vs. {0.4197209743172402`,{V->11.41499238569322`,λ->0.00723894298148576`}}
*)
uncrel[1.,7.221064554765229,1.,1.,0.,100.,{0.,0.,0.},{0.009999959154289148`,0.009999959154289148`}]
uncrel[1.,9.393004827858011,1.,1.,0.,100.,{0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.00781394527921633`,0.009999959154289148`}]
uncrel[1.,10.386440207427597,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]
uncrel[1.,10.91473000281113,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]
uncrel[1.,11.22239888748624,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051,0.007269560988958051,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]
uncrel[1.,11.41449491980886,1.,1.,0.,100.,{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},{0.009999959154289148`,0.00781394527921633`,0.0074508497062710805`,0.007326447541946075`,0.007269560988958051,0.007238948047382406,0.007238948047382406,0.007269560988958051,0.007326447541946075`,0.0074508497062710805`,0.00781394527921633`,0.009999959154289148`}]