Negative electron compressibility in the Hubbard Model

TL;DR

This study investigates negative electron compressibility in the 2D Hubbard model using strong coupling diagrams, revealing phase separation near certain chemical potentials at low temperatures.

Contribution

It introduces a detailed diagrammatic analysis of electron interactions in the Hubbard model, highlighting negative compressibility and phase separation regions.

Findings

Negative compressibility occurs for μ<0 and μ>U at low temperatures.

States with phase separation emerge near these regions at even lower temperatures.

The strong coupling diagram technique effectively captures electron fluctuation effects.

Abstract

Using the strong coupling diagram technique, we study the one-band repulsive Hubbard model on a two-dimensional square lattice in a wide range of chemical potentials $\mu$. Infinite sequences of diagrams describing interactions of electrons with spin and charge fluctuations are taken into account. At low temperatures, regions of the negative electron compressibility are found for $\mu<0$ and $\mu>U$, where $U$ is the Hubbard repulsion. For even lower temperatures, in the neighborhood of these regions, states with the phase separation are revealed.

Equations

$$H=\sum _{\mathbf{l}{\mathbf{l}}^{\prime }\sigma }{t}_{\mathbf{l}{\mathbf{l}}^{\prime }}{a}_{{\mathbf{l}}^{\prime }\sigma }^{†}{a}_{\mathbf{l}\sigma }+\frac{U}{2}\sum _{\mathbf{l}\sigma }{n}_{\mathbf{l}\sigma }{n}_{\mathbf{l},-\sigma },$$

$$G({\mathbf{l}}^{\prime },{\tau }^{\prime };\mathbf{l},\tau )=⟨\mathcal{T}{\bar{a}}_{{\mathbf{l}}^{\prime }\sigma }({\tau }^{\prime })a_{\mathbf{l}\sigma }(\tau )⟩$$

$${\bar{a}}_{\mathbf{l}\sigma }(\tau )={\mathrm{e}}^{\mathcal{H}\tau }{a}_{\mathbf{l}\sigma }^{†}{\mathrm{e}}^{-\mathcal{H}\tau }.$$

$$G(\mathbf{k},j)={\left\{{\left[K(\mathbf{k},j)\right]}^{-1}-t_{\mathbf{k}}\right\}}^{-1},$$

$$\displaystyle C_{1}(\tau^{\prime},\tau)=\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma}(\tau^{\prime})a_{{\bf l}\sigma}(\tau)\big{\rangle}_{0}\,,$$

$$\displaystyle C_{2}(\tau_{1},\sigma_{1};\tau_{2},\sigma_{2};\tau_{3},\sigma_{3};\tau_{4},\sigma_{4})=\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma_{1}}(\tau_{1})a_{{\bf l}\sigma_{2}}(\tau_{2})\bar{a}_{{\bf l}\sigma_{3}}(\tau_{3})a_{{\bf l}\sigma_{4}}(\tau_{4})\big{\rangle}_{0}$$

$$\displaystyle\quad\quad-\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma_{1}}(\tau_{1})a_{{\bf l}\sigma_{2}}(\tau_{2})\big{\rangle}_{0}\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma_{3}}(\tau_{3})a_{{\bf l}\sigma_{4}}(\tau_{4})\big{\rangle}_{0}+\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma_{1}}(\tau_{1})a_{{\bf l}\sigma_{4}}(\tau_{4})\big{\rangle}_{0}\big{\langle}{\cal T}\bar{a}_{{\bf l}\sigma_{3}}(\tau_{3})a_{{\bf l}\sigma_{2}}(\tau_{2})\big{\rangle}_{0}\,,$$

$${\mathcal{H}}_{\mathbf{l}}=\sum _{\sigma }\left[(U/2)n_{\mathbf{l}\sigma }{n}_{\mathbf{l},-\sigma }-\mu n_{\mathbf{l}\sigma }\right].$$

$$\theta (\mathbf{k},j)=t_{\mathbf{k}}+t_{\mathbf{k}}^{2}G(\mathbf{k},j).$$

$$K(\mathbf{k},j)$$

$$V_{s\mathbf{k}}(j+\nu ,\sigma ;j,\sigma ;j^{\prime },-\sigma ;j^{\prime }+\nu ,-\sigma )=C_2(j+\nu ,\sigma ;j,\sigma ;j^{\prime },-\sigma ;j^{\prime }+\nu ,-\sigma )$$

$$+T\sum _{{\nu }^{\prime }}{C}_2(j+\nu ,\sigma ;j+{\nu }^{\prime },\sigma ;j^{\prime }+{\nu }^{\prime },-\sigma ;j^{\prime }+\nu ,-\sigma ){\mathcal{T}}_{\mathbf{k}}(j+{\nu }^{\prime },j^{\prime }+{\nu }^{\prime })$$

$$\times V_{s\mathbf{k}}(j+{\nu }^{\prime },\sigma ;j,\sigma ;j^{\prime },-\sigma ;j^{\prime }+{\nu }^{\prime },-\sigma ),$$

$$V_{g\mathbf{k}}(j+\nu ,{\sigma }^{\prime };j,\sigma ;j^{\prime },\sigma ;j^{\prime }+\nu ,{\sigma }^{\prime })=C_2(j+\nu ,{\sigma }^{\prime };j,\sigma ;j^{\prime },\sigma ;j^{\prime }+\nu ,{\sigma }^{\prime })$$

$$+T\sum _{{\nu }^{\prime }{\sigma }^{\prime \prime }}{C}_2(j+\nu ,{\sigma }^{\prime };j+{\nu }^{\prime },{\sigma }^{\prime \prime };j^{\prime }+{\nu }^{\prime },{\sigma }^{\prime \prime };j^{\prime }+\nu ,{\sigma }^{\prime }){\mathcal{T}}_{\mathbf{k}}(j+{\nu }^{\prime },j^{\prime }+{\nu }^{\prime })$$

$$\times V_{g\mathbf{k}}(j+{\nu }^{\prime },{\sigma }^{\prime \prime };j,\sigma ;j^{\prime },\sigma ;j^{\prime }+{\nu }^{\prime },{\sigma }^{\prime \prime }).$$

$$C_1(j)=Z^{-1}\left[\left({\mathrm{e}}^{-\beta E_1}+{\mathrm{e}}^{-\beta E_0}\right)g_{01}(j)+\left({\mathrm{e}}^{-\beta E_1}+{\mathrm{e}}^{-\beta E_2}\right)g_{12}(j)\right],$$

$$C_2(j+\nu ,{\sigma }^{\prime };j,\sigma ;j^{\prime },\sigma ;j^{\prime }+\nu ,{\sigma }^{\prime })=$$

$$-\left[Z^{-1}{\mathrm{e}}^{-\beta E_1}\left({\delta }_{j{j}^{\prime }}{\delta }_{\sigma {\sigma }^{\prime }}-{\delta }_{\nu 0}\right)+Z^{-2}\left({\mathrm{e}}^{-\beta (E_0+E_2)}-{\mathrm{e}}^{-2\beta E_1}\right)\left({\delta }_{j{j}^{\prime }}-{\delta }_{\nu 0}{\delta }_{\sigma {\sigma }^{\prime }}\right)\right]{\beta }^{\prime }F(j)F(j^{\prime }+\nu )$$

$$+Z^{-1}{\mathrm{e}}^{-\beta E_0}{\delta }_{\sigma ,-{\sigma }^{\prime }}U{g}_{01}(j)g_{01}(j^{\prime }+\nu )g_{02}(j,j^{\prime }+\nu )\left[g_{01}(j+\nu )+g_{01}(j^{\prime })\right]$$

$$+Z^{-1}{\mathrm{e}}^{-\beta E_2}{\delta }_{\sigma ,-{\sigma }^{\prime }}U{g}_{12}(j)g_{12}(j^{\prime }+\nu )g_{02}(j,j^{\prime }+\nu )\left[g_{12}(j+\nu )+g_{12}(j^{\prime })\right]$$

$$\displaystyle\quad-Z^{-1}{\rm e}^{-\beta E_{1}}\delta_{\sigma,-\sigma^{\prime}}\big{\{}F(j^{\prime}+\nu)\big{[}g_{01}(j)g_{01}(j^{\prime})+g_{12}(j)g_{12}(j+\nu)-g_{01}(j^{\prime})g_{12}(j+\nu)\big{]}$$

$$\displaystyle\quad+F(j)\big{[}g_{01}(j+\nu)g_{01}(j^{\prime}+\nu)+g_{12}(j^{\prime})g_{12}(j^{\prime}+\nu)-g_{01}(j+\nu)g_{12}(j^{\prime})\big{]}\big{\}},$$

$$K(j)=C_1(j)-\frac{T}{N}\sum _{\mathbf{k}{j}^{\prime }{\sigma }^{\prime }}{C}_2(j,\sigma ;j,\sigma ;j^{\prime },{\sigma }^{\prime };j^{\prime },{\sigma }^{\prime })\theta (\mathbf{k},j^{\prime }).$$

Full text

0.1 Introduction

The electron compressibility is defined as $\kappa=\bar{n}^{-2}({\rm d}\bar{n}/{\rm d}\mu)$, where $\bar{n}$ is the electron concentration and $\mu$ the chemical potential. Positive values of this quantity correspond to the common situation, in which an increase of the concentration leads to a growth of the chemical potential. However, cases of negative $\kappa$ are also possible, when an increase of $\bar{n}$ is accompanied with a counterintuitive decrease of $\mu$. It is a clear indication of non-rigid bands that is the characteristic property of strong electron correlations. Usually, the negative electron compressibility (NEC) is associated with thermodynamic instability and phase separation (see, e.g., Dezi ). This happens when the dependence $\bar{n}(\mu)$ is multivalued in some range of $\mu$ that is an indication of a mixture of two phases. However, there are also cases, in which NEC does not lead to phase separation (see, e.g., Riley ), and one of these cases will be considered below. Most frequently NEC is observed in low-dimensional dilute electron systems Dezi ; Riley ; Kravchenko ; Eisenstein . Recently NEC was also detected in the quasi-three-dimensional crystal (Sr1-xLax)3Ir2O7 at higher electron concentrations He .

In the one-band repulsive Hubbard model on a two-dimensional (2D) square lattice, NEC is observed already in the Hubbard-I approximation Hubbard63 ; Hubbard64 as soon as $\mu$ leaves the usually considered range $0<\mu<U$ and the temperature $T$ is low enough. Here $U$ is the on-site Coulomb repulsion. The origin of the NEC is the discreteness of the energy levels in the Hubbard atom and their crossing, which occurs when $\mu$ traverses the energies 0 and $U$. Due to occupation-probability factors the first-order cumulant of electron operators, which determines Green’s function in the Hubbard-I approximation, changes drastically in the ranges $-T\lesssim\mu\lesssim T$ and $U-T\lesssim\mu\lesssim U+T$. This alteration is connected with the change of the lowest state – for $\mu\ll-T$ and $\mu\gg U+T$ it is empty and doubly occupied states, while for $T\ll\mu\ll U-T$ this state is singly occupied. For $U\gg t$, where $t$ is the hopping constant, the former two states are associated with weakly correlated bands of the Hubbard model, whereas the latter state is responsible for the case of strong correlations. Dependencies $\bar{n}(\mu)$ in these three cases are markedly different that leads to the NEC in transition regions $\mu\approx 0$ and $\mu\approx U$.

Of course, the Hubbard-I approximation gives an oversimplified picture. It is necessary to prove whether these NEC regions survive if charge and spin fluctuations are taken into account. This question is considered in the present work. To answer it we use the strong coupling diagram technique (SCDT) Vladimir ; Metzner ; Pairault ; Sherman16 ; Sherman18 ; Sherman19a ; Sherman19b . In this approach, Green’s functions are calculated using series expansions in powers of the electron kinetic energy. Terms of the series are products of hopping constants and on-site cumulants of electron operators. The linked-cluster theorem is valid and partial summations are allowed in the SCDT (concise description of the approach can be found in Sherman16 ). In SCDT, the interaction of electrons with charge and spin fluctuations is described by diagrams with ladder inserts. If these ladders are constructed from renormalized hopping lines and second-order cumulants, ladders of all lengths can be summed Sherman18 . As a consequence, fluctuations of all ranges are taken into account in an infinite crystal. In Sherman18 ; Sherman19a , it was shown that the spectral function, magnetic susceptibility, double occupancy and squared site spin calculated in the SCDT are in satisfactory agreement with results of exact diagonalizations, Monte Carlo simulations, numerical linked-cluster expansions and experiments with ultracold fermionic atoms in 2D optical lattices in wide ranges of repulsions, temperatures and concentrations. Below this approach is applied for the consideration of the influence of charge and spin fluctuations on the NEC. We find that the fluctuations do not destroy it and the NEC is observed at low temperatures. To attain even lower $T$, we reduce the set of considered diagrams. In this approximation, regions of bistability or phase separation are found near the NEC domains.

0.2 Model and SCDT method

The fermionic 2D one-band Hubbard model Hubbard63 ; Hubbard64 is used in this work. Its Hamiltonian reads

[TABLE]

where 2D vectors ${\bf l}$ and ${\bf l^{\prime}}$ label sites of a square plane lattice, $\sigma=\uparrow,\downarrow$ is the spin projection, $a^{\dagger}_{{\bf l}\sigma}$ and $a_{{\bf l}\sigma}$ are electron creation and annihilation operators, $t_{\bf ll^{\prime}}$ is the hopping constant and $n_{{\bf l}\sigma}=a^{\dagger}_{{\bf l}\sigma}a_{{\bf l}\sigma}$. In this work, only the nearest neighbor hopping constant $t$ is taken to be nonzero.

In this work, we calculate the dependence $\bar{n}(\mu)$ for the Hamiltonian (1). Here $\bar{n}=2\langle n_{{\bf l}\sigma}\rangle$, angular brackets denote the statistical averaging with the operator ${\cal H}=H-\mu\sum_{{\bf l}\sigma}n_{{\bf l}\sigma}$ and the factor 2 is connected with two spin states. To find this dependence the one-particle Green’s function

[TABLE]

is calculated using the SCDT. Here ${\cal T}$ is the chronological operator arranging operators from right to left in ascending order of times $\tau$ and time dependencies are determined by ${\cal H}$,

[TABLE]

As in the weak coupling diagram technique Mahan , diagrams of the SCDT describing Green’s function (2) can be classed into reducible and irreducible. The latter cannot be divided into two disconnected parts by cutting a hopping line. The sum of all such diagrams is termed the irreducible part $K$. The Fourier transform of Green’s function is expressed through this quantity as

[TABLE]

where k is a 2D wave vector, $j$ is an integer defining the Matsubara frequency $\omega_{j}=(2j-1)\pi T$ and $t_{\bf k}$ is the Fourier transform of $t_{\bf ll^{\prime}}$. In the considered case, diagrams in the particle-hole channel make the main contribution to $K$ Bychkov ; Bickers . These diagrams are shown in Figure 1.

In this figure, circles without letters are cumulants of the first and second orders

[TABLE]

where the subscript 0 near brackets indicates that time dependencies and averages are determined by the Hamiltonian of the Hubbard atom

[TABLE]

Due to the translation symmetry, the cumulants are identical on all lattice sites. The last two terms in $K({\bf k},j)$ correspond to diagrams with ladder inserts and describe the interaction of electrons with charge and spin fluctuations. The vertices $V_{s}$ and $V_{g}$ satisfy the Bethe-Salpeter equations (BSE) shown in Figures 1(b) and (c), in which circles with letters are irreducible four-leg diagrams which cannot be divided into two disconnected parts by cutting a pair of horizontal oppositely directed hopping lines. Thanks to the possibility of the partial summation, bare internal lines $t_{\bf k}$ in Figure 1 can be transformed into dressed ones,

[TABLE]

In this work, the irreducible four-leg diagrams in Figures 1(b) and (c) are approximated by the respective second-order cumulants. As mentioned in the Introduction, this approach gave results in a satisfactory agreement with reliable calculation and experimental methods, at least in the range of chemical potentials $T\ll\mu$, $T\ll U-\mu$. Besides, magnetic susceptibilities obtained in this approximation were compared with those calculated in a more complex approach, in which the longitudinally irreducible vertices in Figures 1(b) and (c) were approximated by full sums of transversal ladders constructed from dressed hopping lines and second-order cumulants. The susceptibilities appeared to be close in these two approaches in the mentioned range of $\mu$. This gives grounds to use the simpler approximation, with second-order cumulants as irreducible vertices, for investigating a wider range of chemical potentials. In this case, $V_{s}$ and $V_{g}$, apart from frequencies, depend only on the transfer momentum. The irreducible part and BSEs read

[TABLE]

[TABLE]

In Equations (6)–(8), cumulants of the first and second orders are given by the equations Vladimir ; Sherman06 ; Sherman08

[TABLE]

where $E_{0}=0$, $E_{1}=-\mu$ and $E_{2}=U-2\mu$ are eigenenergies of the Hamiltonian (4) corresponding to zero, one and two electrons on the Hubbard atom, $Z={\rm e}^{-\beta E_{0}}+2{\rm e}^{-\beta E_{1}}+{\rm e}^{-\beta E_{2}}$ is the respective partition function, $g_{nm}(j)=({\rm i}\omega_{j}+E_{n}-E_{m})^{-1}$, $g_{02}(j,j^{\prime})=({\rm i}\omega_{j}+{\rm i}\omega_{j^{\prime}}+E_{0}-E_{2})^{-1}$, $F(j)=g_{01}(j)-g_{12}(j)$, $\beta=1/T$ and ${\cal T}_{\bf k}(j,j^{\prime})=N^{-1}\sum_{\bf k^{\prime}}\theta({\bf k+k^{\prime}},j)\theta({\bf k^{\prime}},j^{\prime})$ with $N$ the number of sites.

The above equations describe the transition to the long-range antiferromagnetic order Sherman18 . At the derivation of the equation (10) $\beta^{\prime}=\beta$ Vladimir ; Sherman06 ; Sherman08 . However, for such a $\beta^{\prime}$ the antiferromagnetic transition occurs at a small but finite value of $T=\zeta$. In the considered 2D system with the continuous symmetry, this result is in contradiction with the Mermin-Wagner theorem Mermin . Hence the used approximation – with the second-order cumulants as longitudinally irreducible four-leg vertices – overestimates somewhat spin fluctuations. We note that the above-mentioned improvement – the use of the full sums of vertical ladders as such vertices – gives only a slight betterment. To remedy this defect and shift the transition temperature to zero $\beta^{\prime}$ in (10) was set equal to $1/(T+\zeta)$, and the value of $\zeta$ was determined in the case of the strongest antiferromagnetic fluctuations – at half-filling, $\mu=U/2$. Obtained values of $\zeta$ depend on $U/t$ – for $U\gtrsim 4t$ it equals to $0.24t$, while for $U=t$ the parameter $\zeta=0.05t$. It may be interesting to note that the value of $\zeta$ is close to $T_{\rm spin}$, which was associated in Paiva with the maximum in the temperature dependence of the uniform susceptibility. This maximum is a manifestation of the change in the character of long-range spin correlation with temperature.

Expressions (3)–(10) form a closed set of equations, which can be solved by iteration for given values of $U$, $T$, $\mu$ and the considered set of hopping constants $t_{\bf ll^{\prime}}$. As a starting function in this iteration, the Hubbard-I solution was used. This solution is obtained from the above formulas if the irreducible part in (3) is approximated by the first term in the right-hand side of (6) – the first-order cumulant $C_{1}(j)$. No artificial broadening was introduced in these calculations.

Expressions (9) and (10) are considerably simplified for chemical potentials satisfying the conditions $T\ll\mu$ and $T\ll U-\mu$. For $T\ll U$ these conditions correspond to cases of half-filling and moderate doping, which are usually considered in works on the Hubbard model. For this range of $\mu$, occupation-probability factors ${\rm e}^{-\beta E_{0}}$ and ${\rm e}^{-\beta E_{2}}$ are exponentially small in comparison with ${\rm e}^{-\beta E_{1}}$ and terms containing them can be dropped from the formulas. As was shown in Sherman18 , this leads to the considerable simplification of the BSEs (7) and (8) reducing them to small systems of linear equations. In this work, chemical potentials outside of this range are also considered. Therefore, the mentioned simplification is inapplicable. Instead of (8), it is convenient to consider equations for combinations of $V_{g}$, which are symmetrized and antisymmetrized over spin variables. The antisymmetrized combination coincides with $V_{s}$, and both combinations do not depend on the sign of spin projection. BSEs for them are systems of linear equations with matrix indices $\nu$ and $\nu^{\prime}$, which depend parametrically on $j$, $j^{\prime}$ and k. In our calculations, the rank of these systems can range up to 800. Therefore, to decrease the number of numeric operations the iterative solution of these systems was used. The iterations started from the respective values of $C_{2}$. For considered parameters 10–15 iterations were enough to achieve convergence. These iterations were performed on each step of the general iterative solution of equations (3)–(10). In these external iterations, convergence was achieved in 10–40 steps.

0.3 Results and discussion

An example of results obtained using the above-discussed procedure is shown in Figure 2.

The figure demonstrates the temperature and filling evolution of the density of electron states (DOS) $\rho(\omega)=-(N\pi)^{-1}\sum_{\bf k}{\rm Im}G({\bf k},\omega)$. The analytic continuation from the imaginary to real frequency axis was carried out using the maximum entropy method Press ; Jarrell ; Habershon . The electron concentrations were estimated using the equation $\bar{n}=(2T/N)\sum_{{\bf k}j}{\rm e}^{{\rm i}\omega_{j}\eta}G({\bf k},j)$ with $\eta\rightarrow+0$. Using formulas of the previous section, $G({\bf k},j)$ was calculated for $j\leq 100$. Values of Green’s function for larger $j$, which are required for evaluating $\bar{n}$ from this equation, were approximated with the asymptotics Vilk $G({\bf k},j)\rightarrow 1/({\rm i}\omega_{j})$. It was found that calculated Green’s functions settle into this asymptote already at $j\ll 100$. Values of the electron concentration obtained from analytically continued $\rho(\omega)$ differ from these $\bar{n}$ by a few percent.

Due to the particle-hole symmetry of the considered model values of $\mu\leq U/2$ only will be considered below. DOSs calculated for $\mu\gg T$ are similar to those obtained earlier with a somewhat different calculation procedure Sherman18 ; Sherman19a ; Sherman19b . For half-filling, the Mott gap is seen around the Fermi level located at $\omega=0$ (Figures 2(a) and (b)). Suppressions of the DOS near $\omega=\pm U/2$ are connected with the strong reabsorption of electrons on the transfer frequencies of the Hubbard atom $\omega=-\mu$ and $U-\mu$ Sherman16 . With doping, for the higher temperature, a pseudogap appears around the Fermi level (Figures 2(c) and (e)). For the lower temperature, the Fermi-level peak arises inside the pseudogap (Figures 2(d) and (f)). As discussed in Sherman19b , this peak is a manifestation of bound states of correlated electrons and mobile spin excitations. This peak is a 2D analog of the infinite-dimensional resonance peak of the dynamic mean-field approximation Georges and the Abrikosov-Suhl resonance of the Anderson impurity model Hewson . With further doping, the contribution of the upper Hubbard band decreases considerably and the DOS acquires some features of the weakly correlated spectrum (Figures 2(g) and (h)).

The dependence $\bar{n}$ vs. $\mu$ for $U=8t$ and $T=0.32t$ is shown in Figure 3. It has the usual shape – the concentration decreases monotonously with a fall in the chemical potential. However, already at this comparatively high temperature, some change in the dependence can be noticed for $\mu<0$.

The dependence changes radically with temperature lowering.

As seen in Figure 4, there appears a NEC region in the range $-0.4t\lesssim\mu\lesssim 0$. The mechanism of its occurrence can be understood from the expression for $C_{1}$ (9). Due to the discreteness of the energy levels of the Hubbard atom (4), there are three regions of $\mu$, which correspond to its different fillings. In the first region, $\mu\ll-T$, $C_{1}(j)\approx 1/({\rm i}\omega_{j}+\mu)$, the atom is empty, and the solution of Hubbard model describes a gas of weakly correlated electrons. The second region, $T\ll\mu$ and $T\ll U-\mu$, $C_{1}(j)\approx(1/2)[1/({\rm i}\omega_{j}+\mu)+1/({\rm i}\omega_{j}+\mu-U)]$, relates to a lattice with predominantly singly occupied sites and strong electron correlations. The third region, $T\ll\mu-U$, $C_{1}(j)\approx 1/({\rm i}\omega_{j}+\mu-U)$ corresponds to a lattice with predominantly doubly occupied sites, in which single occupancies behave like a gas of weakly correlated species. It is clear that dependencies $\bar{n}(\mu)$ are different in these three situations that leads to NEC in transient regions $-T\lesssim\mu\lesssim T$ and $-T\lesssim U-\mu\lesssim T$. This simplified picture based on the Hubbard-I approximation is not modified considerably by charge and spin fluctuations, as can be seen in Figure 4.

With a further decrease of temperature, the NEC regions become narrower and the amplitude of changes in $\bar{n}$ gets apparently larger. Calculations for $T<0.1t$ need a larger number of frequency points and become too time and memory consuming for the used set of diagrams. However, for parameters corresponding to cuprates, $t\approx 0.2$ eV (see, e.g., Sherman18 ) and $T=0.1t\approx 200$ K. Therefore, temperatures in the range of $0.01t-0.03t$ are closer to the usual experimental conditions. To attain this temperature region the irreducible part was simplified and approximated by the two lowest-order diagrams

[TABLE]

In this approximation, $K$ does not depend on momentum. The dependence $\bar{n}(\mu)$ calculated using (11) looks similar to those shown in Figures 3 and 4 for the same temperatures. However, for lower temperatures, a new feature appears in this dependence – a narrow region of bistability in the vicinity of the NEC.

In this region, there are two values of the electron concentration for each chemical potential. Two branches of the dependence in Figure 5 were obtained by taking a solution for some $\mu$ as a starting one in the iteration procedure, which was carried for a slightly larger or smaller chemical potential. For higher $T$ no bistability is observed.

The coexistence of two solutions with different $\bar{n}$ for a given value of the chemical potential points to the phase separation. Though this result has to be verified in a more elaborate approximation with a larger number of processes in $K$, it indicates the possibility of a purely electronic mechanism for the phase separation in crystals with strong electron correlations. The phase separation is ubiquitous in such compounds (see, e.g., Tranquada ; Vershinin ; Chu ; Leshen ). However, underlying mechanisms of its creation are not well understood and strongly debated (see, e.g., Sherman08 ; Hizhnyakov ; Sboychakov ; Aichhorn ; Otsuki ). In this connection, it is worth noting that in the range of chemical potentials $T\ll\mu$, $T\ll U-\mu$ no indication of the divergence in the charge susceptibility were found using the SCDT Sherman18 ; Sherman19b . This fact excludes a purely electronic mechanism of the phase separation in this range of $\mu$ in the considered model. However, outside of the range, such a mechanism may exist, as seen in Figure 5.

0.4 Concluding remarks

In this work, the dependence of the hole concentration on the chemical potential was studied in the one-band repulsive Hubbard model on the two-dimensional square lattice in the case of strong correlations. Regions $\mu\leq 0$ and $\mu\geq U$ were of particular interest since due to the discreteness of the spectrum of the Hubbard atom and level crossing the negative charge compressibility was expected here. In the calculations of this work, the strong coupling diagram technique was used. Infinite sequences of diagrams describing interactions of electrons with spin and charge fluctuations were taken into account. This was done with the aim to prove the robustness of the NEC with respect to these fluctuations. The NEC regions were really found near $\mu=0$ and $U$ for low temperatures. For even lower temperatures, in the neighborhood of the NEC, there appear the bistability regions, in which two values of the electron concentration correspond to each value of the chemical potential. This is an indication of the phase separation in these regions. However, it should be noted that this result was obtained without regard for charge and spin fluctuations and needs in further verification. Due to more complex expressions for the second-order cumulants, a new technique was used for the solution of Bethe-Salpeter equations for charge and spin vertices. In the range $T\ll\mu$, $T\ll U-\mu$ obtained densities of states were similar to those derived in earlier works with the other method. For strong correlations, half-filling and moderate temperatures the DOS demonstrates the Mott gap near the Fermi level and reabsorption pseudogaps near $\pm U/2$. With doping, a pseudogap appears near the Fermi level at moderate temperatures. With the temperature decrease, the sharp peak develops at or in the nearest vicinity of the Fermi level. This peak is a manifestation of the bound states of correlated electrons and mobile spin excitations. It is the two-dimensional analog of the infinite-dimensional resonance peak of the dynamic mean-field approximation and the Abrikosov-Suhl resonance of the Anderson impurity model. For $\mu\lesssim 0$ and $\mu\gtrsim U$ the DOS acquires some features of a weakly correlated spectrum.

Index