Parquet approximation for molecules: spectrum and optical conductivity of the Pariser-Parr-Pople model

TL;DR

This paper applies the parquet approximation to the Pariser-Parr-Pople model of conjugated molecules, revealing differences from Hubbard models and emphasizing the importance of non-local interactions for optical properties.

Contribution

It demonstrates the effectiveness of the parquet approximation for molecular systems with non-local interactions, providing insights into spectral functions and optical conductivity.

Findings

Weak quasiparticle renormalization compared to Hubbard model

Static self-energy increases the band gap in PPP

Vertex corrections significantly impact optical conductivity

Abstract

We study a simple model system for the conjugated $\pi$-bonds in benzene, the Pariser-Parr-Pople (PPP) model, within the parquet approximation (PA), exemplifying the prospects of the PA for molecules. Advantages of the PA are the polynomial scaling with the number of orbitals, and the natural calculation of one- and two-particle spectral functions as well as of response and correlation functions. We find large differences in the electronic correlations in the PPP model compared to a Hubbard model with only local interactions. The quasiparticle renormalization (or mass enhancement) is much weaker in the PPP than in the Hubbard model, but the static part of the self-energy enhances the band gap of the former. Furthermore, the vertex corrections to the optical conductivity are much more important in the PPP model. Because non-local interactions strongly alter the self-energy, we conclude that the PA is more suitable for calculating conjugated $\pi$-bonds in molecules than single site dynamical mean-field theory.

Tables

Table 1. Table 1: Parameters of the Fermi-liquid-like renormalization, i.e., quasiparticle energy ϵ 𝐤 ∗ subscript superscript italic-ϵ 𝐤 \epsilon^{*}_{\mathbf{k}} and quasiparticle weight Z 𝑍 Z compared to the non-interacting energies ϵ 𝐤 subscript italic-ϵ 𝐤 \epsilon_{\mathbf{k}} and the (predominant) peaks ϵ P ​ A ​ ( E ​ D ) subscript italic-ϵ 𝑃 𝐴 𝐸 𝐷 \epsilon_{PA(ED)} of the PA (ED) from Fig.  5 . For better comparison the non-interacting dispersion relation ϵ 𝐤 subscript italic-ϵ 𝐤 \epsilon_{\mathbf{k}} and the real part of the self-energy at ν 0 = π / β subscript 𝜈 0 𝜋 𝛽 \nu_{0}=\pi/\beta is also listed.

		$U$-only					PPP
$𝐤$	${ϵ}_{𝐤}$	${ϵ}_{𝐤}^{\ast }$	${ϵ}_{\text{PA}}$	${ϵ}_{\text{ED}}$	$\text{Re}\Sigma$	$Z$	${ϵ}_{𝐤}^{\ast }$	${ϵ}_{\text{PA}}$	${ϵ}_{\text{ED}}$	$\text{Re}\Sigma$	$Z$
$0$	-2	-2.0	-1.9	-2.0	-0.38	0.83	-3.3	-3.1	-3.3	-1.47	0.96
$\frac{\pi }{3}$	-1	-1.1	-1.1	-1.3	-0.34	0.80	-2.1	-2.1	-2.2	-1.23	0.94
$\frac{2\pi }{3}$	1	1.1	1.1	1.3	0.34	0.80	2.1	2.1	2.2	1.23	0.94
$\pi$	2	2.0	1.9	2.0	0.38	0.83	3.3	3.1	3.3	1.47	0.96

Equations

$$\mathcal{H}=-t\sum _{i\sigma }{c}_{i\sigma }^{†}{c}_{i+1\sigma }^{}+\frac{1}{2}\sum _{ij,\sigma {\sigma }^{\prime }}{V}_{ij}(n_{i\sigma }-1)(n_{j{\sigma }^{\prime }}-1),$$

$${\Sigma }_k$$

$$U_d^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_m^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_s^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$U_t^{\mathbf{k}{\mathbf{k}}^{\prime }\mathbf{q}}$$

$$F_d^{k{k}^{\prime }q}$$

$${\Phi }_{d/m}^{k{k}^{\prime }q}$$

$$\to$$

$${\Phi }_{s/t}^{k{k}^{\prime }q}$$

$$\to$$

$$\begin{array}{rl}{\lambda }_{d/m}^q{\varphi }_{d/m}^{kq}& =\frac{1}{N\beta }\sum _{k_1}{\Gamma }_{d/m}^{k{k}_{1}q}{G}_{k_1}{G}_{k_1+q}{\varphi }_{d/m}^{k_{1}q}\\ {\lambda }_{pp}^q{\varphi }_{pp}^{kq}& =\frac{1}{2N\beta }\sum _{k_1}{\left[{\Gamma }_s-{\Gamma }_t\right]}^{k{k}_{1}q}{G}_{k_1}{G}_{q-k_1}{\varphi }_{pp}^{k_{1}q}.\end{array}$$

$${\Sigma }_{\mathbf{k}}(\nu )\approx {\text{Re}{\Sigma }_k|}_{\nu \to 0}+i\nu {\frac{\partial \mathrm{Im}{\Sigma }_k}{\partial \nu }|}_{\nu \to 0}.$$

$$Z_{\mathbf{k}}={\left[1-{\frac{\partial \mathrm{Im}{\Sigma }_k}{\partial \nu }|}_{\nu \to 0}\right]}^{-1}.$$

$${ϵ}_{\mathbf{k}}^{\ast }=Z_{\mathbf{k}}\left[{ϵ}_{\mathbf{k}}+{\text{Re}{\Sigma }_k|}_{{\nu }_n\to 0}\right].$$

$$\sigma (\omega )$$

$$\chi ={\int }_0^{\beta }d\tau ⟨j_{\mathbf{q}}(\tau )j_{-\mathbf{q}}(0)⟩e^{i{\omega }_n\tau },$$

$$j_{\mathbf{q}}(\tau )=$$

$$\begin{array}{rl}\chi & =\chi +\chi \\ & =\frac{2}{\beta N}\sum _k{\left[{\gamma }_{\mathbf{k}}^{\mathbf{q}}\right]}^2{G}_{q+k}{G}_k\\ & +\frac{2}{(\beta N{)}^2}\sum _{k,k^{\prime }}{\gamma }_{\mathbf{k}}^{\mathbf{q}}{\gamma }_{{\mathbf{k}}^{\prime }}^{\mathbf{q}}{G}_{k^{\prime }}{G}_{q+k}{F}_d^{k{k}^{\prime }q}{G}_{q+k^{\prime }}{G}_k,\end{array}$$

Full text

I Introduction

The calculation of strongly correlated electron systems is particularly challenging1 since strong correlations imply that many-body perturbation theory (MBPT), e.g., up to first or second order in the interaction $V$, is not a reliable approximation. Even relatively small systems pose a problem. An exact diagonalization (ED)2; 3; 4 of the Hamiltonian, also coined full configuration interaction (CI)5 in chemistry, is still possible for $N={\cal O}(20)$ states (orbitals times spins). Quantum Monte Carlo simulations (QMC)6, on the other hand, suffer from the Fermionic sign problem, restricting this method to $N={\cal O}(100)$.

Against this background, the methods of choice in solid state physics, where (infinitely extended) periodic crystals are typically considered, are density matrix renormalization group (DMRG)7; 8 and dynamical mean field theory (DMFT).9; 10; 11 The former is favorable for one-dimensional systems which can be typically described by matrix product states (MPS) with relatively small matrices. The latter becomes exact for infinite dimensions and is a good approximation for three-dimensional problems, where a local DMFT self-energy often yields a good description. Recently, DMFT has been extended diagrammatically12 to account for non-local correlations beyond the local DMFT self-energy; this is done in the dynamical vertex approximation (D$\Gamma$A),13; 14; 15 the dual fermion (DF)16 and related approaches17; 18; 19; 20; 21; 12. In chemistry, on the other hand, the coupled cluster (CC) method22; 23 has become the standard approach.

The CC method is a resummation of Feynman diagrams but in contrast to MBPT up to infinite orders in $V$. However, only particular diagrams are considered. These can be associated with single, double, or triple excitations from a (typically) Hartree-Fock background and certain copies thereof. Taking single and double excitations into account goes under the name of CCSD, whereas also taking triple excitations into account is coined CCSDT.23

Another resummation of Feynman diagrams is better known in the field of physics: the parquet approximation (PA).24 It is based on an exact set of equations, the parquet equations (PE). The PE require as (only) input the two-particle fully irreducible vertex $\Lambda$ and the non-interacting, one-particle Hamiltonian. If we approximate $\Lambda$ by the bare interaction, $\Lambda=V$, we obtain the PA,24 which includes not only all diagrams up to third order in $V$ (i.e., corrections are $\sim V^{4}$) but in addition some diagrams up to infinite order in $V$ (including all ladder diagrams). The parquet equations have been developed in the 60s 25; 26 and revisited in the 90s,27 at a time when sufficiently powerful computers were not yet available. More recently we have seen a revival in the context of diagrammatic extensions of DMFT, specifically in the context of D$\Gamma$A, and for disordered problems.28; 29; 30; 31; 32; 33; 34

In the (parquet) D$\Gamma$A,31; 32 $\Lambda$ is approximated by a local fully irreducible vertex $\Lambda=\Lambda_{\rm loc}$ (supplemented by the bare non-local Coulomb interaction) which can be calculated numerically by solving an Anderson impurity problem with continuous-time quantum Monte Carlo simulations.35; 36; 37; 38 In this context it has also been recognized that the PA works quite well for weak-to-intermediate Coulomb interactions,32 i.e., when we are still quite far below the Mott transition. Hitherto these studies have been mainly focused on the Hubbard model with only a local Coulomb interaction both in one31 and two32 dimensions. Notwithstanding, a non-local interaction can be added to the Hamiltonian of the system without changing the form of the parquet equations.

In the present paper, we employ the PA for one of the simplest but nonetheless relevant models with local ($U$) and non-local interactions ($V$): the Pariser-Parr-Pople (PPP) model39; 40 for a benzene ring. It consists of six orbitals (one $p_{z}$ orbital on each carbon site) or $N=6\times 2$ states with a nearest-neighbor hopping and both local and non-local Coulomb interactions. It has the advantage that an exact ED solution is still possible. For comparison, we show results for the six site Hubbard model31 ($U$-only model) which has the same Hamiltonian except for neglecting all non-local Coulomb interactions $V$, as well as for the intermediate case, the extended Hubbard model with only the nearest-neighbor interaction ($U+V_{1}$ model).

The paper is organized as follows: In Section II, we introduce the PPP model including the values of the parameters used. In Section III the methods employed, i.e., the PE and the PA, are briefly discussed, as well as the necessary changes to the victory 41 code solving the PE. The results obtained are presented in Section IV: We start in Section IV.1 with the structure of the two-particle vertex function and changes thereof when including non-local interaction in the PPP model. Further we discuss in Section IV.2 the leading instabilities (largest susceptibilities) on the basis of the largest eigenvalues in the corresponding channel. Section IV.3 is devoted to the self-energy and the one-particle spectral function; and Section IV.4 to the optical conductivity. Finally Section V provides a brief summary and outlook.

II Pariser-Parr-Pople model

The Pariser-Parr-Pople (PPP) model 39; 40 was developed as a realistic albeit simple one-band model of the conjugated $\pi$-bonds in carbon-based organic molecules. The model contains a $p_{z}$-$p_{z}$ hopping term $t$ between neighboring carbon sites and a long-range Coulomb interaction $V_{ij}={U}/{\sqrt{1+\alpha|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|^{2}}}$. Here, $U$ is the local Coulomb interaction, $\alpha=[U/(14.397\;\text{eV\AA})]^{2}$, and $d_{ij}=|{\mathbf{r}}_{i}-{\mathbf{r}}_{j}|$ the distance between carbon site $i$ and $j$. The resulting Hamiltonian has the form

[TABLE]

where $n_{i\sigma}=c^{\dagger}_{i\sigma}c_{i\sigma}$, and $c^{\dagger}_{i\sigma}$ and $c_{i\sigma}$ represent creation and annihilation operators for an electron with spin $\sigma$ at carbon site $i$, respectively.

We restrict ourselves to the PPP model for a benzene molecule (cf. left panel of Fig. 1), i.e., Hamiltonian (1) lives on a ring with six sites with periodic boundary conditions.

The second term in Eq. (1) contains the effective density-density interaction, including the chemical potential shift that makes the system particle-hole symmetric with $n=1$ electrons per site. In the geometry of the benzene ring $d_{i(i+1)}=1.4$ Å, $d_{i(i+2)}=\sqrt{3}\times 1.4$ Å, and $d_{i(i+3)}=2\times 1.4$ Å. Within this semiempirical approach, the parameters of the model were obtained in Ref. 42 by fits to experimental data of the low-lying excitation energies of benzene, yielding $t=2.539\,$ eV, $V_{ii}=3.962t$, $V_{i(i+1)}=2.832t$, $V_{i(i+2)}=2.014t$, and $V_{i(i+3)}=1.803t$. In the literature, this model for benzene has been solved within ED 43 and many-body perturbation theory schemes. 44; 45; 46; 47

Besides the PPP model, we also study the corresponding Hubbard model ($U$-only model) by setting $V_{ij}=U\delta_{ij}$, and the extended Hubbard model ($U+V_{1}$ model) with nearest-neighbor interaction with $V_{ii}=3.962t$, $V_{i(i+1)}=2.832t$ only. The interaction potential of these models is likewise displayed in the right panel of Fig. 1. In the following, we set our unit of energy to the hopping element, i.e., $t\equiv 1$.

The Fourier transformation to momentum (${\mathbf{k}}$) space yields the one-particle energies $\epsilon_{{\mathbf{k}}_{i}}=-2t\cos{\mathbf{k}}_{i}$ 48, with ${\mathbf{k}}_{i}\in\{0,\pi/3,2\pi/3,\pi,4\pi/3,5\pi/3\}$; and interactions $V_{\mathbf{q}}=\sum_{l\neq 0}V_{0l}e^{il\mathbf{q}}$. Here we have excluded the local interaction $U=V_{00}$ as it is treated separately in the following. In this work the temperature is set to $T=0.1t$.

The benzene molecule belongs to the symmetry group $D_{6h}$.49 The invariance of benzene under a $C_{6}$ rotation, which is an element of this group, yields a local one-particle Green’s function $G_{ii}$ which is equal for all sites. In momentum space, the reflection symmetry of benzene yields the equivalence $\mathbf{k}\!=\!\frac{\pi}{3}\leftrightarrow\mathbf{k}=\frac{5\pi}{3}$, and likewise, $\mathbf{k}=\frac{2\pi}{3}\leftrightarrow\mathbf{k}=\frac{4\pi}{3}$. For $n\!=\!1$, the invariance under particle-hole transformation yields the symmetry $\mathbf{k}\leftrightarrow\mathbf{k}\pm\pi$. Therewith the number of independent functions reduces. For the one-particle functions, for instance, we can restrict ourselves to consider only two momenta, e.g., $\mathbf{k}=0$ and $\mathbf{k}=\frac{\pi}{3}$, all other properties are obtained by means of the symmetries above.

III Method: extending victory to treat non-local interactions

The physical quantities studied in this paper, such as transport properties or scattering rates, require the determination of both, one- and two-particle, Green’s functions. Hence, a theory that treats one- and two-particle Green’s functions on an equal footing and in a self-consistent way is preferable. To this end, we apply the parquet formalism, 50 which is a set of self-consistently coupled equations, namely the Bethe-Salpeter equation, the parquet equation51; 26; 52 and the Schwinger-Dyson equation. A state-of-the-art implementation is available in the recent victory release. 41 This package is an efficient parquet solver for electron systems in one or two dimensions, which properly takes the high-frequency asymptotics into account via so-called kernel functions. 32; 53

While the parquet equations are an exact set of equations, they require the fully irreducible two-particle vertex $\Lambda$ as an input. But in general $\Lambda$ is not known. It can be approximated by its local contribution (known as parquet-variant of the dynamical vertex approximation, D$\Gamma$A 13; 31; 32; 12) which needs to be calculated by means of other methods, such as ED and QMC simulations. Or one can take the bare Coulomb interaction as $\Lambda$, which is also known as parquet approximation (PA). 27; 24

The parquet equations are numerically quite involved, and require an iterative solution scheme. In particular, since the crossing relations between vertex functions in different channels involve combination of momentum and frequency indices, the range in which the knowledge of the vertex functions is required increases at each iteration of the parquet equation (see below). Due to translational invariance (periodic boundary conditions), the momenta can be restricted to the first Brillouin zone. The frequency, however, must be kept in an, in principle, infinite range. The so-called kernel approximations handle this issue, by taking care of the vertex asymptotics in the frequency space, see Ref. 32. It is beyond the scope of the present article to introduce and derive the parquet equations in full, and we refer the reader to Refs. 24; 12; 54; 32 instead.

Hitherto only a local Hubbard interaction has been treated in victory. Dealing with non-local interactions requires a modification of the self-energy calculation by further contributions, i.e., the one-particle Fock term, a constant Hartree term $(U/2+V_{\mathbf{q}=0})n$ (with the density fixed to $n=1$ electrons per site), and by a two-particle vertex contribution:

[TABLE]

Here, $N$ ($N=6$ for the benzene molecule) corresponds to the total number of points in the discretized momentum-space in the first Brillouin zone, and $\beta=1/T$ is the inverse temperature. The one-particle Green’s function is labeled by $G$, and the two-particle full vertex function by $F$. The vertex function $F$ is directly linked to the connected part of the two-particle Green’s function by cutting (“amputating”) four one-particle Green’s functions from the latter. We employ a combined four-vector notation $k=(\mathbf{k},\nu_{n})$ for the momentum $\mathbf{k}$ and the fermionic Matsubara frequency $\nu_{n}=(2n+1)\pi/\beta$, and $q=(\mathbf{q},\omega_{n})$ for the (one-dimensional) momentum $\mathbf{q}$ and the bosonic Matsubara frequencies as $\omega_{n}=2n\pi/\beta$, with $n\in\mathbb{Z}$. Due to the $SU(2)$-symmetry of the Hamiltonian, the vertex function $F$ can be written in a spin-diagonalized form, with a density ($d$) and magnetic ($m$) channel; cf. Refs. 55; 12.

Beside the modification in the self-energy, the momentum-dependent first-order contribution $V_{\mathbf{q}}$ has to be added explicitly to the fully irreducible vertex function $\Lambda$. This is obvious for the PA, but it is also used as an approximation in the ab initio D$\Gamma$A which only takes into account local vertex corrections besides the bare interaction in $\Lambda$. 56; 21 Both cases are treated in our extension of the victory code, by adding the bare non-local interaction to the local $\Lambda$. In terms of formulas, we replace the bare interaction vertex $U$ in the PA by

[TABLE]

where the subscripts $s$ and $t$ refer to the singlet and triplet channels which are conveniently introduced in the spin-diagonalized notation. 32

The parquet equation read, e.g., in the density channel (cf. Refs. 32; 41 for all channels and further details)

[TABLE]

where $\Phi_{r}$ denotes the reducible vertex function in a given parquet channel $r=d,m,s,t$. The evaluation of the parquet equations requires frequency and momenta outside the box for which $\Phi$ is stored 32, due to frequency-momenta combinations such as $k^{\prime}-k$ and $k+k^{\prime}+q$ in Eq. (4). While the periodic boundary conditions resolve this issue for the momenta, there is no periodicity in frequency space. For the frequencies we hence calculate the asymptotic kernel functions by scanning the surface of the known frequency box as described in Ref. 32. However, since the bare vertex now depends on $\mathbf{k}$, $\mathbf{k}^{\prime}$, and $\mathbf{q}$, so does the high frequency asymptotics of the reducible vertex functions $\Phi$. The reducible vertex function in a given channel $\Phi_{r}$ ($r=d,m,s,t$) is related through the Bethe-Salpeter equation to the full vertex $F$ and the irreducible vertex in channel $r$, $\Gamma_{r}$. The asymptotics of $\Phi$ can thus be obtained 12; 57 by replacing $F$ and $\Gamma$ in the Bethe-Salpeter equation by $U$, yielding

[TABLE]

We utilize this high-frequency behavior when solving the Bethe-Salpeter equations by extending the frequency sum over $\nu_{1}$ (as part of the $k_{1}$ sum) outside the box for which $F$ and $\Gamma$ are stored, i.e., using Eqs. (6) and (8) instead of Eqs. (5) and (7). These terms are actually evaluated as a difference and in imaginary time, see Ref. 41. A similar regularization (extension to high frequencies) is employed for the Dyson-Schwinger equation (2).

With these changes we can employ victory for studying non-local interactions, keeping in mind that a larger number of frequencies $N_{f}$ are needed compared the case with local interactions only, due to the slower convergence to the asymptotic behavior.

IV Results

IV.1 Two-particle vertex function

Let us start by inspecting the two-particle full vertex functions which are obtained in the PA. Fig. 2 shows the full vertex function $F^{kk^{\prime}q}_{d}$ for the benzene ring, comparing the PPP and $U$-only model. Here, we concentrate on the density channel ($d$) since it is the relevant channel for the optical conductivity. The magnetic channel gives a qualitatively similar analysis.

Comparing the upper and lower row of Fig. 2, we see that including a non-local interaction in the PPP model causes a strong momentum selectivity of the full vertex function. The main difference between the full vertex of the PPP and $U$-only model stems from the additional bare non-local interaction, i.e., from $\Lambda_{d}=U+2V_{\mathbf{q}}-V_{\mathbf{k}^{\prime}-\mathbf{k}}$ in Eq. (4). This frequency-independent but non-local ($\mathbf{q}$-dependent) interaction causes quite different background values for the different momentum patches (different colors in Fig. 2). For the $\mathbf{q}=0$ contribution shown in Fig. 2, this leads to a clear main diagonal structure ${\mathbf{k}}={\mathbf{k}}^{\prime}$ in the vertex of the PPP model.

On top of this, in the patches with $q-k\mp k^{\prime}=\textit{const.}$, we observe main and secondary diagonal structures also in the frequency subindex of Fig. 2, for both models. These stem from vertex diagrams where the external legs pairwise connect to the same interaction term, 12 which also determine the asymptotics in the first kernel approximation, i.e., Eqs. (6) and (8). Additionally, when $q\pm k^{(\prime)}=\textit{const.}$, there is a hardly discernible plus structure which stems from vertex diagrams where only one pair of the external legs connect to the same interaction term 12 (they constitute the second kernel function 32).

Except for the constant background given by the bare interaction, it seems that the frequency structure of the vertex for the different models is quite comparable. Picking for instance the momentum patch $(\mathbf{k},\mathbf{k}^{\prime})=(\pi,0)$ for $q_{0}=(0,0)$ in the third column of Fig. 2, where the constant background $U+V_{\mathbf{q}}$ is already subtracted, the frequency structures shown look similar. On a quantitative level, the side diagonal $\nu_{n}=-\nu^{\prime}_{n}$ is somewhat sharper for the $U$-only model than for the PPP model, which instead shows a broader structure. The latter characteristic is the reason why an accurate simulation of the PPP model requires large frequency boxes. This aspect is even more important for the bosonic frequency $\omega_{n}$, as the kernel functions are set to zero outside of the inner bosonic frequency box.

IV.2 Eigenvalues of the BSE

Valuable insights on the physics of the benzene ring can be obtained by analyzing the eigenvalues of the Bethe-Salpeter equation, i.e., the eigenvalues $\lambda$ of

[TABLE]

Here, we have for each $q$, a characteristic equation 111evaluated for the $N_{f}$ frequencies of the inner frequency box for the eigenvalue $\lambda^{q}$, or more precisely three equations for the density and magnetic ($d$/$m$) channels, as well as for the particle-particle ($pp$) channel. For the latter the given combination of singlet and triplet ($s$/$t$) channels allows for a simultaneous study of all superconducting instabilities. In the following, we consider the value of $q_{\text{max}}=(\mathbf{q}_{\text{max}},\omega_{\text{max}})$, corresponding to the leading eigenvalue $\lambda^{q_{\text{max}}}$ in the respective channel. The leading eigenvalues are shown in Fig. 3. For all three models considered, the dominating instabilities are in the $d$ and $m$ channels, which display the largest eigenvalues at $q_{\text{max}}=(\pi,0)$, while the leading eigenvalues in the $pp$ channels, at $q_{\text{max}}=(0,0)$, are much smaller.

Within the $U$-only model, $\lambda_{m}$ is the largest among the leading eigenvalues in the different channels, thus revealing a strong tendency toward spin ordering, with an antiferromagnetic (AF) pattern associated to $\mathbf{q}_{\text{max}}=\pi$. In contrast, within the $U+V_{1}$ model, the system is at the verge of a charge density wave (CDW) instability ($\lambda_{d}\approx 1$). The corresponding wave vector is again $\mathbf{q}_{\text{max}}=\pi$, i.e., there is a strong tendency toward an alternating occupation of sites. Such an ordering is clearly favored by the nearest neighbor repulsive interaction $V_{1}$. The PPP model is somehow in between the previous scenarios, with a similarly large $\lambda$ in both the $m$ and $d$ channels, since the long-range tail of the Coulomb repulsion ”frustrates" both the AF and CDW staggered order. Finally, for all three models, the eigenvalue $\lambda_{pp}$ indicates a pairing tendency at $\mathbf{q}=0$. This is slightly enhanced for the $U+V_{1}$ model, but very far from a superconducting instability for all models.

IV.3 Self-energy and spectra

Let us now turn to the self-energy shown in Fig. 4 for the PPP model and the $U$-only model. Here, the Fock (and Hartree) term is subtracted for a better comparison of the two models. As an exact diagonalization (ED) for a one-band Hamiltonian with six sites is still feasible, the ED results for the self-energy are displayed in Fig. 4 as well.59 Further we show the (${\mathbf{k}}$-independent) DMFT self-energy obtained by continuous time Monte Carlo simulations in the interaction expansion.222Here, we use the interaction expansion continuous time Monte Carlo (CT-QMC) method86; 35 in the implementation of Gang Li.

In the top panel of Fig. 4, we show the imaginary part of the self-energy for two $\mathbf{k}$-points. In the case of the $U$-only model, the DMFT self-energy provides a good description of the system, as it roughly corresponds to a $\mathbf{k}$-average of the ED data. The PA slightly underestimates the self-energy, but it nevertheless provides the correct tendency with respect to the momentum dependence, i.e., the self-energy is (in absolute terms) larger at $\mathbf{k}=(\pi/3)$ than at $\mathbf{k}=(0)$. This confirms previous calculations, 61; 31; 32 showing that the PA gives excellent results in the weak-to-intermediate coupling regime for the Anderson impurity model. Stronger deviations may instead be expected toward strong coupling. As long as the DMFT solution gives a reasonably correct description of the local physics, the parquet results can be improved by using the local vertex $\Lambda_{\rm loc}$ instead of the bare $U$ as the starting point for the parquet equations 61; 31; 32, i.e., using the parquet D$\Gamma$A instead of the PA. However, as seen in Fig. 4, the DMFT self-energy is clearly inaccurate for the PPP model, and a D$\Gamma$A calculation starting from a DMFT solution without self-consistency can therefore not be expected to yield reliable results.

A similar good agreement with the exact ED result is also obtained for the real part of the PA self-energy, shown in the lower panel of Fig. 4. In contrast, the real part of the self-energy is identically zero in DMFT because of particle-hole symmetry, i.e., it is qualitatively and quantitatively very different from the ED and PA result.

When we include non-local interactions and consider the PPP model, the self-energy in Fig. 4 is dramatically suppressed within ED. This aspect is completely missing in the DMFT, which yields exactly the same self-energy for both models, as long as there is no CDW order, since non-local interactions are only included at the Hartree level. 62. Instead, the PA provides an excellent description of the PPP model, and it is able to reproduce not only the overall suppression of the self-energy, but also the change of sign of ${\rm Re}\left[\Sigma_{\mathbf{k}=0}-\Sigma^{HF}_{\mathbf{k}=0}\right]$. One can rationalize the smaller self-energy for the PPP model, by considering the extreme case that the non-local interaction $V_{ij}$ is independent of the distance $|i-j|$. For $V_{ij}=V_{00}=U$ and a fixed number of electrons, all configurations in the occupation space compatible with the Pauli exclusion principle have the same interaction energy. Hence the system behaves more like the non-interacting system. In particular there is neither a renormalization (narrowing) effect (corresponding to the large linear part of ${\rm Im}\Sigma$ for small frequencies) nor a different ${\rm Re}\Sigma$ for each $\mathbf{k}$ points (however occupied/unoccupied states are shifted against each other). In the PPP model we do not yet have this extreme case, but the self-energy is already very much suppressed compared to the $U$-only model with the same $U$. As the PPP model is less strongly correlated, it is maybe not surprising that the PA provides even better results than for the $U$-only model. This suggests that the PA is a good approximation to describe the conjugated $\pi$-bonds in carbon-based molecules where non-local interactions partially compensate the effect of the local interaction.

A separation of the self-energy in a local dynamical part and a non-local static part according to $\Sigma_{(\mathbf{k},\nu)}=\Sigma^{loc}_{\nu}+\Sigma^{{}^{\prime}}_{\mathbf{k}}$ is discussed in Ref. 63. Similarly as in Ref. 64, we observe no clear separation of this kind. The Fock contribution can be assigned straightforwardly to the non-local static part $\Sigma^{{}^{\prime}}_{\mathbf{k}}$. However, comparing the $\mathbf{k}$ and $\nu$ dependence of the self-energy in Fig. 4, one can notice that the separation works better for the $U$-only model. Certainly, this is not the case for the PPP model, as the self-energy at different $\mathbf{k}$-points display a different low-frequency behavior (note the sign change of the dynamical contribution in the real part when comparing the first to the second column of Fig. 4)

The insets of Fig. 4 show the convergence of our results with respect to the size of the frequency box used, now including the Fock term. We see that at $N_{f}=320$ (total number of frequencies: positive and negative) deviations from the extrapolated $N_{f}\rightarrow\infty$ result are very minor. Hence $N_{f}=320$ has been employed for the main panel. Also shown in the inset is the ED result59 (as a triangle at $N_{f}=\infty$), and we see that PA approaches the ED results for large frequency boxes $N_{f}$ and high frequencies ($\nu_{20}$), whereas for smaller frequencies ($\nu_{5}$) there are deviations, as is also visible in the main panel. The right panel of Fig. 4 shows that PA not only correctly reproduces the frequency independent Fock (and Hartree) term but also captures the asymptotic $1/\nu$ behavior of the imaginary part of the self-energy (which is $\mathbf{k}$-independent) and the $1/\nu^{2}$ behavior of the real part of the self-energy.

Reproducing the high-frequency behavior correctly is an important non-trivial aspect of the PA and parquet variants of the D$\Gamma$A. In fact, within the ladder D$\Gamma$A, the correct $1/\nu$ asymptotic is only recovered if a Moriya $\lambda$-correction is included. 15; 12 In dual fermion calculations, 16 the correct asymptotics is obtained if the theory is truncated at the level of the two-particle vertex, but not any longer if three-particle vertex terms are included. 65; 66

Fig. 5 shows the ${\mathbf{k}}$-integrated and ${\mathbf{k}}$-resolved spectral function $A_{\mathbf{k}}(\nu)$, which for the PA and DMFT are obtained by an analytic continuation using the Padé method 333For the Padé interpolation, about $15$ grid points equally spaced with respect to Matsubara frequency are utilized and chosen such that the high-frequency behavior is not included. The correct tail is anyway automatically obtained by the assumed continued fraction.. The discrete peaks of the ED have been broadened by a Lorentzian of width $\delta=0.05t$. Again we show results for the PPP (right) and $U$-only model (left).

For the top panel of Fig. 5, i.e., the ${\mathbf{k}}$-integrated spectral function, we see that the (HOMO-LUMO) gap is considerably larger for the PPP model than for the $U$-only model. To get a better understanding of this difference let us start by looking at the self-energy, which within the band gap is accurately described by its leading order terms

[TABLE]

The static term $\text{Re}\Sigma_{k}(0)$, which here also includes the chemical potential, has opposite sign for the occupied part of the spectrum (${\mathbf{k}}=0,\pi/3$) and the unoccupied one (${\mathbf{k}}=2\pi/3,\pi$). It is added to the non-interacting energies $\epsilon_{\mathbf{k}}$, and contributes therefore directly to the band gap. The linear term, given by the slope of ${\rm Im}\Sigma$ in Fig. 4, contributes instead indirectly to the band gap through the quasiparticle renormalization factor

[TABLE]

The two terms together yield the following positions for the poles $(\epsilon^{*}_{\mathbf{k}})$ of the one-particle Green’s function

[TABLE]

The effect of $Z$ is to make the whole spectrum and hence also the band gap more narrow. Numerically, we determine the slope ${\partial{\rm Im}\Sigma_{k}}/{\partial\nu}$ for a finite $T$ by fitting a second order polynomial in $\nu$ to ${\rm Im}\Sigma_{\mathbf{k},\nu_{n}}$. The quasiparticle excitation energies Eq. (12) nicely fit to the peak positions of the PA and ED data, as shown in the blue shaded fields of Tab. 1.

Now, to understand why the gap increases for the PPP model compared to both, the $U$-only model and the non-interacting spectrum, we can analyze $Z$ and $\text{Re}\Sigma_{k}(0)$ in Tab. 1. We find $Z\approx 0.8$ for the $U$-only model and $Z\approx 0.95$ in the PPP model. The narrowing effect of $Z$ in the $U$-only model is almost completely canceled by the contribution from $\text{Re}\Sigma_{k}(0)$, which brings the band gap very close to its non-interacting value. In the PPP model the effect of $Z$ is instead almost negligible, and ${\rm Re}\Sigma_{k}(0)$ is dominated by a large Fock term, proportional to the non-local interaction strength, which strongly enhances the band gap. Since ${\rm Re}\Sigma_{k}\approx 1.3$ is quite similar for ${\mathbf{k}}=0$ and ${\mathbf{k}}=\pi/3$, the effect is like cutting the non-interacting spectrum at the Fermi level with a pair of scissors and shifting the occupied (unoccupied) down (up). This scissor operator is also a major effect of the $GW$ approximation for semiconductor band gaps.68; 69

We can conclude that, despite the sometimes problematic analytic continuation procedure, the PA provides a very good description of the ED spectrum, with the spectral features being somewhat broadened. In particular, the PA predicts very accurately the spectral gap. In contrast, DMFT slightly underestimates the spectral gap for the $U$-only model, which is not surprising given the previous analysis of the self-energy.

However, the DMFT results are the same for the $U$-only and PPP model (not plotted again in the right panel of Fig. 4), and thus completely miss the enhancement of the spectral gap due to the $\mathbf{k}$-dependence of ${\rm Re}\Sigma_{k}$ in case of the PPP model.

The good description of the PA is even more explicit when we inspect the $\mathbf{k}$-dependent spectral function in the lower of the two plots of Fig. 5. The insets therein show a zoom of the respective main panels, and hence resolve all the spectral features. 444Note that the broadening of the PA spectral function in the top panel is larger than for each $\mathbf{k}$-resolved component, since the Padé fit was performed after the $\mathbf{k}$-summation. We observe that the binding energy of all the spectral features of the ED calculation are excellently reproduced by the $\mathbf{k}$-resolved components within the PA. 555The weight of the spectral features depends slightly on the details of the Padé fit procedure. Please note that, the spectral function at $\mathbf{k}=\pi$ and $\mathbf{k}=\frac{2\pi}{3}$ (above the Fermi energy) can be obtained by mirroring (at the Fermi energy the spectrum of $\mathbf{k}=0$ and at $\mathbf{k}=\frac{\pi}{3}$. Overall, we can record that the PA provides an excellent approximation of the one-particle spectrum.

IV.4 Optical conductivity

We study the linear response of the benzene molecule to an electric field. The response is related to the perturbation by the relation $J(\omega)=\sigma(\omega)E(\omega)$, where the optical conductivity $\sigma(\omega)$ corresponds to the regular part of the current-current correlation function

[TABLE]

Here, the current-current correlation function \raisebox{2.0pt}{\chi}_{jj,q} is defined as

[TABLE]

where

[TABLE]

The current-current correlation function\raisebox{2.0pt}{\chi}_{jj,q} is evaluated in terms of the Green’s function and the full vertex as

[TABLE]

where we have separated the contribution from the bare bubble \raisebox{2.0pt}{\chi}^{\text{bub}}_{jj,q} and that from the vertex corrections \raisebox{2.0pt}{\chi}^{\text{ver}}_{jj,q}. The matrix elements $\gamma_{\mathbf{k}}^{\mathbf{q}}$ (not to be confused with the two-particle vertices of the parquet equation) depend on the applied perturbation. That is, it depends on the direction in which the electric field is applied and in which direction the current response is measured.

In the following, we consider the optical response for two cases: (i) The external field is a magnetic field applied in the direction perpendicular to the benzene ring. It generates a magnetic flux, and is equivalent to a circular electric field applied along the ring. In this case, the coupling with the electrons at small $\mathbf{q}$ is approximately given by the Peierls contribution to the (current) vertex $\gamma_{\mathbf{k}}=\gamma_{\mathbf{k}}^{\mathbf{q}=0}=2t\sin(\mathbf{k})\left(\equiv\partial\epsilon_{\mathbf{k}}/\partial\mathbf{k}\right)$. The respective correlation function is labeled by $P$ (for Peierls), namely \raisebox{2.0pt}{\chi}_{P}\equiv\raisebox{2.0pt}{\chi}_{jj} and likewise the optical conductivity via $\sigma_{P}$. The second case (ii) is the constant $\gamma_{\mathbf{k}}$ case, for which we use the label $C$, i.e. \raisebox{2.0pt}{\chi}_{C} and $\sigma_{C}$. This case corresponds to the dynamic compressibility or an electric field perpendicular to the ring in the limit of infinitesimally small perpendicular hopping.666 Please note that applying an electric field within the plane of the benzene ring results in a much more complicated equation than Eq. 16 since it breaks translational invariance along the ring and momentum is not conserved any longer.

Within victory, Eq. (16) can be easily evaluated. Let us stress that, in order to obtain reasonable results, especially for \raisebox{2.0pt}{\chi}^{\text{bub}}_{jj}, the kernel approximations are extensively exploited. For obtaining the optical conductivity $\sigma(\omega)$ we evaluate the Matsubara current-current correlation function and then perform a numerical analytic continuation \raisebox{2.0pt}{\chi}_{jj}(\omega_{n})\rightarrow\raisebox{2.0pt}{\chi}_{jj}(\omega+i\delta) via the Padé method. The analytic continuation is a delicate procedure, and in order to obtain physical result, only a subset of the bosonic Matsubara frequencies has been used, as shown in the insets of Fig. 6. As the temperature considered ($\beta=10/t$) is low enough, a sparse frequency grid can be used for the interpolation without qualitative changes in the outcome. Note, in fact, that a larger number of input points do not result in actual poles, as they are canceled by the numerator in the continued fraction. 777We further checked that the peaks calculated by the Padé analytical continuation reproduce our data when transformed back to imaginary frequency where each peak gives a Lorentzian with different weight (weight of the peak in real frequencies) and width (energy of the peak).

Let us discuss first the case (ii), i.e., the optical response to a perpendicular electric field, with $\gamma_{\mathbf{k}}\propto const.$ (labeled by $C$). The results for both $\sigma_{C}(\omega)$ and \raisebox{2.0pt}{\chi}_{C,\mathbf{q}=0}(i\omega_{n}) are shown in the main panel and in the inset of Fig. 6, respectively. For \raisebox{2.0pt}{\chi}_{C,\mathbf{q}=0}(i\omega_{n}), we show the full dataset as well as the actual subset used for the analytic continuation. The optical conductivity obtained by the analytic continuation displays a single peak. We analyze the evolution of the position of the peak at frequency $\omega^{\textrm{peak}}$ for the different interaction models and we observe two different effects, as discussed below. Taking into account just the bare bubble contribution \raisebox{2.0pt}{\chi}^{\text{bub}}_{C}, the peak if found at $\omega^{\textrm{peak}}\approx 4.9$ for the $U$-only model. It shifts to a higher frequency $\approx 7.6$ ($U+V_{1}$) and to $\approx 7.4$ (PPP) when increasing the range of the Coulomb interaction. This is in agreement with the larger gap in the spectral function Fig. 5.

The second effect is a reduction of $\omega^{\textrm{peak}}$ due to vertex corrections. While this effect is relatively weak in the $U$-only model and still within the uncertainties of the Padé fit ($\omega^{\textrm{peak}}\approx 4.8$ $\rightarrow$ $4.9$), it becomes significantly stronger both for the $U+V_{1}$ ($7.6\rightarrow 6.7$) and the PPP models ($7.4\rightarrow 5.2$). This behavior is consistent with the formation of excitons due to non-local interactions, which cause a shift of weight in the optical conductivity to lower frequencies. The fact that this effect is more pronounced in the PPP model than in the $U+V_{1}$ one is unexpected. It shows that, although the one-particle self-energy is suppressed by the non-local interactions in the PPP model, the system does not necessarily resemble an effective non-interacting limit. Electronic correlations induced by non-local interactions have a much stronger effect on the two-particle response than on the one-particle properties.

Unfortunately, a quantitative comparison of the results obtained for the PPP model to actual experimental data is difficult. The electron-impact spectra for benzene74; 75; 76 observe a rich structure of resonances. The low-energy excitations detected in the experiments in the range of $4.5-5.5$ eV and $6.2-6.5$ eV have been discussed to occur due to vibrational distortions of the benzene ring, 76 which are not included in our model. Our results for $\sigma_{C}(\omega)$ are instead compatible with the broad experimental spectrum observed at energies above $7.5$ eV. However, since the analytic continuation yields a single broad peak with a relatively large width, we do not have the resolution to identify any substructure.

For the bubble contribution, we can also compare the frequency structure of $\sigma_{C}(\omega)$ to the results obtained from the one-particle Green’s function for the local and the $\mathbf{k}$-resolved spectral function $A_{\mathbf{k}}(\omega)$ in Fig. 5. Optical transitions must occur between occupied and unoccupied parts of $A_{\mathbf{k}}(\omega)$ (i.e., between contributions below and above the Fermi level which is put to zero in Fig. 5). Further the optical transitions must occur between the same momentum $\mathbf{k}$ since the transferred momentum is $\mathbf{q}=0$. The lowest-energy for transitions within the $U$-only and PPP models occur both for $\mathbf{k}=\frac{\pi}{3}$ (and symmetrically for $\mathbf{k}=\frac{2\pi}{3}$), at a transferred energy $\Delta\omega=3.9$ and $\Delta\omega=6.9$, respectively. These values are off by 7% and 20% for the PPP model and $U$-only model, respectively. But given the ambiguity of the Padé method the results are actually quite comparable. Hence, we obtain a consistent picture between the bubble contribution to $\sigma_{C}(\omega)$ and the spectral function.

Finally, we can perform a similar analysis as the one above for the optical response to a perpendicular magnetic field, i.e., case (i) with $\gamma_{\mathbf{k}}=\partial\epsilon_{\mathbf{k}}/\partial\mathbf{k}$. The results are shown in Fig. 7. The peak positions of the bubble contribution to the conductivity nicely reproduce the optical gaps of the spectral function; and thus are not affected by the type of field either magnetic or electric. Also, the effect of vertex corrections on the position of the peak is qualitatively similar to the case of perpendicular electric field, for all models. The main difference is that the response is overall larger [taking a dipole element $\gamma_{\mathbf{k}}=1$ of the same magnitude in case (ii)]—even more so if vertex corrections are taken into account, and that the width of the main peak is substantially reduced.

V Discussion and Conclusion

We have studied the PPP model for the benzene molecule within the PA. A thorough investigation allowed us to gain insight of local and non-local effects on the physics of the system. The results obtained for the PPP model with long-range Coulomb repulsion are compared with those of the $U+V_{1}$ and $U$-only models, where the interaction range is restricted up to nearest-neighbors and on-site only, respectively. We solved these models within the PA, which account for non-local dynamical electronic correlations at all length-scales, as well as within the single-site dynamical mean-field theory, which only retains local quantum fluctuations. To strengthen our analysis, we further compared our results against the exact numerical solution, obtained by the ED of the Hamiltonian. Eventually, this allows us to better understand the effects of non-local interactions and non-local correlations within a unified formalism.

In particular, we find that the vertex functions are strongly enhanced in the PPP model compared to the $U$-only model, simply because the non-local interaction $V$ itself contributes to the vertex. Except for this effect, both (reducible) vertices are quite similar with the vertex for the PPP model being somewhat more smeared out over a larger frequency range than that of the $U$-only model. As for the leading instabilities, we find that antiferromagnetic fluctuations are dominant for the $U$-only (Hubbard) model, as expected. Whereas, charge fluctuations are so strongly enhanced for the $U+V_{1}$ model that this model is actually on the verge of a transition to a charge density wave. For the PPP model these charge fluctuations are again suppressed by second and third nearest neighbors repulsion, so that antiferromagnetic and charge density wave fluctuations turn out to be on par.

Despite the larger vertex, we find that the imaginary part of the self-energy (and with this the quasiparticle renormalization) of the PPP model is actually much smaller than that of the $U$-only model. In this respect electronic correlations are reduced. This can be understood as follows. If we consider two electrons, they are strongly correlated in the case of the Hubbard ($U$-only) model where finding both electrons on the same lattice site is largely suppressed. The non-local interactions, on the other hand, balance the local repulsion and therefore mitigate differences where the two electrons are located, hence reducing electronic correlations. However, the non-local interactions also lead to a large static contribution to the $\mathbf{k}$-dependent self-energy, which in a scissors-like way enhances the splitting of the occupied and unoccupied states. For this reason the long range interactions within the PPP model yields a considerably larger one-particle gap compared to the one of the Hubbard. For the optical conductivity we find that, for the PPP model, vertex corrections result in a much smaller optical gap than the gap in the (one-particle) spectral function. This is not the case for the $U$-only model, suggesting an important effect of non-local electronic correlations in the PPP model.

Since the self-energies of the PPP and $U$-only model are so different, DMFT is certainly not a good approximation for the PPP model (or for the conjugated $\pi$-bonds in carbon based molecules) as it yields the same self-energy (with and without non-local interactions) as long as there is no symmetry breaking. This is mitigated to some extent when using density functional theory (DFT) as a staring point in the so-called DFT+DMFT approach,77; 78; 79; 80; 81 but for getting the scissors-like self-energy we observe even a more involved $GW$+DMFT 82; 83; 84 is needed. This also implies that doing one-shot (ab-initio) D$\Gamma$A 21; 85 calculations with the vertex taken from a DMFT impurity problem is not justified; DMFT is too far away from the correct solution. A self-consistent D$\Gamma$A calculation of the vertex is necessary.

Besides self-consistent D$\Gamma$A, of course also QMC, ED (CI), CC, and (for the special geometry considered) DMRG are viable alternatives for molecules. Our work shows that this list should be supplemented by the PA which has a similar scaling as the CCSD, i.e., the numerical effort scales $\sim(N^{2})^{2.4}$ ($N$: number of states/orbitals; the factor of 2.4 is for an efficient matrix inversion). Employing symmetries and vanishing orbital off-diagonal elements, the effort can often be reduced, e.g., to $\sim N(N)^{2.5}$ in case of the highly symmetric molecules. For comparison the scaling of CCSD is $\sim N_{u}^{2}N^{4}_{o}$ ($N_{u/o}$: number of unoccupied states/occupied states which are both in general $\sim N$)23; and the ED even scales exponentially with $N$. The biggest bottleneck of the PA is the memory required to store the vertex functions, which scales with both the number of Matsubara frequencies and states (or $\mathbf{k}$-points in case of corresponding symmetries).

The advantages of the PA are that all one- and two-particle quantities are naturally calculated without any significant further efforts, including one-particle spectra measured experimentally in photoemission spectroscopy (PES), two-particle spectra measured e.g., through the optical conductivity, and response function as e.g., the magnetic susceptibility. For dynamical quantities (involving excited states) an analytic continuation is necessary which has some broadening effect. This broadening is more relevant for small molecules than for large molecules, where the actual spectrum is already “smeared out” due to the large number of energy levels. Last but not least, an advantage of the PA is that a physical understanding is fostered: we have, among others, the self-energy and vertex corrections from various channels at hand, not only the total energy and the (in general way to complicated) wave function.

Acknowledgment

We would like to thank Felix Hörbinger for helping to set up the ED calculations and Tin Ribic for helpful discussions on the optical conductivity. PP, AK and KH were supported by European Research Council under the European Union’s Seventh Framework Program (FP/2007-2013) through ERC Grant No. 306447, and by the Austrian Science Fund (FWF) through SFB ViCoM F41 and project P 30997. PP has been further supported by the FWF Doctoral School W 1243 “Building Solids for Function”; and AV acknowledges financial support from the Austrian Science Fund (FWF) through the Erwin Schrödinger fellowship J3890-N36, and through project ’LinReTraCe’ P 30213. Calculations have been done in part on the Vienna Scientific Cluster (VSC).