Effect of propagator renormalization on the band gap of insulating solids

TL;DR

This paper investigates how propagator renormalization affects the band gaps of insulating solids using self-consistent second-order perturbation theory, showing improved physical behavior for small-gap systems and providing a framework for advanced embedding methods.

Contribution

It introduces a fully self-consistent, thermodynamically consistent perturbation theory framework applicable to embedding theories like DMFT.

Findings

Large gap systems show reasonable results with both bare and self-consistent methods.

Small gap systems require self-consistent propagator adjustments.

No unphysical band gap behavior observed unlike quasiparticle-based approaches.

Abstract

We present momentum-resolved spectral functions and band gaps from bare and self-consistent second-order perturbation theory for insulating periodic solids. We establish that, for systems with large gap sizes, both bare and self-consistent perturbation theory yield reasonable gaps. However, smaller gap sizes require a self-consistent adjustment of the propagator. In contrast to results obtained within a quasiparticle formalism used on top of bare second-order perturbation theory, no unphysical behavior of the band gap is observed. Our implementation of a fully self-consistent, $Phi$-derivable, and thermodynamically consistent finite-temperature diagrammatic perturbation theory forms a framework on which embedding theories such as the dynamical mean-field theory and self-energy embedding theories can be implemented.

Tables

Table 1. Table 2: HF and GF2 total energies (E) and band gaps ( e g subscript 𝑒 𝑔 e_{g} ) of LiH for different system sizes in pob-TZVP basis

Size	E (HF)	$e_g$ (HF)	E (GF2)	$e_g$ (GF2)
3$\times$3$\times$3	-8.0659	11.82	-8.1097	10.08
4$\times$4$\times$4	-8.0629	11.44	-8.1079	7.73
5$\times$5$\times$5	-8.0618	11.31	-8.1076	6.66
6$\times$6$\times$6	-8.0612	11.24	-8.1075	5.93

Table 2. Table 3: HF and GF2 total energies (E) and band gaps ( e g subscript 𝑒 𝑔 e_{g} ) of LiF for systems of size 3 × \times 3 × \times 3 and 4 × \times 4 × \times 4 in pob-TZVP basis.

Size	E (HF)	$e_g$ (HF)	E (GF2)	$e_g$ (GF2)
3$\times$3$\times$3	-107.092	21.95	-107.325	13.36
4$\times$4$\times$4	-107.089	21.84	-107.321	13.03

Equations

$$V_{ijkl}^{{\mathbf{k}}_1{\mathbf{k}}_2{\mathbf{k}}_3{\mathbf{k}}_4}=\iint \frac{{\varphi }_{i{\mathbf{k}}_1,\mathbf{r}}^{\ast }{\varphi }_{l{\mathbf{k}}_4,\mathbf{r}}{\varphi }_{j{\mathbf{k}}_2,{\mathbf{r}}^{\prime }}^{\ast }{\varphi }_{k{\mathbf{k}}_3,{\mathbf{r}}^{\prime }}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}d\mathbf{r}d{\mathbf{r}}^{\prime },$$

$${\Sigma }_{ij}^{\mathbf{k},(2)}(\tau )$$

$$\times {\tilde{V}}_{i\mathbf{k},p{\mathbf{k}}_1}^Q$$

$$\Omega$$

$$G^{-1}$$

$$G_{HF}^{-1}$$

$$G_0^{-1}$$

$$\mathrm{Tr}\left[\ln \left[-G^{-1}\right]\right]=\mathrm{Tr}\left[\ln \left[{\Sigma }_c-G_{HF}^{-1}\right]\right]=$$

$$\mathrm{Tr}\left[\ln \left[(-G_{HF}^{-1})(-G_{HF}\Sigma +1)\right]\right]=$$

$$\mathrm{Tr}\left[\ln \left[-G_{HF}^{-1}\right]\right]+\mathrm{Tr}\left[\ln \left[1-G_{HF}\Sigma \right]\right]$$

$$\Phi \left[G\right]=\sum _{n=1}^{\infty }{\Phi }^{(n)}\left[G\right]=\sum _{n=1}^{\infty }\frac{1}{2n}\mathrm{Tr}\left[{\Sigma }^{(n)}\left[G\right]G\right],$$

$$\mathrm{Tr}\left[\Sigma \left[G\right]G\right]=\mathrm{Tr}\left[{\Sigma }^{(1)}\gamma \right]+\mathrm{Tr}\left[{\Sigma }^{(2)}G\right].$$

$${\Omega }^{[2]}$$

$$=\frac{1}{\beta }\left\{\frac{1}{2}\mathrm{Tr}\left[{\Sigma }_{\infty }\gamma \right]+\frac{1}{4}\mathrm{Tr}\left[{\Sigma }^{(2)}G\right]-\mathrm{Tr}\left[ln\left[-G_{HF}^{-1}\right]\right]-\mathrm{Tr}\left[ln\left[1-G_{HF}{\Sigma }^{(2)}\right]\right]-\mathrm{Tr}\left[{\Sigma }^{(1)}\gamma \right]-\mathrm{Tr}\left[{\Sigma }^{(2)}G\right]\right\}=$$

$$=\frac{1}{\beta }\left\{-\frac{1}{2}\mathrm{Tr}\left[{\Sigma }_{\infty }\gamma \right]-\frac{3}{4}\mathrm{Tr}\left[{\Sigma }^{(2)}G\right]-\mathrm{Tr}\left[ln\left[-G_{HF}^{-1}\right]\right]-\mathrm{Tr}\left[ln\left[1-G_{HF}{\Sigma }^{(2)}\right]\right]\right\}$$

$$S=-\frac{\partial \Omega }{\partial T},$$

$$C_V=T\frac{\partial S}{\partial T},$$

$$E=T^2\frac{\partial lnZ}{\partial T},$$

$$F=\Omega +\mu N$$

Full text

I Introduction

A truly ab-initio quantitative many-body description of weakly correlated systems beyond density functional theory (DFT) Kohn (1999); Kohn and Sham (1965) remains challenging despite enormous theoretical Hedin (1965) and computational Hybertsen and Louie (1986); Onida et al. (2002); Govoni and Galli (2015); Deslippe et al. (2012); Shishkin and Kresse (2006); Bruneval and Marques (2013) advances in the last thirty years. While such calculations are important for explaining the behavior of insulating materials, they are also necessary as a first step for many strongly correlated computational techniques such as some derivatives Biermann et al. (2005); Nilsson et al. (2017); Kananenka et al. (2015); Zgid and Gull (2017); Lan and Zgid (2017) of the dynamical mean field theory (DMFT) Georges et al. (1996); Kotliar et al. (2006); Held (2007), where the first calculation step consists of evaluating the system by a perturbative weak coupling method.

Three types of perturbation theories exist: bare (non-self-consistent) perturbation theory based on the expansion of the original Hamiltonian in interaction terms, where neither propagators nor interactions are renormalized; self-consistent perturbation theories where propagators but not the interactions are renormalized; and perturbation theories with both propagators and interactions renormalized.

Bare second order perturbation theory is known as Møller - Plesset second order (MP2) Møller and Plesset (1934); Wiberg (1986); Pople et al. (1976); Ayala et al. (2001); Pisani et al. (2008); Usvyat et al. (2011); Maschio et al. (2010); Hirata and Shimazaki (2009); Marsman et al. (2009); Grüneis et al. (2010) when applied to real materials. The self-consistent second order Green’s function perturbation theory (GF2) Phillips and Zgid (2014); Phillips et al. (2015); Rusakov and Zgid (2016); Welden et al. (2016); Kananenka and Zgid (2017); Dahlen and van Leeuwen (2005) renormalizes propagators but not interactions. While fully self-consistent GW Hedin (1965) renormalizes both propagators and interactions, its non-self-consistent variants such as G0W0 Hybertsen and Louie (1986); Aulbur et al. (2000) (partially) renormalize interactions without renormalizing propagators.

The differences in the treatment of propagators and interactions between all these types of perturbation theories are crucial since approaches that do not renormalize the interactions are expected to fail in metallic 3D systems Fetter and Walecka (2003); Grüneis et al. (2010). This breakdown is not expected to occur in insulators.

For weakly correlated materials, most results so far have been obtained within approaches such as MP2 Usvyat et al. (2011); Maschio et al. (2010); Grüneis et al. (2010); Booth et al. (2016) and GW Wilhelm and Hutter (2017); Kaltak et al. (2014); Grumet et al. (2018) at zero temperature. Only recently, finite-temperature results for solids have started to appear Kutepov et al. (2012); Kutepov (2016, 2017a, 2017b). It was demonstrated that MP2 gaps are wildly inaccurate for materials with band gaps smaller than $6$ eV Grüneis et al. (2010), leading to a breakdown of the band gap estimation for silicon and silicon carbide. In contrast, zero-temperature, non-self-consistent GW Hedin (1965) has been very successful in predicting band gaps for semiconductors. This success is usually attributed to the renormalization of the interactions by an infinite series of ‘bubble’ (RPA Hedin (1965)) diagrams – the same diagrams that render this method convergent in the metallic limit.

The renormalization of interactions and propagators is commonly discussed in many-body textbooks Fetter and Walecka (2003); Mattuck (1992), usually at the example of weakly interacting or uniform systems. Due to both the computational cost and implementation difficulties, their effect in realistic solids is difficult to explore. Since there is no obvious reason for the breakdown of the perturbative series in semiconductors and band insulators, it is interesting to explore how the renormalization with self-consistent propagators but unrenormalized interactions affects the band gaps and to compare these results to the ones obtained by MP2. Moreover, since the MP2 band gap is evaluated using approximate band energies for the lowest unoccupied and highest occupied bands Grüneis et al. (2010) reminiscent of the formulas usually employed to solve the quasi-particle (QP) equations in GW Kotani et al. (2007), it is interesting to compare these band gaps to gaps evaluated without solving the QP equations.

In this paper, we focus on quantifying the effects of the renormalization of propagators and self-energies on the value of band gaps in simple 3D solids. This is made possible by the implementation of a fully self-consistent finite temperature second order perturbation theory (GF2) that was so far not available for realistic 3D systems. Comparison between renormalized calculations, unrenormalized calculations, and experimental data also allows us to indirectly infer the contribution of higher order terms excluded from our calculations. Among these are the ‘RPA’-like screening terms included in the GW approximation.

II Method

We investigate the physics of 3D solids by choosing a finite basis set on each atom, and choosing a finite lattice of atoms periodic in all three spatial directions. This yields a periodic electronic structure Hamiltonian best expressed in a Bloch basis in reciprocal space, resulting in the Coulomb integrals

[TABLE]

where $\phi$ are basis functions in reciprocal space. Translational invariance is guaranteed by the momentum conservation $\mathbf{k}_{1}+\mathbf{k}_{2}=\mathbf{k}_{3}+\mathbf{k}_{4}$. These integrals can be decomposed Whitten (1973); Beebe Nelson and Linderberg (1977); Maschio (2011); Varga (2005) into a product of two low-rank tensors, $V_{ijkl}^{\mathbf{k}_{1}\mathbf{k}_{2}\mathbf{k}_{3}\mathbf{k}_{4}}=\sum_{Q}\tilde{V}_{i\mathbf{k}_{1},l\mathbf{k}_{4}}^{Q}\tilde{V}_{j\mathbf{k}_{2},k\mathbf{k}_{3}}^{Q}$, where $Q$ is an auxiliary basis index. For the first order diagram we employ Ewald summation to treat the divergence at zero momentum. For the second-order diagram we excluded this point to avoid the divergence of the exchange diagram.

The second-order self-energy is then evaluated in reciprocal space in an imaginary time formalism,

[TABLE]

with $\tau$ denoting imaginary time $0\leq\tau\leq\beta$ and $\beta=1/k_{B}T$ the inverse of the physical temperature $T$ (assuming Einstein summation over repeated indices).

In bare second-order perturbation theory, the Hartree-Fock Green’s function is employed in Eq. 2 and the Dyson equation is evaluated only once to yield the interacting Green’s function.

In renormalized perturbation theories such as GF2, to achieve self-consistency in both the density matrix (or, correspondingly, the Fock matrix and the frequency independent term of the self-energy $\Sigma_{\infty}$) and the dynamical self-energy $\Sigma(\tau)$, we employ a modification of the iterative procedure described in Ref. Phillips and Zgid, 2014 for molecular systems. First, we adjust the chemical potential to find the correct particle number of the Hartree-Fock (HF) solution. We then obtain the HF propagator, calculate the second-order self-energy, and recompute the interacting Green’s function and density matrix using the Dyson equation, adjusting the chemical potential until the correct density is found. This propagator is then used for the next second-order self-energy evaluation, until convergence is achieved in all quantities.

Energies, entropies, free energies, and specific heats are then computed using standard thermodynamic formulae Luttinger and Ward (1960); Holm and Aryasetiawan (2000); Fetter and Walecka (2003); Dahlen et al. (2006); Welden et al. (2016), for detailed derivations see appendix B.

Quantities obtained in diagrammatic approximations by thermodynamic integration are in general dependent on the integration path Baym and Kadanoff (1961); Baym (1962); i.e. the integration of a quantity such as the energy or the entropy may differ if it is obtained by integration from zero $T$, infinite $T$, infinite chemical potential or via coupling constant integration. So-called ‘$\Phi$-derivable’ Luttinger and Ward (1960), self-consistent methods, such as the self-consistent GF2 method investigated here or the fully self-consistent GW approximation, avoid this problem and are intrinsically thermodynamically consistent.

Standard finite temperature perturbation theories are formulated on the imaginary axis, thus, static quantities such as the density matrix, the energy, the entropy, the static magnetic susceptibility or the specific heat are directly accessible. In contrast, real-frequency-dependent quantities such as the spectral function, the gap, the optical conductivity, or the dynamical magnetic susceptibility require analytical continuation to the real axis in order to be compared to experiment. This analytical continuation is ill conditioned and leads to an amplification of uncertainties, even for data known up to numerical precision, in particular at high temperature and high frequency. The problem is intrinsic to the finite-temperature field-theory formalism on the imaginary axis and can only be overcome by reformulating the method in frequency space or real time. In the present work, we used the maximum entropy method Jarrell and Gubernatis (1996); Levy et al. (2017). Other methods, such as the spectral method Otsuki et al. (2017), Padé Wang and Chakravarty (2018), stochastic analytical continuation Mishchenko et al. (2000); Mishchenko (2012); Fuchs et al. (2010), or Consistent Constraints Goulko et al. (2017) could be explored, as could the continuation of the self-energy to obtain spectra and gaps via the quasi-particle equation Wang et al. (2009).

III Results

We analyze five solids in this paper: Ne, LiF, MgO, LiH, and diamond. The experimental band gaps of these solids are listed in Table 1, along with literature values obtained by other methods. To confirm the thermodynamic consistency of our implementation, we present the evaluation of thermodynamic properties for solid LiH in Fig. 1. Smooth curves have been obtained by Chebyshev interpolation on a temperature grid.

As these insulators have such vastly different band gaps, a different amount of Green’s function renormalization is expected to be necessary, allowing us to examine how the iterative nature of GF2 changes the results in comparison to bare perturbation theory and Hartree-Fock. Moreover, since evaluating band gaps using either the QP equation as used in Ref. Grüneis et al., 2010 or analytical continuation Jarrell and Gubernatis (1996) may give different results, we compare our values for band gaps (obtained using analytical continuation) to the ones available in the literature, where QP equations were used with the bare perturbation theory.

Our GF2 implementation for periodic systems uses a compact Chebyshev polynomial Gull et al. (2018) representation of Green’s functions that converges exponentially (for alternative techniques see Kananenka et al. (2016a, b)) and relies on the open source ALPS library Gaenko et al. (2017). We use periodic density-fitted integrals in Gaussian orbitals, evaluated using the open source pySCF package McClain et al. (2017); Sun et al. (2017). Our calculations result in a set of imaginary time self-energies on discrete $k$-points in the Brillouin zone. In order to obtain smooth spectral functions, we perform a three-dimensional periodic spline interpolation Fritsch and Carlson (1980); Fuchs et al. (2011). The validity of this procedure is assessed by repeating calculations on a set of grid sizes and obtaining convergence of the interpolated quantities. Throughout this work we show data at an inverse temperature of $\beta=100$Ha*-1* ($T\sim 3158$K $\sim 0.27$ eV). This temperature is much lower than other energies in the system, in particular much lower than the gaps of the solids discussed here, so that results can be considered to be close to the ground state.

Solid neon has a large experimental band gap of $21.7$ eV Boursey et al. (1970). Consequently, we do not expect strong ‘screening’ or a substantial renormalization of the Green’s function. The left panel of Fig. 2 shows the momentum-resolved spectral function obtained with bare and self-consistent second order techniques along a standard path in the Brillouin zone (from $\Gamma$ via $X$ and $W$ back to $\Gamma$) for this system, illustrating that self-consistency leads to a negligible change of the Green’s function renormalization. Data is evaluated in the 6-311+G∗ basis of Gaussian orbitals Krishnan et al. (1980); Clark et al. (1983). Thick lines denote points in the Brillouin zone that coincide with our momentum grid. Thin lines denote interpolated values. All spectral functions are plotted as a function of frequency in eV.

In this system, the band gap and the spectral functions are converged with respect to the momentum discretization. This is illustrated in the middle panel of Fig. 2, where we show results on a 3$\times$3$\times$3 and on a 4$\times$4$\times$4 lattice along the $\Delta$ direction in momentum space. Thick lines denote values on the respective momentum grids (black for 4$\times$4$\times$4, green for 3$\times$3$\times$3), thin lines are obtained by interpolation. The $X$ point is absent on the smaller grid. Data has been obtained in the 6-311+G* basis set.

Our calculations also demonstrate that, in the frequency window shown, the spectral function is relatively insensitive to the choice of the basis set, see right panel of Fig. 2. However, it should be stressed that while our results do not show significant differences between aug-cc-pVDZ Dunning (1989) and 6-311+G∗ basis sets, they may not be converged with respect to the basis set size, since both bases are small. Converging our calculations with respect to the basis set size would require a systematic increase of the cardinal number $X$ in the series of aug-cc-pV$X$Z basis sets. This is exceedingly difficult in ordinary solid state calculations since regular Gaussian basis sets such as aug-cc-pV$X$Z become linearly dependent for higher cardinal numbers.

Hartree-Fock gaps can be extracted directly from the eigenvalues. The determination of the gap in correlated methods via the spectral function leaves some arbitrariness, as finite temperature Green’s functions are broadened by self-energy, temperature, and analytical continuation artifacts. For GF2, we chose to define the band gap as the peak-to-peak distance of the k-space peaks closest to the Fermi energy from above and below. For solid neon, GF2 in both the aug-cc-pVDZ and the 6-311+G∗ basis yields a band gap of $20.55$ eV at the $\Gamma$ point. We find that the difference between the GF2 gap ($20.55$ eV) and the experimental band gap ($21.7$ eV) is consistent with the previous bare perturbation theory studiesGrüneis et al. (2010).

Solid LiF also has a wide experimental band gap of $14.2$ eV Piacentini et al. (1976). For the 4$\times$4$\times$4 k-point grid, we observe a band gap of $13.03$ eV in self-consistent GF2, as extracted from the peak-to-peak distance of the spectral function. The difference of the spectral function between the first iteration and self-consistent GF2 is negligible, see left panel of Fig. 3. For k-grid convergence of the unit cell energy and band gap see Tab. 3 and Fig. 5 in appendix A.

Both LiH and diamond have moderate band gaps. Therefore, a significant change of the band gap between the first iteration of GF2 and the fully self-consistent result is expected. Indeed, both the middle and the right panel of Fig. 3 confirm that the band gaps obtained in the initial GF2 iteration are much wider than the self-consistent result, and gradually shrink during the self-consistent iteration progress. Here, to evaluate spectral functions during the first iteration of GF2, we use the analytical continuation of $G(i\omega)=[(i\omega+\mu)1-F_{\text{HF}}-\Sigma^{(2)}]^{-1}$, where $\Sigma^{(2)}$ was obtained in the first iteration of the GF2 method, i.e. using $G_{\text{HF}}(i\omega)=[(i\omega+\mu)1-F_{\text{HF}}]^{-1}$ as the propagator in Eq. 2. The Fock matrix $F_{\text{HF}}$ comes from a preceding HF calculation. Consequently, the first iteration of GF2 lacks two types of renormalization: first, the renormalization coming from the self-consistently updated $\Sigma_{\infty}$ and then consequently updated Fock matrix and, second, the one from the fully self-consistent evaluation of $\Sigma^{(2)}$ which in subsequent iterations is evaluated with renormalized propagators.

For diamond in the pob-TZVP basis, we observe an indirect band gap of $4.8$ eV between $\Gamma$ and the point halfway between the $\Gamma$ and $X$ points. The direct band gap at the $\Gamma$-point is about $6.6$ eV. The positions and values of the direct and indirect band gaps are in good agreement with previous experimental and theoretical results. Phillip and Taft (1964); Saravia and Brust (1968); Clark et al. (1964); Saslow et al. (1966) It is worth emphasizing that the band gap for diamond obtained from self-consistent GF2 is $4.8$ eV, while the MP2 band gap obtained using a QP formalism listed in Ref. Grüneis et al., 2010 is $1.9$ eV. This underestimation of band gaps smaller than $6$ eV is very noticeable (as listed in Ref. Grüneis et al., 2010) for bare second order perturbation theory coupled with band gap evaluation based on the QP formalism. This deficiency seems to be avoided when the fully self-consistent GF2 is employed, as we observe both in the case of LiH and of diamond.

In order to compare results to implementations of MP2, scGW, and QSGW, we also present results for MgO, which has an experimental band gap of $7.8$ eV.

For LiH, we plot the convergence of the band gap and unit cell energy for different k-grids in Tab. 2 and in Fig. 4. The data clearly show that while the unit cell energy converges rather quickly, the convergence of the band gap is less rapid. While the width of the obtained band gap is bigger than the experimentally observed value, we did not achieve convergence with system size. An extrapolation of the gap with inverse system size yields a gap value close to $4.0$ eV.

IV Conclusions

In conclusion, we illustrated the effect of propagator renormalization on the example of gaps and momentum-resolved spectral functions for Ne, LiH, MgO, LiF, and diamond. In all of these cases, we found reasonable agreement of GF2 with experimental values. Wide-gap insulators were found to be fully converged w.r.t the k-point grid, whereas larger momentum grids are needed to converge the band gap (but not the total energy) of LiH. We showed that thermodynamic consistency is obeyed in our calculations, opening the door for systematic thermodynamic calculations of the electronic system of real materials.

Our study reveals three major results. First, the comparison between our results for bare perturbation theory and self-consistent perturbation theory to experiment shows that renormalized propagator diagrams are responsible for most of the difference between bare results and the experiment. This illustrates the importance of propagator renormalization to obtain the reasonable band gaps for insulators. Second, the fact that our results from bare second order perturbation theory yield reasonable gap values that differ substantially from the published MP2 values indicates a breakdown of the QP formalism for gap extraction. The approximations inherent to this formalism allow the expression of the results in a convenient real-frequency ‘band’ picture but, in light of our discrepancies, will need to be revisited. Finally, our results show that controlled self-consistent diagrammatic many-body calculations in standardized Gaussian basis sets are now routinely possible.

Acknowledgements.

S.I. and E.G. acknowledge support by the Simons Foundation via the Simons Collaboration on the Many-Electron Problem. D.Z. and A.R. are supported by NSF-CHE-1453894. E.G and D.Z thank the Simons Foundation for sabbatical support. The Flatiron Institute is a division of the Simons Foundation. S.I. thanks Miguel A. Morales for fruitful discussion. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A System-size convergence for LiF

Tab. 3 shows convergence of the unit cell energy and band gap with momentum grid. Fig. 5 shows the band structure obtained from two momentum grids.

Appendix B Evaluation of thermodynamic properties

The grand potential is defined in terms of Green’s functions, self-energies, and a $\Phi$-functional as Luttinger and Ward (1960)

[TABLE]

In practice, the evaluation direct evaluation of the second term in this form is complicated by the slow decay of Green’s functions as a function of frequency Dahlen et al. (2006). We therefore by defining

[TABLE]

in terms of Matsubara frequencies $i\omega_{n}$, the chemical potential $\mu$, the overlap matrix $S$, the non-interacting $(V=0)$ Hamiltonian $H_{0}$ and the Fock matrix $F$ can evaluate the logarithmic term as Dahlen et al. (2006)

[TABLE]

The $\Phi$-functional is expressed as

[TABLE]

where $\Sigma^{(n)}\left[G\right]$ is the total n-th order self-energy part. Within the second-order theory, the self-energy is approximated as $\Sigma\approx\Sigma^{[2]}=\Sigma^{(1)}+\Sigma^{(2)}$, and using Eq. (6) $\Phi\approx\Phi^{[2]}=\Phi^{(1)}+\Phi^{(2)}=\frac{1}{2}\mathrm{Tr}\left[\Sigma^{(1)}\gamma\right]+\frac{1}{4}\mathrm{Tr}\left[\Sigma^{(2)}G\right]$. Using Eqs. (4) $\Sigma^{(1)}=F-H_{0}$, such that

[TABLE]

Combining Eq. 5-Eq. 7 we obtain the second order approximation of the grand potential as

[TABLE]

Standard textbook relations yield the entropy, specific heat, total energy and free energy as thermodynamic derivatives

[TABLE]

Alternatively, the entropy can be evaluated from the Gibbs-Duhem relation $\Omega=E-TS-\mu N$; the specific heat from its definition $C_{V}=\frac{\partial E}{\partial T}$; and the energy from the Galitskii-Migdal formula Holm and Aryasetiawan (2000).