Dual Boson approach with instantaneous interaction

TL;DR

This paper introduces a simplified version of the Dual Boson approach by enforcing an instantaneous interaction, which preserves key physical laws and accurately captures certain phase transitions in strongly correlated systems.

Contribution

It develops a static interaction variant of the Dual Boson method, maintaining charge conservation and physical consistency, suitable for multiband systems.

Findings

Enforces absence of phase transitions in finite systems

Respects the Mermin-Wagner theorem in 2D systems

Captures metal to charge ordered phase transition accurately

Abstract

The Dual Boson approach to strongly correlated systems generally involves a dynamic (frequency-dependent) interaction in the auxiliary impurity model. In this work, we explore the consequences of forcing this interaction to be instantaneous (frequency-independent) via the use of a self-consistency condition on the instantaneous susceptibility. The result is a substantial simplification of the impurity model, especially with an eye on realistic multiband implementations, while keeping desireable properties of the Dual Boson approach, such as the charge conservation law, intact. We show and illustrate numerically that this condition enforces the absence of phase transitions in finite systems, as should be expected from general physical considerations, and respects the Mermin-Wagner theorem. In particular, the theory does not allow the metal to insulator phase transition associated with the formation of the magnetic order in the two-dimensional system. At the same time, the metal to charge ordered phase transition is allowed, as it is not associated with the spontaneous breaking of a continuous symmetry, and is accurately captured by the introduced approach.

Equations

$$H=$$

$$+\frac{1}{2}\sum _{i,j}{V}_{ij}\left(n_i-⟨n_i⟩\right)\left(n_j-⟨n_j⟩\right).$$

$$\mathcal{S}=$$

$$+\sum _{\mathbf{k},\nu ,\sigma }{t}_{\mathbf{k}}^{}{c}_{\mathbf{k}\nu \sigma }^{\ast }{c}_{\mathbf{k}\nu \sigma }^{}+\frac{1}{2}\sum _{\mathbf{q},\omega }{V}_{\mathbf{q}}{n}_{\mathbf{q}\omega }{n}_{-\mathbf{q}-\omega }$$

$$=$$

$${\mathcal{S}}_{\text{imp}}=$$

$$+U\sum _{\omega }{n}_{\omega \uparrow }{n}_{-\omega \downarrow }+\frac{1}{2}\sum _{\rho ,\omega }{\Lambda }_{\omega }^{\rho }{\rho }_{\omega }{\rho }_{-\omega }.$$

$${\mathcal{S}}_{\text{imp}}=$$

$$X_{\mathbf{q}\omega }^{-1}=$$

$${\mathcal{X}}_{\mathbf{q}\omega }=$$

$$\sum _{\omega }{\chi }_{\omega }^{\text{ch/sz}}=\sum _{\omega }{X}_{\text{loc, }\omega }^{\text{ch/sz}},$$

$${\chi }^{\text{ch/sz}}(\tau =0)=X_{\text{loc}}^{\text{ch/sz}}(\tau =0).$$

$${\Lambda }_{\text{new}}={\Lambda }_{\text{old}}+\xi \frac{X_{\text{loc}}(\tau =0)-\chi (\tau =0)}{X_{\text{loc}}(\tau =0)\chi (\tau =0)}.$$

$$\sum _{\omega }{\chi }_{\omega }^{\text{sz}}=\frac{1}{N}\sum _{\omega }\sum _{\mathbf{q}}{X}_{\mathbf{q}\omega }^{\text{sz}}.$$

$$\int d{\mathbf{q}}^{\prime }|q^{\prime }{|}^{-2}\propto {\int }_{0}d|q^{\prime }||q^{\prime }{|}^{d-1}|q^{\prime }{|}^{-2}={\int }_{0}d{q}^{\prime }|q^{\prime }{|}^{d-3}$$

$${\chi }_{\omega }^{-1}=\chi (U=0{)}_{\omega }^{-1}\mp U^{\prime },$$

$${\mathcal{X}}_{\mathbf{q}\omega }=\mathcal{X}(U=0{)}_{\mathbf{q}\omega }^{-1}\mp U^{\prime },$$

$$X_{\mathbf{q}\omega }^{-1}={\mathcal{X}}_{\mathbf{q}\omega }^{-1}+\Lambda ,$$

$${\left(X_{\mathbf{q}\omega }^{\text{ch}}\right)}^{-1}=$$

$$=$$

$${\left(X_{\mathbf{q}\omega }^{\text{sz}}\right)}^{-1}=$$

$$\sum _{\omega }\frac{\chi (U=0{)}_{\omega }}{1-(U+{\Lambda }^{\text{ch}}-{\Lambda }^{\text{sz}})\chi (U=0{)}_{\omega }}=\frac{1}{N}\sum _{\mathbf{q},\omega }\frac{X(U=0{)}_{\mathbf{q}\omega }}{1-(U-{\Lambda }^{\text{sz}})X(U=0{)}_{\mathbf{q}\omega }},$$

$$\sum _{\omega }\frac{\chi (U=0{)}_{\omega }}{1+(U+{\Lambda }^{\text{ch}}-{\Lambda }^{\text{sz}})\chi (U=0{)}_{\omega }}=\frac{1}{N}\sum _{\mathbf{q},\omega }\frac{X(U=0{)}_{\mathbf{q}\omega }}{1+(U+{\Lambda }^{\text{ch}})X(U=0{)}_{\mathbf{q}\omega }}.$$

$$U^{\prime }\sum _{\omega }{\chi }^2(U=0{)}_{\omega }=\frac{U^{\prime }\pm {\Lambda }^{\text{ch}/\text{sz}}}{N}\sum _{\mathbf{q}\omega }{X}^2(U=0{)}_{\mathbf{q},\omega },$$

$${\Lambda }^{\text{ch}/\text{sz}}=\pm U^{\prime }\cdot C=\pm U\frac{1}{2-C^{-1}},$$

$$C=\frac{{\sum }_{\omega }{\chi }^2(U=0{)}_{\omega }}{\frac{1}{N}{\sum }_{\mathbf{q}\omega }{X}^2(U=0{)}_{\mathbf{q},\omega }}-1.$$

$${\Sigma }_{\nu }^{\text{lat}}\stackrel{\nu \to \infty }{=}\frac{U⟨n⟩}{2}+\frac{U^2{⟨n⟩}^2}{4}\frac{1}{i\nu }+\dots .$$

$${\Sigma }_{\nu }^{\text{imp}}\stackrel{\nu \to \infty }{=}\frac{U^{\prime }⟨n⟩}{2}+\frac{U^{\prime 2}{⟨n⟩}^2}{4}\frac{1}{i\nu }+\dots .$$

$${\tilde{\Sigma }}_{\nu }^{\text{exact}}\stackrel{\nu \to \infty }{=}\frac{(U-U^{\prime })⟨n⟩}{2}+\frac{(U^2-U^{\prime 2}){⟨n⟩}^2}{4}\frac{1}{i\nu }+\dots .$$

$${\tilde{X}}_{\mathbf{q}\omega }^{(0)}$$

$${\tilde{X}}_{\mathbf{q}\omega }^{(0)}/{\chi }_{\omega }^2$$

$$\approx \Lambda =U^{\prime }-U.$$

$${\gamma }_{{\nu }_1{\nu }_2{\omega }_1{\omega }_2}^{(2,2)}\times {\chi }_{{\omega }_1}{\chi }_{{\omega }_2}\sim g_{{\nu }_3}{\delta }_{{\omega }_1{\omega }_2}{\delta }_{{\nu }_1{\nu }_2}{\delta }_{{\nu }_3-{\nu }_1+{\omega }_1}.$$

$${\tilde{\Sigma }}_{\nu }=$$

$$\approx$$

Full text

I Introduction

The Hubbard model Hubbard (1963); Kanamori (1963); Gutzwiller (1963); Hubbard (1964) describes interacting itinerant electrons. The Coulomb interaction between these electrons leads to correlations and makes this a many-body problem. When interaction and itinerancy are almost equally strong, the resulting electronic correlations make the system notoriously hard to study.

The Dynamical Mean-Field Theory (DMFT) Georges et al. (1996) is an important tool to describe correlated electrons. The main idea of the method is to use an auxiliary “impurity” model that consists of only a single site and contains the most important correlations of the full system. This approach is exact in the limit of infinite dimension Metzner and Vollhardt (1989), in finite dimensional systems it still serves as a useful approximation.

Extensions of DMFT have been developed Rohringer et al. (2018) to add the spatial correlations and non-local interactions that DMFT ignores. These extensions add additional correlation effects on top of DMFT, but they can also change the impurity model that is used as a starting point. To account for screening by the non-local Coulomb interaction, Extended DMFT (EDMFT) Sengupta and Georges (1995); Si and Smith (1996); Kajueter (1996); Smith and Si (2000); Chitra and Kotliar (2000, 2001) and its diagrammatic extensions Rohringer et al. (2018), introduce dynamic interactions into the impurity model via an effective frequency dependent bosonic hybridization function. While the EDMFT still considers electronic correlations at the level of the local impurity problem, the Dual Boson (DB) theory Rubtsov et al. (2012); Stepanov et al. (2016) and other EDMFT-based approaches take additional nonlocal correlation effects into account diagrammatically.

The DB approach aims to treat fermionic and collective bosonic degrees of freedom on an equal footing. Nevertheless, until recently the DB calculations were mostly focused on the description of charge degrees of freedom, with the bosonic hybridization function introduced only in the charge channel. This means that the screening by collective spin fluctuations was missing. The introduction of dynamic interactions in the spin channel generally comes at a computational cost Otsuki (2013); Steiner et al. (2015). For single-band systems, a dynamic spin density ($S^{z}$-$S^{z}$) interaction can be incorporated into CT-QMC impurity solvers at moderate cost, in the same way as the dynamic charge interaction Werner and Millis (2007, 2010). However, this immediately breaks the spin-rotational invariance of the local impurity problem. On the other hand, the use of the same frequency dependent bosonic hybridization function for every ($x$, $y$, and $z$) spin channel to obey the rotational invariance immediately leads to a violation of conservation laws Krien et al. (2017).

To address these situations, a simplification of the Dual Boson approach that does not require dynamic interactions in the impurity model has been introduced recently Stepanov et al. (2018). We use the fact that the dynamic interactions in Dual Boson are a priori free parameters in the Hubbard-Stratonovich decoupling that leads to the impurity model. The dynamic interaction is usually determined using a set of self-consistency conditions, one for every frequency. By forcing the interaction to be instantaneous (independent of frequency), the number of free parameters is drastically reduced, with the benefit of having a much simpler impurity model that accounts for the screening by collective charge and spin fluctuations in a spin-rotational invariant form that does not violate local conservation laws. From the physical point of view, the constant form of the hybridization function in the spin channel can be motivated by the fact that collective spin fluctuations are slower and have lower energy than single-particle (electronic) excitations. Thus the interaction between spins and electrons is instantaneous ($\delta$-function in time), which leads to a constant hybridization function.

So far, this method has not been systematically tested. In this work we aim to systematically test its performance in the description of collective electronic instabilities of the two- and three-dimensional extended Hubbard model. We compare this instantaneous interaction Dual Boson approach with the traditional dynamic interaction Dual Boson approach and with the Dynamical Mean-Field which always keeps the interaction of the original Hubbard model. We also discuss the relationship with the two-particle self-consistent approach Vilk et al. (1994); Y.M. Vilk and A.-M.S. Tremblay (1997) and the Moriya correction in D$\Gamma$A Toschi et al. (2007), which are also based on self-consistently renormalizing an effective interaction.

II Model and Method

We consider the half-filled extended Hubbard model on the square and cubic lattice,

[TABLE]

Here $c^{\dagger}_{i\sigma}$ and $c^{\phantom{\dagger}}_{i\sigma}$ are creation and annihilation operators for an electron on site $i$ with spin $\sigma$ and $n_{i\sigma}=c^{\dagger}_{i\sigma}c^{\phantom{\dagger}}_{i\sigma}$ is the corresponding number operator. The total electronic density on site $i$ is equal to $n_{i}=n_{i\uparrow}+n_{i\downarrow}$. The physical parameters of the model are the hopping amplitude $t_{ij}$, the on-site Coulomb interaction $U$ and the intersite Coulomb interaction $V_{ij}$. We consider only nearest-neighbour hopping and interaction, $t_{ij}=t$ and $V_{ij}=V$ when $i$ and $j$ are nearest-neighbors and 0 otherwise.

For $V=0$, this model is simply the Hubbard model. In this work, we consider the (extended) Hubbard model at half-filling, $\left\langle n_{i}\right\rangle=1$. This is obtained by setting $\mu=U/2$. We restrict ourselves to phases without explicit symmetry breaking. We do calculate the susceptibility associated with antiferromagnetic (AF) and charge density wave (CDW) phases to check for instabilities towards ordered phases.

The idea of the Dual Fermion Rubtsov et al. (2008) and Dual Boson Rubtsov et al. (2012) methods is to decouple the extended Hubbard model into two parts. The first is an impurity part that only contains local degrees of freedom. This part can be solved numerically exactly. The impurity part should be chosen in such a way that it contains the most important correlation effects, since all correlations that are present in the impurity model are treated exactly. The second part is the remainder, which is a lattice model just like the original (extended) Hubbard model. The difference with the original model is that the original degrees of freedom are transformed to new dual degrees of freedom when the impurity model is integrated out. These dual degrees of freedom are less correlated, so that the correlations in this dual part can be addressed perturbatively.

Here, we provide an overview of the formulation of the Dual Boson method, for a full derivation and analysis of the method we refer the reader to the original works Rubtsov et al. (2012); van Loon et al. (2014a); Hafermann et al. (2014); Stepanov et al. (2016) and the review Rohringer et al. (2018). Mathematically, the action formalism is the most convenient way to perform the dual decoupling. We write the action in terms of the Fourier transforms of the hopping and interaction, $t_{\mathbf{k}}$ and $V_{\mathbf{q}}$, and in terms of Matsubara frequencies

[TABLE]

This separation is made by introducing fermionic and bosonic hybridization functions $\Delta$, $\Lambda^{\rho}$, where $\rho$ is a bosonic channel, e.g., charge or $S^{z}$.

[TABLE]

These hybridizations $\Delta$, $\Lambda$ can be chosen freely. Below in Section II.1 we discuss how to make this choice. In general, the Dual Boson approach considers $\Lambda$ to be a function of $\omega$. In this work, we restrict ourselves to constant $\Lambda$. In that case, the impurity action can simply be rewritten as

[TABLE]

where $U^{\prime}=U+\Lambda_{\text{ch}}-\Lambda_{\text{sz}}$. For the impurity model to also be at half-filling, we take $\mu^{\prime}=U^{\prime}/2$. This choice of hybridization functions is useful, since the resulting form of the impurity action does not introduce higher-order terms in the Ward identities Krien et al. (2017); Stepanov et al. (2018). As we discuss below in Section II.2, the impurity problem is solved numerically exactly and provides full frequency dependent local quantities such as the one- and two-particle Green’s functions and fermion-fermion and fermion-boson vertices needed for the construction of the dual diagrams.

The spatial correlations in the dual part of the problem are addressed using diagrammatic methods. In this work, we concentrate on the correlation effects in the susceptibility and we do not calculate dual self-energy corrections. The following definition for the susceptibility $X_{\mathbf{q}\omega}=-\left\langle\rho\rho\right\rangle_{\mathbf{q}\omega}$ can be applied both to the charge channel, with $\rho=n-\left\langle n\right\rangle$, and to the spin channel, with $\rho=S^{z}-\left\langle S^{z}\right\rangle$. The Dual Boson expression is given by

[TABLE]

where

[TABLE]

and $\chi_{\omega}=-\left\langle\rho\rho\right\rangle^{\text{imp}}_{\omega}$ is the corresponding susceptibility of the impurity model. The non-local correlation effects are contained in $\tilde{\Pi}$, the dual polarization operator. In this work, we use the ladder equation to calculate $\tilde{\Pi}$ Rubtsov et al. (2012); van Loon et al. (2014a); Hafermann et al. (2014); Stepanov et al. (2016), which describes repeated particle-hole scattering mediated by the fermion-fermion vertex of the impurity model. The ladder DB approach with the frequency dependent bosonic hybridization shows a quantitatively good result for the susceptibility in agreement with Quantum Monte Carlo and Dynamical Cluster Approximation calculations in a broad range of Coulomb interactions and dopings van Loon et al. (2016a); Terletska et al. (2017); Stepanov et al. (2019). The restriction to instantaneous interactions simply leads to the replacement of $\Lambda_{\omega}$ by a constant $\Lambda$ in Eq. (6).

Dynamical Mean-Field Theory usually deals with the Hubbard model, with $V=0$. The DMFT susceptibility is obtained by setting $\Lambda=0$ in Eq. (6). This formulation of the DMFT susceptibility is equivalent Hafermann et al. (2014) to more traditional ways of computing it Georges et al. (1996). The extended Hubbard model’s intersite Coulomb interactions can be included in the DMFT susceptibility via Eq. (6) with $\Lambda=0$, this essentially corresponds to an RPA-like treatment of these additional interactions, in which they do not alter the impurity or the single-particle properties.

II.1 Self-consistency condition

So far, we have not specified the interaction renormalization $\Lambda$, which is determined by a self-consistency condition. We use the “lattice” self-consistency condition proposed for the self-consistent Dual Boson approach Stepanov et al. (2016), except that here we have only two free parameters $\Lambda^{\text{ch}}$, $\Lambda^{\text{sz}}$ instead of two dynamic functions of frequency. This means we only need one self-consistency condition per channel instead of having a self-consistency for every Matsubara frequency.

[TABLE]

where $X_{\text{loc}}$ is the local part of the lattice susceptibility, obtained as the average over momenta of $X_{\mathbf{q}}$. As shown in Ref. Stepanov et al., 2018, this choice of the self-consistency condition follows from the invariance of the initial lattice problem with respect to the variation of the introduced hybridization functions, and it fulfills the Pauli principle.

The sum over frequencies in the self-consistency condition corresponds to taking the equal-time component of the susceptibility111The equal time component is obtained by dividing both sides by $\beta$ to normalize the frequency sum, and the sum runs over both positive and negative Matsubara frequencies., so we can write

[TABLE]

The difference between the local, equal time charge and spin susceptibility determines the double occupancy $D=\left\langle n_{\uparrow}n_{\downarrow}\right\rangle$. The self-consistency condition used here assures that the double occupancy of the lattice susceptibility is equal to the one of the impurity model van Loon et al. (2016b). On the other hand, the potential energies will be different since the lattice model and the impurity model have different interactions $U$ and $U^{\prime}$. This type of inconsistency between impurity and lattice generally occurs in approximate methods based on DMFT Krien et al. (2017). Note that the argument about potential energies applies to the Hubbard model ($V=0$), for the extended Hubbard model there is no direct correspondence between the potential energy of the impurity and the lattice since the latter also contains intersite contributions.

Dual Boson calculations with dynamic interactions in the $S^{z}$ channel break the spin rotational symmetry of the impurity model if they do not have the same dynamic interaction in the $S^{x}$ and $S^{y}$ channels. There is no such rotational symmetry breaking in the current approach, since the impurity model has the rotationally invariant interaction $U^{\prime}$. In the dual part of the calculations, it is sufficient to calculate only the $S^{z}$ channel since the $S^{x}$ and $S^{y}$ channel follow from rotational symmetry.

Appendix A describes analytical results for the self-consistency condition in several simplifying limits. These results carry over from the self-consistent Dual Boson approach with dynamic interactions Stepanov et al. (2016).

The reader might wonder why the self-consistency condition is important in the first place. The relation between the original and the dual theory is exact for any choice of $\Delta$ and $\Lambda$. However, the dual theory is only solved approximately, perturbatively. The choice of $\Lambda$ becomes important in this approximate theory. In this sense, there is some similarity with the Fierz ambiguity Baier et al. (2000); Jaeckel and Wetterich (2003); Ayral et al. (2017), where the decomposition of $U$ into channels determines the computational outcome of approximate methods even though the exact system is unchanged by the choice of decomposition. The present situation can also be seen through this lense, since the impurity part of the calculation depends only on the combination $U^{\prime}=U+\Lambda_{\text{ch}}-\Lambda_{\text{sz}}$, whereas the dual part of the calculation and the transformation between dual and lattice quantities are affected by how $U^{\prime}$ is decomposed into the charge and magnetic channel.

As a corollary, the choice of self-consistency condition should be informed by the approximations made in the dual theory, since only these approximations make the self-consistency condition relevant. Presently, the approximations are the elimination of vertices beyond the two-particle level, the restriction to ladder diagrams and finally the choice of the density and magnetic channels for the ladder. Starting with the last point, the choice of channels for the ladder and subsequently for the susceptibility directly informs the choice of $\Lambda_{\text{ch}}$ and $\Lambda_{\text{sz}}$ as the self-consistent parameters (see, e.g., Ref. Del Re et al., 2018 for a recent discussion on the role of channels in the attractive Hubbard model). Diagrammatic Monte Carlo approaches to the dual theory Iskakov et al. (2016); Gukelberger et al. (2017) do not restrict themselves to ladder diagrams and are therefore not restricted by the second and third point.

II.2 Implementation

Apart from the self-consistency condition and the effective interaction, our method follows the Dual Boson approach Rubtsov et al. (2012) and its implementation van Loon et al. (2014a). For the impurity model, we use a modified version of the open source CT-HYB solver Hafermann et al. (2013); Hafermann (2014) based on the ALPS libraries Bauer et al. (2011).

The implementation of the Dual Boson approach consists of two computationally heavy parts, the impurity solver and the evaluation of the dual diagrams, both implemented in C++ and linked together by a lighter python interface. The self-consistency condition is based entirely on the impurity susceptibility $\chi$ and the local part of the lattice susceptibility $X_{\text{loc}}$. These are part of the usual output of the impurity solver and the dual program. This means that the self-consistency condition can be implemented in the python interface, the renormalized effective interaction $U^{\prime}=U+\Lambda_{\text{ch}}-\Lambda_{\text{sz}}$ enters the impurity solver as a parameter and no other changes are needed.

We use an update formula with damping parameter $\xi$,

[TABLE]

A smaller value of $\xi$ makes the self-consistency procedure more stable at the cost of slower convergence, we take $\xi\in[0.1,1]$. Similar to the self-consistent procedure to determine $U(\omega)$ used in Ref. van Loon et al., 2016b, we find that convergence slows down substantially at larger interaction strengths.

II.3 Phase transitions and the Mermin-Wagner theorem

The goal of extensions of DMFT is to add the physical aspects of finite (low) dimensional physics to the $d=\infty$ solution. Dimensionality is especially important for the understanding of phase transitions to ordered phases. These show up as a divergence in the corresponding susceptibility, at zero frequency and at a specific momentum $q^{\ast}$ which describes the ordering pattern, e.g., $\mathbf{q}_{\text{AF, 2d}}=(\pi,\pi)$ or $\mathbf{q}_{\text{AF, 3d}}=(\pi,\pi,\pi)$ for antiferromagnetic ordering.

There is a natural link between phase transitions and the self-consistency condition when the susceptibility appears in the latter Y.M. Vilk and A.-M.S. Tremblay (1997). We have previously written the self-consistency condition as

[TABLE]

Here, the momentum sum signifies a Brillouin Zone integral in the case of an infinitely large system. A phase transition shows up as a divergence in $X_{\mathbf{q}^{\ast},\omega=0}$ on the right-hand side of Eq. (11). However, the left-hand side is the correlation function of a single site and always stays finite. In two dimensions or in the case of a finite lattice, this would lead to a contradiction 222The sign of $X$ is fixed, so there can be no cancellation of divergences on the right hand side., and we indeed know that there should be no phase transition in a finite system or in a two-dimensional system at finite temperature (Mermin-Wagner). For $d>2$, the integral

[TABLE]

remains finite. This means that a phase transition can occur while still satisfing Eq. (11), but only when taking the limit $N\rightarrow\infty$.

We should point out that many mean-field based methods produce phase transitions in finite systems. In fact, the system size does not even enter traditional Curie-Weiss mean-field theory. The DMFT susceptibility of a finite system can be divergent and non-local extensions like Dual Fermion or Dual Boson do not automatically correct this issue. Here, it is enforced by the self-consistency condition. This situation is in some sense reminiscent of the difference between finite-size and quantum cluster approaches Maier et al. (2005): finite-size simulations only obtain a phase transition after extrapolating to infinite system size, quantum cluster approaches already find the transition at a finite system size. The former situation is technically correct, the latter is convenient in practical situations.

We recognize that the self-consistency condition serves to enforce the Mermin-Wagner theorem Mermin and Wagner (1966) in a two-dimensional system, which forbids spontaneous symmetry breaking of continuous symmetries in two-dimensional systems at finite temperature. For the Hubbard model, the spin rotational symmetry cannot be spontaneously broken. In particular, this means that there can be no antiferromagnetic state at finite temperature. On the other hand, at $T=0$, the system is an antiferromagnetic insulator at any $U>0$ if only nearest-neighbor hopping $t_{ij}$ is considered. Together, this means that at low temperature, the system features very strong and long-ranged antiferromagnetic correlations that are almost truly long-ranged ordered. A similar situation takes place in low-dimensional Heisenberg model, see Refs. Irkhin et al., 1999; Chakravarty et al., 1989, and references therein.

This situation, with long but not infinitely ranged correlations, is challenging to reproduces in (computational) approximations. According to Vilk and Tremblay Y.M. Vilk and A.-M.S. Tremblay (1997), a sufficient criterion for ensuring that an approximation satisfies the Mermin-Wagner theorem is to verify that the double occupancy $D=\left\langle n_{\uparrow}n_{\downarrow}\right\rangle$ obtained from taking the local, equal-time part of the susceptibility stays within the physical range $[0,n^{2}/2]$. The self-consistency condition of the current method ensures that the double occupany of the Hubbard model is equal to that of a reference impurity model that is solved exactly. The latter stays within the physical bounds, so that the former does as well and the method satisfies the Mermin-Wagner theorem.

II.4 Exact properties and approximate solutions

Let us take a step further back to put these developments in the general context of consistency in approximate solutions to many-body problems and in particular to Hubbard-like models. The central point is that (many) exact statements can be made about the true, exact solution of the model, to name just a few: equivalence of thermodynamic quantities and response functions according to the Kubo formula; the Mermin-Wagner theorem; conserved quantities corresponding to symmetries and Goldstone modes arising when these symmetries are broken; sum rules and high-frequency asymptotics derived from commutation relations. 333It is worth observing that these relations occur on very different length and time scales: The Mermin-Wagner theorem and Goldstone modes are long wavelength, low frequency phenomena whereas the high-frequency asymptote of the self-energy is a local, high frequency phenomenon. These concepts are central to our understanding of condensed matter physics.

However, approximate solutions are not guaranteed to satisfy these exact properties. In fact, they will usually not satisfy all these constraints. This issue goes back to the seminal work of Baym and Kadanoff Baym and Kadanoff (1961); Baym (1962), who formulated functional constructions to ensure that diagrammatic approaches satisfy certain conservation laws. The correspondence between correlation and response functions played an important role in the development of the theories of the electron gas Geldart and Taylor (1970) and of magnetism in itinerant electron systems Moriya and Kawabata (1973); Dzyaloshinskii and Kondratenko (1976); Moriya (1985).

In the context of (extensions of) DMFT, the subject of exact properties came up for the thermodynamic consistency of zero- and one-particle properties Aichhorn et al. (2006); Potthoff (2006) and the consistency between one- and two-particle quantities, in particular the Kubo relation between correlation and response Fishman et al. (2005); Otsuki and Kuramoto (2013); van Loon et al. (2015) and the high-frequency asymptote of the susceptibility Krien et al. (2017). Simultaneously, charge conservation and Goldstone modes were investigated Hafermann et al. (2014); van Loon et al. (2014b); Stepanov et al. (2016); Krien et al. (2017); Geffroy et al. (2018).

The conclusion of these investigations is that the DMFT approximate solution to the finite-dimensional Hubbard model satisfies many of the exact properties, but not all. More specifically, the charge response is consistent van Loon et al. (2015), the DMFT susceptibility satisfies global charge conservation Hafermann et al. (2014); van Loon et al. (2014b) and the lowest order terms in the high-frequency asymptotes of the Green’s function, (local) self-energy and the (momentum-resolved) susceptibility are consistent with exact relations. The DMFT functional is conserving Krien et al. (2017) in the sense of Baym and Kadanoff. On the other hand, the DMFT susceptibility violates the Mermin-Wagner theorem 444See, e.g., the review Rohringer et al., 2018., the potential energy/double occupancy is inconsistent between the one- and two-particle levelvan Loon et al. (2016b); Rohringer and Toschi (2016); Krien et al. (2017) and exact relations for the moments of the lower Hubbard band spectral weigh function in the atomic limit are violated Esterling (2018a, b). As we have already seen, DMFT also predicts phase transitions in finite systems, which is inconsistent with fundamental thermodynamic considerations.

Clearly, even the elegant construction of DMFT is not sufficient to recover all exact relations and it is unlikely that any approximate method can. It is possible to enforce specific relations, although usually at a cost. The self-consistency condition employed here can be seen in this way: it enforces the Mermin-Wagner theorem. In this way, it is similar to the $\lambda$ introduced by Moriya Moriya and Kawabata (1973); Moriya (1985).

A similar Moriya-$\lambda$ correction plays a central role in ladder-D$\Gamma$A Toschi et al. (2007); Katanin et al. (2009). The $\lambda$ of D$\Gamma$A is partially related Rohringer et al. (2018) to the $\Lambda_{\omega}$ of DB. Both enter the (inverse) susceptibility as in Eq. 6, but in DB $\Lambda$ also enters into the impurity model, as follows from the exact dual transformation, whereas in D$\Gamma$A $\lambda$ correction is introduced by hand and is not included in the impurity model. We should point out that the instantaneous $\Lambda$ proposed here is even more similar to the instantaneous $\lambda$ of D$\Gamma$A Rohringer and Toschi (2016).

III Square lattice Hubbard model

We now turn to numerical investigation of the instantaneous DB. We study the square lattice Hubbard ($V=0$) model with $t=1$. We use a $32\times 32$ discretization of momentum space.

III.1 Effective interaction

We start by investigating the effective interactions $\Lambda$, since these are the quantities that enter the impurity model. Figure 1 shows $\Lambda$ in the small to moderate coupling regime. We observe that $\Lambda$ depends strongly on the inverse temperature $\beta=1/T$, that $\Lambda^{\text{ch}}$ and $\Lambda^{\text{sz}}$ have opposite sign and that $\Lambda$ is proportional to $U$ at small interaction strengths.

The limit of small $U$ can be understood in terms of perturbation theory. For the calculation of the susceptibility at small $U$, we can neglect the self-energy and treat the interaction in the RPA fashion, as a geometric series. Appendix A.2 gives the details of this approach, which confirm the opposite sign and the proportionality with $U$. For comparison, the bottom panel of Fig. 1 shows the renormalization of the interaction in the two-particle self-consistent approach (TPSC). These calculations have been performed using the TRIQS package Parcollet et al. (2015). We observe that the interaction renormalization in Dual Boson is substantially smaller than that in TPSC. The self-consistency in the present approach starts from the DMFT susceptibility, which already includes local self-energy insertions and dynamical vertex corrections. At high temperatures, it provides a good starting point so that almost no renormalization of the impurity interaction is needed ($\Lambda\approx 0$). As the temperature is lowered, the magnitude of $\Lambda$ increases. At the parameters studied here, the magnitude of $\Lambda$ is an order of magnitude smaller than $U^{\text{ch/sz}}-U$ in TPSC.

One of the main points of the self-consistent Dual Boson approach is that the self-consistency condition enforces the Mermin-Wagner theorem in two dimensions. The absence of antiferromagnetism in TPSC and the self-consistent Dual Boson approaches is visualized in Fig. 2. DMFT is unstable towards antiferromagnetism after $U=6$ when the inverse of the magnetic susceptibility changes sign. The magnetic susceptibility of RPA diverges already at $U\simeq 2.5$. The self-consistency condition enforces a positive value for the DB results and thus pushes the inverse susceptibility away from zero. However, enforcing this condition is not easy, since close to an instability the system is very sensitive to small changes in the effective impurity interaction. This makes calculations at $U\geq 6$ extremely unstable.

We proceed with an in-depth look at $U/t=5$, $\beta t=3$, this is just before the antiferromagnetic susceptibility in DMFT diverges. We find $U^{\prime}/t\approx 5.54$, a more than 10% increase of the effective interaction. This change in interaction is composed of $\Lambda^{\text{ch}}\approx 0.41$, $\Lambda^{\text{sz}}\approx-0.13$.

Looking at the local magnetic susceptibility shown in Fig. 3, we see a large inconsistency in DMFT between the impurity and local lattice susceptibility. In the self-consistent approach, this inconsistency is removed almost completely, even though the self-consistency only enforces equality between the frequency-averaged susceptibilities. Essentially, the problem in DMFT is that $\operatorname{abs}(X^{\text{sz}})$ is too large, so that the double occupancy from X can even turn negative. The self-consistency condition solves this by reducing $\operatorname{abs}(X^{\text{sz}})$. This reduction in the lattice quantity occurs even though $\operatorname{abs}(\chi^{\text{sz}})$ is increased due to the larger effective impurity interaction. The larger effective impurity interaction also increases the self-energy, shown in Fig. 4 (see also Appendix B).

IV Cubic lattice Hubbard model

We now move to a three-dimensional system and consider $U/t=4$ and $\beta t=2.5$. We find $U^{\prime}/t\approx 4.18$, a roughly 5% change in the effective interaction. The channel decomposition of the effective interaction is $\Lambda^{\text{ch}}\approx+0.12$ $\Lambda^{\text{sz}}\approx-0.06$. The difference in magnitude between the charge and spin renormalization shows that we have clearly left the weakly correlated regime.

The calculations shown here have been performed on a $10\times 10\times 10$ cubic lattice. We have verified that using a $20\times 20\times 20$ lattice leads to very similar results. This conforms to the observation that replacing the integral in Eq. (12) by a finite momentum average is a well-behaved operation. All other potential sources of finite size effects are similar to DMFT and not relevant in this parameter regime.

The self-consistency condition enforces equality between the frequency averages of $X^{\text{loc}}$ and $\chi$. In Fig. 5, we show the difference of these two susceptibilities as a function of frequency 555Note that only positive frequencies are shown, so that all finite frequencies have a negative frequency counterpart that is equal by symmetry. There is an essential difference between finite and zero frequency, in both the charge and the spin channel. The contribution at zero frequency is compensated by the finite frequencies. The sign difference between the charge and spin channel in Fig. 5 corresponds to the sign difference between $\Lambda^{\text{ch}}$ and $\Lambda^{\text{sz}}$.

The self-consistency condition on the two-particle level feeds back to the single-particle level via the impurity model. In Fig. 6 we show the self-energy of the impurity model. We see an enhancement of the self-energy in the self-consistent DB approach. This enhancement originates in the larger effective interaction $U^{\prime}>U$.

Moving away from fixed $U$, Figure 7 shows how the effective interaction and the potential and kinetic energy develop as a function of the interaction $U$. The effective interaction becomes more relevant at larger $U$, leading to a reduction in the potential energy and an increase in the kinetic energy Haule (2007) (note the minus sign in $E_{\text{kin}}=-2\left\langle k\right\rangle/\beta$). We also find that the value of the double occupancy obtained within DMFT and DB theory with constant bosonic hybridization is similar to the one of the TPSC approach.

V Extended Hubbard model

We now turn our attention to the extended Hubbard model with finite $V$. We study this model on a cubic lattice with fixed local interaction $U=5$ and fixed temperature $\beta=2.5$. We use a $10\times 10\times 10$ lattice. Figure 8 shows that screening by the nonlocal interaction $V$ reduces the effective impurity interaction $U^{\prime}$, so that $U^{\prime}<U$. As expected, this increases the potential energy and decreases the kinetic energy. Even though it decreases, the impurity interaction is still repulsive for all cases shown here, $U^{\prime}>0$. Looking at the inverse susceptibility in the charge channel, we find that it approaches zero as $V$ is increased. Linear extrapolation 666Note that we find a slope of roughly $X^{-1}\propto-7V$. Without any correlation effects, a slope of $-6V=-V_{q=(\pi,\pi,\pi)}$ would be expected from Eq. (6). predicts the charge order transition to occur at $V\approx 0.99$. The arrow marks $V=U/z=0.83$, the point where charge-order becomes favorable in terms of the potential energy Vonsovsky and Katsnelson (1979). The actual transition occurs later due to the interplay between potential and kinetic energy and entropy.

The extrapolation here is based on data up to $V=0.95$. In principle, no true divergence is expected in a finite system. This has nothing to do with the Mermin-Wagner theorem, since the charge ordering is not associated with the breaking of a continuous symmetry. Instead, however, for finite-size lattices the self-consistency cycle becomes difficult to stabilize when the phase transition region is approached: small changes in the effective interaction $U^{\prime}$ will lead to big changes in the susceptibility $X$.

We can compare our results obtained with the self-consistent instantaneous interaction with those of Ref. van Loon et al. (2014a), where a frequency dependent interaction was used and self-consistency was only done at the EDMFT level. Figure 9 shows the charge susceptibility for these two methods in the cubic and square lattice Hubbard model. The qualitative behavior of both methods is the same, quantitative differences are visible. The self-consistent approach approaches the phase transition already at smaller values of $U$, that is, the inverse susceptibility is smaller close to the phase transition. The origin of this seems to be the enhanced charge susceptibility of the impurity model that originates in the reduced value of the effective interaction $U^{\prime}<U$, as was visible in Fig. 8. For the two-dimensional situation, the inverse susceptibility in the self-consistent DB solution is approximately linear in the regime accessible here. Based on Eq. (11), the self-consistency condition comes into play only when $X^{-1}\approx 1/L^{2}\approx 10^{-3}$ for the $L=32$ system studied here. For a much smaller $L=2$ system (dotted line), it sets in earlier and the dotted line bends upwards at small $V$. Therefore, the true CDW phase boundary can be obtained by extrapolating the results for the inverse charge susceptibility for different system sizes.

VI Conclusions and discussion

We have presented the Dual Boson approach with instantaneous interaction. By construction, this approach produces a susceptibility that satisfies the charge and spin conservation requirement Hafermann et al. (2014); van Loon et al. (2014b); Stepanov et al. (2016); Krien et al. (2017); Stepanov et al. (2018), and the Mermin-Wagner theorem. The instantaneous interaction assumption means that the method does not need an impurity solver that can handle retarded interactions, an important simplification that makes it more amendable to the simulation of multiband systems. We have illustrated the method in two- and three- dimensional systems and have analyzed how the proposed self-consistency condition affects the appearance of ordered phases. In our method the instantaneous interaction of the impurity model is adjusted by the “lattice” self-consistency condition on the bosonic hybridization function. Surprisingly, this constant hybridization function almost completely washes out the difference between the local part of the frequency dependent lattice and impurity susceptibilities. Mean-field based approaches typically overestimate the tendency towards ordered phases, as is clearest in two-dimensional or finite systems where phase transitions are forbidden. Enforcing self-consistency on the two-particle level can cure this deficiency Y.M. Vilk and A.-M.S. Tremblay (1997). Compared to DMFT (in the absence of the nonlocal interaction $V_{\mathbf{q}}=0$), we find that the feedback of collective excitations leads to an enhanced effective interaction and to a more correlated impurity model, as is visible in the double occupancy, kinetic energy and self-energy of both two and three-dimensional systems.

Acknowledgements.

The authors acknowledge useful discussion with Friedrich Krien. L.P., E.G.C.P.v.L., M.I.K. and E.A.S. acknowledge support from ERC Advanced Grant 338957 FEMTO/NANO. M.I.K. and E.A.S. acknowledge support from the Stichting voor Fundamenteel Onderzoek der Materie (FOM), which is financially supported by the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO). A.I.L. acknowledges support from the excellence cluster “The Hamburg Centre for Ultrafast Imaging - Structure, Dynamics and Control of Matter at the Atomic Scale” and North-German Supercomputing Alliance (HLRN) under the Project No. hhp00042.

Appendix A Analytical results

A.1 $U=0$, $t=0$ or $d=\infty$

In these three cases the more restrictive self-consistency condition Stepanov et al. (2016) $\chi_{\omega}=X_{\text{loc},\omega}$ is fulfilled for all $\omega$ by the solution $\Lambda=0$. This means that $\Lambda=0$ is also satisfies the instantaneous self-consistency condition used in this work. This solution corresponds to the exact DMFT solution of the Hubbard model in these three cases.

A.2 small $U$

We consider the Hubbard model, so $V=0$, and we look at the charge channel of the susceptibility. To lowest order in $U$, we can neglect self-energy corrections to the susceptibility and use geometric (RPA-like) equations to describe the impact of the interaction. The impurity model only knows about the impurity interaction $U^{\prime}$ and not about $U$ and $\Lambda$ separately,

[TABLE]

where the sign is $-$ for the charge channel and $+$ for the spin channel. The dual polarization operator $\tilde{\Pi}$ in the ladder approach is completely determined by the expectation values of the impurity model and by $t$ and $\Delta$, since these enter $\tilde{G}$. In particular, $\tilde{\Pi}$ only depends on $U^{\prime}$ and not on $U$ and $\Lambda$ separately, and the same holds for $\mathcal{X}$, which is also given by a geometric expression (It is essentially equal to the DMFT susceptibility of a Hubbard model with interaction $U^{\prime}$.).

[TABLE]

The lattice susceptibility according to DB is then given by Eq. (6).

[TABLE]

where $\Lambda$ is the interaction in the same channel as the susceptibility. Combining these equations, for the charge susceptibility we find

[TABLE]

where we have used that $\mathcal{X}=X$ and is independent of the channel at $U=0$. Similarly, we find

[TABLE]

Thus, the self-consistency conditions at small $U$ read

[TABLE]

These two equations can be solved numerically for $\Lambda^{\text{ch}}$, $\Lambda^{\text{sz}}$. Expanding the denominators and using the condition $\chi(U=0)=1/N\sum_{\mathbf{q}}X(U=0)$, we find

[TABLE]

where

[TABLE]

It is clear that $\Lambda^{\text{ch}}=-\Lambda^{\text{sz}}$ in the small $U$ regime described by these equations. In addition, both components of $\Lambda$ depend linearly on $U$. The magnitude of $\Lambda$ depends on $C$.

The temperature is implicitly contained in $\chi(U=0)$ and $X(U=0)$, which are Lindhardt bubble expressions, and in the sum over Matsubara frequencies. Figure 10 compares the perturbative formula with the numerical results of Fig. 1.

Particle-hole symmetry on hypercubic lattices implies that the change $U\rightarrow-U$ corresponds simply to the interchange of magnetic ($S^{z}$) and density fluctuations, without any difference on the one-particle level. Under the interchange of $\Lambda^{\text{ch}}$ and $\Lambda^{\text{sz}}$, the expressions given here are clearly consistent with this symmetry.

Appendix B Self-energy asymptotics

The general concept that enforcing some exact properties will break others is visible in the present instantaneous DB as well. An example is the high-frequency asymptote of the local self-energy of the paramagnetic Hubbard model,

[TABLE]

For the impurity model, with interaction $U^{\prime}\neq U$, we find

[TABLE]

At high frequencies, where the denominator $1+g\tilde{\Sigma}$ is equal to unity, the relation between dual and lattice self-energy is $\Sigma^{\text{lat}}=\Sigma^{\text{imp}}+\tilde{\Sigma}$. This means that the exact solution of the dual action should have

[TABLE]

Let us see how this exact expression can arise in dual perturbation theory. We consider a situation where the difference $U^{\prime}-U=\Lambda$ is small, so that dual perturbation theory is justified and we also consider $U^{\prime}$ sufficiently small that the vertices can be simplified by doing perturbation theory in the impurity model.

The bare bosonic dual propagator simplifies in the Hubbard model with instantaneous impurity interaction, it is local 777For the Hubbard model, $\tilde{X}^{(0)}$ is always local. The self-consistency condition $\sum_{\mathbf{q}}\tilde{X}^{(0)}=0$ then automatically ensures that $\tilde{X}=0$ everywhere and that Dual Boson reduces to Dual Fermion. The self-consistency condition used in this work does not require $\sum_{\mathbf{q}}\tilde{X}^{(0)}=0$ and $\tilde{X}$ is finite in the Hubbard model. and given by,

[TABLE]

Here, the last line is obtained since we are interested in small $\Lambda$. This shows that the number of bosonic propagators determines the order in $(U^{\prime}-U)$.

We identify three basic diagrams that could be relevant for the asymptote: a “Hartree” and a “Fock” diagram both containing two two-fermion-one-boson vertices and a diagram with a single two-fermion-two-boson vertex $\gamma^{(2,2)}$. The first two of these diagrams vanish since they contain a local Dual Fermion propagator, which is zero by the fermionic self-consistency condition. It is the third diagram, shown in Fig. 11, that we are interested in here. To evaluate it, we need to find a simplified expression for the vertex. In Dual Boson, vertices involving the bosonic degree of freedom typically have a “trivial” contribution, as has been discussed at length for the fermion-boson vertex Rubtsov et al. (2012); van Loon et al. (2014a, 2018a). It essentially originates in the fact that the number of Wick contractions is reduced when going from $c^{\dagger}c$ to $n$. This also applies to $\gamma^{(2,2)}$, which has a trivial contribution

[TABLE]

Evaluating the self-energy with symmetry factor $\frac{1}{2}$, and combining Eqs. (28) and (29) gives

[TABLE]

exactly the desired expression for the first order term.

Higher-order terms should appear in a similar fashion, although they become progressively harder to evaluate: multiple diagrams will contribute and the diagrams become more involved.

The interpretation of this result is two-fold. On the one hand, as long as $U^{\prime}-U$ is not too large, deviations between the exact and numerical asymptotics of the self-energy will also stay small. On the other hand, it shows that a diagram that includes a three-particle vertex is already needed to recover only the lowest order term in the asymptotic expansion of the self-energy. This occurence of three-particle vertices to achieve consistency is reminiscent of Refs. van Loon et al. (2015); Krien et al. (2017); van Loon et al. (2018b).