Dual Boson Diagrammatic Monte Carlo Approach Applied to the Extended Hubbard Model

In this work we introduce the Dual Boson Diagrammatic Monte Carlo technique for strongly interacting electronic systems. This method combines the strength of dynamical mean-filed theory for non-perturbative description of local correlations with the systematic account of non-local corrections in the Dual Boson theory by the diagrammatic Monte Carlo approach. It allows us to get a numerically exact solution of the dual boson theory at the two-particle local vertex level for the extended Hubbard model. We show that it can be efficiently applied to description of single particle observables in a wide range of interaction strengths. We compare our exact results for the self-energy with the ladder Dual Boson approach and determine a physical regime, where description of collective electronic effects requires more accurate consideration beyond the ladder approximation. Additionally, we find that the order-by-order analysis of the perturbative diagrammatic series for the single-particle Green's function allows to estimate the transition point to the charge density wave phase.

In this work we introduce the Dual Boson Diagrammatic Monte Carlo technique for strongly interacting electronic systems. This method combines the strength of dynamical mean-filed theory for non-perturbative description of local correlations with the systematic account of non-local corrections in the Dual Boson theory by the diagrammatic Monte Carlo approach. It allows us to get a numerically exact solution of the dual boson theory at the two-particle local vertex level for the extended Hubbard model. We show that it can be efficiently applied to description of single particle observables in a wide range of interaction strengths. We compare our exact results for the self-energy with the ladder Dual Boson approach and determine a physical regime, where description of collective electronic effects requires more accurate consideration beyond the ladder approximation. Additionally, we find that the order-by-order analysis of the perturbative diagrammatic series for the single-particle Green's function allows to estimate the transition point to the charge density wave phase.

I. INTRODUCTION
Strongly correlated systems represent a formidable challenge in condensed matter physics. For this reason, the study of model systems can allow us to investigate the effects of strong interactions and analyse the effects of different approximations. Among these models, the Hubbard model [1] has been extensively studied in the past decades due to its capacity of successfully describing the emerging physics of some classes of strongly correlated materials, where local interactions are assumed to be much stronger than non-local ones. A major breakthrough in solution of the Hubbard model was made by dynamical mean-field theory (DMFT) [2]. This method becomes exact in the limit of infinite spacial dimensions or connectivity of the lattice [3], and serves as an accurate approximation for single-particle quantities in finite dimensions [4].
At the same time, many real materials exhibit interesting physical effects, such as a charge density wave (CDW) phase, that can not be described by a local Hubbard interaction term alone. In order to consider these phenomena, non-local interactions have to be taken into account. For this aim, in analogy with DMFT, an extended dynamical mean field theory (EDMFT) has been developed [5][6][7][8][9]. However, in this approach the self-energy and polarization operator are local, meaning that they are frequency dependent, but not momentum dependent. Extensions of EDMFT, such as GW+EDMFT [10][11][12][13][14][15][16], Dual Boson (DB) [17][18][19][20], TRILEX [21][22][23], and TRILEX 2 [24] approaches have been developed to cope with this issue. In particular, the Dual Boson and TRILEX 2 techniques are based on the exact transformation that allows to rewrite the initial action of the extended Hubbard model in terms of the local impurity problem, which can be solved numerically exactly [25][26][27][28][29], and a diagram-matic series around the impurity. Since the effective impurity problem already includes the main contribution coming from local correlations, this looks a naturally convenient starting point for perturbation theory and approximate approaches.
So far, calculations in the framework of the Dual Boson approach have been performed only in the ladder approximation [30][31][32][33][34]. This approach is based on the calculation of a specific sub-set of diagrams that, in principle, can be justified by physical considerations only in the regime of developed collective electronic fluctuations [35,36]. However, the ladder DB approximation provides remarkably good results in comparison with other advanced methods, as for instance dynamical cluster approximation [37,38]. An alternative approach that involves solution of parquet equations based on dual theories was recently proposed in Ref. 39. Comparison between various methods based on extensions of (E)DMFT can be found in Ref. 40.
Another route to study strongly correlated systems has recently been attempted: the use of unbiased methods based on the combination of diagrammatic approaches with Markov chain Monte Carlo [41]. In particular, the bare diagrammatic Monte Carlo (DiagMC) method has been successfully applied to the Hubbard model at weak and moderate Coulomb interactions [42]. This method starts from an expansion in terms of the Hubbard coupling U and constructs all Feynman diagrams up to some finite but high order in U. The algorithm allows to sample all possible diagrams without any restriction to specific topologies. Efficient algorithms that express all connected diagrams of the perturbative expansion up to a given order by means of determinants [43] have been developed for various observables and correlation functions [44][45][46][47], significantly reducing the computational cost of the calculation. Approaches based on a small-coupling expansion work very well in the regime of small to moderate couplings, but start to fail when U is of the order of half of the bandwidth [4,48]. These failure is related to the finite convergence radius of the diagrammatic series and can be improved using resummation techniques [44].
To allow for a non-perturbative treatment of strong correlation effects, a diagrammatic Monte Carlo scheme based on Dual Fermion approach [49] was proposed in Refs. 50 and 51. The advantage of this method in comparison with diagrammatic expansions in terms of the bare Coulomb interaction U is that the impurity problem already accounts for the main effects of local correlations, which strongly screen the bare interaction U. The expansion is thus performed in terms of the renormalized local interaction vertex function, which appears to be naturally more convenient at moderate and strongly interacting regime. Additionally, the diagrams are sampled in continuously in the momentum space without the discretization of the Brillouin zone. Hence, the result of the calculation is not influenced by any finite-size effects.
In this paper we generalize this approach to the extended Hubbard model, performing an additional dual transformation that introduces effective bosonic fields related to non-local interactions. Our Dual Boson diagrammatic Monte Carlo (Di-agMC@DB) method combines the advantages of DMFT, because it already accounts for the screened local interaction in the impurity problem, with the capability of sampling all the possible Feynman diagrams without any restriction. The result is an efficient diagrammatic Monte Carlo algorithm that naturally incorporates non-local Coulomb interaction in the original DiagMC@DF approach [51].

II. DUAL BOSON THEORY
Our starting point is the extended Hubbard model in the action formalism Here, c ( * ) kνσ are Grassman variables corresponding to the annihilation (creation) of electrons with momentum k, fermionic Matsubara frequency ν and spin σ. We also introduced the electronic dispersion ε k and chemical potential µ. The model additionally includes an on-site Coulomb interaction of strength U in terms of the electron density n qωσ at momentum q and bosonic Matsubara frequency ω, as well as a non-local interaction V ς q , where the index ς represents charge (ς = ch) or spin (ς = sp = {x, y, z}) degrees of freedom. Variables ρ ς qω = n ς qω − n ς are expressed in terms of composite quantities n ς qω = kν,σσ c * k+q,ν+ω,σ σ ς σσ c q,ω,σ . In the previous expression σ ch = 1, and σ x,y,z is the corresponding Pauli matrix in spin space. The general idea of dual theories is to split the action into two parts: a local impurity problem, that contains the full local interaction, and a non-local part that can be treated perturbatively. Instead of directly applying a perturbation theory to the non-local part, a transformation that introduces new variables is performed. This allows to dress the non-local part with the local impurity quantities. An additional important consideration is that, once the DMFT impurity is chosen, the dual theories represent a diagrammatic expansion around the DMFT solution. This starting point for the perturbation theory looks appealing, since the DMFT already accounts for local many-body effects, which allows to correctly reproduce both the small and large U limits. In order to perform this transformation, we add and subtract an arbitrary fermionic hybridization function ∆ ν , so that we can isolate a local impurity part of the action. With this fermionic hybridization function the action reads S = i S (i) imp + S nonloc , where the impurity part is and the non-local part reads In the following, imp denotes the average with respect to the local action (2). The impurity problem of Eq. (2) can be solved exactly using continuous-time quantum Monte Carlo solvers [25][26][27][28][29]. In the same way we could include a bosonic hybridization function [17][18][19]. However, this step would require an additional discussion of the self-consistency condition needed to determine it. Therefore, we exclude the bosonic hybridization from the current discussion in order to reduce the number of external parameters in the system. The hybridization function can be defined in an arbitrary way, but some choices are more convenient than others. In the rest of the paper we will use ∆ ν obtained from DMFT calculations.
The dual boson transformation amounts to perform a fermionic and a bosonic Hubbard-Stratonovich transformations over the non-local part of the action S nonloc , which introduce new dual fermionic variables f , f * and bosonic φ ς fields. The action obtained after this transformation is quadratic in the electronic operators c ( * ) , so we can integrate them out [32]. The original problem of interacting electrons is then recast into a new problem in terms of the dual degrees of freedom only. Sigle-and two-particle observables of the original electron system can be exactly related to dual correlation functions, as shown for example in Ref. 32 The bare dual propagators are defined as where g ν and w ς ω are the Green's function and renormalized interaction of the auxiliary impurity problem, respectively, and the impurity susceptibility χ ς ω = − ρ ς −ω ρ ς ω imp . Additionally, α ς ω = 1 + U ς χ ς ω with U ch/sp = ±U/2. The choice of the Matsubara frequency space is natural in this case, because in Eq. (5) the ∼ ν −1 part of the tail in G EDMFT kν is exactly canceled by g ν . This means that the dual fermion Green's function decays as fast as ∼ ν −2 , and there are no convergency issues related to summations over Matsubara frequencies. The interaction term truncated at the two-particle level is given bỹ where Λ ς νω and Γ σσ σ σ νν ω are the impurity fermion-boson and fermion-fermion vertex functions, respectively. We use the same definitions as in Ref. 24, that in terms of impurity variables explicitly read In general, the interaction term also contains all the higher order vertices that conserve the number of dual fermions, but we will limit our study to the two-particle interaction terms only. Terms beyond this approximation were shown to lead to very small corrections in many regimes [52]. As a matter of fact, dual theories with only two-particle vertex functions show a rather good agreement with other unbiased methods, and it is still under debate if deviations with other methods are due to higher-order vertices or to different effects [50,51,53]. Additionally, the inclusion of higher order vertices would enormously increase the complexity of the diagrammatic Monte Carlo scheme. In light of all these considerations, we exclude them from our calculations.
In our case, the solution of the impurity problem is obtained using a continuous-time Monte Carlo solver based on hybridization expansion (CT-HYB) [29]. This gives an access to all the impurity observables needed for the construction of the dual boson diagrammatics. In particular, we compute ∆ ν , g ν , and χ ω for the construction of bare propagators, as well as the fermion-fermion vertex Γ σσ σ σ νν ω and the fermion-boson vertex Λ ς νω . Within this approximation, the dual action (4) with the interaction (7) is quadratic in the bosonic fields. This means, that it is possible to integrate dual bosonic degrees of freedom out exactly and obtain a fully fermionic action. The bosonic Hubbard-Stratonovich transformation is necessary for decoupling of the non-local interaction term, that otherwise would prevent the integration of the local impurity action out. Moreover, the introduction of the bosonic variables dresses the interaction in terms of the impurity quantities, so that the bosonic propagator is already partially screened. In order to construct a form of the full fermion-fermion vertex after the FIG. 1. Schematic diagrammatic interpretation of Eq. 12. The full antisymmetrized fermion-fermion interactionΓ (gray box) is a combination of the impurity vertex Γ (white box) and processes involving a boson exchange (wiggly line). The full vertex acquires a momentum dependence due to the presence of the bosonic lines. White and black dots represent incoming and outgoing particles, respectively. Triangles represent Λ νω vertices. The exact dependence on the channel indices and prefactors is shown in Eq. 12.
integration of the bosons, it is useful to decompose the impurity fermion-fermion vertex (9) in channel representation as The result is a modified dual fermion action where ξ is a formal expansion parameter, which must be set to unity (ξ = 1) in the actual calculations and keeps track of the expansion order. Importantly, we introduced a new momentum dependent fermion-fermion vertex that combines the vertex function of the local impurity problem and the nonlocal interaction between fermions mediated by dual bosonic fields Here,M ς,q ν,ν ,ω = Λ ς ν,ωW ς qω Λ ς ν +ω,−ω . Since Eq. (10) holds for both vertices Γ νν ω and Γ kk q νν ω , it is possible to switch easily between the two representations using the relations Additionally, all the other non-zero components can be simply obtained by applying the SU(2) symmetry in Eq. (10) (14) or by exploiting the fact that a simultaneous flipping of all the spins leads to the same result in paramagnetic case. We note that the structure of the new fermion-fermion vertex function (12) shown in Fig. 1 is reminiscent of the TRILEX 2 theory [24], and appears due to the antisymmetrized form of the interaction. Integrating out of bosonic fields is very important for the calculation of diagrams, because it allows to eliminate the bosonic degrees of freedom from the theory analytically and to avoid their sampling in diagrammatic Monte Carlo. In our implementation we compute the dual self-energy. In order to obtain the single-particle observables of the lattice problem, we can use the standard equation that relates the dual self-energy to the lattice self-energy Σ kν , is the self-energy of the impurity problem, as shown for example in Ref. 32. The lattice Green's function can be obtained via the usual Dyson equation from the lattice self-energy or using its relation with the dual Green's function [32].

III. DIAGRAMMATIC MONTE CARLO SCHEME
The algorithm tested in this paper is an extension of the DiagMC@DF method proposed in Refs. 50 and 51. Our DiagMC algorithm computes numerically exactly the coefficients a n (k, ν) in the expansion of the dual self-energỹ for the action (11) up to some maximum order N max . The value of the dual self-energy can be recovered by setting ξ = 1. We will call it diagrammatic Monte Carlo for Dual Bosons (DiagMC@DB). In the same way as the original algorithm, our method is based on bare diagrammatic Monte Carlo approach [54]. This algorithm allows to construct all the Feynman diagrams up to any finite order and to sum over them using Markov chain Monte Carlo. According to Refs. 41 and 55, any correlation function O can be expressed as a sum of diagrams as follows where y is a combined index that contains all the dependence on external points, n indicates the number of vertices that appear in the diagram, C n are the topologies, and O C n is the value of a specific diagram. Additionally x i is shorthand notation for the internal degrees of freedom (k, ν, σ) i corresponding to momentum, Matsubara frequency and spin that originate from the presence of loops of Green's functions. This statement is true provided that the limit in Eq. (17) is well defined and convergent for the chosen parameters as N max → +∞. Divergencies of the diagrammatic series are often related to physical instabilities, as we show in Sec. IV, or to some unphysical behavior of the starting point, for example the antiferromagnetic phase transition of DMFT [51]. The summation over the perturbation order n, topologies and internal degrees of freedom is performed using a Metropolis-Hastings scheme [54], where the function to be sampled is the sgn (O), and the probability distribution is given by he amplitude | O| in order to respect the requirement of positive weight function. This approach automatically satisfies the detailed balance condition for the Markov chain (see Ref. 56), given that the acceptance probability to go from a configuration C n to another configuration G m is constructed as where P C n and P G n are the probabilities of the initial and final configuration respectively. There are no substantial changes from DiagMC@DF in the acceptance-rejection scheme adopted, except for the fact that in our case the bare fermion-fermion vertex function (12) is momentum dependent. Each contribution to the series expansion (17) can be written as a combination of two kinds of diagrammatic elements: fermionic lines that represent dual Green's functions G (called also propagator lines) and verticesΓ described in Eq. (12). Each vertex is attached to four propagator lines, two incoming and two outgoing. The terms at order n in the expansion are represented in terms of Feynman diagrams with n vertices connected by lines in all the possible combinations. These diagrams give an intuitive and efficient picture that allows us to design the updates so that all the contributions to the expansion (17) can be generated by changing how the vertices are connected to each other by mean of the propagator lines. In particular, we use the same worm algorithm described in the Ref. 51 to update the diagram topologies. The aim of the worm algorithm is to enforce momentum conservation, which is a non-local property of the diagram, by means of updates that act locally on few elements of a diagram. The worm algorithm introduces a set of unphysical updates that allow the transition between all the different possible topologies contributing to the dual self-energyΣ kν . This means that we sample all the diagrams with one incoming and one outgoing line that are also irreducible with respect to a cut of a fermionic line. This can be practically implemented by the condition that no internal line can carry the same momentum and frequency dependence of the external lines [51]. All the subtleties and details related to the implementation are discussed in detail in Ref. 51.
Each configuration is identified by an ordered set of vertices, where each vertex is stored together with the incoming and outgoing frequencies, momenta, spins and the connections with the other vertices. The original implementation of Ref. 51 works with unsymmetrized diagrams, in order to avoid topology-dependent prefactors. However, the local vertex Γ itself has an antisymmetric form in spin space, and we find convenient to introduce also the non-local vertex corrections M ς, q ν,ν ω in the antisymmetrized form shown in Eq. (12). The corresponding vertex in the unsymmetrized diagrammatic theory can be obtained by simply dividing this vertex by two [51].
The simultaneous sampling of contributions coming from the fermion-fermion scattering and boson exchange processes efficiently reduces the number of topologies. On the other hand, we do not distinguish between local and non-local di- agrams, so that we can not exclude local diagrams from the beginning simply by a proper choice of the DMFT selfconsistency condition, as it was done in DiagMC@DF calculations. Using the Metropolis-Hastings scheme allows us to compute observables up to a normalization factor. In order to keep track of the normalization, we sample the absolute value of an additional diagram that we can calculate explicitly outside Monte Carlo and we store its value in a suitable accumulator N norm . The chosen diagram is simply a single vertex with unitary value with the upper corners connected by a single bare dual Green's function. Its value is given by which is computed directly from the analytical expression for the bare dual propagatorG. The normalized dual self-energỹ Σ kν is then straightforwardly computed from the normalization accumulator N norm using the following equatioñ IV. RESULTS

A. Computational details
We perform our calculations on a 2D square lattice with the nearest-neighbor hopping amplitude t = 1 that fixes the energy units. We would like to stress, that this is not a limitation of the method, which works with any dispersion and with a general form of the interaction as a function of momentum. We start from the description of the output of the calculation, namely the dual self-energyΣ kν . Since diagrammatic Monte   FIG. 3. Most important self-energy diagrams. Top row shows the only nonzero contribution to the first order diagram, taking into account that we can not connect two slots of the same local vertex with a propagator line due to DMFT self-consistency condition. The middle row shows the second-order diagramsΣ (2) . If V is small compared to U/4, it can be approximated by the second-order dual fermion diagram on the right hand side. The last termΣ corr shows few diagrams that enterΣ in our calculations, but are not included in the ladder DB.
Carlo sums diagrams up to a given order, one has to check that the result is converged with respect to the order. In Fig. 2 we show a converged output of our calculations. In particular, the result for the maximum order of the diagrammatic expansion, order 5, differs from order 4 of ∼ 1%, hence we consider the result converged. Practically, this means that the performed calculation can be considered converged already at order 4. This is the case for most of presented calculations.
Far away from instabilities, it is not possible to observe any improvement in going beyond the 5th order of the diagrammatic expansion, and the computation necessary to achieve convergence at the order 6 increases significantly. Thus, a standard DiagMC@DB computation requires around 12 hours with a hundred parallel runs in order to obtain a converged result at the 5th order. Instead, for a converged result at the 6th order, the required computational time increases to more than 24 hours in order to reach a reasonable accuracy. For this reason, all results presented in this work are calculated up to 5th order of expansion, except the ones that are used for the analysis of the phase transition.
An important remark is that the contribution coming from the non-local interaction can be quite large, even up to values around 3/4 of the bandwidth. Additionally, we observe that the main contribution to the real part of the self-energy comes from two kind of diagrams that are shown in Fig. 3. The first is the single boson diagramΣ (1) that contains only one fac-torM ς, q ν,ν ω , which is the only non-zero contribution at the first order of the diagrammatic expansion in terms of the vertex function. This can be already seen in Fig. 2, where the first order contribution accounts for around 50% of the real part of the dual self-energy. The second important contribution is the second order dual fermion diagramΣ (2) , that contains two fermion-fermion vertices connected to each other. At values of V far away from the CDW instability, other contribution to ReΣ are rather small compared to these ones. On the other hand, the imaginary part of the self-energy is much more sen- sitive to higher order corrections. In Fig. 2 we can see that the second order is way off compared to the third order, accounting for only around 50% of the contributions to ImΣ. Interestingly, the third order already accounts for most of the contributions. We deduce, that the inclusions of third-order diagrams in our expansion that contain multiple fermion-fermion scattering and bosonic exchanges are important for the imaginary part of the self-energy. These diagrams contribute to around 40% of ImΣ at high symmetry points for U = 5, and their impact on dual quantities becomes even more important at larger U. Orders larger than the third typically amount to a correction of less than 10% of ImΣ at high symmetry points.
However, the overall momentum dependence of the lattice self-energy is still dominated by ReΣ in the regimes where U ≤ 4 or U ≥ 8, because of the denominator in Eq. (15), as shown also in Ref. 51. The most important corrections to the lattice self-energy coming from ImΣ appear exactly in the regime between half of the bandwidth and the bandwidth, where it can account for around 40% of the difference with DMFT solution at high symmetry points. Even though second-order is thought to already account for the most important contributions far away from instabilities, as shown in Ref. 4, the inclusion of two-boson exchanges and third-order  In all the figures we show the result of the calculation at the 5th order. In particular, we show results from a quarter of the bandwidth U = 2 up to the bandwidth U = 8. We note that the agreement between these two methods is substantially exact up to a half of the bandwidth for all considered values of the non-local interaction. In fact, in this regime the ladder DB result for the lattice self-energy lies inside the error bars of the DiagMC@DB calculation. For larger values of the on-site Coulomb interaction exceeding the half of the bandwidth, the difference between two theories is more noticeable, especially for small strength of the non-local interaction V. In order to quantify the difference between these two methods, we look at the M = (π, π) point in the momentum space and calculate the following quantity where Σ M,ν 0 is the difference between the self-energy in the specified method and the DMFT self-energy Σ imp ν 0 (15), and ν 0 is the lowest positive Matsubara frequency. We measure differences with respect to DMFT self-energy, because the latter is constant in momentum space and quite large. If we want to resolve relatively small differences in momentum space, we have to exclude its contribution. Additionally, we choose the M point, because it shows the largest difference between the two curves in the Brillouin zone. In this way, we are sure that the δ M parameter contains information only about the maximum mismatch coming from the dual corrections. The reason for taking the real part of this quantity comes from the fact that the imaginary part of the dual self-energy ImΣ shows a systematic shift already at V = 0, i.e. at the Dual Fermion level (see Ref. 51). Here, we aim to assess the behaviour of the self-energy as a function of the non-local V rather than to investigate this aspect. The result for the mismatch parameter δ M is summarised in a tentative phase diagram shown in Fig. 8. We can conclude that the difference between DiagMC@DB and ladder DB approaches is negligible at small U below the half of the bandwidth and further increases with the local interaction. This behavior can be explained considering that for small local Coulomb interaction U the regime is still perturbative in the Dual Boson theory, so we expect all the methods to give quantitatively similar results. This finding is also in agreement with the result of DiagMC@DF calculations [51] obtained for the zero non-local Coulomb interaction. On the other hand, we observe that the mismatch is more severe at V = 0 and decreases as V increases. Indeed, when the non-local Coulomb interaction is large, charge fluctuations in the horizontal channel are expected to give the main contribution to physical observables such as self-energy and susceptibility [36], because the system lies close to the charge density wave (CDW) phase. Ladder DB approach accounts for this kind of fluctuations by construction, and for this reason the mismatch δ M decreases. From Fig. 8 we find, that the values of U at which the largest mismatch occurs lie in the region where the phase transition to the Mott-insulating state was predicted by cluster DMFT [57] and dual fermion [58] calculations at lower temperature. This means, that in this regime contributions not included in the ladder approximation cease to be negligible. These contributions corresponds to bosonic lines in a direction orthogonal to the ladder direction (see, e.g.Σ corr in Fig. 2), which are included in DiagMC@DB. It means, that the correct description of this phase transition, especially at V = 0 would require advanced approaches beyond the ladder approximation. Another consideration that emerges from this analysis as a function of U and V is that up to a half the bandwidth a momentum dependence of the real part of the self-energy at V = U/4 is dominated by the non-local interaction V. It plays a very important role even for U = 6, where we would expect the local interaction to give the most important contribution.
C. Monitoring the CDW phase transition from single particle observables In the current implementation, the DiagMC@DB theory is based on a solution of a single-site impurity problem, which does not allow for a description of broken-symmetry phases. In particular, this results in a divergence of the infinite diagrammatic expansion in terms of bare dual quantities at the phase transition [18,20,32,36,51]. The most interesting phase of the extended Hubbard model that is associated with the presence of the non-local Coulomb interaction is the charge density wave phase, i.e. a checkerboard configuration in the real space with alternating empty sites and doubly occupied sites. The phase transition to this state occurs when V is large enough to overcome the effect of the on-site Coulomb repulsion that favours a single-electron occupation of lattice sites. For the nearest-neighbor Coulomb interaction V this qualitatively happens at V U/z, where z is the number of neighboring lattice sites.
The description of the system inside the CDW state requires an inclusion of symmetry breaking terms in the theory. However, the instability can be identified already in the normal phase studying the charge susceptibility [18,20,32]. In particular, we expect the susceptibility to show a very sharp peak when the instability occurs. This trend can be seen in the upper right panel of Fig. 9, where the inverse of the charge susceptibility linearly decreases. Our ladder Dual Boson calculations predict a transition point V CDW at V = 0.77 for U = 2.5 and V = 1.09 for U = 3.5. However, the critical value V CDW for the CDW phase transition can also be found in a controlled way from the analytic structure of the dual self-energyΣ as a function of the complex expansion parameter ξ. Since the dual action (11) is explicitly constructed for the translationallysymmetric phase, the critical point V = V CDW is marked by a singularity appearing in the functionΣ(ξ) at ξ = 1. When V is increased beyond V CDW in the symmetry-broken phase, this singularity must move along the real axis towards the origin for the physicalΣ(ξ = 1) to remain inaccessible by its power-series expansion (16). The method introduced in Ref. 59 and routinely applied in the context of DiagMC [44] allows to accurately evaluate the specific location of the singularity ξ s . It assumes a generic power-law behavior near the singularity, which is typical for a continuous phase transition, Σ(ξ) ∝ (ξ s − ξ) η for |ξ − ξ s | 1 with some real number η, and extracts ξ s from the behavior of the series coefficients a n in Eq. (16). As shown in Ref. 44, ξ s can be found from a finite number of coefficients {a n } with a reliable error bar that includes both the systematic and statistical (Monte Carlo) error. The result of this procedure for different values of V is shown in the bottom right panel of Fig. 9, where ξ s (V) is obtained from {Re a n (k, ν 0 )}, for n = 1, . . . , N max = 6 projected on the first A 1g -symmetric harmonic ψ s (1,0) (k) = cos(k x ) + cos(k y ) to produce a numerical series from the functional one. The condition ξ s (V) = 1 then gives the critical value V CDW = 0.82 (1).
In combination with this controlled method, we propose an additional empirical way to obtain the instability point. It is important to notice that the checkerboard configuration of electrons is insulating. This means that strong charge fluctuations create a pseudogap in the electronic spectrum, which can be detected calculating the spectral function. In particular, the spectral function at the Fermi energy is directly connected to the local Green's function G(β/2) calculated at imaginary time τ = β/2 by the relation A(E F ) ≈ −β G(β/2)/π (see for example Ref. 18), without the need of analytical continuation from Matsubara to real frequencies. This situation is conceptually similar to the antiferromagnetic pseudogap, but the analysis in the framework of our theory is very different. In fact, the divergence of the dual fermion series is not associated with a true physical instability, as discussed previously, hence the divergence of the diagrammatic series in terms of the local interaction does not have a clear physical interpretation. On the contrary, the non-local interaction V enters only the effective fermion-fermion vertex function of Eq. (12) in a trivial form through the dual boson propagator, which up to a local prefactor is proportional to [32] We consider the local Green's function obtained by replacing the dual self-energy up to order n into the Dyson equation G (n) (β/2) as a function of the non-local interaction V, keeping the local interaction U fixed. If we inspect left panel of Fig. 9, the behaviour of this function resembles a Fermi function where V * n is the critical value of the non-local interaction at which the function shows a steep drop, and α n is a numerical value that defines the broadening of the Fermi function at order n. From these empirical and physical considerations, we expect that V * n → V CDW as the order n → +∞, which means that, if we extrapolate the central point of the Fermi function as a function of n, we expect it to converge to the value V CDW . Results based on this analysis are plotted in the left panel of Fig. 9, and the expected behaviour is clearly visible in the figure. Fitting the value of V * as a function of the order n allows to extrapolate the value of V CDW by letting the order go to infinity. The inset in the left panel of Fig. 9 clearly shows the convergence of the values of V * n to some finite value as the order increases. We compared the extrapolated values with the susceptibility at M = (π, π) calculated from ladder calculations. With this simple analysis we obtained values for V CDW compatible with the ladder results for U = 2.5 and U = 3.5 at β = 4. The specific values are V CDW = 0.77 (2) at U = 2.5 and V CDW = 1.09(2) at U = 3.5. These points are highlighted in red in Fig. 8. The methods presented in this section show how order-by-order or cumulative analysis of the series allows for an accurate extrapolation of results in the limit of infinite order of perturbation expansion. Our results for U = 2.5 show a very good agreement with the value obtained with GW+DMFT in Ref. 16. Additionally there is a good agreement with previous dual boson [18,32,36] and DCA calculations [37,60].

V. CONCLUSIONS
Even though the method presented in the previous sections is not exact, because higher-order impurity vertices are neglected, the DiagMC@DB scheme allows to include contributions coming from all the possible diagrams with no restriction on a particular class of topologies. In other words, Di-agMC@DB is the exact solution of the dual boson action truncated at the level of two-particle scattering. Due to this consideration, the results provided by DiagMC@DB are based on theoretically much more stable grounds than other approximations based on partial resummation of specific diagrams, as in the case of the ladder Dual Boson. Additionally, there are no finite-size effects since we worked in momentum space and Matsubara frequencies directly.
In our calculations we observed a very accurate agreement between DiagMC@DB and ladder calculations for U up to half of the bandwidth. Even above the bandwidth, the ladder Dual Boson seems to capture the main contributions and the difference between the two methods is quite small. In fact, we have never observed a value of the δ M parameter bigger than 8% in the region of the parameter space where series converges. This offers a further validation of the ladder Dual Boson technique over a very wide range of interaction strengths. The presence of strong non-local interaction V enhances nonlocal bosonic excitations in the charge channel that are accounted for in ladder approximation. This can be captured looking at the real part of the self-energy, which coincides for the two methods at large values of V.
The advantage of the DiagMC@DB is that is allows to consider diagrams order by order and investigate the convergence properties of the series in an unbiased and systematic way. In particular, starting from DMFT impurity, the solution is a series is terms of a complicated function of V. Since V doesn't enter the impurity, it is possible to use resummation techniques to estimate the value of V at which the charge order occurs already from the study of single-particle quantities. Different choices of the hybridization functions, obtained for instance from ladder Dual Boson calculations, can extend the convergence radius of the series at lower temperatures (see Ref. 51).
Another strategy that could improve the efficiency of sampling could be the formulation of the series in terms of the semi-bold Green's function, in which some diagrams are included in the bare dual propagator from the very beginning, or the fully bold Green's function, substantially reducing the configuration space. It is expected that both approaches could improve the convergence properties as well, but a systematic study is required. We are currently working in the direction of extending this method to calculate also two-particle observables in two-dimensional heterostructures. In addition, the inclusion of a checkerboard configuration with two nonequivalent sublattices can allow to study the Extended Hubbard model inside the broken symmetry phases, as the CDW phase or the antiferromagnetic phase.