% (c) 2017 Aavishkar A. Patel
% This code solves the imaginary time Dyson equations iteratively for the c-f-SYK model
% on a discrete-time grid, by Fourier transforming back and forth between
% time and frequency representations of the two point functions.
% Run, and hold down return until the plot converges.
% You can choose what to plot where the plot command appears.
% Gt is the time-domain (-\beta to \beta) SYK
% Green's function and Gct is the same for Gc. S is a shorthand for \Sigma,
% the self-energy, and w is a shorthand for \omega.
% Values of couplings, bandwidth, number of grid points,
% and other parameters in the model. Also sets the linear update weight x.
J = 100.0; T = 1.0; mu = 0.0; itermax = 1000; points = 2^20; x = 0.05;
g = 80.0; Lambda = 3.0; muc = 0.0*Lambda; nu0 = 2.0*pi/Lambda; MbyN = 0.5;
% Appropriate time and frequency steps. Set up arrays for time and
% frequency domains. The time domain array is directly -\beta to \beta,
% left to right. The frequency domain array stores positive frequencies
% in the left half, going left to right in increasing magnitude,
% and negative frequencies on the right end, going right to left in increasing
% magnitude.
tstep = 1/(points*T);
wstep = pi/points;
w1 = zeros(1,2*points-1);
w1(1:points)=wstep*(1:points);
w1(points:2*points-1) = wstep*((points:2*points-1)-2*points);
% Initialize green's functions and self energies in both representations.
% When the initialization is nonzero, it contains only fermionic Matsubara
% frequencies.
Gt = 0.5*sign(points+0.5-(1:2*points));
Gw = -conj(0.5*tstep*fft(Gt));
St = zeros(1,2*points);
Sw = zeros(1,2*points);
Gct = 0.5*sign(points+0.5-(1:2*points));
Gcw = -conj(0.5*tstep*fft(Gct));
Sct = zeros(1,2*points);
Scw = zeros(1,2*points);
% Begin iterations
for n=1:itermax
% SYK Part
St = -J^2*(Gt.^2).*fliplr(Gt) - MbyN*g^2*Gt.*Gct.*(fliplr(Gct)); % Flipping the time domain array left to right implements \tau --> -\tau
Sw = -conj(0.5*tstep*fft(St)); % -Conj is used because MATLAB's fft has the exponent convention +i\omega\tau, but we want -i\omega\tau for the forward transform. Also, since the time-domain is doubled (i.e. -\beta to \beta), the - sign is required as the matsubara frequencies are offset by a constant \pi T.
Gw(2:2*points) = 1./(-(exp(-1i*w1)-1)/tstep+mu-Sw(2:2*points)).*mod(1:2*points-1,2); % The i\omega term has to be "lattice-regulated" in order to get the correct values of G close to \tau = 0. The mod restricts to fermion Matsubara frequencies.
Gt = (1-x)*Gt-x*T*2*points*ifft(conj(Gw)); % Another -conj for the same reason, and linear weighted update is required for convergence.
% Itinerant Part. Same strategy.
Sct = -g^2*Gct.*Gt.*fliplr(Gt);
Scw = -conj(0.5*tstep*fft(Sct));
Gcw(2:2*points) = (nu0/(2.0*pi))*(log(Lambda+2*muc-2*(exp(-1i*w1)-1)/tstep-2*Scw(2:2*points))-log(2*muc-Lambda-2*(exp(-1i*w1)-1)/tstep-2*Scw(2:2*points))).*mod(1:2*points-1,2); % This is the finite-bandwidth equation from the paper.
Gct = (1-x)*Gct-x*T*2*points*ifft(conj(Gcw));
plot(real(Gct)); % Plot what you want here
b = waitforbuttonpress; % You'll need to hold down a key to make the iteration go. Release when it converges.
end