Multiloop functional renormalization group for the two-dimensional Hubbard model: Loop convergence of the response functions
We present a functional renormalization group (fRG) study of the two dimensional Hubbard model, performed with an algorithmic implementation which lifts some of the common approximations made in fRG calculations. In particular, in our fRG flow, (i) we take explicitly into account the momentum and the frequency dependence of the vertex functions; (ii) we include the feedback effect of the self-energy; (iii) we implement the recently introduced multiloop extension which allows us to sum up all the diagrams of the parquet approximation with their exact weight. Due to its iterative structure based on successive one-loop computations, the loop convergence of the fRG results can be obtained with an affordable numerical effort. In particular, focusing on the analysis of the physical response functions, we show that the results become independent from the chosen cutoff scheme and from the way the fRG susceptibilities are computed, i.e., either through flowing couplings to external fields, or through a ”post-processing” contraction of the interaction vertex at the end of the flow. The presented substantial refinement of fRG-based computation schemes paves a promising route towards future quantitative fRG analyses of more challenging systems and/or parameter regimes.
pacs:71.10.Fd, 71.10-w, 71.27.+a
Over the last two decades, functional renormalization group (fRG) methods have been broadly used for analyzing two-dimensional (2D) lattice electron systems (for reviews, see Refs. [Metzner et al., 2012; Platt et al., 2013]). The main advantage of the fRG lies in the exploration of the leading low-energy correlations and instabilities towards long-range ordered states, similar to what has been investigated earlier for one-dimensional systemsSólyom (1979); Bourbonnais et al. (2002). However, in one dimension, other methods like Bethe-Ansatz, bosonizationvon Delft and Schoeller (1998); Giamarchi (2004) and DMRGJeckelmann (2002) exist, which are for certain aspects more controlled. Hence, assessing the precision of RG methods in one-dimensional systems was not really in the foreground. The situation evidently changes for two- and three-dimensional systems, where the specific simplifications associated to the peculiar one-dimensional geometry are not applicable. At the same time, spatial correlations in 2D are strong enough to induce qualitative correctionsSchäfer et al. (2015); Rohringer and Toschi (2016) with respect to another class of rigorous many-body approaches, such as the Dynamical Mean-Field Theory (DMFT)Georges et al. (1996); Kotliar et al. (2006); Held (2007) which allows to include all purely local dynamical correlations.
In fact, due to the intrinsic complexity of the many-electron problem in 2D, the development of unbiased quantitative methods applicable to a wide energy range from electronic structures on the scale of a few eV down to, e.g., ground state ordering in the (sub-)meV region, is still on the wishlist. This goal has motivated, in the last decade, the development of several algorithmic schemes for treating electronic correlations in 2D from different perspectivesMaier et al. (2005); Metzner et al. (2012); Rohringer et al. (2018). In this context, the fRG has already unveiled quite promising features: The fRG has the potential of resolving band structures and Fermi surface details and to treat competing orders on low energy scales in a rather unbiased way, since it does not require preliminary assumptions about dominating scattering channels. Recent applications range from studies of cuprate high- superconductorsUebelacker and Honerkamp (2012a); Eberlein (2014) over iron superconductorsPlatt et al. (2013); Lichtenstein et al. (2014) to few-layer graphene systemsScherer et al. (2012a, b), to cite a few.
We also note that, while the current applicability of the fRG is generally restricted to the weak to intermediate coupling regimes, its combinationTaranto et al. (2014); Metzner (2015) with the DMFT might allow, in the future, to access much more strongly correlated parameter regions, including the ones in proximity of the Mott-Hubbard metal-insulator transition. This is achieved by constructing a fRG flow starting from the DMFT solution of the considered lattice problem to the exact solution, i.e., in practice, using the DMFT to determine the initial conditions for the fRG flowTaranto et al. (2014). Similarly to other diagrammatic extensionsRohringer et al. (2018) of DMFT, such as the Dynamical Vertex Approximation (DA)Toschi et al. (2007) or the Dual FermionRubtsov et al. (2008) approach, one might work either with the physical degrees of freedom (as in the so-called DMFRGTaranto et al. (2014)) or in the space of auxiliary (dual) fermionsWentzell et al. (2015), introduced by means of a suitableRubtsov et al. (2008); Rohringer et al. (2018) Hubbard-Stratonovich transformation.
Yet, what is hitherto missing is a thorough analysis of the quantitative reliability of the fRG for a well-defined test case. More precisely this would require to clarify how much the fRG results, going beyond the correct estimation of general physical trends, depend on the approximations inherent in the used fRG scheme. This study within the fRG would then also provide a solid basis for future comparisons with other numerical techniques.
The mentioned approximations can be grouped in three categories:
(i) Momentum/frequency discretization: As the fRG algorithm typically exploits the flow of vertex functions that depend continuously on multiple momenta and frequencies, various approximations are performed to mitigate numerical and memory costs. Early on, -patch discretizations of the momentum dependencies through the Brillouin zone were used. Later, it was noticed that channel-decompositions in conjunction with form factor expansionsHusemann and Salmhofer (2009); Giering and Salmhofer (2012); Wang et al. (2013) lead to physically appealing approximations featuring advantageous momentum resolution and numerical performanceLichtenstein et al. (2017). Clever prescriptions for the treatment of the high-frequency tails of the vertex function have been devisedRohringer et al. (2012); Wentzell et al. (2016) which are also used in this work.
(ii) Self-energy feedback: In many applications of the fRG the self-energy and its feedback on the flow of the -particle () vertex functions has not been accounted for. While there are arguments that the self-energy may be important mainly when the interactions are close to a flow to strong coupling (see Appendix in Ref. [Honerkamp et al., 2001]), more quantitative results should overcome this deficit. In fact, neglecting the self-energy feedback was mainly motivated by the disregarded frequency dependence of the interactions in earlier fRG studies: Within a static treatment the self-energy lacks the effects of quasiparticle degradation, so that its inclusion became less important. Within the current frequency-dependent fRG treatments, the self-energy feedback can be included in a meaningful way. A number of works have already investigated the self-energy effects in the flows to strong coupling in Hubbard-type modelsZanchi, D. (2001); Honerkamp (2001); Honerkamp and Salmhofer (2003); Rohe and Metzner (2005); Katanin et al. (2005); Husemann et al. (2012); Giering and Salmhofer (2012); Uebelacker and Honerkamp (2012b); Eberlein (2015); Vilardi et al. (2017), mainly exploring the quantitative effects, besides signatures of pseudogap openingsKatanin et al. (2005); Rohe and Metzner (2005) and non-Fermi liquid behaviorGiering and Salmhofer (2012) in particular cases.
(iii) Truncation of the flow equation hierarchy: Finally, one should also consider the truncation of the hierarchy of flow equations for the -point one-particle irreducible (PI) vertex functions. This is usually done at “level-II” as defined in Ref. [Metzner et al., 2012], also referred to as one-loop () approximation, i.e. the 1PI six-point vertex is set to zero. Due to this truncation, the final result of an fRG flow might depend - to a certain degree - on the cutoff scheme adopted for the calculation.
In this perspective, it was noticed by KataninKatanin (2004) that replacing the so-called single-scale propagator in the loops on the r.h.s. of the flow equation for the four-point vertex by a scale-derivative of the full Green’s function allows this scheme to become equivalent to one-particle self-consistent (a.k.a. mean-field) theories in reduced models, and then to go beyond such self-consistent approximations in more general models. Another significant comparison can be made with the parquet-based approachesBickers (2004); Janiš (1999), such as the parquet approximation (PA)Yang et al. (2009); Tam et al. (2013); Valli et al. (2015); Li et al. (2016); Wentzell et al. (2016). The latter represents the “lowest order” solution of the parquet equations, where the two-particle irreducible vertex is approximated by the bare interaction. In fact, although the diagrams summed in the truncation of the fRG are topologically the same as in the PA, the way the single contributions are generated during the flow leads in general to differences with respect to the PATaranto (2014); Wentzell et al. (2016). This is due to some internal-line combinations, e.g. in particle-hole corrections to the particle-particle channel, which are suppressed by the cutoff functions attached to the propagators and not fully reconstructed during the flow because of the truncation. A quantitative analysis of this effect has been performed for the single impurity Anderson model in Ref. [Wentzell et al., 2016]. These differences are absent for single-channel summations (e.g. RPA), but could lead to more pronounced quantitative errors in presence of channel coupling, e.g. in the generation of superconducting pairing through spin fluctuations. Furthermore, while the Mermin-Wagner theorem is fulfilled within the PABickers and Scalapino (1992), it is typically violated by fRG calculations. First steps to remedy this shortcoming were undertaken in various worksKatanin (2009); Maier and Honerkamp (2012); Eberlein (2015), but only recently a comprehensive path of how the PA contributions can be recovered in full extent was presented within the multiloop extension of the fRG (mfRG)Kugler and von Delft (2018a); Kugler and von Delft (2018b). The mfRG flow equations incorporate all contributions of the six-point vertex that complement the derivative of diagrams already part of the flow, as organized by their loop structure. A key insight in this approach is that the higher-loop contributions can be generated by computing 1 flows for scale-differentiated vertices, with an effort growing only linearly with the loop order that is fully kept. The multiloop corrections stabilize the flow by enabling full screening of competing two-particle channels, ultimately recovering the self-consistent structure of the PA. As the PA corresponds to a well-defined subset of diagrams, a converged mfRG flow able to reproduce the PA is by construction independent of the adopted cutoff.
In this paper, we present a fRG study of the 2D Hubbard model performed with an algorithm combining the most recent progress on all three approximation levels. We use (i) the so-called “truncated unity” fRGLichtenstein et al. (2017) (TUfRG) formalism to describe the momentum dependence of the vertex and, in addition, keep the full frequency dependence as a function of three independent frequencies. Differently from the approach adopted in Ref. [Vilardi et al., 2017], we employ a refined scheme to treat the high-frequency asymptoticsWentzell et al. (2016) that allows us to reduce the numerical effort considerably. Within this scheme, we can consistently include (ii) the (frequency-dependent) self-energy feedback in our fRG flow equations. Finally, we present (iii) first data for the 2D Hubbard model computed with the multiloop extension proposed by Kugler and von DelftKugler and von Delft (2018a). In this context, we have also generalized the multiloop formalism to compute the flow of the response functions, and illustrated the loop convergence of the fRG results for the 2D Hubbard model. In particular, we show that including up to loops in the fRG flow yields a clear convergence of the data with the loop order and the final results are independent of the cutoff. This represents an important check and illustrates that fRG flows can be brought in quantitative control for 2D problems. Finally, our multiloop analysis of the response functions demonstrates that the two different ways to compute susceptibilities in the fRG, either by tracking the renormalization group flow of the couplings to external fieldsMetzner et al. (2012) or by contracting the final interaction vertex in the corresponding Bethe-Salpeter equation (see, e.g. Ref. Taranto et al., 2014), converge to the same value with increasing loop order. This confirms that the output of this improved fRG scheme can indeed be trusted on a quantitative level.
The paper is organized as follows: The formalism and theory of the linear response functions and their computation by mfRG flow equations are introduced in Section II. In Section III we present the actual implementation scheme for the full momentum- and frequency-dependent fRG. In Section IV we show the results for the 2D Hubbard model, with a detailed analysis of the effects of the different approximation levels and in particular of the convergence with the loop order. A conclusion and outlook is provided in Section V.
Ii Theory and formalism
ii.1 Definitions and formalism
In this section we provide the definitions of the linear response functions to an external field, before describing their computation with the fRG. We focus on correlation functions of fermionic bilinears. In particular, in a time-space translational-invariant system, we consider the charge (density) and spin (magnetic) bilinears, both charge invariant, \cref@addtoresetequationparentequation
and the non-charge invariant pairing (superconducting) bilinears \cref@addtoresetequationparentequation
where and represent the Grassman variables and () a fermionic (bosonic) quadri-momentum (). The integral includes a summation over the Matsubara frequencies (), normalized by the inverse temperature , and an integral over the first Brillouine Zone normalized by its volume . The function , determines the momentum and frequency structure of the bilinears in the different physical channels. In the present case we restrict ourselves to a static external source field, such that the function acquires only a momentum dependence, whose structure is specified by the subscript and explicitly shown in Table 1 (in the present work we will mostly focus on the - as well as -wave momentum structure). Note that by using a different frequency-momentum notation, centered in the center of mass of the scattering process (see “symmetrized” notation in Appendix A), one should account for an additional shift of the momentum dependence by means of the momentum transfer .
We can now define the correlation functions of these bilinears in the three channels \cref@addtoresetequationparentequation
with . In linear response theory, these correlation functions correspond to the physical susceptibilities in the corresponding channels. Divergences in indicate spontaneous ordering tendencies or instabilities of the system. The above definition encodes not only the real-space pattern or wavevector with which the system starts ordering, but also the symmetry of the order parameter associated to the instability. In the D Hubbard model study presented here (see Section IV) we detect various response functions growing considerably towards low , such as the spin-density wave (SDW) response, characterized by the isotropic -wave magnetic susceptibility at , as well as - and -wave pairing response functions at and Pomeranchuk instabilitiesHalboth and Metzner (2000a). Inserting Eqs. (1) or Eq. (2) into Eq. (3), the susceptibilities appear as two-particle Green’s functions. In particular, they can be determined from the two-particle vertex by \cref@addtoresetequationparentequation
where represent the Pauli matrices () and we made use of the spin conservation. Eqs. (4) can be considerably simplified by making use of the SU() symmetry. The “bare bubbles” appearing in Eq. (4), with , read: \cref@addtoresetequationparentequation
By exploiting the SU() symmetry,
we can drop the spin dependencies for the bare bubbles. In presence of the above symmetries, we can introduce the following definitions for (spin-independent) channels of the two-particle vertex \cref@addtoresetequationparentequation
with the Levi-Civita symbol. The resulting spin-independent expression of the physical susceptibilities reads
We conclude this section by recalling the definition of the so-called fermion-boson vertexAyral and Parcollet (2015), which, for the considered symmetries, reads \cref@addtoresetequationparentequation
where, because of the SU() symmetry, we dropped the spin dependence of the fermion-boson vertex \cref@addtoresetequationparentequation
ii.2 Flow equations for the response functions
In this section we derive the mfRGKugler and von Delft (2018a) flow equations of the response functions and discuss the improvement with respect to the versionMetzner et al. (2012). In the following we provide the main steps of the derivation in the 1PI formulationMetzner et al. (2012); Salmhofer and Honerkamp (2001) (see also Ref. [Halboth and Metzner, 2000a] for the Wick-ordered formulation), for the details we refer to Appendix B.
We note that, although appears as a functional dependence in our derivation, it is not an integration variable since our system is fully fermionic (for an fRG formulation of coupled fermion-boson systems, see [Metzner et al., 2012; Kopietz et al., 2010; Schütz et al., 2005; Kugler and von Delft, 2018c]).
By expanding the scale-dependent effective action in powers of the fermionic fields, as well as of the external bosonic source field, we obtain
Note that the index combines the spin index and the fermionic quadrivector (here we disregard additional quantum dependencies, e.g. orbital), while refers to the momentum structure of the coupling to the bilinears, , and to the bosonic quadrivector . In Eq. (13) the first term on the r.h.s. represents the expansion of the effective action in absence of external field (see Section III), while the functional derivatives in the following terms represent the -dependent susceptibility and the fermion-boson vertex in the different channels. Taking the derivative with respect to the scale parameter (see Appendix B), yields the following flow equations for the susceptibility and fermion-boson vertex (assuming SU(2) symmetry and momentum-frequency as well as spin conservation) \cref@addtoresetequationparentequation
and respectively \cref@addtoresetequationparentequation
represents the single-scale propagator. The function , differently from the (fermionic) two-particle vertex , represents a mixed bosonic-fermionic vertex, i.e. with two bosonic and two fermionic legs where we summed over its spin dependences
while the spin-independent form for used in Eqs. (15) reads \cref@addtoresetequationparentequation
The conventional approximationsHalboth and Metzner (2000a); Salmhofer and Honerkamp (2001); Metzner et al. (2012) disregard the first terms on the r.h.s. of Eqs. (14) and (15). This approximation is consistent with the corresponding approximation of (see Appendix C) and justified in the weak-coupling regime. Using the notation of Refs. [Kugler and von Delft, 2018a; Katanin, 2004], one can rewrite the approximation of Eqs. (14) and (15) in a more concise tensor-form \cref@addtoresetequationparentequation
So far we pinpointed two possible ways to compute the susceptibility and fermion-boson vertex from an fRG calculation: (i) Solving Eqs. (14) and (15) alongside the ones for and (at the same level of approximation), and (ii) by means of Eqs. (8) and (10) at the end of the fRG flow, using and , later referred to as ”post-processing”. These two procedures are non-equivalent in the presence of approximations, e.g., if one restricts oneself to the level. This leads to an ambiguity in practical implementations of the fRG. In fact, as shown in Appendix D, the two results deviate at for the case (for a larger number of loops the deviations occur at higher orders in the effective interaction ). In order to solve this ambiguity we note that the exact relations (8) and (10) are fulfilled in the PA. At the same time, the recently introduced multiloop extension allows to sum up all parquet diagrams. Hence, generalizing the multiloop flow to the computation of the response functions recovers the equivalence of the two procedures.
In order to derive the mfRG equations for the response functions, we first recall the channel-decomposition of the two-particle vertex as known from the parquet formalism. The latter divides in the two-particle reducible vertex (all diagrams that can be divided into two separate ones by removing two internal fermionic propagators) and the two-particle irreducible vertex (which can be not be divided). Depending on the direction of the propagation lines the diagrams are reducible onto either parallel, longitudinal antiparallel or transverse antiparallel, corresponding to the particle-particle, particle-hole, and particle-hole crossed channel, respectively. Besides this diagrammatic channel decomposition, there is also a distinct physical channel decomposition that identifies the components and which we will use in the following. Inserting this decomposition into the flow equation for the two-particle vertex, we obtain
While the usual diagrammatic channel decompositionBickers (1991) leads to simple expressions for the two-particle irreducible vertex , the latter assumes a more complicated form in the physical channels \cref@addtoresetequationparentequation
where we approximated the fully two-particle irreducible vertex by its first-order contribution in the interaction , which is known as PA.
We now derive the mfRG flow equations for the response functions, which mimics the effect of the mixed fermion-boson vertices and in the exact flow equations (14) and (15). First, one performs the so-called Katanin substitution Katanin (2004) () in the flow equations (19b), (19a). One observes that all differentiated lines in these flow equations come from . Secondly, differentiated lines from the other channels are contained in the higher-loop terms of the expansion \cref@addtoresetequationparentequation
Using the channel decomposition (21), we can immediately write down the correction to the flow of the fermion-boson vertex. It accounts for the leading-order diagrams in the effective interaction coming from in Eq. (15) (see Appendix B) and reads
On the three- and higher-loop level, we can use in an analogous fashion. Furthermore, one needs to account for vertex corrections to the right of the differentiated lines, such that
Considering the flow equation of the susceptibility (19a), we see that the fermion-boson vertices provide the vertex corrections on both sides of the differentiated lines in . Hence, for all higher-loop corrections we can simply connect to both fermion-boson vertices, thereby raising the loop order by two. We obtain as well as
In total, these equations together with the multiloop flow of the fermionic two-particle vertex (see Section III.2) allow us to sum up all differentiated parquet diagrams of and , so that the aforementioned two ways of computing the response functions within the fRG framework become equivalent. In order to provide a consistent fRG scheme, it is important to adopt the same level of approximation (truncating the sums in Eq. (23a) to a certain finite -loop level) for all flowing quantities.
Iii Numerical implementation
iii.1 Full frequency and momentum parametrization
In order to illustrate the fRG algorithm adopted in the present work, let us start from the flow equations for the one particle irreducible (PI) fermionic vertex in the fRG approximation. In the following, the SU(), spin conserving symmetry will be always assumed. Exploiting this symmetry, the self-energy and two-particle fermionic vertices can be written as
where the fourth argument of is determined by in a momentum and energy conserving system. The spin-independent flow equation for the self-energy reads
where the diagrammatic channel index distinguishes between particle-particle, particle-hole and particle-hole exchange diagrams and the first dependence of the functions refers to the bosonic four-momentum transfer in the internal loop of their corresponding equations \cref@addtoresetequationparentequation
Each of the above equations depends, besides the aforementioned bosonic transfer dependence (, and ), on two fermionic dependencies. Such mixed ‘bosonic-fermonic’ notation, referred to as ‘non-symmetrized’ notation, has been substituted in some work (e.g. in Ref. [Lichtenstein et al., 2017]) by a different notation where the dependencies of the four fermionic propagators involved in the scattering process have been chosen symmetrically with respect to the bosonic four-momentum transfer. This symmetrized notation simplifies the implementation of the symmetries exploited in the fRG code (see Appendix F and Ref. [Lichtenstein et al., 2017]) but leads to less compact flow equations. The equations 31c, 31 and 31a generate the two-particle reducible vertices () of the diagrammatic parquet decomposition
The two-particle fermionic vertex can be reconstructed by using Eq. (32). The use of a mixed ‘bosonic-fermonic’ notation allows us to identify the transfer bosonic four-momentum as the strongest dependence, while the two fermionic dependencies can be treated with controllable approximations. In the following we illustrated two efficient ways to simplify the treatment of both momentum and frequency dependencies.
iii.1.1 Truncated Unity fRG
The approximation for the fermionic momentum dependencies in TUfRGLichtenstein et al. (2017) is done by the expansion of the fermionic momentum dependencies in form factors, illustrated here for the channel
while the expansion of the and leads analogously to the objects and . As proposed in other papers Husemann and Salmhofer (2009); Husemann et al. (2012); Giering and Salmhofer (2012); Maier et al. (2013); de la Peña et al. (2016); Lichtenstein et al. (2017), we choose the form factors such that they correspond to a specific shell of neighbors in the real space lattice. The unity
which is inserted in the flow equations, contains a complete basis set of form factors. As observed also in previous works (Lichtenstein et al., 2017), converged results can be obtained already with a small set of form factors. For a fast convergence it is convenient to include one shell after another, starting from the local constant form factor and increasing the distance of neighbors taken into account. The form factors in this paper are listed in Table 1.
A major difficulty in this approach is the feedback of the different channels into each other. In addition to the dressing of the objects by the form factors, the translation of the notation in momentum and frequency from one to another channel has to be considered. Computationally time consuming integrations in momentum space can be avoided by Fourier transformation and evaluation in real space. Furthermore the expression of the projection in terms of a matrix multiplication allows for the precalculation of the projection matrices which can be found in the Appendix F.
iii.1.2 Dynamical fRG
In frequency space, we adopt the simplifications proposed in Refs. [Wentzell et al., 2016; Li et al., 2016]. For all systems with an instantaneous microscopic interaction one can use diagrammatic arguments to proof that, in the high frequency regime, the two-particle fermionic vertex exhibits a simplified asymptotic structure. In this region one can reduce the three frequencies dependence of using functions with a simplified parametric dependence. It is straightforward that, sending all three frequencies to infinity, reduces to the instantaneous microscopic interaction, which in the present case is represented by the Hubbard on-site interaction . The contribution of the reducible vertices to becomes non-negligible if the bosonic frequency transfer is kept finite, while sending the two secondary fermionic frequencies to infinity. This contribution, depending on a single bosonic frequency transfer on a given channel , has been denoted . For models with a instantaneous and local microscopic interaction, one observes that the momentum dependencies disappear alongside the frequency dependencies when performing such limits. In the limit where just one fermionic frequency is sent to infinity, the vertex can be parametrized by the function . By subtracting the asymptotic functions to the full object we obtain the so-calledWentzell et al. (2016),“rest-function” which decays to zero within a small frequency window. The parametrization given by allows us to reduced the numerical cost of computing and storing the two-particle fermionic vertices. In fact, for any of the three channels, we calculate the fRG flow of the three-frequency dependent function on a small low-frequency region and add the information on the high frequencies by computing the flow of functions and which are numerically less demanding. The full two-particle reducible vertex can be recovered by
where can be obtained from by exploiting the time reversal symmetries (see Appendix A.3).
iii.1.3 Flow equations for TU-dynamical fRG
Finally, applying the aforementioned projection on the form-factor basis one can write matrix-like approximated fRG flow equations for the self-energy, the two-particle fermionic vertex, the fermion-boson vertex and the susceptibility: \cref@addtoresetequationparentequation