Large N Sachdev-Ye-Kitaev model

Example of working with Fourier transforms and iterative algorithms.

SetDirectory@NotebookDirectory[];

Schwinger-Dyson equation solver

The phase change defined below is necessary due to Mathematica’s Fourier/InverseFourier conventions.
One needs to find it by hand, and there is some ambiguity related to the choice of Euclidean-time discretization.
In this case we choose to discretize τ∈[0,β) using (Subscript[τ, j])=(β/(2K))(j-1) with for j=1,2,...

ClearAll@SDSolve;
SDSolve[k_Integer,α_Real,ϵ_Real][β_Real]:=
Block[{phase=β Table[(((-1))^(j-1))(E^(I π(j-1)/(2k))),{j,2k}],fourier,invfourier,ωs,G,Gt,S,St,dist,cur=1},
invfourier=(phase^(-1))InverseFourier[#,FourierParameters->{-1, 1}]&;
fourier=Fourier[# phase,FourierParameters->{-1, 1}]&;
ωs=((2π)/β)(Range[-k+1,k]-(1/2));
Gt[0]=(((-I ωs))^(-1));
G[n_]:=G[n]=invfourier@Gt[n];
S[n_]:=S[n]=((G[n])^3);
St[n_]:=St[n]=fourier@S[n];
Gt[n_]:=Gt[n]=(1-α)Gt[n-1]+α(((-I ωs-St[n-1]))^(-1));
dist[n_]:=dist[n]=Mean@Abs[Gt[n]-Gt[n-1]];
While[dist[cur]>ϵ,cur++];
Gt[cur]
]

Plots

In order to plot the solutions in Euclidean-time, it is convenient to add and subtract the free-fermion solution to avoid Fourier-transforming discontinuous functions

ClearAll@PlotSolution;
PlotSolution[k_Integer,α_Real,ϵ_Real][β_Real][color_]:=
Block[{phase=β Table[(((-1))^(j-1))(E^(I π(j-1)/(2k))),{j,2k}],fourier,invfourier,ts,ωs},
invfourier=(phase^(-1))InverseFourier[#,FourierParameters->{-1, 1}]&;
fourier=Fourier[# phase,FourierParameters->{-1, 1}]&;
ωs=((2π)/β)(Range[-k+1,k]-(1/2));
ts=(β/(2k))Range[0,2k-1];
ListPlot[Transpose@{(ts/β),invfourier[SDSolve[k,α,ϵ][β]-(((-I ωs))^(-1))]+(1/2)},Joined->True,PlotRange->All,AxesLabel->{"τ/β","G(τ)"},AxesStyle->Black,TicksStyle->Medium,ImageSize->1024,PlotStyle->color,PlotLegends->LineLegend[{color},{"β = "<>ToString@β}],AxesStyle->Black,TicksStyle->Medium]
]

Show[{
PlotSolution[(2^16),.5,(10.^(-9))][5.][Blue],
PlotSolution[(2^16),.5,(10.^(-9))][10.][Orange],
PlotSolution[(2^16),.5,(10.^(-9))][100.][Red]
},PlotRange->{0,.5},AxesOrigin->{0,0}
]

Export["SolutionExample.pdf",%]

Thermodynamics

Computing the entropy is straightforward...

SYKEntropy[k_Integer,α_Real,ϵ_Real][β_Real]:=
Block[{ωs,Gt=SDSolve[k,α,ϵ][β]},
ωs=((2π)/β)(Range[-k+1,k]-(1/2));
(1/2)Log@2-Sum[Log@Abs[-I ωs[[i]]Gt[[i]]]+(5/4)(1-Re[-I ωs[[i]]Gt[[i]]]),{i,k+1,2k}]
]

We can extrapolate linearly to T->0, finding good agreement

LinearModelFit[Table[{(β^(-1)),SYKEntropy[(2^16),.5,(10.^(-9))][β]},{β,{100.,110.}}],T,T][0]

(Catalan/(2π))+((Log@2)/8)//N

When computing thermodynamic functions, one should be careful about truncation errors becoming very important at low-temperatures...

SYKEntropy[(2^16),.1,(10.^(-9))][200.]

SYKEntropy[(2^16),.1,(10.^(-9))][500.]

SYKEntropy[(2^16),.1,(10.^(-9))][1000.]

Optional parts

A solver with arbitrary starting point and a mass term

ClearAll@SDSolve2;
SDSolve2[k_Integer,α_Real,ϵ_Real,Gt0_][m_Real,β_Real]:=
Block[{phase=β Table[(((-1))^(j-1))(E^(I π(j-1)/(2k))),{j,2k}],fourier,invfourier,ωs,G,Gt,S,St,dist,cur=1},
invfourier=(phase^(-1))InverseFourier[#,FourierParameters->{-1, 1}]&;
fourier=Fourier[# phase,FourierParameters->{-1, 1}]&;
ωs=((2π)/β)(Range[-k+1,k]-(1/2));
Gt[0]=Gt0;
G[n_]:=G[n]=invfourier@Gt[n];
S[n_]:=S[n]=((G[n])^2)Reverse@G[n];(* ugly but should work! *)
(*S[n_]:=S[n]=-((G[n])^2)Prepend[Reverse@G[n][[2;;]],First@G[n]];*)
St[n_]:=St[n]=fourier@S[n];
Gt[n_]:=Gt[n]=(1-α)Gt[n-1]+α(((m-I ωs-St[n-1]))^(-1));
dist[n_]:=dist[n]=Mean@Abs[Gt[n]-Gt[n-1]];
While[dist[cur]>ϵ,cur++];
Gt[cur]
]

The high-temperature starting point

ClearAll@TreeLevelStart;
TreeLevelStart[k_Integer][m_Real,β_Real]:=With[{ωs=((2π)/β)(Range[-k+1,k]-(1/2))},(((m-I ωs))^(-1))];

Non-vanishing-mass solutions are asymmetric

ClearAll@PlotSolution3;
PlotSolution3[m_Real,β_Real][Gt_][color_]:=
Block[{k=(1/2)Length@Gt,Gt0,phase,fourier,invfourier,ts},
phase=β Table[(((-1))^(j-1))(E^(I π(j-1)/(2k))),{j,2k}];
invfourier=(phase^(-1))InverseFourier[#,FourierParameters->{-1, 1}]&;
fourier=Fourier[# phase,FourierParameters->{-1, 1}]&;
Gt0=TreeLevelStart[k][m,β];
ts=(β/(2k))Range[0,2k-1];
ListPlot[Transpose@{(ts/β),invfourier[Gt-Gt0]+((E^(m(β-ts)))/((E^(β m))+1))},Joined->True,PlotRange->All,AxesLabel->{"τ/β","G(τ)"},AxesStyle->Black,TicksStyle->Medium,ImageSize->1024,PlotStyle->color,PlotLegends->LineLegend[{color},{"β = "<>ToString@β}],AxesStyle->Black,TicksStyle->Medium]
]

With[{k=(2^14),m=0.1,β=20.},
SDSolve2[k,.25,(10.^(-9)),TreeLevelStart[k][m,β]][m,β];
]

PlotSolution3[.1,20.][%][Black]

Thermodynamics

ClearAll@SYKEntropy2;
SYKEntropy2[m_Real,β_Real][Gt_]:=
Block[{k=(1/2)Length@Gt,Gt0},
Gt0=TreeLevelStart[k][m,β];
Log[1+(E^(β m))]-((β m)/((E^(-β m))+1))-Sum[2m Re[Gt[[i]]-Gt0[[i]]]+2Log@Abs[((Gt[[i]])/(Gt0[[i]]))]+(5/2)(1-Re[((Gt[[i]])/(Gt0[[i]]))]),{i,k+1,2k}]
]

With[{k=(2^14),m=0.,β=20.},
SDSolve2[k,.25,(10.^(-9)),TreeLevelStart[k][m,β]][m,β];
]

SYKEntropy2[0.,20.][%]

2SYKEntropy[(2^16),.25,(10.^(-9))][20.]

Dynamic@{m,cur,dist[cur]}

IncMassData=Block[{k=(2^14),β=20.,sol},
sol=TreeLevelStart[k][0.,β];
Table[
sol=SDSolve2[k,.25,(10.^(-9)),sol][m,β];
{m,SYKEntropy2[m,β][sol]}
,{m,0.,.6,.01}]
]

DecMassData=Block[{k=(2^14),β=20.,sol},
sol=TreeLevelStart[k][0.,β];
Table[
sol=SDSolve2[k,.25,(10.^(-9)),sol][m,β];
{m,SYKEntropy2[m,β][sol]}
,{m,.6,0.,-.01}]
]

ListPlot[{IncMassData,DecMassData},Joined->True]