Holstein impurity model coupled to a bath
Synopsis
In this example and the next example, we consider electron-phonon (el-ph) coupled systems. We start with a simple model describing electrons on two sites \(i=0,1\) and a phonon mode coupled to site 0. The Hamiltonian reads
Here, \(\hat{c}^\dagger_{i\sigma}\) is the creation operator of an electron with spin \(\sigma\) at lattice site \(i\), \(\hat{n}_i = \sum_{\sigma}\hat{c}_{i\sigma}^{\dagger}\hat{c}_{i\sigma}\), and \(\hat{a}^\dagger\) creates an Einstein phonon coupled to the 0th site. \(J(t)\) is the hopping parameter of the electrons, \(\epsilon_i\) is the energy level of site \(i\), \(\omega_0\) is the phonon frequency and \(g(t)\) is the el-ph coupling. If we regard the 0th site as an impurity, this model corresponds to a Holstein impurity model with a single bath site. This type of impurity model is self-consistently solved in a dynamical mean-field theory calculation, which we explain in the next example.
In the present example, we consider excitations via a modulation of the hopping parameter or the el-ph coupling. We treat the dynamics within the self-consistent Migdal approximation or the unrenormalized Migdal approximation.
Details and implementation
The source code is organized as follows:
|
main program for the unrenormalized Migdal approximation |
|
main program for the self-consistent Migdal approximation |
|
subroutines for the evaluation of the phonon energy and the phonon displacement |
|
implementation of self-energy approximations |
The corresponding declarations of the routines are in the .hpp files with the same name as the .cpp files. We will first discuss how to solve the problem, and then explain the Migdal approximations.
We first introduce the Green’s functions for the electrons and phonons as
where \(\hat{X} = \hat{a}^\dagger + \hat{a}\) and \(\Delta\hat{X}(t) = \hat{X}(t) - \langle \hat{X}(t) \rangle\). Here, we consider the spin symmetric case.
Since the system is coupled to the phonon only at site 0, one can trace out the contribution from the bath site (site 1) and focus on the Green’s function for site 0, \(G_{00}\). Then the electron and phonon Green’s functions are determined by the Dyson equations
and the phonon displacement, \(X(t)=\langle \hat{X}(t)\rangle\), which can be calculated as
Here the mean-field contribution (\(\Sigma^{\rm MF}(t)\)) corresponds to
and \(\Delta(t,t')\) is the hybridization function, which can be expressed as
with \([i\partial_t -\epsilon_j]G_{0,jj}(t,t')=\delta_\mathcal{C}(t,t')\). \(\Sigma^{\rm corr}(t,\bar{t})\) is the beyond-mean-field contribution to the self-energy, \(D_{0}(t,t')\equiv-i\langle \Delta\hat{X}(t) \Delta\hat{X}(t') \rangle_0\) is the Green’s function for the free phonon system, \(\Pi\) is the phonon self-energy and \(n_0(t)= \langle \hat{n}_0(t) \rangle\).
For a general impurity model, the hybridization function \(\Delta(t,t')\) takes an arbitrary form; however, the structure of equations (Eqs. (49), (50), (51), (52)) stays the same. This type of impurity problem is solved within nonequilibrium dmft, where the hybridization function is self-consistently determined (see Idea of DMFT).
Once \(G_{00}\) and \(\Sigma\) are obtained, the remaining electron Green’s functions can be easily calculated. For example, \(G_{10} = -G_{0,11} * J * G_{00}\). The energies can be also evaluated in a post processing operation. The kinetic energy (the expectation value of the first two terms of the Hamiltonian (48)) is
The interaction energy (the expectation of the third term of Eq. (48)) can be expressed as
where the first term is the mean-field contribution and the second term is the correlation energy. The phonon energy (the expectation value of the last term of Eq. (48)) is
Here \(D_{\rm PP}(t,t')=-i\langle T_\mathcal{C} \Delta\hat{P}(t) \Delta\hat{P}(t')\rangle\) with \(\hat{P}=\frac{1}{i}(\hat{a}-\hat{a}^\dagger)\), \(P(t) \equiv \langle\hat{P}(t)\rangle\) and \(\Delta\hat{P}(t) \equiv \hat{P}(t) - \langle \hat{P}(t)\rangle\).
Self-consistent Migdal approximation as an impurity solver: sMig
The self-consistent Migdal approximation is the lowest order (\(\mathcal{O}(g^2)\)) self-consistent approximation for the self-energies, based on renormalized Green’s functions. Within this approximation, one can treat the dynamics and the renormalization of the phonons induced by the electron-phonon coupling. Even though the phonons can absorb energy from the electrons, the total energy of the system is conserved in this approximation. The impurity self-energies for the electrons and phonons are approximated as
In the sample program, the sMig self-energy is computed by the routine Sigma_Mig (found in programs/Holstein_impurity_impl.cpp). We provide two interfaces for \(0\leq\) tstp \(\leq\) SolverOrder (bootstrapping part) and tstp \(=-1\), tstp \(>\) SolverOrder (Matsubara part and the time-stepping part), respectively. Here, we show the latter as an example:
void Sigma_Mig(int tstp, GREEN &G, GREEN &Sigma, GREEN &D0, GREEN &D, GREEN &Pi, GREEN &D0_Pi, GREEN &Pi_D0, CFUNC &g_el_ph, double beta, double h, int SolverOrder, int MAT_METHOD){
int Norb=G.size1();
int Ntau=G.ntau();
//step1: phonon self-energy gGgG
GREEN_TSTP gGg(tstp,Ntau,Norb,FERMION);
G.get_timestep(tstp,gGg);
gGg.right_multiply(g_el_ph);
gGg.left_multiply(g_el_ph);
GREEN_TSTP Pol(tstp,Ntau,1,BOSON);
Pi.set_timestep_zero(tstp);
Bubble1(tstp,Pol,0,0,gGg,0,0,G,0,0);
Pi.incr_timestep(tstp,Pol,-2.0);
//step2: solve phonon Dyson to evaluate D^corr
step_D(tstp, D0, D, Pi, D0_Pi, Pi_D0, beta, h, SolverOrder, MAT_METHOD);
//step3: evaluate electron self-energy
Bubble2(tstp,Sigma,0,0,D,0,0,gGg,0,0);
}
This routine contains three steps.
The phonon self-energy (Eq. (56)) is evaluated using
cntr::Bubble1.The Dyson equation of the phonon Green’s function (Eq. (50)) is solved using the routine
step_D, see below. HereMAT_METHODspecifies which method is used to solve the Dyson equation of the Matsubara part.The electron self-energy (Eq. (55)) is evaluated using
cntr::Bubble2.
For the bootstrapping part, we carry out these steps for tstp \(\leq\) SolverOrder at once, using start_D instead of step_D.
The routine for solving the Dyson equation for the phonons is implemented using the cntr routines vie2_**. The interface for solving the Matsubara Dyson equation and time stepping (step_D) is given by
void step_D(int tstp,GREEN &D0, GREEN &D, GREEN &Pi, GREEN &D0_Pi, GREEN &Pi_D0,
double beta, double h, int SolverOrder, int MAT_METHOD){
//step1: set D0_Pi=-D0*Pi
cntr::convolution_timestep(tstp,D0_Pi,D0,Pi,beta,h,SolverOrder);
D0_Pi.smul(tstp,-1.0);
cntr::convolution_timestep(tstp,Pi_D0,Pi,D0,beta,h,SolverOrder);
Pi_D0.smul(tstp,-1.0);
//step2: solve [1-D0*Pi]*D=[1+D0_Pi]*D=D0
if(tstp==-1){
cntr::vie2_mat(D,D0_Pi,Pi_D0,D0,beta,SolverOrder,MAT_METHOD);
}else {
cntr::vie2_timestep(tstp,D,D0_Pi,Pi_D0,D0,beta,h,SolverOrder);
}
}
This routine is separated into two steps.
\(-D_0 * \Pi\), which is the kernel of Eq. (50), and its conjugate for \(-\Pi * D_0\) are prepared by
cntr::convolution_timestep.Eq. (50) is solved as a
viewithvie2_**for the corresponding time step.
In start_D, we perform these steps for \(0 \leq\) tstp \(\leq\) SolverOrder at the same time, and use vie2_start.
The free phonon Green’s function, \(D_{0}(t,t')/2\) is computed using the cntr routine
cntr::green_single_pole_XX(D0,Phfreq_w0,beta,h);
Unrenormalized Migdal approximation as an impurity solver: uMig
In the unrenormalized Migdal approximation, the phonons are kept in equilibrium and play the role of a heat bath. Therefore, the total energy is not conserved. The impurity self-energy for the electrons is approximated as
We call this approximation uMig. The implementation is simpler than in the case of sMig and can be found in the routine Sigma_uMig.
Structure of the example program
The example programs for the sMig and uMig approximations are provided in programs/Holstein_impurity_singlebath_Migdal.cpp and programs/Holstein_impurity_single_bath_uMig.cpp. The main body of the example programs consists of two parts. In the first part, we focus on \(G_{00}\) and evaluate it by solving the equations with the corresponding approximations for the self-energies. Once \(G_{00}\) and the self-energies are determined, all other Green’s functions and the energies are evaluated in the second part.
As in the case of the Hubbard chain, the first part consists of three main steps: (i) solving the Matsubara Dyson equation, (ii) bootstrapping (tstp \(\leq\) SolverOrder) and (iii) time propagation for tstp \(>\) SolverOrder. Since the generic structure of each part is similar to that of the Hubbard chain, we only show here the part corresponding to the time propagation.
for(tstp = SolverOrder+1; tstp <= Nt; tstp++){
// Predictor: extrapolation
cntr::extrapolate_timestep(tstp-1,G,integration::I<double>(SolverOrder));
// Corrector
for (int iter=0; iter < CorrectorSteps; iter++){
cdmatrix rho_M(1,1), Xph_tmp(1,1);
cdmatrix g_elph_tmp(1,1), h00_MF_tmp(1,1);
//update correlated self-energies
Hols::Sigma_Mig(tstp, G00, Sigma, D0, D, Pi, D0_Pi, Pi_D0, g_elph_t, beta, h, SolverOrder);
//update phonon displacement
G00.density_matrix(tstp,rho_M);
rho_M *= 2.0;//spin number=2
n0_tot_t.set_value(tstp,rho_M);
Hols::get_phonon_displace(tstp, Xph_t, n0_tot_t, g_elph_t, D0, Phfreq_w0, SolverOrder, h);
//update mean-field
Xph_t.get_value(tstp,Xph_tmp);
g_elph_t.get_value(tstp,g_elph_tmp);
h00_MF_tmp = h00 + Xph_tmp*g_elph_tmp;
h00_MF_t.set_value(tstp,h00_MF_tmp);
//solve Dyson for G00
Hyb_Sig.set_timestep(tstp,Hyb);
Hyb_Sig.incr_timestep(tstp,Sigma,1.0);
cntr::dyson_timestep(tstp, G00, 0.0, h00_MF_t, Hyb_Sig, beta, h, SolverOrder);
}
}
Here the hybridization \(\Delta(t,t')\) is initially fixed via Eq. (53). At the beginning of each time step, we extrapolate \(G_{00}\) to the current time step, which works as a predictor. Next, we iterate Eqs. (49), (50), (51), (52), (55) and (56) (or (57)) to self-consistently determine the self-energies, which works as a corrector. During the iteration, we also update the phonon displacement, \(X(t)\), using Eq. (51). The update of \(X(t)\) is implemented in the routine get_phonon_displace using a cntr routine for the (Gregory) integration, integration::Integrator.
In the second part, we evaluate the rest of the Green’s functions. For the energies, the cntr routines correlation_energy_convolution and convolution_density_matrix are used. In particular, the calculation of the phonon energy (54) is implemented by the routine evaluate_phonon_energy.
Storing Green’s functions
In the example programs, we store Green’s functions with the hdf5 format using routines provided in cntr::herm_matrix, cntr::herm_matrix::write_to_hdf5_slices and cntr::herm_matrix::write_to_hdf5_tavtrel. The corresponding part looks
char fnametmp[1000];
std::sprintf(fnametmp,"%s_green.h5",flout);
hid_t file_id = open_hdf5_file(std::string(fnametmp));
hid_t group_id = create_group(file_id,"parm");
//parameters
store_int_attribute_to_hid(group_id,"Nt",Nt);
store_int_attribute_to_hid(group_id,"Ntau",Ntau);
store_double_attribute_to_hid(group_id,"dt",dt);
store_double_attribute_to_hid(group_id,"beta",beta);
close_group(group_id);
//green's function
//G at each t_n: G^R(t_n,t), G^<(t,t_n), G^tv(t_n,\tau) for t<=t_n
group_id = create_group(file_id, "G00_slice");
G.write_to_hdf5_slices(group_id,OutEvery);
close_group(group_id);
// G(t_rel,t_av) for t_av = [0,OutEvery*dt,2*OutEvery*dt...] for t_rel >= 0
group_id = create_group(file_id, "G00_tavrel");
G.write_to_hdf5_tavtrel(group_id,OutEvery);
close_group(group_id);
It is useful to store the number of time steps (Nt,Ntau), the size of the time step (dt) and the inverse temperature as parm in the same hdf5 file. The routine write_to_hdf5_slices stores the contour function \(F(t,t')\) in the form of \(F^R({\rm tstp},n), F^<(n,{\rm tstp}), F^{\rceil}(n,\tau)\) for \(n\leq {\rm tstp}\) at \({\rm tstp} = m\cdot {\rm OutEvery}\). The routine write_to_hdf5_tavtrel stores the retarded and lesser parts in the form of \(F(t_{\rm rel};t_{\rm av})\), where \(t_{\rm rel}=t-t'\) and \(t_{\rm av}=\frac{t+t'}{2}\), for \(\frac{t_{\rm av}}{dt}\) = \([0,{\rm OutEvery} ,2{\rm OutEvery},...]\) for \(t_{\rm rel}\geq 0\).
Running the example programs
There are two programs for sMig and uMig, respectively: Holstein_impurity_singlebath_Migdal.x and Holstein_impurity_singlebath_uMig.x. These programs allow to treat excitations via a modulation of the hopping parameter and the el-ph coupling. We use \(\epsilon_{0,\mathrm{MF}}\equiv \epsilon_{0} + gX(0)\) as an input parameter instead of \(\epsilon_{0}\). (\(\epsilon_{0}\) is determined in the post processing.) The driver script demo_Holstein_impurity.py located in the utils/ directory provides a simple interface to these programs. After defining the system parameters, calculation scheme (sMig or uMig, time step, convergence parameter, etc.) and excitation parameters, the script creates the corresponding input file and starts the simulation. So simply run
python utils/demo_Holstein_impurity.py
In demo_Holstein_impurity.py, the shape of the excitation is
and one needs to specify \(T_{\rm end}\), \(\Omega\), \(\Delta g\) and \(\Delta J\). At the end of the simulation, the number of particles per spin for each site (\(n_{i,\sigma}(t)\)), the phonon displacement (\(X(t)\)), phonon momentum (\(P(t)\)) and the energies are plotted. In addition, the spectral functions of the electrons and phonons,
for the average time \(t_{\rm av}=\frac{Nt \cdot h}{2}\). Here \(t_{\rm rel}\) indicates the relative time, and \(F_{\rm gauss}(t_{\rm rel})\) is a Gaussian window function, which can also be specified in demo_Holstein_impurity.py. Specifically, the evaluation of the spectral function is done in two steps as
t_rel,G_ret,G_les = read_gf_ret_let_tavtrel(output_file + '_green.h5','Gloc_tavrel',dt,tstp_plot)
A_w = evaluate_spectrum_windowed(G_ret,t_rel,ome,method,Fwindow)
In the first line, with a python routine read_gf_ret_let_tavtrel in python3/SpectrumCNTRhdf5.py, the relative time and the corresponding retarded and lesser Green’s functions are obtained from the .h5 file for the Green’s functions stored in the form of \(G (t_{\rm rel};t_{\rm av})\). Here \(t_{\rm av}\) specified with tstp_plot. In the second line, we operate Eq. (59) with a subroutine evaluate_spectrum_windowed in python3/SpectrumCNTRhdf5.py. Here one can specify the order of the interpolation of Green’s functions for the Fourier transformation, and Fwindow is the window function, which we take a Gauss function in this example.
Discussion
In Fig. 12 (a)–(b), we show the electron spectral function measured at site 0 and the phonon spectrum evaluated within the sMig approximation. For \(J=0\), the site 0 is decoupled from the site 1. The electron spectrum shows a main peak around \(\omega=0.5=\epsilon_{0,\mathrm{MF}}\) and side peaks around \(\epsilon_{0,\mathrm{MF}}\pm\omega_0\). As the hopping parameter is increased, the electron levels at the site 0 and the site 1 hybridize and an additional sharp peak at \(\omega<0\) appears in the electron spectrum, see panel (a). The phonon spectrum shown in panel (b) reveals a shift of the peak position from the bare value \(\omega_0=1.0\) and a finite width of the peak (phonon renormalization).
In Fig. 12 (c)–(e), we show the time evolution of the phonon displacement, the density and the total energy after an excitation via hopping modulation. These results were obtained using the sMig approximation. The phonon displacement shows damped oscillations, see panel (c). Panel (d) demonstrates that after the excitation, the number of electrons at site 0 and site 1 changes with time, while the total number is conserved. Within sMig the total energy is also conserved after the excitation, as is shown in panel (e). On the other hand, when uMig is used, the total energy is not conserved after the excitation, see panel (f). In uMig, the phonons can act as a heat bath and dissipate energy, which is clearly evident when the method is used to treat extended systems, see next section.
Fig. 12 Simulation results for the Holstein impurity model. (a)(b) Spectral functions of the electrons at site 0 (a) and phonon spectrum (b) for various \(J\) obtained using sMig. (c)(d)(e) Time evolution of the phonon displacement, the electron density and the total energy after excitation via hopping modulations with different strengths obtained within sMig. (f) Time evolution of the total energy after excitation via hopping modulations with different strength obtained within uMig. For (a-f), we use \(g=0.5\), \(\omega_0=1.0\), \(\beta=2.0\), \(\epsilon_{0,mf}=0.5\), \(\epsilon_1=-0.5\). For (a-b), a window width \(\sigma_{\rm window}=20.0\) is used for \(F_{\rm gauss}(t_{\rm rel})\). For (c-f), we use \(J_0=0.5\) and a \(\sin^2\) envelope for the hopping modulation with excitation frequency \(\Omega=2.0\) and pulse duration \(T_{\rm end}=12.56\).