Hubbard chain
Synopsis
In this example we consider the one-dimension (1D) Hubbard model with the Hamiltonian
where \(\langle i,j\rangle\) constrains the lattice sites \(i,j\) to the nearest neighbors, while \(\bar\sigma=\uparrow,\downarrow\) for \(\sigma=\downarrow,\uparrow\). We consider \(M\) lattice sites with open boundary conditions. Furthermore, we restrict ourselves to the paramagnetic case with an equal number of spin-up (\(N_\uparrow\)) and spin-down (\(N_\downarrow\)) particles. The number of particles determines the filling factor \(\bar{n}=N_\uparrow/M\).
In this example, the system is excited with an instantaneous quench of the on-site potential of the first lattice site to \(w_0\):
Here, we will treat the dynamics with respect to the Hamiltonian (44) within the second-Born (2B) and the \(GW\) approximation.
Details and implementation
The source code is organized as follows:
|
main program for the 2B approximation |
|
main program for the \(GW\) approximation |
|
implementation of self-energy approximations |
Corresponding declarations of the routines are in the .hpp files with the same name as the .cpp files. We will discuss the implementation of the 2B approximation in detail and then remark on how to adopt it for the \(GW\) approximation.
Second-Born approximation
The 2B approximation corresponds to the second-order expansion in terms of the Hubbard repulsion \(U(t)\), which we treat as time dependent here for generality. Defining the Green’s function with respect to the lattice basis \(i,j\), the 2B is defined by
Due to the paramagnetic constraint, \(G_{ij,\sigma}(t,t^\prime) = G_{ij,\bar\sigma}(t,t^\prime)\); the spin index can be dropped throughout. The 2B self-energy (45) is implemented in two steps.
The (per-spin) polarization \(P_{ij}(t,t^\prime)= -i G_{ij}(t,t^\prime)G_{ji}(t^\prime,t)\) is computed using the routine
cntr::Bubble1and subsequent multiplication by (-1).The self-energy is then given by \(\Sigma_{ij}(t,t^\prime) = i U(t) U(t^\prime) G_{ij}(t,t^\prime) P_{ij}(t,t^\prime)\), which corresponds to a bubble diagram computed by the routine
cntr::Bubble2.
The 2B self-energy is computed by the routine Sigma_2B (found in programs/hubbard_chain_selfen_impl.cpp) as follows:
void Sigma_2B(int tstp, GREEN &G, CFUNC &U, GREEN &Sigma){
int nsites=G.size1();
int ntau=G.ntau();
GREEN_TSTP Pol(tstp,ntau,nsites,BOSON);
Polarization(tstp, G, Pol);
Pol.right_multiply(U);
Pol.left_multiply(U);
for(int i=0; i<nsites; i++){
for(int j=0; j<nsites; j++){
cntr::Bubble2(tstp,Sigma,i,j,G,i,j,Pol,i,j);
}
}
}
First, the polarization Pol, which represents \(P_{ij}(t,t^\prime)\), is defined for the given time step. After computing \(P_{ij}(t,t^\prime)\) by the function
void Polarization(int tstp, GREEN &G, GREEN_TSTP &Pol){
int nsites=G.size1();
for(int i=0; i<nsites; i++){
for(int j=0; j<nsites; j++){
cntr::Bubble1(tstp,Pol,i,j,G,i,j,G,i,j);
}
}
Pol.smul(-1.0);
}
the lines
Pol.right_multiply(U);
Pol.left_multiply(U);
perform the operation \(P_{ij}(t,t^\prime) \rightarrow P_{ij}(t,t^\prime) U(t^\prime)\) and \(P_{ij}(t,t^\prime) \rightarrow U(t)P_{ij}(t,t^\prime)\), respectively. Finally, cntr::Bubble2 computes \(\Sigma_{ij}(t,t^\prime)\).
Mean-field Hamiltonian and onsite quench
The total self-energy also includes the Hartree-Fock (HF) contribution, which we incorporate into the mean-field Hamiltonian \(\epsilon^{\mathrm{MF}}_{ij}(t) = \epsilon^{0}_{ij}(t) + U (n_i-\bar{n})\) with the occupation (per spin) \(n_i(t)= \mathrm{Im}[G^<_{ii}(t,t)]\). The shift of chemical potential \(-U \bar{n}\) is a convention to fix the chemical potential at half filling at \(\mu=0\). In the example program, the mean-field Hamiltonian is represented by the cntr::function eps_mf. Updating eps_mf is accomplished by computing the density matrix using the class routine cntr::herm_matrix::density_matrix.
The general procedure to implement a quench of some parameter \(\lambda\) at \(t=0\) is to represent \(\lambda\) by a contour function \(\lambda_n\): \(\lambda_{-1}\) corresponds to the pre-quench value which determines the thermal equilibrium, while \(\lambda_n\) with \(n\ge 0\) governs the time evolution. In the example program, we simplify this procedure by redefining \(\epsilon^{0}_{ij} \rightarrow \epsilon^{0}_{ij} + w_0 \delta_{i,1}\delta_{j,1}\) after the Matsubara Dyson equation has been solved.
Structure of the example program
The program is structured similar to the previous examples. After reading variables from file and initializing the variables and classes, the Matsubara Dyson equation is solved in a self-consistent fashion. The example below illustrates this procedure for the 2B approximation.
tstp=-1;
gtemp = GREEN(SolverOrder,Ntau,Nsites,FERMION);
gtemp.set_timestep(tstp,G);
for(int iter=0;iter<=MatsMaxIter;iter++){
// update mean field
hubb::Ham_MF(tstp, G, Ut, eps0, eps_mf);
// update self-energy
hubb::Sigma_2B(tstp, G, Ut, Sigma);
// solve Dyson equation
cntr::dyson_mat(G, MuChem, eps_mf, Sigma, beta, SolveOrder);
// self-consistency check
err = cntr::distance_norm2(tstp,G,gtemp);
if(err<MatsMaxErr){
break;
}
gtemp.set_timestep(tstp,G);
}
Updating the mean-field Hamiltonian (hubb::Ham_MF), the self-energy (hubb::Sigma_2B) and solving the corresponding Dyson equation (cntr::dyson_mat) is repeated until self-consistency has been reached, which in practice means that the deviation between the previous and updated Green’s function is smaller than MatsMaxErr. For other self-energy approximations, the steps described above (updating auxiliary quantities) have to be performed before the self-energy can be updated.
Once the Matsubara Dyson equation has been solved up to the required convergence threshold, the starting algorithm for time steps \(n=0,\dots,k\) can be applied. To reach self-consistency for the first few time steps, we employ the bootstrapping loop:
for (int iter = 0; iter <= BootstrapMaxIter; iter++) {
// update mean field
for(tstp=0; tstp<=SolveOrder; tstp++){
hubb::Ham_MF(tstp, G, Ut, eps0, eps_mf);
}
// update self-energy
for(tstp=0; tstp<=SolveOrder; tstp++){
hubb::Sigma_2B(tstp, G, Ut, Sigma);
}
// solve Dyson equation
cntr::dyson_start(G, MuChem, hmf, Sigma, beta, h, SolveOrder);
// self-consistency check
err=0.0;
for(tstp=0; tstp<=SolverOrder; tstp++) {
err += cntr::distance_norm2(tstp,G,gtemp);
}
if(err<BootstrapMaxErr && iter>2){
break;
}
for(tstp=0; tstp<=SolverOrder; tstp++) {
gtemp.set_timestep(tstp,G);
}
}
Finally, after the bootstrapping iteration has converged, the time propagation for time steps n>SolveOrder is launched. The self-consistency at each time step is accomplished by iterating the update of the mean-field Hamiltonian, Green’s function and self-energy over a fixed number of CorrectorSteps. As an initial guess, we employ a polynomial extrapolation of the Green’s function from time step \(n-1\) to \(n\), as implemented in the routine cntr::extrapolate_timestep. Thus, the time propagation loop takes the form
for(tstp = SolverOrder+1; tstp <= Nt; tstp++){
// Predictor: extrapolation
cntr::extrapolate_timestep(tstp-1,G,SolveOrder);
// Corrector
for (int iter=0; iter < CorrectorSteps; iter++){
// update mean field
hubb::Ham_MF(tstp, G, Ut, eps0, hmf);
// update self-energy
hubb::Sigma_2B(tstp, G, Ut, Sigma);
// solve Dyson equation
cntr::dyson_timestep(tstp,G,MuChem,eps_mf,Sigma,beta,h,SolveOrder);
}
}
After the Green’s function has been computed for all required time steps, we compute the observables. In particular, the conservation of the total energy provides a good criterion to assess the accuracy of the calculation. The total energy for the Hubbard model (43) is given in terms of the Galitskii-Migdal formula
The last term, known as the correlation energy, is most conveniently computed by the routine:
Ecorr = cntr::correlation_energy(tstp, G, Sigma, beta, h);
GW approximation
Next we consider the \(GW\) approximation. We remark that we formally treat the Hubbard interaction as spin-independent, while the spin-summation in the polarization \(P\) (which is forbidden by the Pauli principle) is excluded by the corresponding prefactor.
Within the same setup as above, the \(GW\) approximation is defined by
where \(\delta W_{ij}(t,t^\prime)\) denotes the dynamical part of the screened interaction \(W_{ij}(t,t^\prime) = U \delta_{ij}\delta_\mathcal{C}(t,t^\prime) + \delta W_{ij}(t,t^\prime)\). We compute \(\delta W_{ij}(t,t^\prime)\) from the charge susceptibility \(\chi_{ij}(t,t^\prime)\) by \(\delta W_{ij}(t,t^\prime) = U(t) \chi_{ij}(t,t^\prime) U(t^\prime)\). In turn, the susceptibility obeys the Dyson equation
The strategy to compute the \(GW\) self-energy consists of three steps:
Computing the polarization \(P_{ij}(t,t^\prime)\) by
cntr::Bubble1(as in Second-Born approximation).Solving the Dyson equation (47) as VIE. By defining the kernel \(K_{ij}=-P_{ij}(t,t^\prime)U(t^\prime)\) and its hermitian conjugate, Eq. (47) amounts to \([1+K]\ast \chi=P\), which is solved using
cntr::vie2.Computing the self-energy by
cntr::Bubble2.
The implementation of step 1 was already discussed above. For step 2, we distinguish between the equilibrium and time stepping on the one hand, and the starting phase on the other hand. For the former, we have defined the routine
void GenChi(int tstp, double h, double beta, GREEN &Pol,
CFUNC &U, GREEN &PxU, GREEN &UxP, GREEN &Chi, int order){
PxU.set_timestep(tstp, Pol);
UxP.set_timestep(tstp, Pol);
PxU.right_multiply(tstp, U);
UxP.left_multiply(tstp, U);
PxU.smul(tstp,-1.0);
UxP.smul(tstp,-1.0);
if(tstp==-1){
cntr::vie2_mat(Chi,PxU,UxP,Pol,beta,CINTEG(order));
} else{
cntr::vie2_timestep(tstp,Chi,PxU,UxP,Pol,CINTEG(order),beta,h);
}
}
Here, PxU and UxP correspond to the kernel \(K_{ij}\) and its hermitian conjugate, respectively. Analogously, the starting routine is implemented as
void GenChi(double h, double beta, GREEN &Pol, CFUNC &U,
GREEN &PxU, GREEN &UxP, GREEN &Chi, int order){
for(int n = 0; n <= order; n++){
PxU.set_timestep(n, Pol);
UxP.set_timestep(n, Pol);
PxU.right_multiply(n, U);
UxP.left_multiply(n, U);
PxU.smul(n,-1.0);
UxP.smul(n,-1.0);
}
cntr::vie2_start(Chi,PxU,UxP,Pol,CINTEG(order),beta,h);
}
Finally, the self-energy is computed by
void Sigma_GW(int tstp, GREEN &G, CFUNC &U, GREEN &Chi, GREEN &Sigma){
int nsites=G.size1();
int ntau=G.ntau();
GREEN_TSTP deltaW(tstp,ntau,nsites,CNTR_BOSON);
Chi.get_timestep(tstp,deltaW);
deltaW.left_multiply(U);
deltaW.right_multiply(U);
for(int i=0; i<nsites; i++){
for(int j=0; j<nsites; j++){
cntr::Bubble2(tstp,Sigma,i,j,G,i,j,deltaW,i,j);
}
}
}
The above implementation of the \(GW\) self-energy can be found in programs/hubbard_chain_selfen_impl.cpp. The structure of the program (programs/hubbard_chain_gw.cpp) is analogous to Structure of the example program.
Running the example programs
There are two programs for the 2B and the \(GW\) approximation, respectively: hubbard_chain_2b.x, hubbard_chain_gw.x. The driver script demo_hubbard_chain.py located in the utils/ directory provides a simple interface to these programs. Simply run
python utils/demo_hubbard_chain.py
After defining the parameters and convergence parameters, the script creates the corresponding input file and launches all three programs in a loop. The occupation of the first lattice site \(n_1(t)\) and the kinetic and total energy are then read from the output files and plotted. The script demo_hubbard_chain.py also allows to pass reference data as an optional argument, which can be used to compare to exact results. For some parameters, results of the exact treatment are available in the data/ directory.
Discussion
As an example, we consider the Hubbard dimer (\(M=2\)) for \(U=1\) at half filling (\(\mu=0\), \(\bar{n}=1/2\)). The system is excited by a strong on-site quench \(w_0=5\). The results are obtained by running
python utils/demo_hubbard_chain.py data/hubbardchain_exact_M2_U1_n1_w5.dat
This strong excitation leads to the phenomenon of artificial damping: although the density \(n_1(t)\) exhibits an oscillatory behavior for all times in an exact treatment, the approximate NEGF treatment — with either self-energy approximation considered here — shows damping until an unphysical steady state is reached.
Fig. 11 Dynamics of the Hubbard chain: Here, we consider the Hubbard dimer with M=2 and U=1 at half filling. We used Nt=400 and Ntau=400 at an inverse temperature of beta=20.0 and a time step of h=0.025.