# Nearly universal crossing point of the specific heat curves of Hubbard models

## I Introduction

Recently attention was drawn to the fact that in various strongly correlated systems the curves of the specific heat $`C(T,X)`$ vs. temperature $`T`$ cross once or twice when plotted for different values of a second thermodynamic variable $`X`$. For example, crossing points are observed for different pressures ($`X=P`$) in normalfluid <sup>3</sup>He (Ref. ) and heavy-fermion systems such as CeAl<sub>3</sub> and UBe<sub>3</sub>. By changing the magnetic field ($`X=B`$) the same feature is seen in heavy-fermion compounds such as CeCu<sub>6-x</sub>Al<sub>x</sub> and Nd<sub>2-x</sub>Ce<sub>x</sub>CuO<sub>4</sub>. Crossings of the specific heat curves are also observed in the simplest lattice model for correlated electrons, the Hubbard model,

$$\widehat{H}=\underset{𝐤\sigma }{}(ϵ_𝐤\mu )\widehat{a}_{𝐤\sigma }^+\widehat{a}_{𝐤\sigma }^{}+U\underset{i}{}\widehat{n}_i^{}\widehat{n}_i^{},$$
(1)   
where $`ϵ_𝐤`$ is the dispersion of a single electronic band, $`\mu `$ the chemical potential, and $`U`$ the local interaction. At half filling the curves $`C(T,U)`$ vs. $`T`$ always cross at two temperatures. This is observed, for example, in the case of the model with nearest-neighbor hopping in $`d=1`$,  $`d=2`$, and $`d=\mathrm{}`$, as well as for long-range hopping in $`d=1`$; for the latter two systems the specific heat is shown in Fig. 1. Furthermore, crossing is found in $`d=1`$ when a magnetic field $`B`$ is changed at constant $`U`$. The fact that these crossing points may be very sharp was analyzed in Ref. , and was traced to the properties of certain generalized susceptibilities of the system.

In the following, we will consider only the crossing of the specific heat curves occuring for $`X=U`$ in the paramagnetic phase of the Hubbard model with a symmetric half-filled band ($`n=1`$). We will investigate yet another observation, namely that for small $`U`$ the specific heat at the high-temperature crossing point has practically the same value of approximately $`0.34k_B`$ for all dimensions $`d`$ and dispersions $`ϵ_𝐤`$, which can be seen also in Fig. 1. This is surprising, because the temperatures at which this crossing occurs are very different for different dispersions and dimensions, and because the maximum value of $`C(T,U)`$ and its value at the low-temperature crossing point vary strongly, too. It should be noted that the specific heat, like the entropy, is a dimensionless quantity when expressed in units of $`k_B`$.

We denote by $`T_+`$ the temperature at which the curves $`C(T,U)`$ vs. $`T`$ cross for different values of $`U`$. Then the specific heat is independent of $`U`$ at the crossing temperature $`T_+(U)`$, defined by

$$\frac{C}{U}|_{T_+(U)}=0$$
(2)   
Since we are not concerned with the dependence of $`T_+(U)`$ on $`U`$, but rather with the crossing point value of $`C(T,U)`$ for different lattice systems, we consider only the limit of small $`U`$ and define

$`T_+`$ $``$ $`\underset{U0^+}{lim}T_+(U),`$ (3)   
$`C_+`$ $``$ $`\underset{U0^+}{lim}C(T_+(U),U).`$ (4)   
This weak-coupling crossing point can be calculated without approximation using second order perturbation theory in $`U`$. Higher orders in perturbation theory would be necessary to determine the dependence of $`C(T_+(U),U)`$ on $`U`$.

It is the purpose of this paper to show that the value of $`C_+`$ at the high temperature crossing point is almost universal, and to analyze the origin of this peculiar feature. We calculate $`C_+`$ for a half-filled band with a symmetric density of states (DOS) $`N(ϵ)=𝑑𝐤\delta (ϵϵ_𝐤)`$. We also show that the weak dependence of $`C_+`$ on lattice properties can be understood by starting from the limit $`d=\mathrm{}`$ and using expansions in terms of two small parameters.

This paper is organized as follows. In Sec. II we review the origin of crossing points in the Hubbard model and show how to calculate $`T_+`$ and $`C_+`$ in second order perturbation theory, the details of which are contained in Appendix A. Various non-interacting lattice systems are listed in Sec. III, and the values of $`C_+`$ at the high-temperature crossing point for these systems are presented in Sec. IV. Sec. V contains expansions of $`C_+`$ that reveal the influence of the density of states and the lattice dimension. We close with a conclusion in Sec. VI.

## II Crossing points in the specific heat of the Hubbard model

The entropy per lattice site $`S(T,U)`$ is given by

$$S(T,U)=\underset{0}{\overset{T}{}}𝑑T^{}\frac{C(T^{},U)}{T^{}}.$$
(5)   
For the Hubbard model $`S(T,U)`$ approaches a constant when $`T\mathrm{}`$. Taking the derivative of Eq. (5) with respect to $`U`$ we find

$$0=\underset{0}{\overset{\mathrm{}}{}}\frac{dT}{T}\frac{C(T,U)}{U}.$$
(6)   
Since $`C/U`$ is not identically zero, there must exist temperature regions where it has positive and negative values. We assume that no phase transitions occur, so that $`C/U`$ is a continuous function of $`U`$. Then there exist temperatures where $`C/U`$ changes sign; at these temperatures the curves $`C(T,U)`$ vs. $`T`$ cross \[see Eq. (2)\].

There are two such crossing points in the paramagnetic phase of the half-filled Hubbard model, as can be seen from the sign of $`C/U`$ at very low and very high temperatures. For intermediate values of $`U`$ the specific heat of this model shows the following general features: $`C(T,U)`$ starts linear in $`T`$ at low temperatures and develops a two-peak structure, one at temperatures $`T4t^2/U`$ due to spin excitations ($`t`$ is the hopping amplitude), and one at temperatures $`TUW`$ due to charge excitations, where $`W`$ is the bandwidth. For low temperatures the spin exitations become stronger for increasing $`U`$, thus $`C/U>0`$. At high temperatures an increase in $`U`$ pushes out the charge peak and thus increases $`C(T,U)`$, which tends to $`(a+bU^2+O(U^4))/T^2`$ with $`a,b>0`$. Hence $`C/U>0`$ for both high and low temperatures, so that the sum rule \[Eq. (6)\] yields an intermediate region with $`C/U<0`$. There are thus two sign changes corresponding to two crossing points.

To determine the location of the crossing points defined by Eqs. (3) and (4), we calculate the internal energy per lattice site in perturbation theory in $`U`$,

$$E(T,U)=E^{(0)}(T)+\frac{1}{4}U+U^2E^{(2)}(T)+O(U^3).$$
(7)   
Here ($`k_B1`$, $`\beta =1/T`$)

$$E^{(0)}(T)=2𝑑ϵ\frac{N(ϵ)ϵ}{1+\mathrm{exp}(\beta ϵ)},$$
(8)   
which is the internal energy for the non-interacting system, $`U/4`$ is the Hartree contribution, and the second-order correlation energy is given by (see Appendix A for details)

$$E^{(2)}(T)=\frac{}{\beta }\frac{\beta ^2}{32}\underset{0}{\overset{\text{ }1}{}}𝑑x𝑑𝐤𝑑𝐩𝑑𝐪\frac{\mathrm{cosh}[\frac{1}{2}x\beta (ϵ_𝐤+ϵ_𝐩+ϵ_{𝐩+𝐪}+ϵ_{𝐤+𝐪})]}{\mathrm{cosh}(\frac{1}{2}\beta ϵ_𝐤)\mathrm{cosh}(\frac{1}{2}\beta ϵ_𝐩)\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤+𝐪})\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐩+𝐪})},$$
(9)   
where the integrations, e. g. $`𝑑𝐤d^dk/(2\pi )^d`$, run over the first Brillouin zone.

In the limit of infinite spatial dimensions this expression can be simplified further. In this case, momentum conservation at vertices becomes irrelevant, so that the integrals factorize (see Appendix A)

$`E^{(2)}(T)`$ $`=`$ $`{\displaystyle \frac{}{\beta }}{\displaystyle \frac{\beta ^2}{32}}{\displaystyle \underset{0}{\overset{\text{ }1}{}}}𝑑x\left[{\displaystyle 𝑑ϵN(ϵ)\frac{\mathrm{cosh}(\frac{1}{2}x\beta ϵ)}{\mathrm{cosh}(\frac{1}{2}\beta ϵ)}}\right]^4.`$ (10)   
Note that as usual in infinite dimensions the dispersion $`ϵ_𝐤`$ enters into one-particle quantities only via the DOS $`N(ϵ)`$. Therefore this expression is much easier to evaluate numerically than Eq. (9).

The specific heat $`C(T,U)=E/T`$ has the expansion

$$C(T,U)=C^{(0)}(T)+U^2C^{(2)}(T)+O(U^4),$$
(11)   
where

$`C^{(0)}(T)`$ $`=`$ $`{\displaystyle \frac{\beta ^2}{2}}{\displaystyle 𝑑ϵ\frac{N(ϵ)ϵ^2}{\mathrm{cosh}^2(\frac{1}{2}\beta ϵ)}},`$ (12)   
and the function $`C^{(2)}(T)`$ can be written as (see Appendix A)

$`C^{(2)}(T)`$ $`=`$ $`{\displaystyle \frac{\beta ^2}{32}}{\displaystyle \frac{^2}{\beta ^2}}\beta ^2{\displaystyle \underset{0}{\overset{\text{ }1}{}}}𝑑x{\displaystyle \underset{m}{}}\left[f_m(x,\beta )\right]^4.`$ (13)   
Here the sum runs over lattice sites $`𝐑_m`$ and the functions $`f_m(x,\beta )`$ are given by

$`f_m(x,\beta )`$ $`=`$ $`{\displaystyle 𝑑𝐤\frac{\mathrm{exp}\left(i𝐤𝐑_m+\frac{1}{2}x\beta ϵ_𝐤\right)}{\mathrm{cosh}(\frac{1}{2}\beta ϵ_𝐤)}}.`$ (14)   
Comparison with Eq. (10) shows that in $`d=\mathrm{}`$ only the local term with $`𝐑_m=0`$, i. e. $`f_0(x,\beta )`$ $`=`$ $`𝑑ϵ`$ $`N(ϵ)`$ $`\mathrm{cosh}(\frac{1}{2}x\beta ϵ)`$ $`/`$ $`\mathrm{cosh}(\frac{1}{2}\beta ϵ)`$, contributes to the sum in Eq. (13).

The crossing point in the specific heat occurs at the temperature $`T_+(U)`$ for which $`C(T,U)`$ is independent of $`U`$. In view of Eqs. (2), (3), and (11), the crossing temperature $`T_+`$ in the limit $`U0`$ is given by the root of the equation

$$C^{(2)}(T_+)=0,$$
(15)   
and the specific heat at the crossing point, Eq. (4), is

$$C_+=C^{(0)}(T_+).$$
(16)   
These equations will be evaluated for several lattices and dimensions that are described in the next section.

## III Momentum dispersion and density of states

We consider only one-band systems at half-filling with a symmetric density of states on lattices in finite and infinite spatial dimensions. For systems with finite bandwidth we set $`W/21`$, where $`W`$ is the bandwidth, while for infinite bandwidth we use a unit second moment of the density of states, i. e. $`𝑑ϵN(ϵ)ϵ^21`$.

1. Finite Dimensions. For the linear chain, square lattice, and simple cubic lattice, i. e., the hypercubic lattices in $`d=1,2,3`$, we use the tight-binding dispersion $`ϵ_𝐤=2t_{i=1}^d\mathrm{cos}k_i,`$ $`|k_i|\pi `$, which describes nearest neighbor (NN) hopping with amplitude $`t1/2d`$. Furthermore we study the body-centered cubic (bcc) lattice in $`d=3`$ with NN hopping, which can be regarded as a subset of the simple cubic lattice with hopping across the space-diagonal, so that $`ϵ_𝐤=8t\mathrm{cos}(k_x)\mathrm{cos}(k_y)\mathrm{cos}(k_z)`$ with $`t\frac{1}{8}`$. For these systems we use perturbation theory as described in Sec. II. Finally, for one-dimensional long-range hopping $`t(r)1/r`$ the known interacting dispersion can be used instead of perturbation theory. The free dispersion is $`ϵ_k=tk`$, $`t1/\pi `$, with a constant density of states. Fig. 2 shows the various densities of states used in $`d=1,2,3`$.

2. Infinite dimensions. We consider first the hypercubic lattice and generalized honeycomb lattice with NN hopping. For the hypercubic lattice, the hopping must be scaled as $`t=1/\sqrt{2d}`$ to obtain a nontrivial limit for $`d\mathrm{}`$. In this case the density of states becomes a Gaussian with unit variance, $`N(ϵ)=\mathrm{exp}(ϵ^2/2)/\sqrt{2\pi }`$, whereas for the generalized honeycomb lattice the same scaling leads to $`N(ϵ)=|ϵ|\mathrm{exp}(ϵ^2)`$. We also study the Bethe lattice with infinite connectivity and semicircular density of states, $`N(ϵ)=\frac{2}{\pi }\sqrt{1ϵ^2}`$, $`|ϵ|1`$.

Furthermore we can take advantage of the fact that in $`d=\mathrm{}`$ only the density of states of the non-interacting system appears in Eq. (10). Hence it may be chosen at will even when no corresponding dispersion is known. We consider three such functions each containing a tunable real parameter $`\alpha >0`$. This allows us to study the behavior of $`C_+`$ for a wide range of DOS shapes. The first is the “metallic” density of states

$$N(ϵ)=\frac{1+\alpha }{2\alpha }(1|ϵ|^\alpha ),|ϵ|1,$$
(17)   
bearing this name due to its finite value at the Fermi energy, while the following “semi-metallic” density of states

$$N(ϵ)=\frac{1+\alpha }{2}|ϵ|^\alpha ,|ϵ|1,$$
(18)   
vanishes at the Fermi energy. Furthermore we employ a semi-metallic density of states “with tails” having infinite bandwidth and unit variance

$$N(ϵ)=\frac{c^c}{\mathrm{\Gamma }(c)}|ϵ|^\alpha \mathrm{exp}(cϵ^2),c\frac{1+\alpha }{2}.$$
(19)   
It reduces to the DOS for the generalized honeycomb lattice with NN hopping in the case of $`\alpha =1`$. The important special case of a constant rectangular density of states, $`N(ϵ)=\frac{1}{2}`$ for $`|ϵ|1`$, is contained in Eq. (17) for $`\alpha \mathrm{}`$ and in Eq. (18) for $`\alpha 0`$. Note that in the limit $`\alpha \mathrm{}`$, the DOS in Eq. (18) approaches two $`\delta `$-peaks, $`N(ϵ)=\frac{1}{2}(\delta (1+ϵ)+\delta (1ϵ))`$. This particular case is of interest only because in this case the integrals in Eqs. (8) and (10) can be calculated analytically. For general $`\alpha `$, on the other hand, the functions in Eqs. (17) and (18) model typical DOS shapes for finite dimensions. Several densities of states used in $`d=\mathrm{}`$ are depicted in Fig. 3.

## IV Results for the specific heat at the high-temperature crossing point

In this section we present results for the specific heat $`C_+`$ at the high-temperature crossing point for the density of states discussed in Sec. III. To calculate $`C_+`$ according to Eqs. (15) and (16), the integrals appearing in Eqs. (12) and (13) have to be calculated numerically. We determine them to high precision (typically $`10^8`$) by either Monte-Carlo integration (using the VEGAS algorithm) or by a high-temperature expansion (described in Appendix A). For several cases both methods were applied and yielded the same results within numerical accuracy. Typical results for the functions $`E^{(2)}(T)`$ and $`C^{(2)}(T)`$ are shown over a wide temperature range in Fig. 4 for the linear chain with nearest neighbor hopping and for the constant DOS in infinite dimensions. There is a maximum at lower and a minimum at higher temperatures in $`E^{(2)}(T)`$, corresponding to the two zeros of $`C^{(2)}(T)`$. These are the temperatures where the specific heat curves cross for $`U0`$.

Numerical values for $`T_+`$ and $`C_+`$ are listed in Table I. In view of the drastically different DOS shapes (see Fig. 2 and 3) it is quite remarkable that the values of $`C_+`$ for these systems are very similar, ranging from 0.331 to 0.358. For the DOS with tunable parameter $`C_+`$ is plotted vs. $`\alpha `$ in Fig. 5. Again, $`C_+`$ varies only by a small amount although the shapes of the DOS change strongly. In particular, for the semi-metallic DOS \[Eq. (18)\] $`C_+`$ lies between the values for the constant and two $`\delta `$-peak DOS, since these are the limits of Eq. (18) for $`\alpha 0,\mathrm{}`$.

These results raise the following questions which will be addressed in the next section: (i) Why is $`C_+`$ at the high-temperature crossing point so insensitive against changes of the DOS of the non-interacting electrons and the spatial dimension? (ii) What determines the (small) spread in the values of $`C_+`$? (iii) Why is the value of $`C_+`$ at the low-temperature crossing point varying so much stronger?

## V Expansions of $`C_+`$

We will see that the influence of the DOS and of the crystal lattice on $`C_+`$ can be understood by expanding around the limit of $`d=\mathrm{}`$. In addition, in $`d=\mathrm{}`$ the dispersion $`ϵ_𝐤`$ enters only via the DOS $`N(ϵ)`$, so that the effect of its form on $`C_+`$ can be studied by expanding in terms of the difference between $`N(ϵ)`$ and a reference DOS $`\overline{N}(ϵ)`$.

### A Influence of the DOS in $`d`$ $`=`$ $`\mathrm{}`$

First we consider infinite-dimensional systems with an arbitrary symmetric DOS $`N(ϵ)`$ with finite bandwidth. The DOS is compared to a rectangular-shaped DOS with the same bandwidth, $`\overline{N}(ϵ)=\frac{1}{2}`$ for $`|ϵ|<1`$. Their difference is characterized by the quantity

$$\delta N=\underset{1}{\overset{1}{}}𝑑ϵ|N(ϵ)\overline{N}(ϵ)|,$$
(20)   
which will serve as an expansion parameter. We will expand $`C_+`$ to lowest order in $`\delta N`$, making use of the known functions $`\overline{f}_0(x,\beta )`$, $`\overline{C}^{(2)}(\beta )`$, $`\overline{C}^{(0)}(\beta )`$ pertaining to $`\overline{N}(ϵ)`$, as well as the known crossing points for the rectangular DOS, $`\overline{T}_+=1/\overline{\beta }_+=0.13801`$ and $`\overline{C}_+=0.44046`$ at low temperatures, while $`\overline{T}_+=1/\overline{\beta }_+=0.57190`$ and $`\overline{C}_+=0.33935`$ at the high-temperature crossing point.

We begin by expanding $`C^{(2)}(\beta )`$ for small $`\delta N`$:

$$C^{(2)}(\beta )=\overline{C}^{(2)}(\beta )+4\left[\frac{\beta ^2}{32}\frac{^2}{\beta ^2}\beta ^2\underset{0}{\overset{\text{ }1}{}}𝑑x[f_0(x,\beta )\overline{f}_0(x,\beta )]\overline{f}_0(x,\beta )^3\right]+O((\delta N)^2).$$
(21)   
Using standard inequalities it can be shown that the neglected terms indeed vanish like $`(\delta N)^2`$ (for fixed $`\beta <\mathrm{}`$). Now $`\overline{C}^{(2)}(\beta )`$ is expanded around the known crossing point $`\overline{\beta }_+`$ where $`\overline{C}^{(2)}(\overline{\beta }_+)=0`$. Then the new crossing point temperature is determined from the condition

$$0\stackrel{!}{=}C^{(2)}(\beta _+)=(\beta _+\overline{\beta }_+)\overline{C}^{(2)}(\overline{\beta }_+)+4\left[\frac{\beta ^2}{32}\frac{^2}{\beta ^2}\beta ^2\underset{0}{\overset{\text{ }1}{}}𝑑xf_0(x,\beta )\overline{f}_0(x,\beta )^3\right]_{\beta =\overline{\beta }_+}+O((\delta N)^2).$$
(22)   
To lowest order the shift in $`\beta _+`$ is hence given by

$$\delta \beta _+=\beta _+\overline{\beta }_+=\underset{1}{\overset{1}{}}𝑑ϵN(ϵ)\mathrm{\Delta }\beta _+(ϵ)+O((\delta N)^2),$$
(23)   
where

$$\mathrm{\Delta }\beta _+(ϵ)=\frac{\overline{\beta }_+^2}{8\overline{C}^{(2)}(\overline{\beta }_+)}\left[\frac{^2}{\beta ^2}\frac{\beta ^2}{\mathrm{cosh}(\frac{1}{2}\beta ϵ)}\underset{0}{\overset{\text{ }1}{}}𝑑x\overline{f}_0(x,\beta )^3\mathrm{cosh}(\frac{1}{2}x\beta ϵ)\right]_{\beta =\overline{\beta }_+}.$$
(24)   
Finally we expand $`C^{(0)}(\beta )`$ in $`\beta \overline{\beta }_+`$ to first order and evaluate it for $`\beta =\overline{\beta }_++\delta \beta _+`$. Thus the specific heat at the crossing point $`C_+`$ is obtained to lowest order in $`\delta N`$ as

$$C_+=\overline{C}_++\underset{1}{\overset{1}{}}𝑑ϵN(ϵ)\mathrm{\Delta }C_+(ϵ)+O((\delta N)^2),$$
(25)   
where

$`\mathrm{\Delta }C_+(ϵ)`$ $`=`$ $`{\displaystyle \frac{\overline{\beta }_+^2ϵ^2}{2\mathrm{cosh}^2(\frac{1}{2}\overline{\beta }_+ϵ)}}`$ (27)   
$`+\overline{C}^{(0)}(\overline{\beta }_+)\mathrm{\Delta }\beta _+(ϵ)\overline{C}^{(0)}(\overline{\beta }_+).`$   
The last two equations show how the first order effect on $`C_+`$ of a deviation of $`N(ϵ)`$ from a rectangular shape can be determined by a single integration.

The numerical evaluations of the function $`\mathrm{\Delta }C_+(ϵ)`$ are plotted in Fig. 6 for both crossing points. The amplitudes of the function $`\mathrm{\Delta }C_+(ϵ)`$ corresponding to the low- and high-temperature crossing points are seen to differ greatly, i. e. by a factor of about 40. This implies a much greater sensitivity towards changes of the DOS and the dimension of the value of $`C_+`$ at the low-temperature crossing point.

It is also clear that $`C_+`$ is not entirely universal at the high-temperature crossing point, since for a general density of states $`N(ϵ)`$ the integral in Eq. (25) does not vanish. We can in fact estimate the maximum value of the shift $`\delta C_+=C_+\overline{C}_+`$ for arbitrary $`N(ϵ)`$ with finite bandwidth, using simple integral inequalities:

$`|\delta C_+|`$ $``$ $`\mathrm{min}(a_1,a_2\underset{0ϵ1}{\mathrm{max}}\left[N(ϵ)\right])+O((\delta N)^2),`$ (28)   
with $`a_1\mathrm{max}_{0ϵ1}|\mathrm{\Delta }C_+(ϵ)|`$, $`a_2_1^1𝑑ϵ|\mathrm{\Delta }C_+(ϵ)|`$. At the high-temperature crossing point we find $`a_1=0.02268`$ and $`a_2=0.02200`$, i. e. to order $`O(\delta N)`$ we have $`|\delta C_+|\overline{C}_+`$, which is the reason for the insensitivity of $`C_+`$ to changes in $`N(ϵ)`$. Furthermore the predicted range of values $`0.339\pm 0.023`$ corresponds well to the observed range of $`0.331`$-$`\mathrm{\hspace{0.17em}0.358}`$; see Table I. On the other hand, at the low-temperature crossing point we find $`a_1=0.9330`$ and $`a_2=0.7727`$, so that $`|\delta C_+|\overline{C}_+`$. Hence $`C_+`$ is indeed not confined to a small interval in this case.

To check the validity of this expansion we applied Eq. (25) to several infinite-dimensional systems at the high-temperature crossing point. Results for the Bethe lattice, as well as for the metallic DOS \[Eq. (17)\] and the semi-metallic DOS, \[Eq. (18)\] at several values of the parameter $`\alpha `$ are given in Table II. As expected we find that the lowest-order approximation in Eq. (25) describes the behavior of $`C_+`$ very well if the deviation from the rectangular DOS is not too large. The difference between the exact value of $`C_+`$ and the approximate value $`\overline{C}_++\delta C_+`$ is due to corrections of order $`(\delta N)^2`$ and hence is typically an order of magnitude smaller than the first order correction $`\delta C_+`$.

Hence we have shown that for a small change $`\delta N`$ in the density of states the variation of $`C_+`$ at the low-temperature crossing point is large, while at the high-temperature crossing point it is small and well-described by the first order correction in $`\delta N`$. This gives quantitative answers to the questions posed at the end of Sec. IV.

### B Lattice effects in dimension $`d`$ $`=`$ 1, 2, 3

To determine the influence of the lattice dimension on the value of $`C_+`$ we study the hypercubic lattice with NN hopping $`t`$ and the scaling $`t=1/\sqrt{2d}`$ and expand $`C_+`$ in $`1/d`$. Since the large variation of $`C_+`$ at the low-temperature crossing point can already be understood from the results of the previous subsection, we will perform this calculation only for the high-temperature crossing point.

As discussed in Sec. II, for $`d=\mathrm{}`$ only the local term involving $`f_0(x,\beta )`$ remains in the lattice sum in Eq. (13). The $`1/d`$-corrections are given by the contribution from the $`2d`$ nearest neighbors vectors such as $`𝐑_1=(1,0,0,\mathrm{})`$,

$`{\displaystyle \underset{m}{}}f_m(x,\beta )^4=f_0(x,\beta )^4+2df_1(x,\beta )^4+\mathrm{}`$ (30)   
$`=f_0(x,\beta )^4+{\displaystyle \frac{1}{2d}}\left[{\displaystyle \frac{2}{\beta }}{\displaystyle \frac{}{x}}f_0(x,\beta )\right]^4+O(d^2),`$   
where a partial derivative w.r.t. $`x`$ was employed in the second line to remove the factor $`\mathrm{cos}(𝐤𝐑_1)`$ from $`f_1(x,\beta )`$. Hence the corrections to $`C_+`$ in order $`1/d`$ can be calculated entirely from $`f_0(x,\beta )`$. This function contains an integral over the DOS $`N(ϵ)`$ only, which in turn must be expanded in $`1/d`$. Then the finite-$`d`$ correction to $`T_+`$ and $`C_+`$ can be calculated by expanding Eqs. (15) and (16) to first order in $`1/d`$, with the following result for the high-temperature crossing point:

$`C_+`$ $`=`$ $`0.343630+0.013599{\displaystyle \frac{1}{d}}+O(d^2),`$ (31)   
$`T_+`$ $`=`$ $`\left[0.847667+0.082650{\displaystyle \frac{1}{d}}+O(d^2)\right]{\displaystyle \frac{1}{\sqrt{2d}}}.`$ (32)   
Here the temperature scale has been reset to our previous choice of unit half bandwidth in finite dimensions (i. e., the extra factor $`1/\sqrt{2d}`$ must be omitted in order to recover $`T_+`$ in $`d=\mathrm{}`$). Again this expansion compares very well with the numerical results; see Table III. Note that the coefficient of $`d^1`$ in Eq. (31) is already much smaller than $`C_+`$ in $`d=\mathrm{}`$, which is the reason for the insensitivity of $`C_+`$ to lattice effects, such as Brillouin zone shape and momentum conservation. We expect that for longer-range hopping the smallness of the deviations in the value of $`C_+`$ may equally be traced to $`1/d`$-corrections.

## VI Conclusion

For many correlated electronic systems the curves of the specific heat vs. temperature obtained for different values of a second thermodynamic parameter $`X`$ are known to intersect. For the Hubbard model at half-filling the specific heat curves for different values of the Hubbard interaction $`U`$ cross twice, the crossing point at high-temperatures being remarkably sharp up to intermediate values of $`U`$. We observed that $`C_+`$, the value of the specific heat at this crossing point in the weak-coupling limit, is practically the same for a several different lattice systems. We analyzed the origin of this conspicuous feature by calculating $`C_+`$ in perturbation theory in $`U`$. We found the values of $`C_+`$ at the high-temperature crossing point to occur in very small interval, i. e. $`C_+0.34`$ is indeed almost independent of dimensionality, crystal lattice, and energy dispersion. This is not the case for the crossing point at low temperatures, where $`C_+`$ varies on a much larger scale.

Qualitatively, the reason for this difference can be traced to the relevant energy scales on which $`C(T,U)`$ varies. At high temperatures, the energy scale for $`T`$ is essentially determined by the bandwidth, i. e. by the hopping amplitude $`t`$ in the dispersion $`ϵ_𝐤`$. At low temperatures, on the other hand, the generation of low-energy excitations (which are responsible for the strong enhancement of the low-temperature specific heat and $`C/U>0`$) leads to a renormalized energy scale $`tt_{\text{eff}}t`$. The first maximum in $`C(T,U)`$ occurs at a temperature that is of the order of $`t_{\text{eff}}`$ (see Fig. 4). As a consequence, the first sign change in $`C/U`$ is also linked to $`t_{\text{eff}}`$, so that the intersection of $`C(T,U)`$ and $`C(T,0)`$ at low temperatures does not occur at any predetermined value. In contrast, the second sign change in $`C/U`$ is determined only by energy scales that also appear in the non-interacting system, leading to a nearly universal value for the high-temperature intersection of $`C(T,U)`$ and $`C(T,0)`$.

To gain a more quantitative understanding we identified two small parameters which determine the crossing point values $`C_+`$. The starting point for expansions in these small parameters is the limit of infinite dimensions ($`d=\mathrm{}`$). (i) For $`d=\mathrm{}`$ the dependence of $`C_+`$ on the shape of a DOS $`N(ϵ)`$ with finite bandwidth is well described by the first-order correction in the parameter $`\delta N=𝑑ϵ|N(ϵ)\frac{1}{2}|`$. This parameter is a measure of the difference between $`N(ϵ)`$ and a constant rectangular DOS. It turns out that at the high-temperature crossing point this correction is small for almost all DOS, while it is large at the low-temperature crossing point. (ii) For hypercubic lattices in dimensions $`1d<\mathrm{}`$ the value of $`C_+`$ may be obtained by an expansion around $`d=\mathrm{}`$ in powers of $`1/d`$. At the high-temperature crossing point the value of $`C_+`$ of the $`d`$-dimensional system is already accurately determined by the first-order correction in $`1/d`$ even for $`d`$ as low as $`d=1`$, due to the smallness of the prefactor of this term. These expansions show in detail why $`C_+`$ has an almost universal value at the high-temperature crossing point.

## A Internal energy at weak coupling

The internal energy per lattice site $`E`$ is given by

$$E(T,U)=T𝑑𝐤\underset{\omega _n,\sigma }{}e^{i\omega _n\eta }\frac{ϵ_𝐤+\frac{1}{2}\mathrm{\Sigma }_\sigma (𝐤,i\omega _n)}{G_0(i\omega _n,𝐤)^1\mathrm{\Sigma }_\sigma (𝐤,i\omega _n)},$$
(A1)   
where $`\omega _n`$ denotes fermionic Matsubara frequencies, $`\eta 0^+`$, and $`\mathrm{}1`$. The non-interacting Green function is $`G_0(i\omega _n,𝐤)^1=i\omega _n(ϵ_𝐤\mu )`$, and $`\mathrm{\Sigma }_\sigma (𝐤,i\omega _n)`$ is the self-energy for spin $`\sigma `$. In the paramagnetic phase the sum over spins just gives a factor of two, and the spin indices on the self-energy can be dropped.

For a symmetric DOS the chemical potential $`\mu `$ at half filling is given by $`U/2`$ for all temperatures due to particle-hole symmetry. It is useful to define shifted functions $`\widehat{G}_0^1=G_0^1U/2`$ and $`\widehat{\mathrm{\Sigma }}=\mathrm{\Sigma }U/2`$, with the new chemical potential fixed at 0. Up to second order in $`U`$ $`\widehat{\mathrm{\Sigma }}`$ is given by only one Feynman diagram where the lines now represent $`\widehat{G_0}`$. We write $`\widehat{\mathrm{\Sigma }}(𝐤,i\omega _n)=U^2\widehat{\sigma }(𝐤,i\omega _n)+O(U^3)`$, with

$`\widehat{\sigma }(𝐤_1,i\omega _n)=\text{}`$ $`=`$ $`T^2{\displaystyle \underset{\omega _l,\nu _m}{}}{\displaystyle 𝑑𝐤_2𝑑𝐤_2𝑑𝐤_3}`$ (3)   
$`\times \widehat{G}_0(𝐤_2,i\omega _l)\widehat{G}_0(𝐤_3,i\omega _l+i\nu _m)\widehat{G}_0(𝐤_4,i\omega _n+i\nu _m){\displaystyle \underset{𝐊}{}}\delta (𝐤_1𝐤_2+𝐤_3𝐤_4𝐊).`$   
where $`\nu _m`$ denote bosonic Matsubara frequencies and the sum is over reciprocal lattice vectors $`𝐊`$. Expansion of the internal energy \[Eq. (A1)\] in powers of $`U`$ yields Eqs. (7), (8), and

$$E^{(2)}(T)=T𝑑𝐤\underset{\omega _n}{}\frac{e^{i\omega _n\eta }}{i\omega _nϵ_𝐤}\left[\widehat{\sigma }(𝐤,i\omega _n)+2ϵ_𝐤\frac{\widehat{\sigma }(𝐤,i\omega _n)}{i\omega _nϵ_𝐤}\right].$$
(4)   
The frequency summations are carried out as usual, with the result

$$T\underset{n}{}\frac{\widehat{\sigma }(𝐤_1,i\omega _n)}{i\omega _nϵ_{𝐤_1}}=𝑑𝐤_2𝑑𝐤_3𝑑𝐤_4\frac{\mathrm{sinh}[\frac{1}{2}\beta (ϵ_{𝐤_1}ϵ_{𝐤_2}+ϵ_{𝐤_3}ϵ_{𝐤_4})]}{8[ϵ_{𝐤_1}ϵ_{𝐤_2}+ϵ_{𝐤_3}ϵ_{𝐤_4}]_{i=1}^4\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤_i})}.\underset{𝐊}{}\delta (𝐤_1𝐤_2+𝐤_3𝐤_4𝐊),$$
(5)   
Here we eliminate the energy denominator using the identity $`\mathrm{sinh}(\frac{1}{2}\beta y)=\frac{1}{2}\beta _0^1𝑑x\mathrm{cosh}(\frac{1}{2}x\beta y)`$. Next we rewrite the second term in square brackets in Eq. (4) as a formal derivative $`2ϵ_{𝐤_1}/ϵ_{𝐤_1}`$ of Eq. (5). We make the integrand symmetric in all $`ϵ_{𝐤_i}`$ by shifting $`𝐤_2`$ and $`𝐤_4`$ by $`𝐐`$, where $`ϵ_{𝐤+𝐐}=ϵ_𝐤`$. Then the derivative can be replaced by $`2\beta /\beta `$ because only products $`\beta ϵ_{𝐤_i}`$ appear. Taking $`\beta `$ inside the derivative we finally arrive at

$`E^{(2)}(T)`$ $`=`$ $`{\displaystyle \frac{}{\beta }}{\displaystyle \frac{\beta ^2}{32}}{\displaystyle \underset{0}{\overset{\text{ }1}{}}}𝑑x{\displaystyle 𝑑𝐤_1𝑑𝐤_2𝑑𝐤_3𝑑𝐤_4}`$ (7)   
$`\times {\displaystyle \frac{\mathrm{cosh}[\frac{1}{2}x\beta (ϵ_{𝐤_1}+ϵ_{𝐤_2}+ϵ_{𝐤_3}+ϵ_{𝐤_4})]}{\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤_1})\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤_2})\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤_3})\mathrm{cosh}(\frac{1}{2}\beta ϵ_{𝐤_4})}}{\displaystyle \underset{𝐊}{}}\delta (𝐤_1𝐤_2+𝐤_3𝐤_4𝐊).`$   
In $`d=1,2,3`$ this simplifies to Eq. (9), while for $`d=\mathrm{}`$ the $`\delta `$-function can be omitted and the numerator can be replaced by $`_{i=1}^4\mathrm{cosh}(\frac{1}{2}x\beta ϵ_i)`$, leading to Eq. (10).

The integrals in Eq. (7) also factorize if we express momentum conservation as a sum over lattice vectors $`𝐑_m`$,

$`{\displaystyle \underset{𝐊}{}}\delta (𝐤_1𝐤_2+𝐤_3𝐤_4𝐊)={\displaystyle \underset{m}{}}\mathrm{exp}\left[i(𝐤_1𝐤_2+𝐤_3𝐤_4𝐊)𝐑_m\right],`$ (8)   
leading to Eqs. (13) and (14). Since the functions $`f_m(x,\beta )`$ in Eq. (14) cannot be calculated in closed form for general $`ϵ_𝐤`$ we employ a high-temperature expansion, which yields

$`f_m(x,\beta )`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{\beta ^n}{n!}}E_n(\frac{1+x}{2})I_n(𝐑_m),`$ (9)   
$`I_n(𝐑_m)`$ $`=`$ $`{\displaystyle 𝑑𝐤(ϵ_𝐤)^n\mathrm{exp}\left(i𝐤𝐑_m\right)},`$ (10)   
where $`E_n(x)`$ are Euler polynomials. The remaining Brillouin zone integrals are calculated as follows. For the $`d`$-dimensional hypercubic lattice the lattice vectors are $`𝐑_m=(m_1,\mathrm{}m_d)`$ with integer $`m_i`$, so that for NN hopping $`t`$ we have

$$I_n^{\text{hc}}(𝐑_m)=(t)^n\underset{n_1+\mathrm{}n_d=n}{}\left(\genfrac{}{}{0pt}{}{n}{n_1,\mathrm{}n_d}\right)\underset{i=1}{\overset{d}{}}\left[\underset{\pi }{\overset{\pi }{}}\frac{dk_i}{2\pi }(2\mathrm{cos}k_i)^{n_i}\mathrm{cos}(k_i|m_i|)\right].$$
(11)   
The integral in square brackets equals $`\left(\genfrac{}{}{0pt}{}{n_i}{r_i}\right)`$ if $`0n_i|m_i|2r_i`$ with integer $`r_i`$, and zero otherwise. Hence

$$I_n^{\text{hc}}(𝐑_m)=(t)^nn!\underset{r_1,\mathrm{},r_d}{}^{}\underset{i=1}{\overset{d}{}}\frac{1}{r_i!(r_i+|m_i|)!},$$
(12)   
where the sum is restricted to $`2r_1+\mathrm{}+2r_d=nM^{\text{hc}}`$, with $`M^{\text{hc}}_{i=1}^d|m_i|`$. Note that $`I_n^{\text{hc}}`$ vanishes if $`M^{\text{hc}}>n`$ or if $`n+M^{\text{hc}}`$ is odd. For the bcc lattice in $`d=3`$ we use the hypercubic lattice basis $`𝐑_m`$ and the dispersion given in Sec. III. We obtain

$$I_n^{\text{bcc}}(𝐑_m)=(t)^n\underset{i=1}{\overset{3}{}}\left[\underset{\pi }{\overset{\pi }{}}\frac{dk_i}{2\pi }(2\mathrm{cos}k_i)^n\mathrm{cos}(k_i|m_i|)\right]=(t)^n\underset{i=1}{\overset{3}{}}\left(\genfrac{}{}{0pt}{}{n}{r_i}\right),$$
(13)   
with $`r_i`$ defined as above. In particular $`I_n^{\text{bcc}}`$ vanishes if $`M^{\text{bcc}}\text{max}_i|m_i|>n`$ or if $`n+|m_i|`$ is odd.

The advantage of the present high temperature expansion is that terms with $`M>n`$ vanish when expanding to order $`\beta ^n`$, so that the lattice sum terminates at $`m_i=n`$. In addition, lattice symmetries can be used to further reduce the computational effort. This makes the calculation feasible even in dimension $`d=3`$ where numerical Brillouin zone integration is usually difficult. At the high-temperature crossing point $`O(\beta ^{80})`$ is typically sufficient to obtain $`T_+`$ and $`C_+`$ to an accuracy of $`10^8`$.