Schwinger-Dyson Equation Method to Calculate Total Energy in Periodic Anderson Model with Electron-Phonon Interactions

TL;DR

This paper introduces a novel method to calculate the total energy in the periodic Anderson model with electron-phonon interactions using a path integral approach, enabling broader applicability beyond Hamiltonian-based systems.

Contribution

A new total energy calculation method derived from path integrals, applicable to systems lacking a Hamiltonian description, demonstrated on the periodic Anderson model.

Findings

Successfully derived total energy formula from path integral formalism.

Applied method to the periodic Anderson model with electron-phonon interactions.

Results match traditional Hamiltonian-based formulas.

Abstract

We have recently employed periodic Anderson model with electron-phonon interactions to describe Cerium volume collapse ($\gamma \rightarrow \alpha$ transition) under pressure. To describe the volume collapse transition in Cerium, we have tried to plot the pressure versus volume curves and see when a kink structure emerges. One way to obtain the pressure versus volume curve is to calculate the total energy of the system at different temperatures. In order to solve the periodic Anderson model with electron-phonon interactions, we have used the continuous time quantum Monte Carlo algorithm by integrating out the phonons to obtain a retarded electron-electron interaction. Monte Carlo simulation results give us the model's electronic Green function and self energy, from which we can calculate the system's total energy. A well known formula for calculating total energy through the knowledge of electron's self energy and Green function is derived from the model's Hamiltonian, which is not available for our case. Here, we have devised a new method to derive the total energy formula purely from a path integral description of the system, without resort to its Hamiltonian formulation. Since not all systems can be described by a Hamiltonian, thus our method has a wider applicability. It is noteworthy that the total energy formula that we derived here takes an identical form with that obtained from a model's Hamiltonian.

Equations

$$E$$

$$+\frac{1}{\beta }\sum _{i{\omega }_n}\sum _{\mathbf{k},\sigma }{\Sigma }_{\mathbf{k},\sigma }(i{\omega }_n)G_{\mathbf{k},\sigma }(i{\omega }_n),$$

$$\hat{H}=ϵ{\hat{c}}^{†}\hat{c}+\omega {\hat{a}}^{†}\hat{a}+g{\hat{c}}^{†}\hat{c}(\hat{a}+{\hat{a}}^{†})$$

$$Z$$

$$\displaystyle S=\int_{0}^{\beta}d\tau\Big{(}\bar{c}(\partial_{\tau}+\epsilon)c+\bar{a}(\partial_{\tau}+\omega)a+g\bar{c}c(a+\bar{a})\Big{)}$$

$$Z$$

$$=\frac{1}{1-e^{-\beta \omega }}\int \mathcal{D}[\bar{c},c]e^{-S_{eff}},$$

$$S_{eff}={\int }_0^{\beta }d\tau \bar{c}(\tau )({\partial }_{\tau }+ϵ)c(\tau )+\frac{g^2}{2}\iint d\tau d{\tau }^{\prime }\bar{c}(\tau )c(\tau )D_0(\tau -{\tau }^{\prime },\omega )\bar{c}({\tau }^{\prime })c({\tau }^{\prime })$$

$$D_0(\tau ,\omega )=-\frac{e^{\omega |\tau |}}{e^{\beta \omega }-1}-\frac{e^{-\omega |\tau |}}{1-e^{-\beta \omega }},-\beta <\tau <\beta$$

$$G({\tau }^{\prime }-\tau )=-⟨T\hat{c}({\tau }^{\prime }){\hat{c}}^{†}(\tau )⟩$$

$$G({\tau }^{\prime }-\tau )=-\frac{1}{Z}\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]c({\tau }^{\prime })\bar{c}(\tau )e^{-S}$$

$$\frac{\partial G({\tau }^{\prime },\tau )}{\partial {\tau }^{\prime }}=-\delta ({\tau }^{\prime }-\tau )+ϵ⟨T\hat{c}({\tau }^{\prime }){\hat{c}}^{†}(\tau )⟩+g⟨T\hat{c}({\tau }^{\prime })\hat{\varphi }({\tau }^{\prime }){\hat{c}}^{†}(\tau )⟩$$

$$({\partial }_{{\tau }^{\prime }}+ϵ)G({\tau }^{\prime },\tau )=-\delta ({\tau }^{\prime }-\tau )+g⟨T\hat{c}({\tau }^{\prime })\hat{\varphi }({\tau }^{\prime }){\hat{c}}^{†}(\tau )⟩$$

$$c\to c^{\prime }=c+\delta c,$$

$$\bar{c}\to {\bar{c}}^{\prime }=\bar{c}+\delta \bar{c},$$

$$\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]\bar{c}(\tau )e^{-S[\bar{c},c;\bar{a},a]}=\int \mathcal{D}[{\bar{c}}^{\prime },c^{\prime }]\mathcal{D}[\bar{a},a]{\bar{c}}^{\prime }(\tau )e^{-S[{\bar{c}}^{\prime },c^{\prime };\bar{a},a]}$$

$$S[{\bar{c}}^{\prime },c^{\prime };\bar{a},a]$$

$$\displaystyle+\int_{0}^{\beta}d\tau\Big{(}\delta\bar{c}(\partial_{\tau}+\epsilon)c+\bar{c}(\partial_{\tau}+\epsilon)\delta c+g\delta\bar{c}c(a+\bar{a})+g\bar{c}\delta c(a+\bar{a})\Big{)}$$

$$=:S[\bar{c},c;\bar{a},a]+\delta S$$

$$\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]\bar{c}(\tau )e^{-S[\bar{c},c;\bar{a},a]}$$

$$=\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a](\bar{c}(\tau )+\delta \bar{c}(\tau ))e^{-S}(1-\delta S)$$

$$=\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]\bar{c}(\tau )e^{-S[\bar{c},c;\bar{a},a]}+\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]e^{-S}(\delta \bar{c}(\tau )-\bar{c}(\tau )\delta S)$$

$$0=\int \mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]e^{-S}(\delta \bar{c}(\tau )-\bar{c}(\tau )\delta S)$$

$$\delta \bar{c}(\tau )-\bar{c}(\tau )\delta S$$

$$\displaystyle\int_{0}^{\beta}d\tau^{\prime}\Bigg{[}\delta(\tau-\tau^{\prime})-(\partial_{\tau^{\prime}}+\epsilon)\frac{1}{Z}\int\mathcal{D}[\bar{c},c]\mathcal{D}[\bar{a},a]e^{-S}c(\tau^{\prime})\bar{c}(\tau)$$

$$({\partial }_{{\tau }^{\prime }}+ϵ)G({\tau }^{\prime },\tau )-g⟨T\hat{c}({\tau }^{\prime })\hat{\varphi }({\tau }^{\prime }){\hat{c}}^{†}(\tau )⟩=-\delta (\tau -{\tau }^{\prime })$$

$$S={\int }_0^{\beta }d\tau \bar{c}(\tau )({\partial }_{\tau }+ϵ)c(\tau )+S_{\varphi },$$

$$\displaystyle S_{\phi}=\int_{0}^{\beta}d\tau\Big{(}\bar{a}(\partial_{\tau}+\omega)a+g\bar{c}c(a+\bar{a})\Big{)}$$

$$Z$$

$$=Z_{\varphi }\int \mathcal{D}[\bar{c},c]e^{-S_{eff}}$$

$$\int \mathcal{D}[\bar{c},c]e^{-S_{eff}[\bar{c},c]}\bar{c}(\tau )=\int \mathcal{D}[{\bar{c}}^{\prime },c^{\prime }]e^{-S_{eff}[{\bar{c}}^{\prime },c^{\prime }]}{\bar{c}}^{\prime }(\tau )$$

$$0=\int \mathcal{D}[\bar{c},c]e^{-S_{eff}[\bar{c},c]}(-\delta S_{eff}\bar{c}(\tau )+\delta \bar{c}(\tau ))$$

$$\delta S_{eff}={\int }_0^{\beta }d\tau \delta \bar{c}({\partial }_{\tau }+ϵ)c+g^2{\int }_0^{\beta }d{\tau }_1{\int }_0^{\beta }d{\tau }_2\delta \bar{c}({\tau }_1)c({\tau }_1)D^0({\tau }_1-{\tau }_2)\bar{c}({\tau }_2)c({\tau }_2)$$

$$\delta \bar{c}(\tau )-\delta S_{eff}\bar{c}(\tau )$$

Full text

††preprint: APS/123-QED

I Introduction

Recently, we employ periodic Anderson model with electron-phonon interactions to study the behavior of $4f$ electrons in Cerium to gain further understanding about the volume collapse that Cerium experiences at high pressurePeng_2013 . We have introduced electron-phonon interactions to the periodic Anderson model to capture the influence of lattice dynamics upon the Cerium volume collapse. It is argued that the Cerium volume collapse should be due to the emergence of a first order phase transition. A good method to characterize this first order phase transition is to obtain pressure versus volume curves at different temperatures and see when the pressure versus volume curve develops a kinkPhysRevB.99.155147 . We can calculate the pressure by taking partial derivative of free energy with respect to volume, whereas free energy can be calculated once we have obtained energy as a function of temperaturePhysRevB.80.140505 . Thus, a central issue we encountered is to calculate the total energy for the periodic Anderson model with electron-phonon interactions. It is well know that once we have a model’s self energy and Green function, we can calculate its total energy via this formula(see, for example, Ref. fetter2003quantum ):

[TABLE]

where $\epsilon_{\bm{k}}$ is electron’s dispersion relation in a lattice, $\Sigma_{\bm{k},\sigma}(i\omega_{n})$ is electron’s self energy, and $G_{\bm{k},\sigma}(i\omega_{n})$ is electron’s Green function, with $\omega_{n}=(2n+1)\pi/\beta$ being the Matsubara frequency. This formula is a natural corollary of equation of motion for electron’s Green function, which can be conveniently derived once we have the model’s Hamiltonian. However, not all physical systems are susceptible to a Hamiltonian formulation, and a notable family of such systems are the ones that are derived from the principle of least actionfeynman2005feynman . It is universally known that the quantum analogue of the principle of least action is Feynman’s path integral method. Therefore, for the systems that can only be properly described by a path integral, the validity of Eq. [1] which is derived with the implicit assumption for the existence of a Hamiltonian is at stake, and a justification or refutation of that formula is on its way.

For our case, in order to solve the periodic Anderson model with electron-phonon interactions, i.e., obtain the electron’s Green function and self energy, we have employed the diagrammatic Monte Carlo method as detailed in Ref. CTQMC_2007 , where the authors have integrated out the phonons to obtain a time-retarded electron-electron interaction. Once the phonons have been integrated out, we are left with an effective action, from which no explicit Hamiltonian can be obtained. Because a Hamiltonian is directly related with the total energy of a system, it becomes tricky to calculate a system’s total energy without an explicit Hamiltonian formulation. It it however noteworthy that not all systems with retarded interactions are bereft of a Hamiltonian formulation, see, for example, Ref. PhysRevB.96.075155 , where the authors can still obtain the total energy of a system with retarded interactions using operator method due to the availability of a Hamiltonian formulation for their system. In our case for which no explicit Hamiltonian formulation is available, we are going to devise a new method to calculate the total energy using effective action, rather than Hamiltonian, as our starting point. We will show that the formula (Eq. [1]) for calculating the total energy which was originally derived from the Hamiltonian is still applicable here when the system can only be described by an effective action, although the mathematical details for the formula derivation are completely different.

The problem of dealing quantum mechanically with a system that is derived from the principle of least action is encountered by Feynman, who designed the path integral approach to quantum mechanics just to address such an issueRevModPhys.20.367 ; feynman2005feynman ; feynman1965quantum . By integrating out the phonons in our model, we also encounter a similar situation, where we are left with an effective system with no Hamiltonian formulation. In order to derive the total energy formula for such a system, we need to employ the Schwinger-Dyson equation method, which furnishes us with an equation motion for Green function which is purely derived from path integrals. Once we have the equation of motion for electron’s Green function, we can easily calculate the total energy of the system. Therefore, the central topic of this paper is the equation of motion for Green function, which we shall derive using several different methods. The equivalence of these methods is also a major topic of this paper.

The organization of the paper is as follows. In section II, we will show the equivalence of Schwinger-Dyson equation method and the operator method in the derivation of equations of motion in a toy model that nevertheless captures the essence of electron-phonon interactions. In this section, we will derive the equation of motion using three methods: the operator method, the path integral method with phonons present and the path integral method with phonons integrated out. All of these methods yield identical results. Having established the validity of Schwinger-Dyson equation method in this toy model, we will continue to show in section III that we can derive the equation of motion, and thus the total energy, of periodic Anderson model with electron-phonon interactions (PAM + phonons). Since we have used the continuous time quantum Monte Carlo algorithm to simulate PAM + phonons model, and during the Monte Carlo simulation process it is advantageous to integrate out phonons to obtain an effective model that contains retarded electron-electron interactions, we thus need to calculate the total energy of a system without an explicit Hamiltonian formulation. It is this situation that compels us to develop the Schwinger-Dyson equation schwinger1951green ; peskin2018introduction method to derive the equation of motion as described in this paper. Finally, we make a conclusion in section IV. In the appendix, we give more details about how to integrate out phonons in periodic Anderson with electron-phonon interactions.

II Schwinger-Dyson equation method to derive equation of motion for a simplified Holstein model

In this section, we are going to present the derivation of the equation of motion for Green function in a simplified Holstein model using three different methods: the operator method, the path integral method with phonons present, and the path integral method with phonons integrated out. The operator method, although most straightforward, is restricted to the case where we have a Hamiltonian formulation of the system, a formulation that is not always available. We will next show that we can get the same results as that from operator method using path integral formulation. Once we have formulated our system using path integrals, we can integrate out the phonons and obtain an effective action that only involves electrons. We shall demonstrate that the equation of motion of the electron’s Green function, and thus the total energy of the system, remains invariant regardless whether the phonons are retained or integrated out. In the Continuous time quantum Monte Carlo simulation process, it is advantageous to integrate out the phonons and obtain a retarded potential between electrons. The results obtained in this section indicate that we can use the same equation to calculate the total energy even if the phonons are integrated out and a Hamiltonian formulation of the system is unavailable. It is noteworthy that the Schwinger-Dyson equation method applies to any model that admits a path integral formulation, especially models with electron-phonon interactions.

II.1 Model description

Our toy model can be described by this Hamiltonian:

[TABLE]

This model contains one spinless electron with energy $\epsilon$, one dispersionless phonon with energy $\omega$, and an interaction term between the electron and the phonon. The partition function for this Hamiltonian is

[TABLE]

where the action of the system iscoleman2015introduction

[TABLE]

Since the complex variables $\bar{a},a$ are quadratic in the exponent, we can perform the Gaussian integral in a closed form as

[TABLE]

where the effective action $S_{eff}$ is

[TABLE]

In the above equation for $S_{eff}$, we have defined the bare phonon propagator:

[TABLE]

Electron’s Green function for Hamiltonian that does not depend on $\tau$ explicitly is defined as (using operator notation)

[TABLE]

It can be also be equivalently defined as (using path integral formulation)

[TABLE]

Now we are going to derive the equation of motion for the electron’s Green function using both operator notation and path integral formulation, and show that these methods yield identical result.

II.2 Operator method to derive equation of motion for Green function

Derivation of the equation of motion for Green function using operator notation method is pretty straightforward. We only need to apply time derivative operator $\frac{\partial}{\partial\tau}$ to Equ. [8] and obtain this result:

[TABLE]

That is,

[TABLE]

Here, we have introduced the phonon displacement operator $\hat{\phi}=\hat{a}+\hat{a}^{\dagger}$. This method is easy to understand and straightforward to follow. However, in order to use this method, we must have a Hamiltonian formulation of the system, which is not always available. Feynman originally designed path integral formulation of quantum mechanics just to cope with the situation where the system is derived from the principle of least action and cannot be described using a Hamiltonian. We also encounter the same situation while employing CT-QMC to simulate models with electron-phonon interactions. Next, we will show that we can equally well derive the equation of motion for Green function even with the absence of an explicit Hamiltonian. The way to circumvent Hamiltonian is to use the path integral method, which we will describe in detail next.

II.3 Path integral method to derive equation of motion for Green function: with phonons present

As we can see, the equation of motion for Green’s function is easy to obtain as long as we have an operator formulation of the system. We can equivalently derive the equation of motion using path integral formulation. For this, we are to employ the Schwinger-Dyson equation. Schwinger-Dyson equation is not well known in condensed matter community, and thus we shall give this equation a brief derivation that should be easy to follow for anyone who is familiar with variational calculus. First, consider the following infinitesimal translation of the Grassmann field variables:

[TABLE]

By definition, the following quantity should be invariant before and after the field translation:

[TABLE]

However, the action $S$ does change under this translation. In the notation of the translated field variables, the action is

[TABLE]

Taking advantage of the fact that $\mathcal{D}[\bar{c},c]=\mathcal{D}[\bar{c}^{\prime},c^{\prime}]$ since the Jacobian determinant for translation is unity, we can express Eq. [13] equivalently as

[TABLE]

After cancellation of equal terms in the above equation, we finally arrive at the Schwinger-Dyson equation:

[TABLE]

This equation should hold for any field variations $\delta c,\delta\bar{c}$. In order to get the equation of motion for Green function, we set $\delta c=0$, and the integrand in the above equation simplifies into this form:

[TABLE]

Plug this back into Equ. [16], we have

[TABLE]

Since this equation holds for any $\delta\bar{c}$, we thus, from the definition of Green function in path integral formulation, have

[TABLE]

This is identical to the equation of motion that was derived using operator notation.

II.4 Path integral method to derive equation of motion for Green function: with phonons integrated out

The advantage of Schwinger-Dyson equation method for derivation of equation of motion for Green function is that it applies to any system that admits a path integral formulation. Following Feynman, we can reformulate a large family of classical Hamiltonian systems into its quantum counterpart using path integrals, however, as Feynman noted in his thesis, there are systems that can only be represented using classical action for which no Hamiltonian expression can be found. In this case, the Schwinger-Dyson equation is the only method available for the derivation of equation of motion. One such case when the system can only be represented in path integrals is the Holstein model with phonons integrated out. This is the standard approach to simulate electron-phonon interactions using quantum Monte Carlo. In order to study the thermodynamical properties of the system, we need to known the model’s entropy and free energy, and the numerical computation of these two quantities requires the calculation of energy. One method to derive model’s energy is to use the equation of motion for Green function. Thus, it is an imperative task to know how to derive the equation of motion for Green function without an explicit presence of Hamiltonian that we usually take for granted. For our simplified Holstein model, the action is

[TABLE]

where $S_{\phi}$ is the part of the action that contains phonons:

[TABLE]

By integrating out phonons, we can recast the partition function into this form:

[TABLE]

Here, Z_{\phi}=\Big{(}1-e^{-\beta\omega}\Big{)}^{-1} is the partition function for bare phonons, and $S_{eff}$ whose explicit form is already given in Eq. [6] is the effective action obtained by integrating out the phonons. Similarly, define a quantity that is invariant under the transformation $c\rightarrow c^{\prime}=c+\delta c,\bar{c}\rightarrow\bar{c}^{\prime}+\delta\bar{c}^{\prime}$, that is,

[TABLE]

Using the techniques that we employed for the derivation of Eq. [16], we arrive at the Schwinger-Dyson equation for this model as

[TABLE]

If we set $\delta c=0$ and notice that $D^{0}(\tau_{1}-\tau_{2})=D^{0}(\tau_{2}-\tau_{1})$, then we obtain the variation of the action as

[TABLE]

Thus the integrand in Equ. [24] is

[TABLE]

We demand that Eq. [24] should hold for any field variation $\delta\bar{c}$. From the principles of variational calculus, we obtain the equation of motion for Green function:

[TABLE]

The equivalence of Eq. [11], Eq. [19] and Eq. [27] is easy to see. In Eq. [27], if we set $\tau=0,\tau^{\prime}=0^{-}$, then we have the expectation value of the interaction energy as

[TABLE]

Similarly, for Eq. [11] and [19], we can also set $\tau^{\prime}=0^{-},\tau=0$ to obtain the electron-phonon interaction energy as

[TABLE]

The equivalence of Eq. [28] and Eq. [29] is obvious. Note that the total energy can be calculated by adding electronic kinetic energy which is trivial to obtain to the electron-phonon interaction energy. Thus far, we have derived the total energy formula via three distinct methods: the operator method, the path integral method with phonon present, and the path integral method with phonons integrated out. All three methods yield identical result.

In this section, we have demonstrated the applicability of Schwinger-Dyson equation method for the derivation of total energy formula for a toy Holstein model. In the next section, we will employ Schwinger-Dyson equation to derive the total energy formula for periodic Anderson model with electron-phonon interactions following the procedures that are delineated in this section.

III Schwinger-Dyson method to derive equation of motion for periodic Anderson model with electron-phonon interactions

III.1 Periodic Anderson model with electron-phonon interactions: A Brief description

Our model Hamiltonian for periodic Anderson model with electron-phonon interaction is

[TABLE]

Here, $c_{i,\sigma}^{\dagger},c_{i,\sigma}(f_{i,\sigma}^{\dagger},f_{i,\sigma})$ creates and destroys a $c(f)$ electron of spin $\sigma$ at lattice site $i$, respectively. $\hat{a}_{i}$ ($\hat{a}^{\dagger}_{i}$) is the annihilation (creation) operator for a phonon located at lattice site $i$. $U$ is the Hubbard repulsion between localized $f$-electrons, and $V$ characterizes the hybridization between $c$ and $f$ electrons. This Hamiltonian in momentum space is

[TABLE]

We are going to use the path integral method to integrate out the phonons and obtain an effective action that describes the retarded density-density interactions between electrons. With path integrals, the partition function for this Hamiltonian can be written as

[TABLE]

Integrating out the phonons, we can factorize the partition function into the form $Z=Z_{0}Z_{eff}$, with $Z_{0}=(1-e^{-\beta\omega})^{-1}$, and $Z_{eff}=\int\mathcal{D}[\bar{c},c]\mathcal{D}[\bar{f},c]e^{-S_{eff}}$. The effective action is

[TABLE]

See the appendix for a detailed derivation of the effective action. With the effective action, we can employ the Schwinger-Dyson equation to derive the equation of motion for electron’s Green function, and further calculate the total energy of the Hamiltonian. Next, we will show how to derive the equation of motion for Green function.

III.2 Equation of motion

Now that we have already obtained the effective action for our model, we shall begin to derive the equation of motion for Green function using Schwinger-Dyson equation method. Consider the variations of field variables,

[TABLE]

With this transformation, the effective action becomes $S_{eff}[\bar{c},c;\bar{f},f]\rightarrow S_{eff}[\bar{c}^{\prime},c^{\prime};\bar{f}^{\prime},f^{\prime}]=S_{eff}[\bar{c},c;\bar{f},f]+\delta S_{eff}$. Define a quantity that is invariant under this transformation,

[TABLE]

After expanding to first order and cancelling identical terms on both sides, we have

[TABLE]

Here, the variation of action is

[TABLE]

From Eq. [47], we can obtain two independent equations, which are

[TABLE]

For the first equation, we can set $\delta c=0,\delta f=0,\delta\bar{f}=0$ and obtain

[TABLE]

Plug this into the path integral, we have the equation

[TABLE]

For the second equation, we can set $\delta f=0,\delta c=0,\delta\bar{c}=0$, and then we have

[TABLE]

Plug this into the path integral, we have

[TABLE]

Combining Eq. [III.2] and [III.2] together, we have this matrix equation ($\phi$ and $\psi$ are just padding elements for the matrix equation, and their explicit forms are of no concern to us):

[TABLE]

Here, in the above equation, we have defined two quantities

[TABLE]

If we set $k=k^{\prime},\sigma=\sigma^{\prime},\tau=0,\tau^{\prime}=0^{-}$ in Equation [68], then we have

[TABLE]

Therefore, the total interaction energy is

[TABLE]

The kinetic energy is trivial to calculate, which is

[TABLE]

Combination of the kinetic energy and the interaction energy gives us the total energy of the system.

IV Conclusion

In this paper, we have demonstrated how to derive the Green function’s equation of motion without resort to an explicit Hamiltonian. The equation of motion for Green function is straightforward to get when we have Hamiltonian formulation of the system. However, Hamiltonian formulation of a system is not always available. One such case where we are bereft of an explicit Hamiltonian formulation is when we integrate out phonons to obtain a retarded electron-electron interaction in Holstein model or its generalizations. We have shown that for a simplified toy Holstein model, Green function’s equation of motion can be derived via three distinct yet equivalent methods, one of which is to integrate out the phonons and use Schwinger-Dyson equation method. We further show that for periodic Anderson model with electron-phonon interactions, we can equally well derive the equation of motion via Schwinger-Dyson equation method. Our ultimate goal is to calculate the total energy of the system, a goal that can be accomplished with some mathematical manipulations once we have already obtained the equation of motion for electronics Green functions. We showed that the well known formula for calculating the total energy of a system is still applicable here, although for our case that formula should be derived through the Schwinger-Dyson equation method.

Appendix A Integrate out phonons in periodic Anderson model with electron-phonon interaction

In this appendix, we will show how to integrate out the phonons in Periodic Anderson model with electron-phonon interactions. Once the phonons are integrated out, we are left with an effective action that only involves the electrons. In the new effective action, electrons experience a retarded attractive potential. We can interpret this as one electron is attracting another electron by emitting or absorbing a virtual phonon. The advantage of integrating out phonons is that we are left with an effective system that involves only electrons, an effective system that can be conveniently simulated using the continuous time quantum Monte Carlo (CT-QMC) method. We have been using CT-QMC to simulate the effective system that results from integrating out phonons, and it is the study of this system that motivates us to develop the Schwinger-Dyson equation method to derive the total energy formula as shown in this paper.

The partition function for our original Hamiltonian is

[TABLE]

Here, $S$ is the action which is

[TABLE]

The phononic contribution to the action is

[TABLE]

Integration with respect to phonons yields

[TABLE]

Here, the partition function for bare phonons is

[TABLE]

and the current $J_{\bm{q},m}$ is defined as

[TABLE]

Its conjugate is

[TABLE]

After integrating out the phonons, we have an effective action which is

[TABLE]

An important property of $J_{\bm{q},m}$ is that $J_{-\bm{q},-m}=\bar{J}_{\bm{q},m}$. Thus, the effective action can be rewritten as

[TABLE]

Fourier transforming back to $\tau$ space, we have

[TABLE]

Here, we have introduced the bare phonon propagator in the effective action,

[TABLE]

The effective action $S^{\phi}_{eff}$ represents the time-retarded density-density interaction between two $c$ electrons mediated by the exchange of a virtual phonon. This interaction differs from the instantaneous Hubbard interaction not only in that it is time-retarded, but also in that the spins of the two electrons do not have to be opposite, as required for on-site Hubbard interactions.

With this observation, the original partition function can be recast into the form

[TABLE]

where, the effective action is

[TABLE]

In our Monte Carlo simulation, we start with the effective action, and obtain Green functions and self-energies for retarded conduction electrons and instantaneous $f$ electrons. Now that we are bereft of a Hamiltonian formulation of the system, we cannot directly use the formula $E_{V}=\frac{1}{\beta}\sum_{\bm{k},\sigma,i\omega_{n}}\Sigma_{\bm{k,\sigma}}(i\omega_{n})G_{\bm{k},\sigma}(i\omega_{n})$ which is derived based on the assumption that the system can be described by a Hamiltonian. Here, we propose an alternative method to derive the total energy with the effective action as our starting point. The derivation of the total energy demands the knowledge of the equation of motion for the Green functions. Schwinger-Dyson equation was originally designed to obtain the equation of motion for the Green function based on the path integral principle, and thus can be readily used here to derive the total energy for our system. In the next section, we will give a brief review of the Schwinger-Dyson equation.