{'fraction_non_alphanumeric': 0.05560500304492434, 'fraction_numerical': 0.013561624584250715, 'mean_word_length': 4.40101201771031, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Linear response within the projection-based renormalization method: Many-body corrections beyond the random phase approximation

The explicit evaluation of linear response coefficients for interacting many-particle systems still poses a considerable challenge to theoreticians. In this work we use a novel many-particle renormalization technique, the so-called projector-based renormalization method, to show how such coefficients can systematically be evaluated. To demonstrate the prospects and power of our approach we consider the dynamical wave-vector dependent spin susceptibility of the two-dimensional Hubbard model and also determine the subsequent magnetic phase diagram close to half-filling. We show that the superior treatment of (Coulomb) correlation and fluctuation effects within the projector-based renormalization method significantly improves the standard random phase approximation results.


Introduction
The most popular approach to evaluate linear response coefficients for many-body systems is probably the standard random phase approximation (RPA), which is used in quite different fields in physics. Besides its physical merits including the fulfillment of conservation laws the popularity of the RPA results from is conceptual simplicity in its derivation as well as from its numerical practicability. Attempts to go beyond the RPA have turned out to be extremely demanding. One of the options to improve the RPA are methods based on conserving approximations (see, e.g., Ref. [1]), thereby following an approach which was introduced by Baym and Kadanoff [2,3]. Conserving approximations are consistent with microscopic conservation laws for particle number, energy or momentum. Other work is based on the time-dependent Hartree-Fock approximation which uses a frequency-dependent local field factor in a modified RPA expression [4]. However, this method turns out to be rather complex and physically unsatisfactory. Attempts to find a numerical solution of the basic integral equation could not be reached without further approximations (see discussion in Ref. [4] and references therein). Recently, Vilk and Tremblay extended the RPA by including vertex corrections, taken into account correlation and exchange effects [5]. Com-paring their ansatz for double occupancies in the Hubbard model, the authors found good quantitative agreement with results from Monte Carlo simulations for singleparticle and two-particle properties [6,7]. Another way to improve the RPA is the so-called self-consistent RPA [8], which is based on a non-perturbative variational scheme. This approach has been adopted to the investigation of various nontrivial models but is however limited to small systems [8,9,10].
For these reasons it is important to develop new manyparticle techniques having the ability to include correlation effects. One approach that overcomes some of the shortcomings of the RPA is the projection-based renormalization method (PRM) [11,12]. In the recent past the PRM has been successfully applied to several physical problems such as superconductivity [13], quantum phase transitions in coupled electron-phonon systems [14,15,16,17], exciton and plasmaron formation [16,18], BCS-BEC transition [19], electronic phase separation [20], valence transitions [21], or the Kondo lattice problem [22]. In the present work, adding time-and wave-vector-dependent external fields, we demonstrate how the PRM can be combined with linear response theory in order to calculate response functions for generic correlation models. In particular we derive an explicit analytical expression for the dynamical spin susceptibility of the Hubbard model. For the two-dimensional (2D) case the PRM phase boundaries between the paramagnetic and antiferromagnetic respectively ferromagnetic phases are determined for weak-tointermediate Hubbard interactions. An elaborate weakcoupling approach is of particular importance in low spatial dimensions since in 1D and 2D also weakly interacting systems tend to be strongly correlated.
The paper is organized as follows. In the Sect. 2, we recapitulate the RPA to the Hubbard model. Section 3 introduces the PRM approach, which is applied to the Hubbard model in Sect. 4, focusing on the response to an external magnetic field. Thereby the renormalization equations for the model parameters, the transformations of the operators and various expectations values are derived. Details can be found in Appendices A and B. Section 4.4 provides our main analytical result: the explicit expression for the dynamical spin susceptibility. Selected numerical results for the 2D Hubbard model can be found in Sec. 5, in particular the ground-state phase diagram in the U -n plane and the wave-vector-and frequency-dependence of the magnetic susceptibility. We conclude in Sect. 6.

Standard RPA approach to the Hubbard model
The Hubbard model is a paradigmatic model for the study of correlation effects in itinerant electron systems. Independently proposed by Gutzwiller [23], Hubbard [24], and Kanamori [25] in 1963, it was originally designed to describe the ferromagnetism of transition metals. Successively, the model has been studied in the context of antiferromagnetism, metal-insulator transition, and high temperature superconductivity. The Hubbard Hamiltonian is given by Here c † iσ (c iσ ) is a fermionic creation (annihilation) operator of a spin σ (=↑, ↓) electron, and n iσ = c † iσ c iσ . U denotes the on-site Coulomb interaction andt are the electron transfer matrix elements between nearest-neighbor Wannier sites i and j. The physics of the model is governed by the competition between itinerancy (t; delocalization, kinetic energy) and short-range Coulomb repulsion (U ; localization, magnetic order), where the fermionic nature of the charge carriers is of great importance (Pauli exclusion principle). Besides the parameter ratio U/t, the particle density n, the temperature T , and the spatial dimension D (geometry of the lattice) are crucial.
Although a tremendous amount of work has been devoted to the Hubbard model, in order to determine its ground-state, spectral and thermodynamic properties, exact results are rare and only a few special cases and limits are ultimately understood. In 1D, the algebraic and thermodynamic Bethe ansatz enables an exact treatment of the model [26,27]. However the Bethe ansatz technique does not provide a complete framework since it generally does not allow the evaluation of the response functions. For D > 1 approximations are unavoidable anyway. There usually the weak-(U/W ≪ 1) and strong-coupling (U/W ≫ 1) limits of the model were studied, with uncertain extrapolations to the region U/W ∼ 1. Here W is the bare electronic bandwidth. For a D-dimensional hypercubic lattice we have W = 4Dt.
In consideration of the magnetic behavior of a Hubbard model system the response to an applied external field is of particular importance. Adding a small magnetic field that periodically oscillates in space and time the Hamiltonian takes the form where is the kinetic energy of electrons in momentum space. In the numerical evaluation below the Hubbard model is considered on a square, where the dispersion is given by with chemical potential µ. The last two terms in Eq. (2) respectively read and Note that the wave vector q is imposed by the external field. s z q is the component of the spin operator in field direction: Then the linear response of the spin expectation value s z −q (t) with respect to h(t) ∼ Re e −iωt is given by where where the expectation value · · · (t) on the right-hand side is formed with H(t) and therefore becomes time dependent. It can be simplified by help of the operator identity Here p is either p = ±q or p = 0. We have also introduced the k-resolved spin operator The expectation value c † kσ c kσ has no linear contribution in h(t) and can be considered as time independent.
Using Eqs. (10) and (11), the Hamiltonian H(t) can be rewritten as where In Eq. (15) we have introduced an effective fieldĥ q (t) that contains an internal field proportional to U . Also the kinetic energy H 0 has acquired a Hartree shift proportional to U . Finally, the part H f (t) reads where the t-dependence enters via the fluctuation operators.
The standard RPA expression for χ(q, ω) is obtained by neglecting the fluctuating part H f (t) completely, i.e., H(t) reduces to The linear response of s z −q (t) to the effective fieldĥ q (t) becomes where the expectation value is now formed with the unperturbed Hamiltonian H 0 . By help of relation (15) one arrives at Here χ 0 (q, ω) is the dynamical susceptibility of the unperturbed system H 0 with η = 0 + . Note that Eqs. (18)- (20) are the usual RPA equations.Thus, the dynamical RPA susceptibility is defined by the prefactor in Eq. (19):

PRM formalism
Our aim is to evaluate the dynamical susceptibility χ(q, ω) beyond the standard RPA by including fluctuation processes, which are induced by the fluctuation part H f (t) of the Coulomb interaction. To this end, we combine linear response theory with the PRM. The PRM in the original version without time-dependent external field starts by separating a given many-particle Hamiltonian into an unperturbed part H 0 and a time-independent perturbation H f . Since H f and H 0 do not commute, the perturbation induces transitions between the eigenstates of H 0 . The basic idea of the PRM is to eliminate successively all transitions due to H f so that finally only the unperturbed, yet renormalized Hamiltonian (now calledH 0 ) remains.
In the present case the many-particle Hamiltonian (13) is time dependent, is time-dependent due to the external field. As before, our aim is to evaluate the response of the expectation value s −q (t) up to linear order in the external field. However, the fluctuation term H f (t) of the Coulomb interaction should now be taken into account. Since both H 0 and H h (t) do not commute with H f (t), the latter Hamiltonian will henceforth be considered as perturbation. In particular, H f (t) again induces transitions between the eigenstates of H 0 . In the PRM these transitions will be eliminated by a sequence of unitary transformations, which are performed in small steps ∆λ by proceeding from large to small transition energies. Let H λ be the Hamiltonian after all transitions with energies larger than some cutoff λ have already been integrated out. The transformation from cutoff λ to a somewhat reduced cutoff λ − ∆λ formally reads Here X λ,∆λ = −X † λ,∆λ is the generator of the unitary transformation from λ to λ − ∆λ, whereas H λ (t), represents the renormalized Hamiltonian after all transitions (in the eigenbasis of H 0,λ ) with energies larger than λ have been eliminated from H f (t). Similarly, H λ−∆λ (t) denotes the Hamiltonian with the somewhat reduced cutoff λ − ∆λ. Due to transformation (22) the parameters of H λ become renormalized, and also new terms can in principle be generated. For the generator X λ,∆λ we chose which agrees with the lowest order result for X λ,∆λ in the absence of an external field [11]. Here L 0,λ is the Liouville operator of the 'unperturbed' Hamiltonian H 0,λ . It is defined by the commutator L 0,λ A = [H 0,λ , A], applied to any operator variable A. Moreover, Q λ−∆λ is a generalized projection operator. It projects on all transition operators in H f,λ (in the basis of H 0,λ ) with transition energies larger than λ − ∆λ. Note that X λ,∆λ also depends on time t since H f,λ (t) depends on t via the external field. The elimination procedure starts from the original Hamiltonain H(t) (where the largest cutoff energy is called λ = Λ) and proceeds in steps ∆λ until λ = 0 is reached. This limit provides the desired effective Hamiltoniañ H(t) := H λ→0 (t) with H λ→0 (t) = H 0,λ→0 + H h,λ→0 (t). Note that the elimination of the transitions leads to renormalized parameters inH(t). Thus, after all transitions from H f (t) have been used up, the final HamiltonianH(t) is diagonal or at least quasi-diagonal and allows to evaluate expectation values. As a matter of course the parameters inH(t) depend on their values in the original model H(t).
Having in mind small renormalization steps ∆λ, the transformation (22) from λ to λ − ∆λ can be restricted to an expansion up to second order in U (and linear order in h(t)). Then H λ−∆λ (t) reads where relation (24) was used and the explicit t-dependence is suppressed. P λ−∆λ = 1 − Q λ−∆λ is the projector on all low-energy transitions with energies smaller than λ − ∆λ.
The commutators in Eq. (25) give rise to renormalization contributions to H λ−∆λ (t). Having in mind an application of linear response theory, Eq. (25) has to be evaluated up to linear order in the external field. Finally, comparing the result of the evaluated right-hand side with the generic form of H λ one is led to renormalization equations, which relate the parameters of the Hamiltonian at cutoff λ with those at cutoff λ − ∆λ.
As is discussed below one also has to evaluate expectation values, which are formed with Hamiltonian H. Exploiting the unitary invariance of operators below a trace we can write where the same unitary transformation as before is applied to operator A, i.e. A(λ) = e X λ Ae −X λ . Here X λ is generator of the unitary transformation between cutoff Λ and λ, andÃ = A(λ → 0). Thus, additional renormalization equations for A(λ) are required.
Let us mention that Wegner and coworkers [28,29] have introduced a theoretical approach related to the PRM. This approach is based on the application of continuous unitary transformations instead of discrete ones as in the present case. To our knowledge it was not applied up to now to the investigation of many-body corrections beyond the random phase approximation discussed in the present work. However, correlation and fluctuation processes can be discussed in this method as well, compare for instance references [30,31,32], or [33]. The relationship between the continuous method and the PRM for the case without time-dependent field is studied in references [19] and [34]. There it was shown that the continuous method can be derived within the PRM framework in the limit of small ∆λ → 0 using a particular choice for the complement part P λ−∆λ X λ,∆λ of generator (24).

PRM for the Hubbard model
We are now in the position to apply the general formalism of Section 3 to the Hubbard model. Thereby, the influence of the fluctuation term H f (t) will be investigated. Following the ideas of the PRM, we have to start from an ansatz for the renormalized Hamiltonian H λ (t). A perturbative evaluation of transformation (22) suggests the use of the Due to the renormalization all coefficients in H 0,λ and H h,λ (t) depend on λ. Moreover, inĤ h,λ (t) new operator contributions are generated. The coefficients u k,k−q,λ (t) and v k,k+p−q;k ′ ,k ′ −p,λ (t) are expected to depend linearly on the external field and are therefore explicitly time-dependent. From hermiticity follows that they obey the relations: Finally P λ in H f,λ (t) projects on the low-energy excitations smaller than λ. One finds which shows that the operators on the right-hand side have different transition energies. Here, we have defined two Θ-functions: They guarantee that only transitions with excitation energies smaller than λ are kept in H f,λ (t). Finally we use relation (10) to regroup H f,λ (t): with , H α f,λ (t) only depends on fluctuation operators, whereas H β f,λ (t) is linear to the external field via the expectation value s z k ′ ,−q . Here we have defined Note that also higher-order fluctuation contributions both to H f,λ (t) and X λ,∆λ (t) could be considered. Their inclusion would extend the range of validity of the present approach, which is restricted to the range from small to intermediate coupling U/W 1, to larger values. However, they would further complicate the evaluation of transformation (25) and will be neglected. Similarly, from Eq. (24) one finds the following expression for the generator X λ,∆λ : The coefficients in Eqs. (40), (41) are defined by In the products of Θ-functions assure that only excitations between λ and λ−∆λ are eliminated by the unitary transformation (22).

Renormalization equations
Integrating out all transitions induced by H f,λ (t), the parameters of H λ (t) will be renormalized. Only the Coulomb coupling U remains λ-independent (apart from the Θfunctions in Eq. (33)). The λ-dependence of the parameters will be derived with transformation (25) for an additional step from λ to λ − ∆λ. The result of the explicit evaluation has to be compared with the generic expression for H λ−∆λ (t), which is obtained by replacing λ in H λ [Eqs. (27)- (33)]. In this way one obtains the desired renormalization equations, which connect the λ-dependent parameters of H λ with those at cutoff λ − ∆λ. According to Appendix B we find The renormalization contributions on the right-hand sides of these equations are of order U and U 2 . They are given in Appendix B: δε k,λ by Eqs. (104) , and (117). In order to reduce the operator structure of H λ−∆λ (t) to operators which appear in H λ (t) an additional factorization of higher operator terms has been performed. Therefore, the expectation values c † kσ c kσ and c † kσ c k±q,σ enter the renormalization contributions.
The renormalization equations have to be solved numerically, starting from the initial parameters of the original model H(t), i.e., Suppose, the expectation values on the right-hand side of Eqs. (47)-(49) are known, the renormalization procedure from Λ to λ = 0 leads to the fully renormalized Hamilto- The tilde symbols denote the fully renormalized quantities at λ = 0. All excitations from H f,λ (t) have been eliminated, leading to the renormalization of H 0,λ and H h,λ (t). The final HamiltonianH(t) describes a system of free renormalized conduction electrons in a renormalized effective field. Thereby, the quantitiesũ k,k−q (t) and v k,k+p−q;k ′ ,k ′ −p (t) depend linearly on the external field and are time-dependent. This follows from the renormalization equations (47)-(49), using the expressions for the renormalization contributions from Appendix A, together with the initial conditions (50). Therefore, relying on linear response theory with respect to the effective field, any expectation value can be evaluated.

Occupation numbers n kσ
The yet unknown expectation values on the right-hand side of the renormalization equations (47)-(49) can be evaluated self-consistently as follows. Let us first consider the averaged occupation number n kσ for fixed spin σ. Using Eq. (26), n kσ can be rewritten as In principle, the last expectation value has to be formed with the time-dependent HamiltonianH(t). However, restricting ourselves to first order in h(t),H(t) can be replaced by the time-independent HamiltonianH 0 .c † kσ is the fully renormalized creation operatorc † where the operator structure of (56) is again taken over from the lowest order expansion in X λ of the unitary transformation.
In analogy to Eq. (22) renormalization equations for the coefficients x k,λ and y kpk ′ ,λ can be derived by evaluating the renormalization step from λ to λ − ∆λ. One finds This equation connects y kpk ′ ,λ−∆λ at cutoff λ − ∆λ with the coefficients x k,λ , y kpk ′ ,λ at cutoff λ. Similarly, also an renormalization equation for x k,λ can be found. Alternatively one may start from the anti-commutator relation of c † kσ (λ) and c kσ (λ). Taking the expectation value with H λ , we get or Here, a factorization approximation for was used. Eq. (59) connects the coefficients x k,λ and y kpk ′ ,λ for any value of λ. The equation for x k,λ−∆λ is found from the sum rule (59), when λ is replaced by λ − ∆λ, i.e., Together with Eq. (57) this equation relates x k,λ−∆λ with x k,λ and y kpk ′ ,λ . Thus Eqs. (61), (57) connect the parameter values at λ with those at λ−∆λ. Integrating Eqs. (57), (61) between Λ and λ = 0 (thereby using x k,Λ = 1 and y kpk ′ ,Λ = 0), we obtaiñ where the tildes again denote the fully renormalized quantities. Thus, for n kσ the final result is which is independent of σ in the paramagnetic state. f (ε k ) is the Fermi function.

Transformation of spin operators
To evaluate the dynamical spin susceptibility we need the transformed spin operator,s z q = s z q,λ→0 , where s z q,λ = e X λ s z q e −X λ . For the λ-dependence we use the following ansatz, which corresponds to ansatz (29) forĤ h,λ . According to Appendix A we write Their initial values at λ = Λ are The renormalization equations for the coefficients are found from the transformation step s z q,λ−∆λ = e X λ,∆λ s z q,λ e −X λ,∆λ , where in linear response theory X λ,∆λ can be replaced by its part X (α) λ,∆λ . A closer inspection shows that the renormalization equations can be taken over from the equations for u k,k−q,λ and v k,k+p−q,k ′ ,k ′ −p,λ . We get where and δβ (1) k,k+p−q,k ′ ,k ′ −p,λ = (68) Note that in Eq. (66) no equivalent to δu (2) k,k−q enters, since the latter contribution was caused by the commutator [X λ,∆λ , H f,λ ] (compare Appendix A). Furthermore, renormalizations from the second part in ansatz (64), being proportional to ∼ β k,k+p−q,k ′ ,k ′ −p,λ , have been neglected.
Let us finally stress that the renormalization contributions to α k,k−q,λ and β k,k+p−q,k ′ ,k ′ −pλ vanish for wave vector q = 0. This is in accord with the rotational invariance of the system in spin space.

Dynamical magnetic susceptibility
The dynamical q-and ω-dependent magnetic susceptibility χ(q, ω) is defined by the linear response of the averaged spin s z −q (t) to the small external field h(t). SinceĤ h,λ=0 itself is proportional to h(t): Here we have again used the unitary invariance (26) of operator expressions under a trace. The expectation value in Eq. (69) is formed with the renormalized one-particle andĤ h,λ=0 (t) are the fully renormalized quantities. They are given bỹ Thus, Eq. (69) reduces to where we have introduced k-resolved susceptibilities and with Note that in both susceptibilities (73) and (74) renormalized energiesε k enter and not the unrenormalized energies ε k as in Eq. (20). According to Sect. 4.2 and Appendix B the quantities u k,k−q (t) andṽ k,k+p−q;k ′ ,k ′ −p (t) depend linearly on the external field h(t) via the k-resolved spin operator expectation values s z k,−q (t). In order to simplify the further calculation we shall trace back these quantities to the expectation value s z −q (t) of the full spin operator. This is done by assuming that the coefficients B Thus, from u k,k−q,λ (t) as well as from v k,k+p−q;k ′ ,k ′ −p,λ (t) a common factor s z −q (t)/N can be extracted: where we introduced time-independent quantities u 0 k,k−q,λ and v 0 k,k+p−q;k ′ ,k ′ −p,λ . They obey the following renormalization equations: with the time-independent renormalization contributions given in Appendix C. Due to Eq. (50) the initial conditions for u 0 k,k−q,λ and v 0 k,k+p−q;k ′ ,k ′ −p,λ at cutoff Λ are Having Eq. (80), (81) we are in a position to rewrite relation (72) in a time-independent form: where, on the right-hand side, we have used the relation The coefficientsũ 0 k,k−q andṽ 0 k,k+p−q;k ′ ,k ′ −p are again fully renormalized. They are found by solving the renormalization equations (80), (81) due to the initial conditions (82). Eq. (83) is an implicit equation for the dynamical susceptibility χ(q, ω) which also enters the quantities δu 0(1) k,k−q,λ and δv 0(1) k,k+p−q,k ′ ,k ′ −p,λ (see Appendix C). Solving for χ(q, ω) we obtain our final analytical result: with The PRM result (85) for the dynamical magnetic susceptibility of the Hubbard model represents an extension of the standard RPA expression. Besides the new coefficientsα k,k−q in the numerator and denominator an extra term ∆(q, ω) occurs in the denominator, which generalizes the overall shape of an RPA expression. Noteworthy the quantitiesũ 0 k,k−q andṽ 0 k,k+p−q;k ′ ,k ′ −p in ∆(q, ω) themselves depend on χ(q, ω) (Appendix C). Thus, χ(q, ω) has to be solved self-consistently from Eqs. (85), (86). Expression (85) reduces to the standard RPA results when all renormalization effects are disregarded. Then,α k,k−q andβ k,k+p−q,k ′ ,k ′ −p keep their original values 1 and 0, whereasũ 0 k,k−q andṽ 0 k,k+p−q;k ′ ,k ′ −p remain zero and also ∆(q, ω) vanishes. Then χ(q, ω) becomes χ RP A (q, ω) given by Eq. (21) with χ 0 (q, ω) being the unrenormalized dynamical magnetic susceptibility of free electrons [Eq. (20)].
In the next section χ(q, ω) will be evaluated numerically. Before, let us study the special case q = 0. Since the total spin s z q=0 commutes with the total Hamiltonian, i.e. [H, s z q ] = 0, the dynamical susceptibility vanishes, which can immediately be seen from the general expression (9) for χ(q = 0, ω) (ω finite). The same conclusion can also be drawn from the PRM formalism. Since the total spin s z q=0 commutes with the Hamiltonian it also commutes with the generator X λ,∆λ . Therefore, the coefficients in representation (70) fors z q,λ will not be renormalized for q = 0. Thus, we haveα k,k = 1 andβ k,k+p,k ′ ,k ′ −p = 0. Similarly, the coefficientsũ 0 k,k−q andṽ 0 k,k+p−q;k ′ ,k ′ −p in expression (71) forĤ h,λ=0 (t) also vanish for q = 0. Thus, the quantity ∆(q, ω) in the denominator of Eq. (85) vanishes and the susceptibility χ(q = 0, ω) takes the standard RPA form Here, according to Eq. (73), the susceptibility which contains the renormalized energiesε k and not the unrenormalized energies ε k . From Eqs. (87) and (88) one immediately concludes that for q = 0 the real and imaginary parts of χ(q, ω) vanish (for any finite ω). Compare also Section 5.3 below. On the other hand, when the limit ω → 0 is taken first, the imaginary part of χ(q, ω) vanishes for any q. This follows from the analytical properties of Imχ(q, ω) or from Eqs. (85) and (88). In contrast, the real part Reχ(q, ω = 0) stays finite and reduces to the static q-dependent susceptibility In the limit q → 0, χ(q) gives the uniform static susceptibility. Last but not least, keep in mind that the renormalization equations derived so far, exclusively apply to the paramagnetic phase.

Numerical Results and Discussion
The set of self-consistency Eqs. (47)-(52), (85), and (86) has to be solved numerically in momentum space. Due to additional internal k-sums in the renormalization contributions, we restrict ourselves to a square lattice with N = 24 × 24 sites using periodic boundary conditions. In contrast, for the non-interacting susceptibility χ 0 k (q, ω) a larger mesh in momentum space of 2000 × 2000 points is used close to the Fermi surface due to the large variations in this quantity. Choosing reasonable initial values for the various expectation values, the renormalization starts from the cutoff Λ of the original model and proceeds in energy steps ∆λ = 0.5t until λ = 0 is reached. Then the expectation values are recalculated. Convergence is assumed to be achieved if all quantities are determined with a relative error less than 10 −5 . We have convinced ourselves that a larger lattice size as well as a smaller value of ∆λ will not modify the presented results. In what follows we measure all energies in units oft.

Band renormalization
Correlation effects, which are included in the PRM scheme, lead to a momentum-dependent renormalization ε k →ε k of the band structure in the paramagnetic phase. This differs from standard Hartree-Fock [35,36], Gutzwiller [23,37], or slave-boson treatments [38,39], where band renormalization either not at all or in terms of a momentum independent band narrowing takes place. To illustrate the PRM renormalization of the quasiparticle band, in Fig. 1 the difference δε k = ε k −ε k is shown for a square lattice Brillouin zone for the particle density n = 0.9 and three different U values. The overall bandwidth is reduced by a factor of 0.22, 0.34 and 0.45 for U = 2, 3 and 4, respectively. Beyond that, we find that the momentum dependence of δε k increases with U and is largest near the center k = (0, 0) and at the corners k = (π, π) of the Brillouin zone. This should have strong impact on the uniform and staggered (static) spin susceptibilities which are studied in the next subsection.

Temperature dependence
Let us first consider the half-filled band case n = 1 and track the behavior of the static susceptibility χ(q; T ) as the temperature T is lowered. The PRM results for χ(q; T ) are shown in Fig. 2 for two fixed wave vectors q = 0 (uniform susceptibility, panel (a)) and q = Q = (π, π) (staggered susceptibility, panel (b)). As one can see, the logarithmic singularity in the density of states at the band center, ρ(E) ∝ ln(ε/4t) for ε → 0, leads to a divergence of the noninteracting susceptibilities χ 0 (0; T ) ∝ − ln(T /t) and χ 0 (Q; T ) ∝ −[ln(T /t)] 2 . This indicates a magnetic instability of the corresponding PRM susceptibilities at some finite T for any U . Therefore, since the divergence of χ 0 (Q; T ) is stronger than that of χ 0 (0; T ), the PRM predicts a transition to a magnetic phase with strong antiferromagnetic fluctuations, which sets in at a higher temperature. The larger the U -values, the higher are the transition temperatures. For very small U the RPA and PRM results are nearly identical. As follows from the preceding section, the PRM renormalization of the uniform susceptibility χ(0; T ) is solely caused by the one-particle energies ε k (compare Eq. (87)), which are barely changed in this limit; cf. Fig. 1. On the other hand, the renormalization of χ(Q; T ) is affected by the coefficientsα k,k−q as well as by the second contribution ∆(Q) in the denominator. This term shifts the zero of the denominator in Eq. (85), with the result that antiferromagnetic transition temperature is reduced, as it should be if the correlations/fluctuations of the Hubbard system are treated better. Acceptably the antiferromagnetic critical temperature stays larger than the ferromagnetic one.
Of course, in 2D the occurrence of a finite transition temperature is an artefact of the approximations in the PRM, which is also known from the standard RPA. According to the Mermin-Wagner theorem, for a 2D model with continuous symmetry long-range order can only occur for T = 0 [40,36]. Indeed, from unbiased numerical approaches [36] it was shown that long-range (antiferromagnetic) order is expected at T = 0 only for half-filling. Clearly, the PRM in the present version does not overcome this shortcoming. However, as seen from Fig. 2, it gives the right tendency. One may expect that higher order fluctuation terms, not included at present in ansatz (27)- (29) for H λ (t) and (39)-(41) for X λ,∆λ (t) improve the situation further. Figure 3 shows how (a) the uniform and (b) the staggered spin susceptibilities vary in the paramagnetic phase as the Hubbard interaction U is enhanced at zero temperature. Now we are off half-filling, n = 0.9. We again find divergencies in both susceptibilities, signaling a tremendous increase of ferromagnetic and antiferromagnetic correlations. If compared to the RPA results (dashed lines), the PRM (solid lines) shifts the critical value of U to lower (higher) values in the former (latter) case. This is easily understood since the renormalization of χ(q) at q = 0 comes (solely) from the PRM band narrowing yielding a higher density of states (χ 0 (0)) and consequently a lower U c than within an RPA treatment. On the other hand, for χ(Q) the term ∆(Q) > 0 is more important, which enhances U c , i.e., suppresses the range where the state with long-range antiferromagnetic correlations pops up.

Magnetic phase diagram
Tracing the divergencies of the uniform and staggered susceptibilities in the n-U model-parameter plane at T = 0, a ground-state phase diagram of the 2D Hubbard model can be derived. Thereby, only paramagnetic states with increasingly strong ferromagnetic and antiferromagnetic correlations can be detected. The such kind determined phase diagrams agree with those obtained from the corresponding order parameter self-consistency equations [41,42,43]. Of course, the paramagnetic-ferromagnetic and paramagnetic-antiferromagnetic phase boundaries have to be determined separately. This has been undertaken in Fig. 4. Since our PRM treatment of the Hubbard model is a weakto-intermediate coupling approach, the calculations were restricted to densities not too far away from half filling (the instabilities appear at large values of U otherwise).
In the whole density range studied (0.7 ≤ n ≤ 1), the antiferromagnetic instability sets in first, i.e. an antiferromagnetic state is established before ferromagnetic order can be established. This corroborates previous Hartree-Fock, RPA, and slave-boson results [35,38,39]. Quantitative deviations from the RPA phase boundaries exist however (cf., in Fig. 4, the corresponding transition lines). A tricritical point, where the ferromagnetic and antiferromagnetic instabilities intersect, is expected to appear at a density slightly smaller than n = 0.7. It is not our aim, however, to map out the phase diagram in more detail; simply because it is now commonly accepted-owing to numerical studies but not rigorously proven-that longrange ordered phases will not be stable in the 2D Hubbard model away from half filling. Thus, as noticed long-time ago, it appears that approximative solutions to the simple 2D Hubbard model might do better in describing the magnetic features of real quasi-2D materials than the (still not available) exact solution [36]. In this respect primarily meaningful will be the improved treatment of correlations by the PRM in the paramagnetic phase.

Dynamic spin susceptibility
The phase diagram Fig. 4 from Sec. 5.2 shows that transitions from the PM to the AFM and from the PM to the FM states approach each other, when the density n gets close to half-filling (n = 1). Thereby the respective critical U values go to zero. In Fig. 5 the PRM results for the imaginary and real parts of χ(q, ω) are displayed for the density n = 0.985 very close to 1 and a small coupling U = 2. Curves are shown for different q values along the diagonal direction q = (q x , q x ) in the Brillouin zone. The steps between subsequent q x curves are chosen as π/12, where the lowest q x value is π/12. Note that Imχ(q, ω) for q x = 0 (q = 0) vanishes due to rotational symmetry of the total spin density s z q=0 . A strong paramagnon peak structure is found in Imχ(q, ω) at ω = 0 around the antiferromagnetic wave vector Q = (π, π). Also in the real part Reχ(q, ω) a strong peak structure appears for the same q and ω values. However, no such strong structure is found at small q. Since Reχ(q, ω) agrees for ω = 0 with the static susceptibility χ(q) = χ(q, ω = 0) one concludes that antiferromagnetic fluctuations at Q dominate ferromagnetic fluctuations (with q ≪ 1) already for small deviations from half-filling.
Next, let us discuss the circumstances at the density n = 0.7, which is slighter above the density where the  . The density is n = 0.985, i.e. close to half-filling, and U = 2. In both parts of χ(q, ω) a strong paramagnon excitation around ω = 0 evolves at the antiferromagnetic wave vector Q = (π, π). Since no enhancement is found in Reχ(q, ω) at the ferromagnetic wave vector q = 0 one concludes that antiferromagnetic are stronger than ferromagnetic correlations.   Fig. 5 for the density n = 0.7 and U = 11.5. This U is expected to be slightly below the U value at the tricritical point where the phase boundaries of the PM-AFM and PM-FM transitions intersect. Again, in both parts of χ(q, ω) a paramagnon structure is found for the antiferromagnetic wave vector Q for ω ≈ 0. Note that the contribution to Reχ(q, ω) at ω = 0 from the smallest q is now enhanced compared to that in Fig. 5. The evaluation of χ(q → 0) = Reχ(q → 0, ω = 0) shows that this quantity rapidly increases by slightly increasing U from 11.5 to the value U P M −F M crit ≃ 12.5, where χ(q = 0) would diverge. The red dashed curves result from the standard RPA for an extremely small q = (0.01, 0.01).
ferromagnetic and antiferromagnetic instabilities are expected to intersect (compare Fig. 4). Fig. 6 shows again the real and imaginary part of χ(q, ω), now for U = 11.5, which is approximately the critical U value for the PM-AFM transition. Indeed, for the antiferromagnetic wave vector Q there is again a paramagnon structure at ω ≈ 0 in χ(q, ω). However, as shown in the real part of χ(q, ω), also ferromagnetic fluctuations (for q ≪ 1 and ω = 0) are considerably enhanced compared to the case of Fig. 5. From the U dependence of the uniform static susceptibility χ(q = 0) (not shown), one finds that it tremendously increases for slightly increasing U and diverges at the critical value U P M−F M crit ≈ 12.5. This divergence would correspond to a transition to a ferromagnetic phase, if it had not been before the transition to the antiferromagnetic phase. The red dashed curves in both panels of Fig. 6 result from a standard RPA calculation for an extremely small q value, q = (0.01, 0.01), which is expected to almost agree with PRM results. Due to the finite lattice of 24 × 24 sites used in the PRM calculation the PRM curves in the figures are restricted to not too small q values.
Finally, let us compare the PRM with the standard RPA. As already discussed in section 2 the RPA arises when all renormalization effects are neglected. Panels (a)-(c) of Fig. 7 show the imaginary part Imχ(q, ω) as a function of ω for an intermediate density n = 0.8 and three different U values, (a) U = 0, (b) U = 2, and (c) U = 6. When q → Q the RPA curves (red dashed curves) for U = 6 exhibit a relatively narrow peak at low frequencies, which is again interpreted as paramagnon peak. In contrast, the PRM curves at the same q and U are much less pronounced. The different behavior is easily understood from the phase diagram in Fig. 4, since U = 6 is much closer to the critical RPA value U RP A crit ≃ 6.75 than to the critical value U P RM crit ≃ 8.8 from the PRM approach. The corresponding PRM result for U = 8.8 is shown in panel (d) which clearly shows a pronounced paramagnon peak as expected.
In all panels (a) to (d) the curves at Q = (π, π) are more pronounced around ω = 0 than those for q values close to the center of the Brillouin zone. From this feature one should not draw the conclusion that antiferromagnetic fluctuations are always more important than ferromagnetic ones. As was already mentioned, for ω finite Imχ(q, ω) always vanishes at q = 0, i.e. at the ferromagnetic wave vector. Therefore, as was done in section 4, a comparison of antiferromagnetic and ferromagnetic fluctuations can only be drawn from the values of the real part Reχ(q, ω) of the dynamical susceptibility at ω = 0, which is equivalent to the static q dependent susceptibility χ(q) = Reχ(q, ω = 0).

Summary
Combining linear response theory with the projector-based renormalization method, we presented a theoretical approach for the evaluations of the susceptibilities that generalizes the standard RPA scheme. In this way important  many-body correlations beyond the RPA level were included. To exemplify the advancement the theory was applied to the two-dimensional paradigmatic Hubbard model, for which an analytical expression for the dynamical spin susceptibility χ(q, ω) was derived that improves the RPA result. While the uniform spin susceptibility, where q = 0, still exhibits the standard RPA form, renormalized quasiparticle energies enter the band susceptibility contributions χ 0 (q, ω). At any finite wave-vector, however, the shape of χ(q, ω) changes. Besides momentum-dependent prefactors in χ 0 (q, ω), an extra term occurs in the denominator of χ(q, ω), which has to be determined self-consistently. This term particularly changes the pole structure of χ(q, ω) in the limit ω → 0. As a result, the magnetic phase diagram, which can be derived from the instabilities of the static spin susceptibilities at wave-vectors q = 0 or q = Q, is modified quantitatively. The same holds for the paramagnon spectrum obtained from the dynamical response. Most notably, the better (PRM) treatment of the Coulomb interaction effects severely reduce the exaggeration of the paramagnons in comparison to the RPA, i.e., the tendency of the Hubbard system to develop long-range ferromagnetic or antiferromagnetic order at certain band fillings is weakened.
where the prefactors u k,k−q,λ (t) and v k,k+p−q;k ′ ,k ′ −p,λ (t) depend on λ and also on the external field h(t).
Similarly, for an ansatz of the transformed spin s z q,λ one starts from s z q,Λ−∆λ = e X Λ,∆λ s z q,Λ e −X Λ,∆λ = s z q,Λ + [X Λ,∆λ , s z q,Λ ] + · · · , where s z q,Λ = s z q . In lowest order perturbation theory one obtains which fulfills s z q,λ = (s z −q,λ ) † . Thus, one arrives at where again the one-particle and two-particle contributions were extracted. This result agrees with Eq. (64).

B Evaluation of commutators
In this appendix we evaluate the commutators from transformation (25), which are responsible for the renormalization of H λ (t).
Since H h,λ (t) is linear in the external field, the generator X λ,∆λ (t) can be limited to the part X α λ,∆λ (t) in Eq. (39). For the part ofĤ h,λ (t) proportional to u , we find Here the factors δu (1) k,k−q,λ , and δv (1) k,k+p−q,k ′ ,k ′ −p,λ are given by δu (1) k,k−q,λ = (96) The second contribution to the above commutator arises from the term in H h,λ (t) linear in v Here, the evaluation has to be followed up by a decomposition into one-particle and two-particle contributions. This leads to renormalization contributions to u k,k−q,λ and v k,k+p−q,k ′ ,k ′ −p,λ . However, they turn out to be small and can be neglected. They are at least by a factor of order O(U/∆ε) smaller than the contributions (96), (97), where ∆ε denotes an energy difference of the order of the conduction electron band width.
Possible contributions to u k,k−q,λ and v k,k+p−q;k ′ ,k ′ −p,λ cancel. The commutator part of (120) with four fermion operators leads to higher order contributions to H f,λ , which have been neglected. The remaining 'mixed' contributions (αβ), (βα) to the commutator [X λ,∆λ , Q λ−∆λ H f,λ ]], which are of order h(t), are negligible as well, since the energies of the two Θ-function products do not coincide. Finally we consider the last commutator in transformation (25), [X λ,∆λ , [X λ,∆λ ,Ĥ h,λ ]]. Here only the four fermion part X α λ,∆λ contributes in first order in h(t). Renormalization contributions to u k,k−q,λ and v k,k+p−q;k ′ ,k ′ −p,λ are again expected to be small due to the appearance of the two Θ-function products from the two generators X α λ,∆λ .