% (c) 2018 Aavishkar A. Patel
% This code solves the imaginary-time Dyson equations iteratively for the RM+U(1) gauge-field model
% on a discrete-time grid, by Fourier transforming back and forth between
% time and frequency representations of the two point functions.
% Gt is the time-domain (-\beta to \beta) fermion
% Green's function and Dt is the same for the boson. S is a shorthand for \Sigma,
% the self-energy, and w is a shorthand for \omega.
% Values of couplings, number of grid points,
% and other parameters in the model. Also sets the linear update weights xf and xb.
t = 1.0; T = 1.0e-5; mu = 0.0; itermax = 600; points = 2^22;
xf = 0.05; xb = 0.02; g = 1.0; MbyN = 0.5; mB = 1.0e-6;
% 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 (bosonic) Matsubara
% frequencies.
skipinit = false; % Set true if you want to reuse a previous solution as initial condition.
if(~skipinit)
Gt = 0.5*sign(points+0.5-(1:2*points)); % Free fermion Green's function.
Gw = -conj(0.5*tstep*fft(Gt));
St = zeros(1,2*points);
Sw = zeros(1,2*points);
Pit = 2.0*t^2*MbyN*Gt.*fliplr(Gt); % Flipping the time domain array left to right implements \tau --> -\tau.
Piw = conj(0.5*tstep*fft(Pit));
Dw = zeros(1,2*points);
Dw(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2+Piw(2:2*points)-Piw(1)).*(1-mod(1:2*points-1,2));
Dt = T*2*points*ifft(conj(Dw));
Dthold = Dt;
end
% Begin iterations
for n=1:itermax
% Fermion Part
St = t^2*Gt.*(1+Dt-Dt(1));
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 for fermions as the matsubara frequencies are offset by a constant \pi T.
Gw(2:2*points) = (1./(1i*w1/tstep+mu-Sw(2:2*points))-tstep./(1i*w1)).*mod(1:2*points-1,2); % The mod restricts to fermion Matsubara frequencies.
% Boson Part. Same strategy.
Pit = 2.0*t^2*MbyN*Gt.*fliplr(Gt);
Piw = conj(0.5*tstep*fft(Pit));
Dw(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2+Piw(2:2*points)-Piw(1)).*(1-mod(1:2*points-1,2));
% Update Green's functions
Gt = (1-xf)*Gt+xf*(0.5*sign(points+0.5-(1:2*points))-T*2*points*ifft(conj(Gw))); % Another -conj for the same reason, and linear weighted update is required for convergence.
Dthold = (1-xb)*Dt+xb*T*2*points*ifft(conj(Dw));
% Overshoot protection mechanism to avoid unstable regime and ensure convergence
if(1-Dthold(1)<0)
xb = 0.9*xb; % If boson overshoot detected, block the boson update and slow down its future updates, allowing the fermion to catch up.
else
Dt = Dthold; % If no overshoot detected, proceed as usual.
end
plot(-imag(Sw(2:2:100))); % Plot negative imaginary part of Matsubara fermion self-energy.
pause(1/60); % Stop by using ctrl+c when converged.
disp([100*(1-Dt(1)) n]);
end
Gw(2:2*points) = (1./(1i*w1/tstep+mu-Sw(2:2*points))).*mod(1:2*points-1,2);
Gw0 = Gw;
Gw0(2:2*points) = (1./(1i*w1/tstep+mu)).*mod(1:2*points-1,2);
Q = Gt(points); % Filling
FF = T*sum(log(Gw(2:2:2*points)./Gw0(2:2:2*points))) - T*log(1+exp(mu/T)) + mu*Q - T*sum(Sw.*Gw) + 0.5*t^2*T*sum(Gw.^2); % Fermion free energy.
Dw0 = Dw;
Dw0(2:2*points) = 1./((w1/tstep).^2/g^2+mB^2).*(1-mod(1:2*points-1,2));
FB = 0.5*T*sum(log(Dw0(3:2:2*points)./Dw(3:2:2*points)))/(2.0*MbyN) + T*log(1-exp(-mB*g/T))/(2.0*MbyN); % Boson free energy.
F = FF + FB; % Total free energy.