fibo[n_]:=
If[n==1,1,
If[n==2,2,fibo[n-1]+fibo[n-2]]];

lfibo[n_]:=Table[{fibo[k],0},{k,1,n}];

Tight-Binding model, nearest neighbors, coupling param following Fibo sequence

Hamiltonian

order=9;
va=1.;
vb=-0.5;
taille = fibo[order]

lcoupl[n_]:=
If[n==1,{va},
If[n==2,{va,vb},Flatten[{lcoupl[n-1],lcoupl[n-2]}]]];

coupl[i_]=(lcoupl[order])[[i]];

lc=lcoupl[order];
HamFiboTB=Reap[
Sow[{{1,taille},lc[[taille]]}];
Sow[{{taille,1},lc[[taille]]}];
Do[{
Sow[{{i,i+1},lc[[i]]}];
Sow[{{i+1,i},lc[[i]]}];

},{i,1,taille-1}]][[2,1]];

tmp=Split[Sort[HamFiboTB],#1[[1]]==#2[[1]]&];
HamFiboTB=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];

lc=lcoupl[order];
HamFiboTBOBC=Reap[
Do[{
Sow[{{i,i+1},lc[[i]]}];
Sow[{{i+1,i},lc[[i]]}];

},{i,1,taille-1}]][[2,1]];

tmp=Split[Sort[HamFiboTBOBC],#1[[1]]==#2[[1]]&];
HamFiboTBOBC=Table[tmp[[i,1,1]]->Tr[Transpose[tmp[[i]]][[2]]],{i,1,Length[tmp]}];

spectre=Sort[Eigenvalues[SparseArray[HamFiboTB]]];

spectreOBC=Sort[Eigenvalues[SparseArray[HamFiboTBOBC]]];

Show[ListPlot[Sort[spectre]],ListPlot[Sort[spectreOBC],PlotStyle->Red]]

IDOS

heav[x_]:=If[x<0,0,1];

idos[x_]:=Sum[heav[x-spectre[[i]]],{i,1,taille}];

Plot[idos[x],{x,-2,2}]

Wave packet

eigvect=Sort[Eigenvectors[SparseArray[HamFiboTB]]];

eigvectOBC=Sort[Eigenvectors[SparseArray[HamFiboTBOBC]]];

ListPlot[Psi0]

Initially spatially localized WP

Psi=Table[Chop[Sum[eigvect[[n]][[ki]]*eigvect[[n]][[k]]*Exp[-I*spectre[[n]]*t],{n,1,taille}]],{k,1,taille}];

Show[ListPlot[Chop[Abs[Psi/.t-> 0.0]^2]],ListPlot[Chop[Abs[Psi/.t-> 0.1]^2],PlotStyle->Red],ListPlot[Chop[Abs[Psi/.t-> 0.2]^2],PlotStyle->Green],ListPlot[Chop[Abs[Psi/.t->1000000000.]^2],PlotStyle->Black]]

ListPlot[Chop[Abs[Psi/.{ki->10,t-> 0.0}]^2]]

(Psi/.{ki->45,t-> 0.0})[[45]]

Initially uniform. spectr. distrib. WP

PsiUni=Table[Chop[Sum[eigvect[[n]][[k]]*Exp[-I*spectre[[n]]*t],{n,1,taille}]],{k,1,taille}];

ListPlot[Chop[PsiUni/.t-> 0.0]]

RMS

Initially localized WP

rmssq=Sum[(k-ki)^2*Chop[Abs[Psi[[k]]]^2],{k,1,taille}];

rmssq115=rmssq/.ki->115;

DiscrListRms115=Table[{n/10.,Chop[rmssq115/.t->n/10.]},{n,1,1000}];

DiscrListRms45v2=Table[{n/100.,Chop[rmssq45/.t->n/100.]},{n,1,1000}];

ListPlot[Transpose[{Transpose[DiscrListRms115][[1]],Transpose[DiscrListRms115][[2]]}]]

ListLogPlot[Transpose[{Transpose[DiscrListRms45v2][[1]],Transpose[DiscrListRms45v2][[2]]}]]

Spatial averaging for init. loc. WP

rms=Reap[
Do[{psikii=Psi/.ki->kii;Sow[Sum[(k-kii)^2*Chop[Abs[psikii[[k]]]^2],{k,1,taille}]]},{kii,1,taille}]][[2,1]];

rmsav=1/taille*Sum[rms[[k]],{k,1,taille}];

rmsav/.t->0.

Chop[rmsav/.t->0.6]

DiscrListRmsAv=Table[{n/10.,Chop[rmsav/.t->n/10.]},{n,1,100}];

DiscrListRmsAvV2=Table[{10+n/10.,Chop[rmsav/.t->10+n/10.]},{n,1,100}];

DiscrListRmsAvV3=Table[{n/10.,Chop[rmsav/.t->n/10.]},{n,1,500}];

Export["DiscrListRmsAvV3.dat",DiscrListRmsAvV3];

ListPlot[DiscrListRmsAv]

lrms=Table[{k,rms[[k]]},{k,1,taille}];

ListPlot[DiscrListRmsAvV2]

ListPlot[DiscrListRmsAvV3]

For uniformally spectr. distrib. init. WP

rmssqUni=Sum[(k-ki)^2*Chop[Abs[PsiUni[[k]]]^2],{k,1,taille}];

rmssqUni115=rmssqUni/.ki->115;

Clear[DiscrListRmsUni45]

DiscrListRmsUni115=Table[{n/10.,Chop[rmssqUni115/.t->n/10.]},{n,1,1000}];

DiscrListRmsUni115V2=Table[{n/1000.,Chop[rmssqUni115/.t->n/10.]},{n,1,1000}];

ListPlot[Transpose[{Transpose[DiscrListRmsUni115][[1]],Transpose[DiscrListRmsUni115][[2]]}]]

Show[ListPlot[Transpose[{Transpose[DiscrListRms115][[1]],Transpose[DiscrListRms115][[2]]}],PlotStyle->Red],ListPlot[Transpose[{Transpose[DiscrListRmsUni115][[1]],10^(-3)*Transpose[DiscrListRmsUni115][[2]]}]]]

ListLogLogPlot[Transpose[{Transpose[DiscrListRmsUni115][[1]],Transpose[DiscrListRmsUni115][[2]]}]]

ListPlot[Transpose[{Transpose[DiscrListRmsUni115V2][[1]],Transpose[DiscrListRmsUni115V2][[2]]}]]

Time-averaging, init. spectr. unif. distrib

mom2=Table[{n,nn,Sum[(k-ki)^2*eigvect[[n]][[k]]*eigvect[[nn]][[k]],{k,1,taille}]},{n,1,taille},{nn,1,taille}];

mom245=mom2/.ki->45;

ListPlot[Table[{n,mom245[[1]][[n]][[3]]},{n,1,taille}]]

ListPlot[Table[{n,mom245[[45]][[n]][[3]]},{n,1,taille}]]

rmsUniTA=Sum[Sinc[(spectre[[n]]-spectre[[nn]])t]*mom2[[n]][[nn]][[3]],{n,1,taille},{nn,1,taille}];

rmsUniTA45=rmsUniTA/.ki->45;

DiscrListRmsUniTA45=Table[{n/10.,Chop[rmsUniTA45/.t->n/10.]},{n,1,1000}];

DiscrListRmsUniTAPow45=Table[{n/10.,(n/10.)^(Log[2]/Log[3])Chop[rmsUniTA45/.t->n/10.]},{n,1,1000}];

DiscrListRmsUniTAPow45v2=Table[{n/10.,(n/10.)^(Log[2]/(2*Log[3]))Chop[rmsUniTA45/.t->n/10.]},{n,1,1000}];

ListLogLogPlot[Transpose[{Transpose[DiscrListRmsUniTA45][[1]],Transpose[DiscrListRmsUniTA45][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrListRmsUniTAPow45v2][[1]],Transpose[DiscrListRmsUniTAPow45v2][[2]]}]]

Time-averaging, init. loc

mom2loc=Table[{n,nn,Sum[(k-ki)^2*eigvect[[n]][[ki]]*eigvect[[nn]][[ki]]*eigvect[[n]][[k]]*eigvect[[nn]][[k]],{k,1,taille}]},{n,1,taille},{nn,1,taille}];

distribeig=Table[{n,eigvect[[n]][[ki]]},{n,1,taille}];

distribeigOBC=Table[{n,eigvectOBC[[n]][[ki]]},{n,1,taille}];

distribeigV2=Table[eigvect[[n]][[ki]],{n,1,taille}];

distribeigV2OBC=Table[eigvectOBC[[n]][[ki]],{n,1,taille}];

ldistribeig=Reap[Do[Sow[distribeig/.ki->n],{n,1,taille}]][[2,1]];

ldistribeigOBC=Reap[Do[Sow[distribeigOBC/.ki->n],{n,1,taille}]][[2,1]];

Show[Table[ListPlot[ldistribeig[[n]]],{n,1,taille}]]

Show[Table[ListPlot[ldistribeigOBC[[n]]],{n,1,taille}]]

ListPlot[Sort[ldistribeigOBC[[55]]]]

ldistribeig[[1]]

Sort[distribeigV2/.ki->1]

Sort[distribeigV2OBC/.ki->25]

mom2loc25=mom2loc/.ki->25;

rmsTA=Sum[Sinc[(spectre[[n]]-spectre[[nn]])t]*mom2loc[[n]][[nn]][[3]],{n,1,taille},{nn,1,taille}];

rmsTA25=rmsTA/.ki->25;

DiscrListRmsTA25=Table[{n/10.,Chop[rmsTA25/.t->n/10.]},{n,1,1000}];

DiscrListRmsTAPow25=Table[{n/10.,(n/10.)^(-Log[2]/Log[3])Chop[rmsTA25/.t->n/10.]},{n,1,1000}];

ListLogLogPlot[Transpose[{Transpose[DiscrListRmsTA25][[1]],Transpose[DiscrListRmsTA25][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrListRmsTAPow25][[1]],Transpose[DiscrListRmsTAPow25][[2]]}]]

Time-average return probability

Init. uniform. spectr. distrib. WP

RetProbUni=1/taille^2*Sum[Sinc[(spectre[[k]]-spectre[[j]])t],{k,1,taille},{j,1,taille}];

DiscrRetProbUni=Table[{n/10.,Chop[N[RetProbUni/.t->n/10.]]},{n,1,1000}];

DiscrRetProbUniT3=Table[{n/10.,Chop[N[RetProbUni/.t->3*n/10.]]},{n,1,1000}];

DiscrRetProbUniPow=Table[{n/10.,(n/10.)^N[Log[2]/Log[3]]*Chop[N[RetProbUni/.t->n/10.]]},{n,1,1000}];

DiscrRetProbUniPowV2=Table[{n/10.,(n/10.)^N[Log[2]/Log[3]]*Chop[N[RetProbUni/.t->n/10.]]},{n,1,5000}];

DiscrRetProbUniPowV3=Table[{n/10.,(n/10.)^N[Log[2]/Log[1+Sqrt[5]]]*Chop[N[RetProbUni/.t->n/10.]]},{n,1,1000}];

N[Log[2]/Log[3]]

N[Log[2]/Log[1+Sqrt[5]]]

Show[ListLogLogPlot[Transpose[{Transpose[DiscrRetProbUni][[1]],Transpose[DiscrRetProbUni][[2]]}]],ListLogLogPlot[Transpose[{Transpose[DiscrRetProbUniT3][[1]],Transpose[DiscrRetProbUniT3][[2]]}],PlotStyle->Red]]

Show[ListPlot[Transpose[{Transpose[DiscrRetProbUni][[1]],Transpose[DiscrRetProbUni][[2]]}]],ListPlot[Transpose[{Transpose[DiscrRetProbUniT3][[1]],Transpose[DiscrRetProbUniT3][[2]]}],PlotStyle->Red]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProbUniPowV2][[1]],Transpose[DiscrRetProbUniPowV2][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProbUniPowV3][[1]],Transpose[DiscrRetProbUniPowV3][[2]]}]]

Init. loc. WP

Clear[DiscrRetProb]

RetProb=Sum[Abs[eigvect[[k]][[k0]]]^2*Abs[eigvect[[j]][[k0]]]^2*Sinc[(spectre[[k]]-spectre[[j]])t],{k,1,taille},{j,1,taille}];

RetProb45=RetProb/.k0->45;

RetProb1=RetProb/.k0->1;

RetProb89=RetProb/.k0->89;

DiscrRetProb45=Table[{n/10.,Chop[N[RetProb45/.t->n/10.]]},{n,1,1000}];

DiscrRetProbPow45=Table[{n/10.,(n/10.)^N[Log[2]/Log[3]]*Chop[N[RetProb45/.t->n/10.]]},{n,1,1000}];

DiscrRetProb1=Table[{n/10.,Chop[N[RetProb1/.t->n/10.]]},{n,1,1000}];

DiscrRetProbPow1=Table[{n/10.,(n/10.)^N[Log[2]/Log[3]]*Chop[N[RetProb1/.t->n/10.]]},{n,1,1000}];

DiscrRetProb89=Table[{n/10.,Chop[N[RetProb89/.t->n/10.]]},{n,1,1000}];

ListLogLogPlot[Transpose[{Transpose[DiscrRetProb45][[1]],Transpose[DiscrRetProb45][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProbPow45][[1]],Transpose[DiscrRetProbPow45][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProb1][[1]],Transpose[DiscrRetProb1][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProbPow1][[1]],Transpose[DiscrRetProbPow1][[2]]}]]

ListLogLogPlot[Transpose[{Transpose[DiscrRetProb89][[1]],Transpose[DiscrRetProb89][[2]]}]]

Participation ratio

Init. loc. WP

ipr=Sum[Abs[Psi[[k]]]^4,{k,1,taille}];

Sum[Abs[(Psi/.{ki->1,t->0.})[[k]]]^2,{k,1,taille}]

ipr22=ipr/.ki->22;

DiscrIpr22=Table[{n/10.,Chop[ipr22/.t->(n/10.)]},{n,1,1000}];

ListLinePlot[Transpose[{Transpose[DiscrIpr22][[1]],Transpose[DiscrIpr22][[2]]}],PlotRange->{{0,100},{0,1}}]

ipr22/.t->0.

Init. uniform. spectr. distrib. WP

PsiUniNorm=PsiUni/Sqrt[55];

iprUni=Sum[Abs[PsiUniNorm[[k]]]^4,{k,1,taille}];

Sum[Abs[(PsiUniNorm/.t->0.)[[k]]]^2,{k,1,taille}]

DiscrIprUni=Table[{n/10.,Chop[iprUni/.t->(n/10.)]},{n,1,1000}];

ListLinePlot[DiscrIprUni,PlotRange->{{0,100},{0,1}}]

For eigenvectors

ipreig=Table[{n,Sum[Abs[eigvect[[n]][[k]]]^4,{k,1,taille}]},{n,1,taille}];

Show[ListPlot[ipreig],ListPlot[lfibo[9],PlotStyle->Red]]