{'fraction_non_alphanumeric': 0.03717769748727213, 'fraction_numerical': 0.016546230908195105, 'mean_word_length': 4.367232260907889, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 15, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Thermodynamics of the 2D Hubbard model

A theoretical analysis of the thermodynamic response functions of the 2D single-band Hubbard model is realized by means of the composite operator method. It is shown that all the features of these quantities can be explained by looking at the dependence of the thermodynamic variables on their conjugate ones. It is found that the numerical data from quantum Monte Carlo techniques for the internal energy and electronic specific heat are well reproduced. The anomalous normal state properties in hole-doped cuprate high Tc superconductors are also well described. Actually, the results for the linear coefficient of the electronic specific heat are in agreement with those obtained by using a pure fermionic theoretical scheme. Finally, we obtain several characteristic crossing points for the response functions when reported as functions of some thermodynamic variables. These peculiar features indicate the existence of more than one energy scale competing with thermal excitations in the 2D Hubbard model.


I. INTRODUCTION
It is believed, on both experimental and theoretical grounds, that superconductivity and charge transport in high T c cuprates are mostly confined to the CuO 2 planes [1,2]; and hence the attention of many physicists has been dedicated to 2D models which contain as an essential feature a competition between the band picture and highly correlated many body effects. Of course, some features of the phase diagram, like the existence of a finite Néel temperature, can only be explained by adding a coupling between the planes.
The bonding combination of Cu and O orbitals turns out to be quite deep below the Fermi level, so that no dynamical freedom is left to treat d and p orbitals separately [3] (there are some strong experimental evidences, mostly based on the study of the Knight shift, that in the CuO 2 plane one spin degree of freedom is observed [4]). Through the Pauli principle, the energy of the p electron excitation is, for example, largely modified by the change of charge and spin states of the neighboring Cu ions. A p electron and charge and spin fluctuations on neighboring Cu ions are simultaneously excited so that electronic excitations are formed on a CuO 2 cluster as a whole. Then, the resulting complex can be described by a single-band Hubbard model [5].
In the simplest form, the Hubbard model, first introduced to describe the correlations of electrons in a narrow d-band of transition metals, contains a kinetic term which describes the motion of the electrons among the sites of the Bravais lattice and an interaction term between electrons of opposite spin on the same lattice site. By varying the model parameters, it is believed that the Hubbard model is capable of describing many properties of strongly correlated fermion systems. Among different examples, the Hubbard Hamiltonian is applicable to describe the metal-insulator transition in a series of transition metal oxides such as Sr 1−x La x T iO 3 [6,7] and V 2 O 3 [8][9][10][11]. The applicability of the model to the superconducting copper-oxides is related to the fact that upon doping most of these compounds exhibit a metal-insulator Mott transition; the superconducting state is near the Néel state and there are many experimental results [12][13][14][15] which show a close relation between the antiferromagnetic correlations in the Cu − O planes and the occurrence of the superconducting phase. However, it is important to stress that an appropriate description of a bad metal with large energy scale spin fluctuations by means of a purely electrostatic Hamiltonian should preserve the symmetry expressed by the Pauli principle that codifies the correct interplay between charge and magnetic configurations [16][17][18].
Although considerable attention has been devoted to the Hubbard model and significant progress was achieved in understanding ground state properties, particularly at half-filling, static and dynamical spin correlations, the optical conductivity and other observables, a clear comprehension of the low-lying excitations is still lacking [2]. The difficulty is not to be found only in the absence of any obvious small parameter in the strong coupling regime. More deeply, it is due to the difficulty of handling simultaneously itinerant aspects (spatial correlations) and atomic aspects (pronounced on-site quantum fluctuations) [19].
In recent years we have been developing a method of calculation, denominated Composite Operator Method (COM) [17,18,[20][21][22][23][24][25][26][27][28][29][30][31][32], that has been revealed to be a powerful tool for the description of local and itinerant excitations in strongly correlated systems. In previous papers, we considered the Hubbard electronic operators for the determination of fundamental excitations. A fully self-consistent calculation of the electronic propagator has been realized by means of a constraint with the physical content of the Pauli principle [17,[23][24][25]. We calculated local quantities, as the double occupancy and the magnetic moment [17,23], the energy per site [24], the chemical potential [24], the magnetic susceptibility [25,29], the density of states and the quasi-particle spectra [28,32]. In all the cases, the results show a good agreement with those obtained by numerical simulation. In particular, the results obtained for the magnetic properties can reproduce the unusual characteristics observed in high T c superconducting materials [25,29,33]. Therefore, the agreement strengthens the idea that a microscopic single-band model contains the essential physical features of the new class of materials.
In this paper we investigate the electronic specific heat and the entropy per site of the 2D Hubbard model for a paramagnetic ground state. It will be shown that all the features of these quantities can be understood by looking at the dependence of the chemical potential and double occupancy on their conjugate thermodynamic variables, that is, the particle concentration and the on-site Coulomb repulsion, respectively. A comprehensive comparison among different methods to compute the specific heat will shed a new light on the approximation used. It will emerge that in our theoretical scheme, even if dynamical effects in the self-energy are neglected promoting unstable collective asymptotic modes to the role of well-defined quasi-particle excitations, extended spin modes can be captured by properly combining symmetry requirements and extended operatorial basis. Indeed, the presence of a low temperature peak that appears when the low-lying spin states are excited will appear as an important feature shared with the quantum Monte Carlo data [34]. In addition, an alternative fermionic scenario for the computation of the electronic specific heat together with the results for the Wilson ratio will confirm cuprates as dominated by conventional fermionic excitations. An extensive study of the thermodynamic response functions will reveal the existence of critical lines which separate different energy scales created by the interplay between charge and spin modes. In other words a study of thermodynamics quantities, such as the double occupancy, the entropy, the chemical potential, the specific heat, indicates lines in the U − T plane which separate a highly-correlated behavior, dominated by spin and charge fluctuations and a non-interacting behavior, dominated by thermal fluctuations. In particular, there emerges a region of filling where the entropy reduces by increasing the filling signalling the set up of an ordered phase. For T → 0 there is a well-defined marginal concentration where a quantum phase transition occurs. A detailed comparison with the non-interacting case will be also presented throughout the paper.
The plan of the article is as follows. In the next Section we present the 2D Hubbard model and the electron propagator in the COM. In Sec. III we review experimental data for some thermodynamic properties. The results for the electronic specific heat are presented in Sec IV, where a theoretical understanding of the different ways to compute the specific heat is also presented. Section V is devoted to a discussion of double occupancy. In Sec. VI the results for the chemical potential versus temperature are discussed. The entropy is analyzed in Sec VII. Some concluding remarks are presented at the end.

II. ELECTRON PROPAGATOR IN THE HUBBARD MODEL
The Hubbard model is defined by The notation is the following. The variable i stands for the lattice vector R i . c(i), c † (i) are annihilation and creation operators of c−electrons at site i, in the spinor notation: t ij denotes the transfer integral and describes hopping between different sites; the U term is the Hubbard interaction between two c−electrons at the same site with being the charge-density operator per spin σ. n(i) is the total charge-density operator. µ is the chemical potential. In the nearest neighbor approximation, for a two-dimensional cubic lattice with lattice constant a, we write the hopping matrix t ij as The scale of the energy has been fixed in such a way that t ii = 0. It should be noted that since the interactions are restricted to the same site, the dimensionality of the system comes in only when a specific form for α(k) is taken. In other words, the stabilization of eventual cooperative phenomena is uniquely governed by the band dispersion.
The point of view adopted in the COM is that the Heisenberg operators {c(i), c † (i)} are not good candidates as a basis for calculations. Because of strong correlations the c−electrons loose their identity and new fields, whose properties are self-consistently determined by the dynamics and by the symmetries of the model, together with the boundary conditions, might be more appropriate as a starting point for the physical description of the system. Due to the on-site Coulomb interaction, it is known that two sharp features develop in the band structure which correspond to the Hubbard subbands and describe interatomic excitations mainly restricted to subsets of the occupancy number. Indeed, a first natural choice for composite fields is given by the Hubbard constrained electronic operators describing the transitions (n = 0) ⇐⇒ (n = 1) and (n = 1) ⇐⇒ (n = 2), respectively. The two-point retarded thermal Green's function is defined as where Ψ(i) is the doublet composite operator The bracket . . . indicates the thermal average and R is the usual retarded operator.
In previous papers, we have shown that the determination of the single-particle Green's function (2.8) can be realized in a fully self-consistent way once a unique approximation is made [17,[23][24][25]. This approximation consists in neglecting the dynamical part in the self-energy and corresponds to a pole expansion of the spectral intensities. As shown in Ref. 23-25, by considering time translational invariance and no magnetic order, the single-particle electronic propagator has the following expression Ω and Ω B being the volume of the unit cell in the direct and reciprocal space, respectively. f F (ω) is the Fermi distribution function. E 1 (k) and E 2 (k) are the energy spectra A 1 (k) and A 2 (k) are the spectral intensities with the definitions where n = c † c is the particle density. The parameters ∆ and p are static intersite correlation functions defined as The notation c α (i) stands to indicate the field c on the first neighbor sites: We are also using the following notation A crucial point in the calculation of the propagator (2.10) is the implementation of the Pauli principle by means of constraining equations to determine the internal parameters µ, ∆ and p. Use of this principle leads to the selfconsistent equations [23,25] The recovery of the Pauli principle, very often violated by other approximations, assures a dynamics bounded to the Hilbert space capable of describing in a correct way the interplay between the charge and the magnetic configurations. Furthermore, we have shown [35] that in the two-pole approximation [36] the set of self-consistent equations (2.18) is the only one which restores the particle-hole symmetry and the Pauli principle, which are intimately connected. It is possible to go beyond the two-pole approximation by enlarging the set of asymptotic fields [18,[20][21][22] or by taking into account the dynamical corrections to the self-energy [31,32].

III. THERMODYNAMICS AS REVEALED BY EXPERIMENTS
The electronic specific heat C(T ) of cuprate high-T c superconductors has been measured. In particular C(T ) of La 2−x Sr x CuO 4 [37,38] has been studied for 0.03 < x < 0.45 in the range of temperatures between 1.5 and 300 K, and of Y Ba 2 Cu 3 O 6+y [39] for 0.16 ≤ y ≤ 0.97 between 1.8 and 300 K. From these experiments the following behavior has been observed for the coefficient γ = C/T of the normal state specific heat: a) for fixed temperature, γ(x, T ) increases with doping; a1) in the case of La 2−x Sr x CuO 4 , γ(x, T ) exhibits a rather sharp maximum at x ≈ 0.25 (near the doping where superconductivity disappears), then starts to decrease; the same behavior for γ(x, T ) has been estimated in Ref. 40, but with a peak located around x ≈ 0.18, close to the optimal doping; for La 2−x Ba x CuO 4 [41] a maximum has been observed at x ≈ 0.22; a2) in the case of Y Ba 2 Cu 3 O 6+y , γ(x, T ) increases smoothly to a plateau or two broad maxima, situated at y ≈ 0.6 and y ≈ 0.9, respectively; b) for fixed doping, γ(x, T ) as a function of temperature exhibits a broad peak moving to lower temperatures with increasing the dopant concentration; c) further increasing y, the T-dependence weakens and in the region of high doping no increase is observed. For Y Ba 2 Cu 3 O 6+y , no substantial increase is observed for y > 0.8.
As noticed by Vollhardt [42], there is a peculiar feature of the specific heat observed in a large variety of systems. The specific heat curves versus T, when plotted for different, not too large, values of some thermodynamic variable, intersect at one or even two well defined temperatures. In 3 He the specific heat C(T, P ) curves versus T at different pressures P intersect at a well defined temperature [43,44]; in heavy fermions CeAl 3 [45] and U Be 13 [46] upon change of P , U P t 3−x P d x [47] and CeP t 3 Si 1−x Ge x [48] upon change of x, CeCu 6−x Au x when either P [49] or the magnetic field B [50] is varied; in semi-metal, Eu 0.5 Sr 0.5 As 3 [51] upon change of B.
The following properties have been observed for the entropy S [38, 39, 52]: a) for a given temperature, S increases with doping; a1) in the case of La 2−x Sr x CuO 4 [38], S(x, T ) reaches a maximum in the vicinity of x ≈ 0.25, then decreases; a2) in the case of Y Ba 2 Cu 3 O 6+y [39,52], S reaches a maximum in the vicinity of y ≈ 0.97; b) for a given dopant concentration, S exhibits a superlinear dependence on the temperature; c) the normal state entropy as a function of T extrapolates to a negative value at T = 0K; d) there is a striking numerical correlation between S/T and aχ 0 , where χ 0 is the bulk susceptibility and a is the Wilson ratio.

A. GENERAL FORMULAS
The specific heat C(T ) is defined as where E is the internal energy density, given by the thermal average of the Hamiltonian N being the number of sites. Calculation of internal energy by means of Eq. (4.2) will generally require the calculation of two-particle Green's functions. An alternative way to calculate the internal energy is the following. By introducing the Helmholtz free energy per site where S is the entropy per site, from the thermodynamics we have Then, it is straightforward to obtain the following formulas from which the specific heat turns out to be In this scheme the thermodynamic quantities are all expressed through the chemical potential, whose determination requires the knowledge of the single-particle Green's function. Summarizing, we have two distinct ways to calculate the internal energy, based on the use of Eqs. (4.2) and (4.7). In principle these equations are equivalent and lead to the same result when an exact solution is available. However, the situation drastically changes when approximations are involved and different results can be obtained. Indeed an open problem in Condensed Matter Physics is to find a unique consistent scheme of approximation capable of treating on an equal footing, both one-and two-particle Green's functions.
Sometimes in the literature the internal energy for electronic systems is calculated by means of the expression where N (ω) is the density of states. However, Eq (4.9) is based on the assumption that the system admits a description in terms of fermionic quasi-particles. It is well known that the correlation among the original electrons can generate extended bosonic modes, and use of expression (4.9) requires special attention. By general argument it can be shown that for interacting systems the correct expression for the internal energy, calculated by means of the density of states, is given by where H I is the non-quadratic part of the Hamiltonian in the chosen canonical representation.

B. NON-INTERACTING CASE
To discuss the specific heat it is useful at first to consider the non-interacting [i.e. U=0] Hubbard model. In this case the thermal retarded Green's function can be exactly calculated and has the expression where the energy spectrum has the expression The chemical potential is determined as a function of n and T by means of the equation where we put As we discussed above, we have different ways of calculating the free energy. In the non-interacting case, where an exact solution is available, all different procedures must give the same result. Since this point will acquire some relevance in the interacting case, we shall examine this in detail.
In the non-interacting case the density of states is given by By substituting this expression into Eq. (4.10) we have where a factor 2 has been included in order to take into account the spin degree of freedom. The same result is obtained by taking the thermal average of the Hamiltonian where use has been made of the expression (4.10) for the single-particle Green's function.
The proof that use of Eq. (4.7) leads to the same result requires some work. By taking the derivative with respect to T of Eq. (4.12) we have By considering that the derivative with respect to n of Eq. (4.12) gives substitution of (4.17) into (4.8) leads to This concludes the proof that in the non-interacting case all expressions (4.2), (4.3) and (4.8) give the same result. We also note that the specific heat can be calculated by means of the following expression In Figs. 1 and 2 the linear coefficient of the electronic specific heat γ(T, n) is plotted as a function of temperature and filling, respectively.
a For a fixed filling γ(T, n) first increases as a function of T, exhibits a maximum at a certain temperature T m , and then decreases. The value of T m decreases by increasing the filling. At half-filling T m is zero and γ(T, n) diverges as T → 0; this is an effect of the van Hove singularity (vHs). The temperature behavior of γ(T, n) is similar to the one exhibited by the static uniform spin magnetic susceptibility χ 0 [12,25], whereas the doping dependence is different.
b1. At T = 0, γ(T, n) increases by increasing the filling, diverges at n = 1 and then decreases. b2. At finite temperature the peak splits in two peaks, symmetric with respect to n = 1. The distance between the two peaks increases by increasing T. This is shown in Fig. 2, where γ(T, n) is plotted as a function of the filling at various temperatures. The shift of the two peaks with respect to n = 1 increases by increasing T.

C. INTERACTING CASE
In the interacting case the expressions (4.2) and (4.7) for the internal energy give different results. By recalling (2.10), Eq. (4.2) gives where D is the double occupancy. We use E H to indicate that the internal energy has been calculated by the average of the Hamiltonian. E F is the part of the internal energy calculated on the assumption of only fermionic elementary excitations Alternatively, we have seen that the internal energy and the specific heat can be calculated by means of Eqs. (4.8) and (4.9). This procedure requires a knowledge of the first and second temperature derivatives of the chemical potential. Let us define By taking the first derivative with respect to T of the self-consistent equations (2.18), we obtain (1) 1 (4.26) Explicit calculation of the derivatives defined in Eq. (4.26) shows that the equations (4.25) provide a set of linear algebraic equations for the three parameters µ 1 , ∆ 1 , p 1 , as functions of the parameters µ, ∆, p. In the same way, by taking the second derivative with respect to T of Eqs. (2.18) we obtain (2) 1 (4.28) These equations provide a set of linear algebraic equations for the three parameters µ 2 , ∆ 2 , p 2 as functions of the parameters µ, ∆, p, µ 1 , ∆ 1 , p 1 . Once the self-consistent calculation of the three parameters µ, ∆, p has been performed by means of the set (2.18), then the calculation of the first and second derivatives of the chemical potential reduces to the solution of simple linear equations.
In the interacting case, because of the approximation used, the different procedures to calculate the internal energy give different results. At first we shall compare our theoretical results with the data obtained by numerical analysis.
The specific heat C(T ) of the 2D Hubbard model has been recently calculated in Ref. 34 by using quantum Monte Carlo techniques. In particular, the Monte Carlo data for the energy per site E = H /N have been fitted by polynomials and the specific heat has been calculated by taking derivatives from these polynomials analytically. Different polynomials have been chosen in different regions of temperature. In the calculation of E , it has been found that finite size effects are strong at weak coupling but become negligible for U ≥ 8.  In Fig. 3 we present the internal energy versus temperature in the range 0 ≤ T ≤ 5 for several values of doping and U = 4. The results are compared with QMC data of Ref. 34. We are using E T to indicate the solution that comes from Eq. (4.8). As a general feature we observe that E T corresponds to the lowest energy solution and agrees very well with the QMC data, in the entire region of temperature and for all studied dopant concentrations. When U increases, the theoretical solution deviates from QMC in the region of low temperatures, indicating that the antiferromagnetic (AF) correlations are not properly taken into account in the strong coupling regime. This is clearly seen in Fig. 4 where the internal energy is plotted versus U for half-filling. We now consider the specific heat. There are two important features in the Monte Carlo calculations: 1) a low temperature peak that appears when the low-lying spin states are excited, and 2) a high temperature peak which appears when states in the upper Hubbard band are excited. In the weak coupling regime the low temperature peak moves to slightly higher temperature as U increases, reaching a turning point at U = 7 where the peak is at T = 0.3. For U > 7 the peak slowly moves to lower temperatures, as U grows. This indicates the beginning of the strong coupling regime. The broad high temperature peak moves to higher temperatures as U increases, as expected since its presence corresponds to the excitation of states across the gap that grows with U . In the strong regime [U > 7] the position of the charge peak increases linearly with U: T charge ≈ 0.24U . In addition, all the curves C versus T for several values of U intersect at T ≈ 1.6. The QMC results for the specific heat of the 2D Hubbard model qualitatively agree with the half-filled 1D Hubbard model [53,54].
The specific heat C T = dE T /dT , calculated by means of Eq. (4.9), is compared with QMC data in Fig. 5 for U = 4 and various dopant concentrations. At low density the agreement is generally good, in both the weak and strong coupling cases [U=8]. At higher densities the QMC data show a double peak structure, which is enhanced at half filling, but also present at n = 0.75 for U = 8. In the discussed range of values of U our results do not show a double peak structure. The presence of two peaks in the specific heat has been attributed to the spin and charge excitations. When U is weak the two peaks overlap and there is no resolution. By increasing U the position of the charge peak moves to higher temperatures and we expect to be able to distinguish the two contributions. A study of the specific heat in the strong coupling regime is given in Figs. 6 and 7, where the two expressions C H = dE H /dT and C T = dE T /dT are plotted, respectively, as functions of T at half-filling. A peak appears at low temperatures when U is rather large [say U ≥ 8 for C H and U ≥ 15 for C T ]. This behavior qualitatively reproduces the QMC results. In order to isolate the spin excitations for lower values of U we can subtract from C T the contribution coming from the pure fermionic part, C F = dE F /dT , responsible in large extent for the charge peak structure. This is done in Fig. 8 where we report C T − C F as a function of T for various values of n and U . A low-energy excitation appears and its contribution generally increases by increasing n and U . In Ref. 34 the position of the charge peak T charge has been calculated for different values of U. In our analysis when U is large T charge can be calculated from C T and the results lie on the line T charge = 0.24U . When U is lower, it is not possible to resolve the two peaks. Since the charge excitation mainly comes from the interband fermionic excitations, we can estimate T charge from the calculation of C F . This study is shown in Fig. 9, where T charge is plotted versus U. The solid line is the position of the charge peak, calculated by means of C F ; the triangles denote T charge calculated by means of C T ; the squares are the QMC results. The linear coefficient of the specific heat as a function of the particle density has been studied by QMC for U = 8 and T varying from 0.5 to 3. Unfortunately, due to the sign problem it is not easy to study low temperatures in QMC. A study of γ(x, T ), by means of C T , shows a good agreement for high temperatures, but not at lower temperatures, where QMC results exhibit a strong downward deviation in the region of filling where we should expect an increase of γ(x, T ) due to the effect of the van Hove singularity.
Up to this point we have performed a detailed and ample comparison with QMC, generally finding a quite reasonable agreement in a large region of values for the model parameters. For intermediate values of U [U=4] the agreement is quite satisfactory. On this basis we are confident that the approximation used is adequate for the 2D Hubbard model and we can pass to examine the next question as to which extent the physics of real systems is retained in the model.
One feature present in a large variety of systems is a characteristic crossing point in the specific heat curves versus T [42][43][44][45][46][47][48][49][50][51]. This behavior has been also found in 1D models [53][54][55] and in the Hubbard model in infinite dimension [19], where a crossing temperature T = 0.59 has been observed in the range 0.5 ≤ U ≤ 2.5. For the 2D Hubbard model the same behavior, as predicted by Vollhardt [42] has been observed by means of quantum Monte Carlo calculations [34], where for the case of half-filling a crossing temperature T = 1.6 ± 0.2 has been observed in the range 2 ≤ U ≤ 12. In Fig. 10 the specific heat C T is given versus T for half-filling and various values of U. When n = 1 the curves cross at the same temperature T ≈ 2.0; when doping is considered the region of crossing spreads out and moves to higher temperatures. From the thermodynamic relations, this crossing temperature T U corresponds to a turning point of the double occupancy as a function of T , that is (∂ 2 D/∂T 2 ) TU = 0. A study of the function (∂ 2 D/∂T 2 ) will be done in the next Section. As mentioned in Section III, the linear coefficient of the specific heat of cuprates exhibits an anomalous behavior in the normal state. To investigate this, in Fig. 11 we present the linear coefficient γ(x, T ) as a function of the doping x = 1 − n for various temperatures. As a general behavior we see that by increasing the doping γ(x, T ) increases up to a certain doping and then decreases. The nature of the peak is due to the fact that the Fermi energy crosses the vHs for a certain critical value x c . The value of x c depends on the ratio U/t and varies between 0 and 1/3, as U increases from zero to infinite. For U/t = 4 it is found x c = 0.27, very close to the experimental value observed in La 2−x Sr x CuO 4 [37,38]. At half-filling the Fermi energy (ǫ F ) is at the center of the two Hubbard bands; by varying the dopant concentration some weight is transferred from the upper to the lower band, ǫ F moves to lower energies and crosses the vHs for a critical value of the doping; increasing x, further moves ǫ F away from the vHs. A study of the Fermi surface shows that for x > x c we have a closed surface which becomes nested at x = x c and opens for x < x c . An enlarged Fermi surface with a volume larger than the non-interacting one has been reported by QMC calculations [56,57] and by other theoretical works [58,59].  The peak position of γ(x, T ) depends on the temperature. In the limit of zero temperature a sharp peak is exactly located at x = x c . By increasing the temperature the peak moves away from x c and broadens into two peaks. The situation is similar to what we have calculated for the non-interacting case [cfr. Fig. 2]; we find that the role played by the interaction manifests itself through the shift of the vHs and the band structure which creates an asymmetry between the two peaks. The behavior described in Fig. 11 well reproduces the experimental situation.
In the case of La 2−x Sr x CuO 4 the experimental data [37,38] for γ(x, T ) are reported in Fig. 12 We also mention that the specific heat of La 2−x Sr x CuO 4 has been estimated in Ref. 40 from the data for the heat capacity anomaly at the superconducting transition temperature by assuming a BCS-type relation. Under this assumption the authors find the same behavior for γ(x, T ), but with a peak located around x ≈ 0.18, close to the optimal doping. In Figs. 13a and 13b we present the linear coefficient γ(x, T ) as a function of the temperature for values of the filling x > x c and x < x c , respectively. At x = x c we see that γ(x, T ) diverges as T → 0; this is an effect of the vHs. When x = x c the Fermi energy moves away from the vHs and the peak exhibited by γ(x, T ) moves away from T=0. As shown in Figs. 13a and 13b, γ(x, T ) as a function of temperature has different behaviors in the two regions x > x c and x < x c . In the overdoped region γ(x, T ) firstly increases as a function of T, exhibits a maximum at a certain temperature T m and then decreases. This behavior is similar to the one exhibited by χ 0 (T ) [12,25]. As shown in Fig.  13a, when the doping decreases the value of T m moves to lower temperatures. This behavior qualitatively agrees with the non interacting case showing that for x > x c the AF correlations are weak. A different situation is observed in the underdoped region, where γ(x, T ) is always a decreasing function of T. When we look at the experimental results for La 2−x Sr x CuO 4 [37,38] and for Y Ba 2 Cu 3 O 6+y [39,52] we find that the behavior of γ(x, T ) as a function of T in the underdoped region is more similar to that for the non-interacting case, showing that cuprate materials could be conveniently described in terms of a Fermi liquid picture, dominated by the Van Hove singularity. Then, if we want to use the Hubbard model to describe high-T c cuprates, we might assume a scheme for non-interacting -like particles and use Eq. (4.23) to calculate the specific heat. When we do this, we find the results reported in Figs. 14, 15a and 15b.    As a function of the doping γ(x, T ) has a similar behavior for that part dominated by the Van Hove singularity, but the band structure plays a different role. By comparing Figs. 11 and 14 it is worthwhile to stress that the presence of a maximum in γ(x, T ) does not necessarily indicate a pure fermionic excitation spectrum; it is a consequence of the enhancement in the density of states when the Fermi energy crosses the vHs. This becomes more evident in Fig. 15a  and 15b where the temperature dependence is studied at fixed doping. Both in the underdoped and overdoped regions γ(x, T ) exhibits a maximum at a certain temperature. The behavior given in these figures qualitatively reproduces the experimental situation reported in Ref. 37,38,39,52. The fact that for Y Ba 2 Cu 3 O 6+y γ(x, T ) is always a decreasing function of T when y > 0.8 is understood because by approaching the critical doping, T m is shifted to low temperatures, below the critical superconducting temperature.

V. THE DOUBLE OCCUPANCY
As a simple thermodynamic quantity indicating the degree of correlation of the system, in this Section we study the double occupancy D, defined as the fraction of doubly occupied sites This quantity can be calculated by means of the expression Then, the first and second temperature derivatives of D can be analytically calculated as once the self-consistent equations have been solved. We display in Figs. 16a and 16b the temperature dependence of D for various values of U at particle concentrations n = 0.7 and n = 0.8. In Fig. 16a the data display a characteristic low temperature behavior: D initially decreases with temperature, reaches a minimum, and increases again. In other words, the curve indicates the presence of a T region where the formation of local magnetic moments is enhanced with increasing T [the double occupancy D determines the local spin-spin correlation function S 2 through the equation S 2 = 3(n − 2D)/4]. This behavior is characteristic of incipient localization effects in a strongly correlated Fermi liquid in a regime dominated by spin fluctuations. Starting from the low temperature Fermi liquid regime, when the temperature increases, the system can gain free energy by localizing the particles (i.e. decreasing D) in order to take advantage of a larger spin entropy [19,53]. In the absence of spin excitations one would observe decreasing values with increasing T.  Fig. 22b indicates the temperature TU where the curve changes curvature.
In Fig. 16b D is a monotonic increasing function of temperature. In this case the values of n and U are large enough to inhibit localization effects due to the increase of temperature. To study this behavior in more detail this behavior, the derivative with respect to the temperature of the double occupancy has been analyzed. The results show that for a given T there exists a critical value of U, say U D (T ), such that The function U D (T ), defined by ∂D/∂T = 0, is given in Fig. 17. We note that at T=0 U D (0) coincides with U c (n, 0), defined as the critical strength of the on-site Coulomb interaction for which the Fermi energy crosses the vHs at fixed dopant concentration. U D (T ) goes to zero for some temperature T D . For n = 0.7 we find that T D = 0.581. When T > T D we have ∂D/∂T > 0 for all values of U. The behavior of T D as a function of n is reported in the Fig.  18. The fact that (∂D/∂T ) U=UD(T ) = 0 implies that at U = U D (T ) the double occupancy does not depend on T. However, the curves of D as function of U for different values of T will not cross in a single point because U D (T ) changes with the temperature in a significant way. At very high temperature T ≫ U , larger than T U where there is a change in the concavity, D asymptotically tends to the atomic value n 2 /4, as expected. We have seen in the previous section that the specific heat curves versus T for different values of U cross almost at the same point T U , determined by the equation A study of this equation by means of formula (5.3) gives the results plotted in Fig. 19, where T U is plotted versus U for n = 0.75.  The function n µ (T ) is presented in Fig. 20. It can be shown that at T = 0 n µ coincides with n c , the critical value where the Fermi level crosses the van Hove singularity. But, the temperature dependence of n µ (T ) is remarkably different from the one of n c (T ). Thus, only for T → 0 we may relate the transition ∂µ/∂T < 0 =⇒ ∂µ∂T > 0 to the reversal of the sign of the derivative for the density of states at the Fermi level. Furthermore, we observe that n µ (T ) reaches the value of 1 for some temperature T µ (for U = 4 we find T µ = 0.843). When T > T µ we have ∂µ/∂T < 0 for all values of n . The fact that (∂µ/∂T ) n=nµ (T ) = 0 implies that at n = n µ (T ) the chemical potential does not depend on T. Therefore, the curves of µ as function of n for various values of T , reported in Fig. 29, will not cross exactly in the same point as claimed in Ref. 63.

VII. THE ENTROPY
The entropy S(T, n), connected to the total number of spin and charge excitations at temperature T and filling n, is a bulk thermodynamic quantity uniquely determined by the spectrum of excitations, whose magnitude and temperature dependence provide an important test for proposed theories. Theoretical works available so far are the following. Bipolaron models propose preformed boson charge carriers at T c and behaving classically at higher temperatures. Apart from some inconsistency related to the magnitude of the entropy, these theories have to resort to the existence of thermally excited triplet bipolarons in explaining the deviation of S from a linear T-dependence in underdoped samples [64]. Theoretical studies of the entropy in the strongly correlated electron systems have also been performed in the framework of the statistical spin liquid [65][66][67][68]. This scheme is based on the assumption that in the strongly correlated metals the doubly occupied single-spin configurations must be excluded not only in the real space representation, but also in the reciprocal space. By means of the spin liquid statistics, the entropy of localized moments is reproduced when the Mott insulator limit is reached for half-filling. Nevertheless, this over imposed statistics freezes the system in a wrong Hilbert space whenever different choices of the parameters modify the interplay between thermal excitations and electronic interactions. Some theories [3] predict decoupled holon (boson) and spinon (fermion) excitations. In these approaches it is difficult to reconcile the experimentally observed magnitude for the entropy with its partition between statistically independent excitations. Moreover, the striking numerical correlation between S/T and aχ 0 is expected for weakly interacting fermions but not if the dominant excitations are those of spinless bosons. In Ref. 63 exact diagonalization studies of the t − J model have been performed; for several thermodynamic quantities a critical doping concentration that marks a change of the Fermi surface character is found.
By means of the relation (4.7), we have calculated the entropy per site S(T, n). Recalling the behavior of ∂µ/∂T , we see that the entropy must have the following behavior (i) for T < T µ , S(T, n) increases with increasing particle concentration, reaches a maximum for n = n µ , then decreases (see Fig. 22a); (ii) for T > T µ , S(T, n) always increases with increasing particle concentration (for U = 4 we find T µ = 0.843) as is shown in Fig. 22b. Again the peak structure reflects a Fermi level crossing the vHs at the critical doping. This behavior is in agreement with the experimental data from Refs. 38, 39 and 52. Indeed, the experiments show a well defined peak structure in a large region of temperature (from 40K to 320K); furthermore, the position of the peak slightly changes with temperature. In the theoretical analysis the position of the peak as a function of temperature is governed by n µ (T ), reported in Fig. 20, that shows a smooth variation in the region of physical relevance (T < 0.05). In Fig. 23 the entropy versus the particle concentration is reported for various values of the Coulomb interaction. The maximum of the entropy shifts to lower values of n by increasing U, varying between 1 and 2/3 when U varies from zero to ∞. In Fig. 24, we see that for T = 0.4 all entropy curves for different U cross at the same particle concentration. By comparison with figure 21 and by numerical analysis, we see that the crossing point is exactly the critical concentration n µ (T ). Recalling the Maxwell relation the behavior shown in Fig. 22 implies (∂D/∂T ) n=nµ = 0. If we remember that the double occupancy D determines the local spin-spin correlation function, it is clear that a sign reversal of its derivative with respect to the temperature represents a crossover from a regime dominated by spin fluctuations, where S is a decreasing function of U, to another regime favouring charge fluctuations (electronic delocalization), where the entropy is an increasing function of U. In Figs. 25 and 26 we report the temperature dependence of S and S/T for several dopant concentrations. The curves have a qualitative agreement with the experimental ones [38,39,52]. In Fig. 25a, where x > x c , S is a decreasing function of x at a fixed temperature; the opposite behavior is observed in Fig. 25b.  In the limit of zero temperature the entropy goes to zero by a linear law. When T increases the entropy deviates from the linear behavior. In the region 0.01 ≤ T ≤ 0.1 the temperature dependence is well described by the law where the coefficients S 0 , S 1 , S 2 are strongly dependent on the filling.  It is worth noticing that in the limit of large temperatures (a weak point of the statistical spin liquid [65][66][67][68]) our results for the entropy asymptotically agree with the exact expression where χ 0 is the bulk susceptibility and a is the Wilson ratio a = π 2 3 (7.5) In the case of La 2−x Sr x CuO 4 [38] and of Y Ba 2 Cu 3 O 6+y [39,52,69], there is a striking numerical correlation between S/T and aχ 0 . As noticed by Loram et al., this resemblance shows that the total spin+charge spectrum over all moments k (from S) and the k = 0 spin spectrum (from χ 0 ) have a similar energy dependence. Also, experimental evidences suggest that the low-energy excitations are predominantly those of conventional fermions, and that the substantial T dependencies of S/T and χ 0 are primarily determined by the energy dependence of the single-particle density of states in the vicinity of the Fermi level. In Fig. 27 we present S/T χ 0 as a function of the temperature for various values of doping (i.e., 0.03 ≤ x ≤ 0.27). For all values of dopant concentration S/T χ 0 is almost constant over a wide range of temperatures.
In addition, the value of S/T χ 0 is less than the non-interacting one. This is due to the fact that by introducing interaction the number of microscopic states accessible to the same macroscopic state is reduced (i.e., the entropy per site) whereas the susceptibility is increased by incipient localization effects.
In order to obtain a better understanding of how the thermodynamics of an electronic liquid is modified by the interaction, we have performed a study of non-interacting Hubbard model (i.e. U = 0). What we learn from the study of this model is that the critical value (i.e. half-filling) above which the entropy looks a decreasing function of the filling is uniquely fixed by the statistics. The temperature has no role. On the contrary, in the interacting case the energy scale of charge configurations has a crucial role in the region of particle concentration between n µ and half-filling. For n µ < n < 1 there is a critical temperature T µ , depending on the filling, above which the behavior is similar to the non-interacting case [i.e. where ∂S/∂n becomes positive, or where ∂µ/∂T becomes negative].
We now consider the physical origin of these results. In the non-interacting case the combinatorics dictated by the Fermi statistics governs the behavior of the entropy. This can be understood if we think that for the single-site problem the entropy has the values 0, ln 2 and 0 for occupancy 0, 1 and 2, respectively. In the interacting case it is natural to look for a critical value of the filling above which the number of permutations satisfying the restrictions of the boundary conditions starts to decrease. For n µ < n < 1, because of the Coulomb interaction, by increasing the particle density the number of microscopic realizations, accessible to the same observable macroscopic state, decreases [the Pauli principle is obviously crucial for a correct counting]. This is true only if the thermal excitations do not exceed the energy scale fixed by the interaction. Definitely, for 0 < n < n µ we have a sort of disordered non-interacting state with ∂S/∂n > 0, whereas for n µ < n < 1 the low-lying excitations characterize a far from random spatial pattern [i.e. ∂S/∂n < 0]. In the range n µ < n < 1 incommensurate magnetism and superconductivity are experimentally observed [12,70].

VIII. CONCLUDING REMARKS
The 2D single-band Hubbard model has been studied by means of the composite operator method. By considering the Hubbard operators as basic set of fields, which describe interatomic excitations restricted to subsets of the occupancy number, the single-particle electronic propagator has been computed in a fully self-consistent way by means of a quasi-particle scheme capable of coherently integrating dynamics, boundary conditions and symmetry principles.
The paper was devoted to the study of the electronic specific heat and entropy per site in the paramagnetic phase. We analyzed these quantities by looking at the dependence of the thermodynamic variables on their conjugate ones, that is, for example, the relation between entropy and temperature, chemical potential and particle concentration, double occupancy and on-site Coulomb repulsion. Once the self-consistent equations for the single-particle propagator have been solved, we have determined the temperature derivatives of the internal parameters by means of exact linear systems of algebraic equations. The determination of the first and second temperature derivatives of the chemical potential has been revealed crucial in determining the thermodynamic response functions under investigation. For the electronic specific heat and internal energy we have presented three different schemes of calculation. All of them allowed the possibility to obtain a deep theoretical understanding of how and to which extent collective excitations can be retained in the description of thermal response functions. We have obtained a good agreement with the data by quantum Monte Carlo techniques for the electronic specific heat and the internal energy [34]. Further on, although Monte Carlo data shared common features with the results from the calculations through the T-derivative of the chemical potential, the experimental data for cuprates, as revealed by the Wilson ratio and linear coefficient of the electronic specific heat, have shown that in such systems the dominant excitations are those of conventional non-interacting fermions [38].
We obtained several characteristic crossing points for the response functions when reported as functions of some thermodynamic variables. These peculiar features, already evidenced by Vollhardt [42], marked turning points where different response functions evolve from a non-interacting behaviour (i) the entropy is an increasing function of U; (ii) the entropy is an increasing function of n; (iii) the double occupancy is a decreasing function of T; (iv) the T-derivative of the chemical potential is a decreasing function of n; (v) the linear coefficient of the specific heat is an increasing function of n; to an unconventional dependence on the conjugate variables (vi) the entropy is a decreasing function of U; (vii) the entropy is a decreasing function of n; (viii) the double occupancy is an increasing function of T; (ix) the T-derivative of the chemical potential is an increasing function of n; (x) the linear coefficient of the specific heat is a decreasing function of n.
Before closing we would like to mention that the region of filling, where (vi)-(x) hold, coincides with that where incommensurate magnetism and superconductivity are experimentally observed in LSCO cuprates family.