Probing Phase Coherence Via Density of States for Strongly Correlated Excitons

Abstract

We present the calculation of the coherent spectral functions and density of states (DOS) for excitonic systems in the frame of the three-dimensional extended Falicov–Kimball model. Using gage-invariant U(1) transformation to the usual fermions, we represent the electron operator as a fermion attached to the U(1) phase-flux tube. The emergent bosonic gage field, related to the phase variables, is crucial for the Bose–Einstein condensation (BEC) of excitons. Employing the path-integral formalism, we manipulate the bosonic and fermionic degrees of freedom to obtain the effective actions related to fermionic and bosonic sectors. Considering the normal and anomalous excitonic Green functions, we calculate the spectral functions, which have the forms of convolutions in the reciprocal space between bosonic and fermionic counterparts. For the fermionic incoherent part of the DOS, we have found the strong evidence of the hybridization gap in DOS spectra. Furthermore, considering Bogoliubov coherence mechanism, we calculate the coherent DOS spectra. For the coherent normal fermionic DOS, there is no hybridization gap found in the system due to strong coherence effects and phase stiffness. The similar behavior is observed also for the condensate part of the anomalous excitonic DOS spectra. We show that for small values of the Coulomb interaction the fermionic DOS exhibits a Bardeen–Cooper–Schrieffer (BCS)-like double-peak structure. In the BEC region of the BCS–BEC crossover, the double-peak structure disappears totally for both: coherent and incoherent DOS spectra. We discuss also the temperature dependence of DOS functions.

1 Introduction

The attractive Coulomb interaction between conduction band electrons and valence band holes plays a crucial role in forming the bosonic electron–hole (e–h) bound pairs [1] (excitons). The Bose–Einstein condensation (BEC) of the excitons represents a very interesting subject in solid-state physics. Mostly, the e–h systems are realized in photoexcited semiconducting materials, and the properties of such systems strongly depend on the e–h density, temperature of the system, pressure, and other important physical quantities. There is a number of works, both experimental and theoretical, where all those effects are investigated intensively [1–6]. A low density system of excitons behaves like the usual Bose gas, and, at cryogenic temperatures, the BEC transition can be expected [4, 5]. On the other hand, the high density system of bound e–h pairs behaves like the system of usual Cooper pairs in superconductors [7]. In this limit, we have the Bardeen–Cooper–Schrieffer (BCS) state of e–h pairs [8, 9]. Despite many experimental investigations to observe the coherent exciton condensates [6, 10–14], there is not yet definitive evidence for such states. Therefore, an expected BCS–BEC crossover represents actually a fascinating problem typical to the excitonic systems [15–20]. Especially, it is interesting from the viewpoint of the difference from the similar crossover in superconductors or trapped atomic Fermi gases [21–24]. The transition to the e–h pair condensed phase, in the limit of weak coupling, is related to the relative motion between electrons and holes [15], implying the BCS regime and is in contrast to the case of strong coupling regime, when the BEC state is related to the motion of the center of mass of excitons [15]. In the whole BCS–BEC transition region, the e–h mass difference leads to a large suppression of the BEC transition temperature, which is proved to not be same as the temperature of excitonic pair (EP) formation, and hence the excitonic insulator (EI) phase [15]. Recently, it has been shown theoretically that the excitonic insulator and the excitonic condensate are not exactly the same states of the matter [15, 25, 26]. The author in Refs. [25] and [26] shows from general considerations that in the low density limit of the excitonic pairs, the critical temperature of excitonic condensation should be much smaller than the temperature of EP formation, in contrast to the previous treatments, [16–20] where the EI state is ad hoc associated with the BEC state of excitons as the same. Similarly, in Ref. [15], it is shown that the EI state is an excitonium state, where incoherent e–h bound pairs are formed and, furthermore, at the lower temperatures, the BEC of excitons appears in consequence of the reconfiguration and coherent condensation of the preformed excitonic pairs. Obviously, in the low density limit, the gas of free excitons undergoes the BEC phase transition at very low temperatures, and in general, the BEC temperature transition line is not coinciding with that of EP formation. The Bose condensation of excitonic pairs is possible only when the macroscopic phase coherence is attended by the system [25]. However, the experimental evidence of the existence of two distinct transition temperatures, for the general case of a three-dimensional (3D) bulk system, is yet lacking in the literature. In fact, the question, whether a true “coherent” BEC transition is present in the system of excitons (in the case of high exciton density limit), is still ambiguous. The experimental proof of it is a very cumbersome problem, because of the dominant role of quantum fluctuations at low temperatures $k_\mathrm{B}T\lesssim E^\mathrm{eh}_\mathrm{Ry}$ (with $E^\mathrm{eh}_\mathrm{Ry}$ being the binding energy of a Mott-Wannier exciton), when very large zero-point oscillations are present.

The continuing growing interest, to the problem of the coherent excitonic condensates, motivates us to calculate excitonic spectral functions and density of states (DOS), as a direct probe mechanism, to compare the results with the high-resolution studies of angle-resolved photoemission spectroscopy (ARPES) [27–30], and with spectral weight measurements on the excitonic materials at the very low temperatures.

The general strategy of our calculations is based on the effective actions method. First, we transform the initial total action of the system to a gage-invariant form, by applying the U(1) gage transformation to the fermion operators. As a result, the electron appears in the theory like a composite object of that of the fermion with the attached U(1) phase “flux-tube.” The electron factorization in terms of two variables has an unprecedented impact on the whole theory. Then we integrate out the phase variables and we get the effective fermionic action in the theory, and we derive a set of self-consistent equations for the EI state. Furthermore, we discuss shortly the results of the quantum rotor model obtained after the integration of fermionic degrees of freedom.

The path-integral formalism elaborated here permits to calculate the correlation functions in the system and, as a result, we obtain the expressions of normal (both incoherent and coherent) and anomalous excitonic spectral functions, and shapes of DOS corresponding. A special attention is paid, when calculating the bosonic phase-stiffness DOS function, which is negative. Furthermore, it is crucial for calculation of the coherent normal and anomalous excitonic DOS functions. Namely, considering the expressions of the bosonic and fermionic DOS functions, we calculate the total, phase-coherent DOS functions, as convolutions from bosonic and fermionic counterparts.

For the incoherent partial normal fermionic DOS functions at $T=0$, we obtain a hybridization gap in the excitation spectra as a direct consequence of the presence of the Hartree-type gap in the single-particle energy scales. For the case of coherent normal and anomalous excitonic DOS functions, this hybridization gap is absent totally for all frequency modes and for all values of the Coulomb interaction. This is due to strong coherence effects in the strongly correlated fermion system at low temperatures and at low densities. For the anomalous excitonic DOS function, we found that the hybridization gap is absent for the case of small and intermediate values of the Coulomb interaction parameter, but there is a finite constant small gap that is opening in the strong coupling regime, signaling the passage to the SC (BEC) side of the SC-SM phase transition (the BCS–BEC crossover).

The paper is organized as follows: in Sect. 2, we introduce the model Hamiltonian and we discuss the main calculation schemes. In Sect. 3, we describe the EI state, with a short discussion about important energy scales in the system. In Sect. 4, we present the analytical calculation of spectral functions and DOS functions of the system in consideration. At the end of Sect. 4, we show the results of numerical evaluations for calculated DOS functions.

2 The Method

2.1 The EFKM Hamiltonian

As the model for study the excitonic condensation at low temperatures, we have chosen the two-band extended Falicov–Kimball model (EFKM) [16, 31–34], due to its large applicability for treatment of the electronic correlations [37, 38]. The Hamiltonian of the EFKM model is

$$\begin{aligned} \mathcal{{H}}&= -\,t_{c}\sum _{\left\langle {\mathbf {r}},{\mathbf {r}}' \right\rangle }\left[ \bar{c}({{\mathbf {r}}})c({{\mathbf {r}}}')+h.c.\right] -\bar{\mu }\sum _{{\mathbf {r}}}n({\mathbf {r}}) \nonumber \\&-\,t_{f}\sum _{\left\langle {\mathbf {r}},{\mathbf {r}}' \right\rangle }\left[ \bar{f}({{\mathbf {r}}})f({{\mathbf {r}}}')+h.c.\right] +\frac{\epsilon _{c}-\epsilon _{f}}{2}\sum _{{\mathbf {r}}}\tilde{n}({\mathbf {r}}) \nonumber \\&+\,U\sum _{{\mathbf {r}}}\frac{1}{4}\left[ n^{2}({\mathbf {r}})-\tilde{n}^{2}({\mathbf {r}})\right] . \end{aligned}$$

(1)

Here, $\bar{f}({{\mathbf {r}}})$ ($\bar{c}({{\mathbf {r}}})$) is the $f$ ($c$) electron creation operator at the lattice position ${\mathbf {r}}$, the summation $\left\langle {\mathbf {r}}, {\mathbf {r}}' \right\rangle$ is over nearest neighbors (n.n) sites on the 3D cubic lattice. The short- hand notations are introduced $n({\mathbf {r}})=n_{c}({\mathbf {r}})+n_{f}({\mathbf {r}})$ and $\tilde{n}({\mathbf {r}})=n_{c}({\mathbf {r}})-n_{f}({\mathbf {r}})$ in order to simplify the calculations. Next, $t_{c}$ is the hopping amplitude for $c$-electrons and $\epsilon _{c}$ is the corresponding on-site energy level parameter. Similarly, $t_{f}$ is the hopping amplitude for $f$-electrons and $\epsilon _{f}$ is the on-site energy level parameter for $f$-orbital. For $t_{c}t_{f}<0$ ($t_{c}t_{f}>0$), we have a direct (indirect) band gap semiconductor. The on-site (local) Coulomb interaction $U$ in the last term of the Hamiltonian in Eq. (1) plays the coupling role between two bands. As we will see later on, the strength of the local Coulomb interaction will tune the semi-metal (SM)–semiconductor (SC) transition in the system. The chemical potential $\bar{\mu }$ is $\bar{\mu }=\mu -\bar{\epsilon }$, where $\bar{\epsilon }=\left( \epsilon _{c}+\epsilon _{f}\right) /2$. Naturally, we adjust the chemical potentials $\mu _{f}$ and $\mu _{c}$ in order to maintain the number of electrons in $f$ and $c$ orbitals separately. Then, the equilibrium value of the chemical potential $\mu \equiv \mu _{f}=\mu _{c}$ in  Eq. (1) will be determined from the half-filling condition, i.e., we suppose that $\left\langle n_{c}({\mathbf {r}})\right\rangle +\left\langle n_{f}({\mathbf {r}})\right\rangle =1$. We will use $t_{c}=1$ as the unit of energy, and we fix the band parameter values $\epsilon _{c}=0$ and $\epsilon _{f}=-1$. The Fermi energy level is assumed to be situated at the level of the $c$-band, thus $\epsilon _{F}=0$. For the $f$-band hopping amplitude $t_{f}$, we consider the values $t_{f}=-0.3$ and $t_{f}=-0.1$ corresponding to the heavy hole and light hole $f$-bands. Throughout the paper, we set $k_\mathrm{B}=1$, $\hbar =1$, and lattice constant $a=1$. Also, we keep the frequency symbol $\nu$ for fermions and $\omega$ for bosons, throughout the paper. In the case of degenerated $f$ and $c$-bands, i.e., when $\epsilon _{f}=\epsilon _{c}$ and $t_{f}=t_{c}$, the EFKM model reduces to the standard Hubbard model [35]. The principal advantage of the EFKM, in comparison with the genuine Falicov–Kimball model (FKM) [36], is that it is taking into account the direct nearest neighbors $f$ electron hoppings [37, 38] ($t_{f}$) and it can be shown [35] that the EI state is unstable when the pure FKM is approached.

In fact, the EFKM Hamiltonian in Eq. (1) is equivalent to the asymmetric Hubbard model, if we associate to the orbitals $c$ and $f$ the spin variables, by replacing the fermionic Hilbert space with the pseudo-fermionic one and then linearizing the interaction term via the bosonic states (see in Ref. [17]).

2.1.1 The Partition Function

The Hamiltonian in Eq. (1) is containing two separate quadratic terms and is suitable for decoupling by functional path integration method [39, 40]. We employ the imaginary-time, fermionic path-integral method, and we introduce the fermionic Grassmann variables [39] ${f}({{\mathbf {r}}}\tau )$ and ${c}({{\mathbf {r}}}\tau )$ at each site ${\mathbf {r}}$ and each time $\tau$ varying in the interval $0\le \tau \le \beta$, where $\beta =1/T$ with $T$ being the thermodynamic temperature [41]. The time-dependent variables ${f}({{\mathbf {r}}}\tau )$ and ${c}({{\mathbf {r}}}\tau )$ are satisfying the anti-periodic boundary conditions ${x}({{\mathbf {r}}}\tau )=-{x}({{\mathbf {r}}}\tau +\beta )$, where $x=f$ or $c$.

The grand canonical partition function of system of fermions written as a functional integral over the Grassmann fields is

$$\begin{aligned} Z_\mathrm{GC}=\int \left[ {D}\bar{c}{D}c\right] \left[ {D}\bar{f}{D}f\right] e^{-\mathcal{{S}}[\bar{c},c, \bar{f},f]}, \end{aligned}$$

(2)

where the action in the exponent is given in the path-integral formulation in the form

$$\begin{aligned} \mathcal{{S}}[\bar{c},c, \bar{f},f]=\sum _{x=f,c}\mathcal{{S}}_\mathrm{B}[\bar{x},x]+\int ^{\beta }_{0}\mathrm{d}\tau {{H}}(\tau ). \end{aligned}$$

(3)

Here, $\mathcal{{S}}_\mathrm{B}[\bar{f},f]$ and $\mathcal{{S}}_\mathrm{B}[\bar{c},c]$ are fermionic Berry terms. They are defined as follows:

$$\begin{aligned} \mathcal{{S}}_\mathrm{B}[\bar{x},x]=\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \bar{x}({\mathbf {r}}\tau )\dot{x}({\mathbf {r}}\tau ), \end{aligned}$$

(4)

and $\dot{x}({\mathbf {r}}\tau )=\partial _{\tau }x({\mathbf {r}}\tau )$ is the time derivative.

2.1.2 Decoupling of Term Proportional to $n^{2}$

We will decouple the quadratic density terms in Eq. (1) using the Hubbard–Stratonovich linearisation procedure [40] and by introducing new variables $V({\mathbf {r}}\tau )$ and $\mathcal{{ \varrho }}({\mathbf {r}}\tau )$, conjugated, respectively, to the density terms $n({\mathbf {r}}\tau )$ and $\tilde{n}({\mathbf {r}}\tau )$. For the quadratic term, proportional to $n^{2}({\mathbf {r}}\tau )$ in the exponent of the partition function in Eq. (2), we have

$$\begin{aligned}&\exp \left[ {-\frac{U}{4}\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau n^{2}\left( {\mathbf {r}}\tau \right) }\right] \nonumber \\&\quad =\int \left[ {D}V\right] e^{-\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ \frac{V^{2}({\mathbf {r}}\tau )}{U}-iV({\mathbf {r}}\tau )n({\mathbf {r}}\tau )\right] }. \end{aligned}$$

(5)

Furthermore, we combine the exponent in Eq. (5) with the effective chemical potential term, linear in total electron density $n({\mathbf {r}})$ (see the second term, in Eq. 1). Then, we decompose the variables $V({\mathbf {r}}\tau )$ into a static and the periodic part

$$\begin{aligned} V({\mathbf {r}}\tau )=V_{0}({\mathbf {r}})+\tilde{V}({\mathbf {r}}\tau ), \end{aligned}$$

(6)

where, from time-periodicity of $\tilde{V}({\mathbf {r}}\tau )$, it follows the relation $\int ^{\beta }_{0}\mathrm{d}\tau \tilde{V}({\mathbf {r}}\tau )=0$. As a result, the integration over variables $V({\mathbf {r}}\tau )$ becomes the integration over the scalar static variables $V_{0}({\mathbf {r}})$ and the integration over the periodic field $\tilde{V}({\mathbf {r}}\tau )$:

$$\begin{aligned} \int \left[ {D}V\right] \cdots =\int \left[ {D}V_{0}\right] \int \left[ {D}\tilde{V}\right] \cdots \end{aligned}$$

(7)

For the periodic part in Eq. (6), we introduce U$(1)$ phase-field variables $\phi ({\mathbf {r}}\tau )$ using a Faraday-type relation [42]

$$\begin{aligned} \tilde{V}({\mathbf {r}}\tau )=\frac{\partial {\phi }({{\mathbf {r}}}\tau )}{\partial {\tau }}\equiv \dot{\phi }({\mathbf {r}}\tau ). \end{aligned}$$

(8)

Thus, for the dynamic part, we transform the integration over the gage variables $\tilde{{V}}({\mathbf {r}}\tau )$, into an integration over generic phase variables $\phi ({\mathbf {r}}\tau )$

$$\begin{aligned} \int \left[ \mathcal{{D}}\tilde{V}\right] \cdots \rightarrow \int \left[ \mathcal{{D}}\phi \right] \cdots \end{aligned}$$

(9)

The periodicity of $\tilde{V}\left( {\mathbf {r}}\tau \right)$ implies that $\phi \left( {\mathbf {r}}\beta \right) =\phi \left( {\mathbf {r}}0\right)$. The integration measure, in Eq. (9), over the variables $\phi$ is defined as

$$\begin{aligned} \int \left[ \mathcal{{D}}\phi \right] \cdots&\equiv \int ^{\infty }_{-\infty }\prod _{{\mathbf {r}}}\mathrm{d}\phi _{0}({\mathbf {r}}) \nonumber \\&\times \int ^{\phi _{f}=\phi ({\mathbf {r}}\beta )}_{\phi _{i}=\phi _{0}({\mathbf {r}})}\prod _{{\mathbf {r}}}\mathrm{d}\phi ({\mathbf {r}}\tau )\ldots , \end{aligned}$$

(10)

where the notations $\phi _{i}$ and $\phi _{f}$ mean the initial and final paths. The path integral in Eq. (10) could be transformed into path integration over the compact U(1) group manifold, since the electromagnetic group U(1), governing the phase field, is compact, i.e., $\phi ({\mathbf {r}}\tau )$ has the topology of a circle ($S_1$); thus we have a non-homotopic mapping of the configuration space onto the U(1) gage group $S_{1}\rightarrow U(1)$. The paths, which loop around a circle in different number of times, are in different homotopy classes, and they cannot be continuously deformed into one another. All these paths can be characterized by their proper winding numbers [43] $m\left( {\mathbf {r}}\right)$. Any two paths which have different winding numbers cannot be continuously transformed one to another, and in order to include all possible phase path contributions, we have to sum over all topologically inequivalent phase configurations described by their winding numbers [43]. Accordingly, the path-integral in Eq. (10) is transformed as

$$\begin{aligned} \int \left[ \mathcal{{D}}\phi \right] \cdots =\int \left[ \mathcal{{D}}\varphi \right] \ldots , \end{aligned}$$

(11)

where the integration measure is now

$$\begin{aligned} \int \left[ \mathcal{{D}}\varphi \right] \cdots&\equiv \sum _{\left\{ m({\mathbf {r}})\right\} }\int ^{2\pi }_{0}\prod _{{\mathbf {r}}}\mathrm{d}\varphi _{0}({\mathbf {r}}) \nonumber \\&\times \int ^{\varphi ({\mathbf {r}}\beta )=\varphi _{0}\left( {\mathbf {r}}\right) +2{\pi }m({\mathbf {r}})}_{\varphi \left( {\mathbf {r}}0\right) =\varphi _{0}\left( {\mathbf {r}}\right) }\prod _{{\mathbf {r}}}\mathrm{d}\varphi ({\mathbf {r}}\tau )\ldots . \end{aligned}$$

(12)

In performing the integration over the phase field, one should take into account that the field configurations satisfy the boundary conditions [43]

$$\begin{aligned} \varphi ({\mathbf {r}}\beta )-\varphi ({\mathbf {r}}0)=2\pi {m({\mathbf {r}})}. \end{aligned}$$

(13)

Thus, integration over all phases $\phi ({\mathbf {r}}\tau )$ amounts to the integration over the $\beta$-periodic field $\varphi ({\mathbf {r}}\tau )$ and the summation over a set of U(1) integer winding numbers $m({\mathbf {r}})$. For the scalar static part $V_{0}({\mathbf {r}})$, we get the following functional integral:

$$\begin{aligned}&\int \left[ {D}V_{0}\right] e^{\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau -\frac{V^{2}_{0}({\mathbf {r}})}{U}+iV_{0}({\mathbf {r}})\left[ n({\mathbf {r}}\tau )-\frac{2\bar{\mu }}{U}\right] }. \end{aligned}$$

(14)

The saddle-point value of ${V}_{0}({{\mathbf {r}}})$ is given by

$$\begin{aligned} {{V}_{0}}=i\frac{Un}{2}-i\bar{\mu }, \end{aligned}$$

(15)

where $n$ is total average particle density $n=n_{c}+n_{f}$ (furthermore, we will fix $n$ as equal to $1$, corresponding to the case of half-filling).

Thereby, after decoupling of the quadratic term proportional to $n^{2}$ in the Hamiltonian in Eq. (1), we get a contribution to the partition function in Eq. (2) in the form

$$\begin{aligned} \exp \left[ {-\mathcal{{S}}\left[ \varphi \right] -\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau {\mu }_{n}}n({\mathbf {r}}\tau )\right] , \end{aligned}$$

(16)

where the emergent phase-only action $\mathcal{{S}}[\varphi ]$ is given as

$$\begin{aligned} \mathcal{{S}}[\varphi ]=\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ \frac{\dot{\varphi }^{2}({\mathbf {r}}\tau )}{U}-\frac{2\bar{\mu }}{iU}\dot{\varphi }({\mathbf {r}}\tau )-i\dot{\varphi }({\mathbf {r}}\tau )n({\mathbf {r}}\tau )\right] \end{aligned}$$

(17)

and the effective chemical potential ${\mu }_{n}$, attached to the total density operator in Eq. (16), is given in the form ${\mu }_{n}=\frac{Un}{2}-\bar{\mu }$.

2.1.3 Decoupling of Term Proportional to $\tilde{n}^{2}$

The decoupling of the quadratic term proportional to $\tilde{n}^{2}({\mathbf {r}}\tau )$ in the exponent of the partition function in Eq. (2) is also straightforward. We obtain

$$\begin{aligned}&\exp \left[ {\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \frac{U}{4}\tilde{n}^{2}({\mathbf {r}}\tau )}\right] \nonumber \\&\quad =\int \left[ {D}{\varrho }\right] e^{-\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ \frac{\varrho ^{2}({\mathbf {r}}\tau )}{U}-\varrho ({\mathbf {r}}\tau )\tilde{n}({\mathbf {r}}\tau )\right] }. \end{aligned}$$

(18)

After combining the expression in the exponent in Eq. (18) with the similar linear term in the expression of the Hamiltonian in Eq. (1) (see the fourth term in Eq. 1), we have

$$\begin{aligned}&\int \left[ {D}{\varrho }\right] e^{\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau -\frac{\varrho ^{2}({\mathbf {r}}\tau )}{4U}+\varrho ({\mathbf {r}}\tau )\left[ \tilde{n}({\mathbf {r}}\tau )-\frac{\epsilon _{c}-\epsilon _{f}}{2U}\right] }. \end{aligned}$$

(19)

The saddle-point evaluation for $\varrho$ gives

$$\begin{aligned} \varrho _{0}=\frac{U\tilde{n}}{2}-\frac{\epsilon _{c}-\epsilon _{f}}{2}, \end{aligned}$$

(20)

where $\tilde{n}=\left\langle \tilde{n}({\mathbf {r}}\tau )\right\rangle$ is the average of the particle density difference function. As a result of the decoupling, we obtain a “Zeeman”-like contribution to the partition function

$$\begin{aligned} \exp \left[ {-\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \mu _{\tilde{n}}\tilde{n}({\mathbf {r}}\tau )}\right] \end{aligned}$$

(21)

with the attached effective chemical potential $\mu _{\tilde{n}}=\frac{\epsilon _{c}-\epsilon _{f}}{2}-\frac{U\tilde{n}}{2}$.

2.1.4 Linearized Action with Phase-Field Contribution

To summarize, the grand canonical partition function of the system, after of both procedures of decoupling, is

$$\begin{aligned} {Z}_{GC}=\int \left[ {D}\bar{c}{D}c\right] \left[ {D}\bar{f}{D}f\right] \left[ {D}\varphi \right] e^{-\mathcal{{S}}[\bar{c},c,{\bar{f}},f,\varphi ]}, \end{aligned}$$

(22)

where the action $\mathcal{{S}}[\bar{c},c,{\bar{f}},f,\varphi ]$ in the exponent is given by

$$\begin{aligned} \mathcal{{S}}[\bar{c},c,{\bar{f}},f,\varphi ]&= \mathcal{{S}}\left[ \varphi \right] +\sum _{x=f,c}\mathcal{{S}}_\mathrm{B}\left[ \bar{x},x\right] \nonumber \\&-\,t_{c}\sum _{\left\langle {\mathbf {r}}, {\mathbf {r}}' \right\rangle }\int ^{\beta }_{0}\mathrm{d}\tau \left[ \bar{c}({{\mathbf {r}}}\tau )c({{\mathbf {r}}}'\tau )+h.c.\right] \nonumber \\&-\,{t}_{f}\sum _{\left\langle {\mathbf {r}}, {\mathbf {r}}' \right\rangle }\int ^{\beta }_{0}\mathrm{d}\tau \left[ \bar{f}({{\mathbf {r}}}\tau )f({{\mathbf {r}}}'\tau )+h.c.\right] \nonumber \\&+\,\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ {\mu }_{n}n({\mathbf {r}}\tau )+\mu _{\tilde{n}}\tilde{n}({\mathbf {r}}\tau )\right] . \end{aligned}$$

(23)

After the Hubbard–Stratanovich linearisation, we got the total action of the system that is linear in terms of fermion densities and contains in addition a phase-dependent term $\mathcal{{S}}\left[ \varphi \right]$, and also the terms, proportional to the effective chemical potentials $\mu _{n}$ and $\mu _{\tilde{n}}$.

2.2 The U$(1)$ Transformation

In the perspective to treat the local and non-local correlations in our excitonic system, it is important to separate the U(1) gage degrees of freedom related to the phase sector. To this end, we perform the local gage transformation to new fermion Grassmann variables $\tilde{f}({\mathbf {r}}\tau )$ and $\tilde{c}({\mathbf {r}}\tau )$. Meanwhile, this procedure will automatically eliminate also the last imaginary term, appearing in the expression of the phase action in Eq. (17).

For the electrons of $f$ and $c$-orbitals, the U$(1)$ transformation is

$$\begin{aligned} \left[ \begin{array}{cc} x({\mathbf {r}}\tau ) \\ \bar{x}({\mathbf {r}}\tau ) \end{array} \right] =\hat{{\mathcal {U}}}(\varphi )\cdot \left[ \begin{array}{cc} \tilde{x}({\mathbf {r}}\tau ) \\ \bar{\tilde{x}}({\mathbf {r}}\tau ) \end{array} \right] , \end{aligned}$$

(24)

where ${\hat{\mathcal {U}}}(\varphi )$ is the U(1) transformation matrix ${\hat{\mathcal {U}}}(\varphi )=\hat{I}\cdot \cos \varphi ({\mathbf {r}}\tau )+i\hat{\sigma }_{z}\cdot \sin \varphi ({\mathbf {r}}\tau )$ with the unit matrix $\hat{I}$ and $\hat{\sigma }_{z}$ being the Pauli matrix. The variables $\tilde{x}=\tilde{f}$, $\tilde{c}$, and we used the bosonic phase variables $\varphi$ introduced in Eqs. (11) and (12). After transformations given in Eq. (24), we obtain the total action of the system in the U(1) gage-invariant form (for comparison, see the action in Eq. 23 before transformations)

$$\begin{aligned} \mathcal{{S}}[\bar{\tilde{c}},\tilde{c},{\bar{\tilde{f}}}, \tilde{f},\varphi ]&= \mathcal{{S}}_{0}[\varphi ]+\sum _{x=\tilde{f},\tilde{c}}\mathcal{{S}}_\mathrm{B}\left[ \bar{\tilde{x}},\tilde{x}\right] \nonumber \\&-\,t\sum _{\left\langle {\mathbf {r}},{\mathbf {r}}' \right\rangle }\int ^{\beta }_{0}\mathrm{d}\tau \left[ \bar{\tilde{c}}({{\mathbf {r}}}\tau )\tilde{c}({{\mathbf {r}}}'\tau )e^{-i\left[ \varphi ({{\mathbf {r}}}\tau )-\varphi ({{\mathbf {r}}}'\tau )\right] }+h.c.\right] \nonumber \\&-\,\tilde{t}\sum _{\left\langle {\mathbf {r}},{\mathbf {r}}' \right\rangle }\int ^{\beta }_{0}\mathrm{d}\tau \left[ \bar{\tilde{f}}({{\mathbf {r}}}\tau )\tilde{f}({{\mathbf {r}}}'\tau )e^{-i\left[ \varphi ({{\mathbf {r}}}\tau )-\varphi ({{\mathbf {r}}}'\tau )\right] }+h.c.\right] \nonumber \\&+\,\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ {\mu }_{n}n({\mathbf {r}}\tau )+\mu _{\tilde{n}}\tilde{n}({\mathbf {r}}\tau )\right] . \end{aligned}$$

(25)

Now, $\tilde{t}$ and $t$ in Eq. (25) are, respectively, $\tilde{f}$-band and $\tilde{c}$-band fermion transfer integrals. We got in Eq. (25) also, a new, emergent, quadratic phase action $\mathcal{{S}}_{0}[\varphi ]$

$$\begin{aligned} \mathcal{{S}}_{0}[\varphi ]=\sum _{{\mathbf {r}}}\int ^{\beta }_{0}\mathrm{d}\tau \left[ \frac{\dot{\varphi }^{2}({\mathbf {r}}\tau )}{U}-\frac{2\bar{\mu }}{iU}\dot{\varphi }({\mathbf {r}}\tau )\right] . \end{aligned}$$

(26)

The partition function of the system in the new variables $\tilde{f}$ and $\tilde{c}$ is

$$\begin{aligned} {Z}_{GC}=\int \left[ D\bar{\tilde{c}}D\tilde{c}\right] \left[ D\bar{\tilde{f}}D\tilde{f}\right] \left[ D\varphi \right] e^{-\mathcal{{S}}[\bar{\tilde{c}},\tilde{c},{\bar{\tilde{f}}}, \tilde{f},\varphi ]}. \end{aligned}$$

(27)

From this form of the partition function, we will generate the effective actions for fermions and for bosonic phase sector (see the general procedure presented in Fig. 1).

Fig. 1

Functional integration procedure. The final actions are obtained after phase (for the EI state) and fermion averaging (for the phase stiffness and excitonic condensate) (Color figure online)

3 Effective Action for Fermions

By following the left-lowest root, presented in Fig. 1, we will integrate out the phase variables. We obtain

$$\begin{aligned} \mathcal{{Z}}=\int \left[ \mathcal{{D}}\bar{\tilde{c}}\mathcal{{D}}\tilde{c}\right] \left[ \mathcal{{D}}\bar{\tilde{f}}\mathcal{{D}}\tilde{f}\right] e^{-\mathcal{{S}}_\mathrm{eff}[\bar{\tilde{c}},\tilde{c},{\bar{\tilde{f}}},f]} \ , \end{aligned}$$

(28)

where the effective, phase-averaged fermionic action in the exponent is

$$\begin{aligned} \mathcal{{S}}_\mathrm{eff}[\bar{\tilde{c}},\tilde{c},{\bar{\tilde{f}}},\tilde{f}]=-\!\ln {\int \left[ \mathcal{{D}}\varphi \right] }e^{-\mathcal{{S}}[\bar{\tilde{c}},\tilde{c},{\bar{\tilde{f}}},\tilde{f},\varphi ]}. \end{aligned}$$

(29)

The Fourier transformation of fermionic variables $\tilde{f}$ and $\tilde{c}$ is

$$\begin{aligned} x({\mathbf {r}}\tau )=\frac{1}{{\beta {N}}}\sum _{{\mathbf {k}},\nu _{n}}x_{{\mathbf {k}}}(\nu _{n})e^{i({\mathbf {k}}{\mathbf {r}}-\nu _{n}\tau )} \end{aligned}$$

(30)

where $N$ is the number of lattice sites, and $\nu _{n}={\pi (2n+1)}/{\beta }$ are the Fermi–Matsubara frequencies [41] with $n=0,\pm 1,\pm 2,\ldots$.

Then, the effective phase-averaged fermionic action of the system in the Fourier space takes the following form:

$$\begin{aligned} \mathcal{{S}}_\mathrm{eff}\left[ \tilde{\bar{c}},\tilde{c},\tilde{\bar{f}},\tilde{f}\right]&= \frac{1}{\beta {N}}\sum _{{\mathbf {k}}\nu _{n}}\left[ \bar{\tilde{c}}_{\mathbf {k}}(\nu _{n}),\bar{\tilde{f}}_{\mathbf {k}}({\nu _{n}})\right] \nonumber \\&\times \,\mathcal{{G}}^{-1}({\mathbf {k}},\nu _{n})\left[ \begin{array}{cc} {\tilde{c}}_{\mathbf {k}}(\nu _{n})\\ {\tilde{f}}_{\mathbf {k}}(\nu _{n}) \end{array} \right] . \end{aligned}$$

(31)

Here, $\mathcal{{G}}^{-1}({\mathbf {k}},\nu _{n})$ is the inverse of the Green function matrix

$$\begin{aligned} \mathcal{{G}}^{-1}({\mathbf {k}},\nu _{n})= \left( \begin{array}{cc} {E}^{\tilde{c}}_{{\mathbf {k}}}(\nu _{n}) &{} \quad -\bar{\Delta } \\ -\Delta &{} \quad {E}^{\tilde{f}}_{{\mathbf {k}}}(\nu _{n}) \end{array} \right) , \end{aligned}$$

(32)

where we have introduced the local excitonic order parameter $\Delta =\left\langle \bar{\tilde{c}}({\mathbf {r}}\tau )\tilde{f}({\mathbf {r}}\tau )\right\rangle$. The single-particle Bogoliubov quasienergies ${E}^{\tilde{f}}_{{\mathbf {k}}}(\nu _{n})$ and ${E}^{\tilde{c}}_{{\mathbf {k}}}(\nu _{n})$ are ${E}^{\tilde{f}}_{{\mathbf {k}}}(\nu _{n})=\bar{\epsilon }_{\tilde{f}}-i\nu _{n}-\tilde{t}_{{\mathbf {k}}}$ and ${E}^{\tilde{c}}_{{\mathbf {k}}}(\nu _{n})=\bar{\epsilon }_{\tilde{c}}-i\nu _{n}-{t}_{{\mathbf {k}}}$. Next, $\tilde{t}_{{\mathbf {k}}}$ and ${t}_{{\mathbf {k}}}$ are band-renormalized hopping amplitudes $\tilde{t}_{{\mathbf {k}}}=2\tilde{t}{\mathrm {g}}_\mathrm{B}\gamma _{{\mathbf {k}}}$ and ${t}_{{\mathbf {k}}}=2t{\mathrm {g}}_\mathrm{B}\gamma _{{\mathbf {k}}}$, where ${\mathrm {g}}_\mathrm{B}$ is the 3D bandwidth renormalization factor

$$\begin{aligned} \mathrm {g}_\mathrm{B}=\left. \left\langle e^{-i[\varphi ({{\mathbf {r}}}\tau )-\varphi ({{\mathbf {r}}}'\tau )]} \right\rangle \right| _{|{\mathbf {r}}-{\mathbf {r}}'|={{d}}} \ , \end{aligned}$$

(33)

and $\gamma _{{\mathbf {k}}}$ is the 3D lattice dispersion $\gamma _{{\mathbf {k}}}=\cos (k_{x}d_{x})+\cos (k_{y}d_{y})+\cos (k_{z}d_{z})$, with $d_{\alpha }$ ($\alpha =x,y,z$), being the components of the lattice spacing vector ${\mathbf {d}}={\mathbf {r}}-{\mathbf {r}}'$ with ${\mathbf {r}}$ and ${\mathbf {r}}'$, being nearest neighbors site positions. For the simple cubic lattice, we have $d_{\alpha }\equiv a=1$.

The quasiparticle energies $\bar{\epsilon }_{\tilde{f}}$ and $\bar{\epsilon }_{\tilde{c}}$ are of Hartree type, and they are defined in the theory by relation $\bar{\epsilon }_{\tilde{x}}=\epsilon _{{x}}-\mu +Un_{\tilde{y}}+i\left\langle \dot{\varphi }({{\mathbf {r}}}\tau )\right\rangle$, where $\tilde{y}$ means orbital, opposite to $\tilde{x}$. The EI low-temperature phase is characterized by the local excitonic order parameter $\Delta$. Without any loss of generality, we can suppose the case of the EI state, with uniform real gap parameter $\Delta =\bar{\Delta }$. The EI state develops from local on-site electron–hole correlations. The expectation value, given in the expression of local EI order parameter, could be calculated in the frame of path-integral method, as well as the fermion density averages of the respective band levels $n_{\tilde{x}}=\left\langle \bar{\tilde{x}}({\mathbf {r}}\tau )\tilde{x}({\mathbf {r}}\tau )\right\rangle$. We get a set of the coupled self-consistent equations for the EI order parameter $\Delta$, single-particle fermion densities $n_{\tilde{x}}$, and EI chemical potential $\mu$

$$\begin{aligned}&\frac{1}{N}\sum _{{\mathbf {k}}}\left[ f(E^{+}_{{\mathbf {k}}})+f(E^{-}_{{\mathbf {k}}})\right] =1, \end{aligned}$$

(34)

$$\begin{aligned}&\tilde{n}=\frac{1}{N}\sum _{{\mathbf {k}}}\xi _{{\mathbf {k}}}\cdot \frac{f(E^{+}_{{\mathbf {k}}})-f(E^{-}_{{\mathbf {k}}})}{\sqrt{\xi ^{2}_{{\mathbf {k}}}+4\Delta ^{2}}}, \end{aligned}$$

(35)

$$\begin{aligned}&\Delta =-\frac{U\Delta }{N}\sum _{{\mathbf {k}}}\frac{f(E^{+}_{{\mathbf {k}}})-f(E^{-}_{{\mathbf {k}}})}{\sqrt{\xi ^{2}_{{\mathbf {k}}}+4\Delta ^{2}}}. \end{aligned}$$

(36)

Here, $f(\epsilon )=1/\left( e^{\beta \epsilon }+1\right)$ is the Fermi–Dirac distribution function, $\xi _{{\mathbf {k}}}=-{t}_{{\mathbf {k}}}+\bar{\epsilon }_{\tilde{c}}+\tilde{t}_{{\mathbf {k}}}-\bar{\epsilon }_{\tilde{f}}$ is the quasiparticle dispersion, and the energy parameters $E^{+}_{{\mathbf {k}}}$ and $E^{-}_{{\mathbf {k}}}$ are defined as

$$\begin{aligned} E^{\pm }_{{\mathbf {k}}}=\frac{1}{2}\left( -{t}_{{\mathbf {k}}}+\bar{\epsilon }_{\tilde{c}}-\tilde{t}_{{\mathbf {k}}}+\bar{\epsilon }_{\tilde{f}}\pm {\sqrt{\xi ^{2}_{{\mathbf {k}}}+4\Delta ^{2}}}\right) . \end{aligned}$$

(37)

In fact, the difference between Eqs. (34)–(36) and the Hartree–Fock results [20] is in the presence of bandwidth renormalization factor $g_\mathrm{B}$. The calculation of the factor $g_\mathrm{B}\left( {\mathbf {r}}-{\mathbf {r}}'\right)$ could be done effectively within the self-consistent harmonic approximation (SCHA) method [44–46]. In this approximation, the quantum rotor description is reduced to classical Hamiltonian one by the Feynmann-Kleinert minimization procedure [44, 45]. Our results for SCHA show that the factor $g_\mathrm{B}$ is equal identically to $1$ at $T=0$ as in the two-dimensional (2D) case (see in Ref. [44] for details). For higher temperatures, it differs from unity, but not much.

The difference between the energy parameters in Eq. (37) defines the charge transfer gap in the system $\Delta _{c}=E^{+}_{{\mathbf {k}}}-E^{-}_{{\mathbf {k}}}=\sqrt{\xi ^{2}_{{\mathbf {k}}}+4\Delta ^{2}}$ (see in Ref. [47] for details).

The numerical solution of the system of equation Eqs. (34)–(36) is discussed in detail in Ref. [47], where finite-difference approximation method is used in numerical evaluations. Different values of $\tilde{t}$ hopping amplitude are considered there. The results are coinciding well with the previous Hartree–Fock (HF), improved slave boson, and 2D constrained path Monte Carlo investigations [15–20, 32, 33]. This good correspondence is ascribed to a rather weak band renormalization at $T=0$ (in our case $g_\mathrm{B}=1$). At finite temperatures, the particle number fluctuations are important and the band renormalization becomes necessary, especially when approaching from the band insulator (BI) high-temperature side. Indeed, as the numerical evaluations show, the transition temperature $T_\mathrm{EI}$ of the e–h pair formation is not vanishing for the case of the vanishing narrow band hopping $\tilde{t}=0$.

The exact solutions for the chemical potential could be obtained from Eqs. (34)–(36), both at the boundary of EI transition (i.e., when $\Delta \left( T_\mathrm{EI}, U\right) =0$) and in the EI state region ($\Delta \left( T < T_\mathrm{EI}, U\right) \ne 0$). The results are discussed in Ref. [47].

Meanwhile, it is also shown in Ref. [47] that the excitonic BEC transition critical temperature $T_{c}$ is much smaller than the critical temperature $T_\mathrm{EI}$ of the excitonic pair formation, in good agreement with previous theoretical predictions [15, 23–26, 44, 48]. The self-consistent numerical solutions for $T_\mathrm{EI}$ and $T_{c}$ are shown in Fig. 2 for $\tilde{t}=-0.3$. It is clear that, for the case of intermediate and strong interaction limits, the coherent BEC transition critical temperature $T_{c}$ is smaller of about two orders of magnitude than temperature $T_\mathrm{EI}$. Contrary, in the very small interaction case, we have the coincidence of both transition temperatures. This result is expected also from general considerations, because at small $U/t$, we are in the BCS limit, which means that the pairing and condensation occur simultaneously.

Fig. 2

EI transition critical temperature $T_\mathrm{EI}/t$, (upper curve) and excitonic BEC transition critical temperature $T_{c}/t$ (lower curve) for the $f$-band hopping $\tilde{t}=-0.3$. Different phase separation regions are shown in the figure. The data for the lower curve were taken from Ref. [47] (Color figure online)

A similar effect is found recently in Ref. [44], considering 2D excitonic systems, where the exciton-superfluid transition critical temperature is found much smaller than the critical temperature $T_\mathrm{EI}$ of exciton pair formation. Such a reduction of the transition temperature, due to the coherent pairing scheme, is given also in Ref. [48], where the BCS–Bose crossover is studied in the 2D attractive Hubbard model. Especially, a HF pair formation temperature is estimated, which corresponds to the regime of the incoherent or local pairs (for a comparison, see the correspondence to our $T_\mathrm{EI}$, and order parameter $\Delta$, which is also local), in difference with superconducting pairing temperature, at which the coherent Cooper pairs start to be formed. The general idea used in Ref. [47] is based on the non-local n.n. excitonic exchange correlation mechanism and Bogoliubov’s mean-field self-consistency assumption [49] for the effective phase action obtained after fermion integration procedure in the action in Eq. (27) (this is represented by the right-lowest root in Fig. 1, given in Sect. 2.2). The obtained Hamiltonian for the phase sector is very similar to the classical Hamiltonian one, [44–47] with an effective phase-stiffness parameter $J$ that emerges after the fermion Wick averaging procedure [41]. The corresponding phase action is [47]

$$\begin{aligned} \mathcal{{S}}_{J}\left[ \varphi \right] =-2J{\mathrm {g}_\mathrm{B}}\cdot \int ^{\beta }_{0}\mathrm{d}\tau \sum _{\left\langle {\mathbf {r}},{\mathbf {r}}'\right\rangle }\cos {\left[ \varphi ({\mathbf {r}}\tau )-\varphi ({\mathbf {r}}'\tau )\right] }. \end{aligned}$$

(38)

It appears that the non-zero value of quantity $J$ is directly related with the pairing gap $\Delta$, since it is shown [47] that $J$ vanishes, when $\Delta =0$. Here, we present only the final analytical result for $J$ (for details see in Ref. [47])

$$\begin{aligned} J&= \frac{16\Delta ^{2}t\tilde{t}}{z^{2}{N^{2}}}\sum _{{\mathbf {k}},{\mathbf {k}}'}\frac{\epsilon \left( {{\mathbf {k}}}\right) \epsilon \left( {{\mathbf {k}}}'\right) }{{\sqrt{\xi ^{2}_{{\mathbf {k}}}+4\Delta ^{2}}}}\left[ \Lambda _{1}({\mathbf {k}},{\mathbf {k}}')\tanh \left( \frac{\beta E^{+}_{{\mathbf {k}}}}{2}\right) \right. \nonumber \\&\left. -\,\Lambda _{2}({\mathbf {k}},{\mathbf {k}}')\tanh \left( \frac{\beta E^{-}_{{\mathbf {k}}}}{2}\right) \right] , \end{aligned}$$

(39)

where $z$ is the number of n.n. sites on the 3D lattice. The parameters $\Lambda _{1}({\mathbf {k}},{\mathbf {k}}')$ and $\Lambda _{2}({\mathbf {k}},{\mathbf {k}}')$ entering in Eq. (38) are given by $\Lambda _{1}({\mathbf {k}},{\mathbf {k}}')=\left( E^{+}({\mathbf {k}}) - E^{+}({\mathbf {k}}')\right) ^{-1}\cdot \left( {E^{+}({\mathbf {k}}) - E^{-}({\mathbf {k}}')}\right) ^{-1}$, $\Lambda _{2}({\mathbf {k}},{\mathbf {k}}')=\left( {E^{-}({\mathbf {k}}) - E^{-}({\mathbf {k}}')}\right) ^{-1}\cdot \left( {E^{-}({\mathbf {k}}) - E^{+}({\mathbf {k}}')}\right) ^{-1}$. Also, $J$ is strictly positive for all the regions of the normalized Coulomb interaction parameter $U/t$. In addition, it follows from the analytical form of $J$ that the macroscopic phase coherence in the system is characterized by an energy scale $J\sim t_{e}\cdot t_{h}/(t_{e}+t_{h})$ for all values of the Coulomb interaction parameter $U$. This is related to the motion of the center of mass of e–h composed of quasiparticle [15], because $t_{e}\cdot t_{h}/(t_{e}+t_{h}) \sim (m_{e}+m_{h})^{-1}$. For the strong interaction case, we are converging with the hard core Boson model, with the kinetic energy scale $\Delta t_{e}\cdot t_{h}/U$ (with $\Delta$ being the local excitonic order parameter). Thereby, it is shown in Ref. [47] that non-local correlations between the electrons and holes of different n.n. excitonic pairs are related with the excitonic BEC condensation. Furthermore, in the frame of the quantum rotor model, the excitonic BEC transition probability function is derived and its temperature dependence is found [47].

As a systematic continuation of the theoretical study given in Ref. [47], we elaborate here on the analytical forms of normal, $f$ and $c$-band, Green functions (incoherent and coherent) and coherent, anomalous excitonic Green function. Also, we will discuss in detail the spectral functions and DOS corresponding, in the next. The numerical evaluations of calculated DOS functions are shown in Sect. 4.4.

4 Single-Particle DOS

4.1 Normal and Excitonic DOS Functions

The spectral density functions of the system of interacting exciton gas will determine the excitonic center-of-mass distribution related to the condensation in the low-temperature limit. Therefore, the calculation of these functions represents an important task. Within our theoretical approach, we can access a variety of correlation functions in the system. We will concentrate now on the $c$- and $f$-band normal spectral functions and also the excitonic anomalous spectral function, represented in terms of the initial operators. The $c$- and $f$-band normal excitonic Green functions are

$$\begin{aligned} G_\mathrm{x,x}\big ({\mathbf {r}}\tau ,{\mathbf {r}}'\tau '\big )=-\big \langle {x}({\mathbf {r}}\tau )\bar{x}({\mathbf {r}}'\tau ')\big \rangle , \end{aligned}$$

(40)

and the anomalous excitonic Green function is defined as

$$\begin{aligned} G_\mathrm{c,f}\big ({\mathbf {r}}\tau ,{\mathbf {r}}'\tau '\big )=\big \langle \bar{c}({\mathbf {r}}\tau )f({\mathbf {r}}'\tau ')\big \rangle . \end{aligned}$$

(41)

After introducing the U(1) transformations, defined in Eq. (24), we will have the Green function’s decomposition as

$$\begin{aligned} G_\mathrm{x,x}({\mathbf {r}}\tau ,{\mathbf {r}}'\tau ')=-\Big \langle \tilde{x}({\mathbf {r}}\tau )\bar{\tilde{x}}\Big ({\mathbf {r}}'\tau '\Big )\Big \rangle \cdot \Big \langle e^{-i\left[ \varphi ({\mathbf {r}}\tau )-\varphi ({\mathbf {r}}'\tau ')\right] }\Big \rangle , \end{aligned}$$

(42)

$$\begin{aligned} G_\mathrm{c,f}({\mathbf {r}}\tau ,{\mathbf {r}}'\tau ')=\Big \langle \bar{\tilde{c}}({\mathbf {r}}\tau )\tilde{f}\Big ({\mathbf {r}}'\tau '\Big )\Big \rangle \cdot \Big \langle e^{-i\left[ \varphi ({\mathbf {r}}\tau )-\varphi ({\mathbf {r}}'\tau ')\right] }\Big \rangle . \end{aligned}$$

(43)

The phase factors, appearing in the definitions of Green functions, define in fact a charge-bosonic propagator $G_{z}({\mathbf {r}}\tau , {\mathbf {r}}'\tau ')$ as follows:

$$\begin{aligned} G_{z}\Big ({\mathbf {r}}\tau , {\mathbf {r}}'\tau '\Big )=\Big \langle e^{-i\left[ \varphi ({\mathbf {r}}\tau )-\varphi ({\mathbf {r}}'\tau ')\right] }\Big \rangle . \end{aligned}$$

(44)

Then, we pass to the Fourier space representation for the U(1) - transformed Green functions

$$\begin{aligned} \tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\Big ({\mathbf {r}}\tau , {\mathbf {r}}'\tau '\Big )&= \frac{1}{\beta {N}}\sum _{{\mathbf {k}},\nu _{n}}\tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n}) \nonumber \\&\times \, e^{i\left[ {\mathbf {k}}\big ({\mathbf {r}}-{\mathbf {r}}'\big )-\nu _{n}\big (\tau -\tau '\big )\right] }, \end{aligned}$$

(45)

$$\begin{aligned} \tilde{G}_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}\Big ({\mathbf {r}}\tau , {\mathbf {r}}'\tau '\Big )&= \frac{1}{\beta {N}}\sum _{{\mathbf {k}},\nu _{n}}\tilde{G}_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu _{n}) \nonumber \\&\times \, e^{i\left[ {\mathbf {k}}\big ({\mathbf {r}}-{\mathbf {r}}'\big )-\nu _{n}\big (\tau -\tau '\big )\right] } \end{aligned}$$

(46)

and

$$\begin{aligned} G_{z}({\mathbf {r}}\tau , {\mathbf {r}}'\tau ')=\frac{1}{\beta {N}}\sum _{{\mathbf {k}},\omega _{n}}G_{z}({\mathbf {k}},\omega _{n})e^{i\left[ {\mathbf {k}}\big ({\mathbf {r}}-{\mathbf {r}}'\big )-\omega _{n}\big (\tau -\tau '\big )\right] }. \end{aligned}$$

(47)

Furthermore, the Fourier transformation of the functions in Eqs. (42)–(43) will be rewritten as a convolution in the reciprocal ${\mathbf {k}}$-space

$$\begin{aligned} G_\mathrm{x,x}({\mathbf {k}},\nu _{n})=\frac{1}{\beta {N}}\sum _{{\mathbf {q}},\omega _{n}}G_{z}\big ({\mathbf {q}},\omega _{n}\big )\cdot \tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\big ({\mathbf {k}}-{\mathbf {q}},\nu _{n}-\omega _{n}\big ) \end{aligned}$$

(48)

and

$$\begin{aligned} G_\mathrm{c,f}({\mathbf {k}},\nu _{n})=\frac{1}{\beta {N}}\sum _{{\mathbf {q}},\omega _{n}}G_{z}\big ({\mathbf {q}},\omega _{n}\big )\cdot \tilde{G}_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}\big ({\mathbf {k}}-{\mathbf {q}},\nu _{n}-\omega _{n}\big ). \end{aligned}$$

(49)

It is worth to mention that frequency summations in Eqs. (48) and (49) are over Bose–Matsubara frequencies $\omega _{n}={2\pi {n}/\beta }$. The Fermionic Green functions will be calculated using the formalism discussed in Sects. 2 and 3 and also, functional derivation techniques [39]. Particularly, for the $f$- and $c$-band Green functions $\tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n})$, we get

$$\begin{aligned} \tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( {\mathbf {k}},i\nu _{n}\right)&= \big \langle \tilde{ x}\big ({\mathbf {k}},\nu _{n}\big )\bar{\tilde{x}}\big ({\mathbf {k}},\nu _{n}\big )\big \rangle \nonumber \\&= \frac{E^{\tilde{y}}_{{\mathbf {k}}}\left( \nu _{n}\right) }{E^{\tilde{x}}_{{\mathbf {k}}}\left( \nu _{n}\right) E^{\tilde{y}}_{{\mathbf {k}}}\left( \nu _{n}\right) -\Delta ^{2}}. \end{aligned}$$

(50)

In general, the experimental observation of hybridization between the valence band and conduction band could be done by examining the ARPES spectra, which measure the spectral intensities just above and below the temperature $T_\mathrm{EI}$ of excitonic pair formation. In ARPES experiments, one observes the imaginary part of the real retarded Green function; therefore, the calculation of it represents a remarkable importance. Indeed, the single-particle DOS is related with the imaginary part of the retarded Green functions, and thus we need to calculate real retarded function, which corresponds to the normal Matsubara Green function $\tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( {\mathbf {k}},i\nu _{n}\right)$. This could be done by the analytical continuation into the upper-half complex semi-plane ($\nu _{n}>0$) of frequency modes $i \nu _{n}$

$$\begin{aligned} \tilde{G}^{R}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu )=\tilde{G}_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\big ({\mathbf {k}},i\nu _{n}\big )|_{i\nu _{n}\rightarrow \nu +i\eta }. \end{aligned}$$

(51)

The single-particle DOS is defined then as

$$\begin{aligned} \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( {\mathbf {k}},\nu \right)&= -\frac{1}{\pi }\mathrm{Im }\tilde{G}^{R}_\mathrm{x,x}({\mathbf {k}},\nu )=\left( \bar{\epsilon }_{\tilde{y}}-\tilde{t}_{{\mathbf {k}}}-\nu \right) ^{2} \nonumber \\&\times \,\delta \left[ \left( \nu ^{2}+A_{\mathbf {k}}\nu +B_{{\mathbf {k}}}\right) \cdot \left( \bar{\epsilon }_{\tilde{y}}-\tilde{t}_{{\mathbf {k}}}-\nu \right) \right] , \end{aligned}$$

(52)

where

$$\begin{aligned}&A_{\mathbf {k}}=t_{{\mathbf {k}}}+\tilde{t}_{\mathbf {k}}-\bar{\epsilon }_{\tilde{x}}-\bar{\epsilon }_{\tilde{y}} \end{aligned}$$

(53)

and

$$\begin{aligned} B_{\mathbf {k}}&= \bar{\epsilon }_{\tilde{x}}\bar{\epsilon }_{\tilde{y}}+4t\tilde{t}\left[ \cos (k_x)+\cos (k_y)+\cos (k_z)\right] ^{2} \nonumber \\&-\,2\tilde{t}\bar{\epsilon }_{\tilde{x}}\left[ \cos (k_x)+\cos (k_y)+\cos (k_z)\right] \nonumber \\&-\,2t\bar{\epsilon }_{\tilde{y}}\left[ \cos (k_x)+\cos (k_y)+\cos (k_z)\right] -\Delta ^{2}. \end{aligned}$$

(54)

${\mathbf {k}}$-summed DOS will be

$$\begin{aligned} \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( \nu \right) =\frac{1}{N}\sum _{\mathbf {k}}\rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( {\mathbf {k}},\nu \right) . \end{aligned}$$

(55)

The summations over the wave vectors in Eq. (55) can be simplified by introducing the appropriate DOS function for the 3D cubic lattice $\rho _{3D}(x)=\frac{1}{N}\sum _{{\mathbf {k}}}\delta (x-\gamma _{\mathbf {k}})$. Then, it is not difficult to show that

$$\begin{aligned} \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( \nu \right)&= \int ^{+3.0}_{-3.0}\mathrm{d}x \rho _{ 3D}(x)\frac{\left[ \bar{\epsilon }_{\tilde{y}}-\tilde{t}(x)-\nu \right] ^{2}}{\sqrt{\xi ^{2}(x)+4\Delta ^{2}}}\nonumber \\&\cdot \left\{ \frac{\delta \left[ \nu -E^{+}(x)\right] }{|\bar{\epsilon }_{\tilde{y}}-\tilde{t}(x)-E^{+}(x)|}+\frac{\delta \left[ \nu -E^{-}(x)\right] }{|\bar{\epsilon }_{\tilde{y}}-\tilde{t}(x)-E^{-}(x)|}\right\} . \end{aligned}$$

(56)

Here, $\tilde{t}(x)=2\tilde{t}x$ and $t(x)=2tx$, and the energy parameters $E^{\pm }(x)$ are continuous versions of parameters defined in Eq. (37). The DOS $\rho _{3D}(x)$, for the simple cubic 3D lattice, is given by

$$\begin{aligned} \rho _{3D}(x)&= \frac{1}{\pi ^{3}}\int ^{\min (1,2-x)}_{\max (-1,-2-x)}\mathrm{d}y \frac{\Theta \left( 1-\frac{|x|}{3}\right) }{\sqrt{1-y^{2}}} \nonumber \\&\times \,{{\mathbf {k}}\left[ \sqrt{1-\left( \frac{y}{2}+\frac{x}{2}\right) ^{2}}\right] }, \end{aligned}$$

(57)

where $\Theta (x)$ is the Heaviside step function, and ${\mathbf {k}}(x)$ is the elliptic function of the first kind: ${\mathbf {k}}(x)=\int ^{\pi /2}_{0}dt 1/\sqrt{1-x^{2}\sin ^{2}t}$.

For the anomalous excitonic Green function, we have

$$\begin{aligned} \tilde{G}_{\tilde{c},\tilde{f}}\left( {\mathbf {k}},i\nu _{n}\right) =\Big \langle \bar{\tilde{c}}({\mathbf {k}},\nu _{n})\tilde{f}({\mathbf {k}},\nu _{n})\Big \rangle =\frac{\Delta }{E^{\tilde{x}}_{{\mathbf {k}}}\left( \nu _{n}\right) E^{\tilde{y}}_{{\mathbf {k}}}\left( \nu _{n}\right) -\Delta ^{2}}. \end{aligned}$$

(58)

The retarded function, which corresponds to it, is then [41]

$$\begin{aligned} \tilde{G}^{R}_{\tilde{c},\tilde{f}}({\mathbf {k}},\nu )=\tilde{G}_{\tilde{c},\tilde{f}}\left( {\mathbf {k}},i\nu _{n}\right) |_{i\nu _{n}\rightarrow \nu +i\eta }. \end{aligned}$$

(59)

The single-particle DOS is forthcoming then as

$$\begin{aligned} \rho _{\tilde{c},\tilde{f}}\left( {\mathbf {k}},\nu \right) =-\frac{1}{\pi }\mathrm{Im } \tilde{G}^{R}_{\tilde{c},\tilde{f}}({\mathbf {k}},\nu ) =\Delta \cdot \delta \left( \nu ^{2}+A_{\mathbf {k}}\nu +B_{{\mathbf {k}}}\right) . \end{aligned}$$

(60)

Obviously, it has more simple form than the function in Eq. (52). The ${\mathbf {k}}$-summed DOS for excitons will be

$$\begin{aligned} \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )=\Delta \cdot \left\{ \frac{\rho _{3D}\left[ \Lambda _{1}(\nu )\right] }{|\chi _{1}\left[ \Lambda _{1}(\nu )\right] |}+\frac{\rho _{3D}\left[ \Lambda _{2}(\nu )\right] }{|\chi _{2}\left[ \Lambda _{2}(\nu )\right] |}\right\} , \end{aligned}$$

(61)

where the dimensionless parameters $\Lambda _{1,2}(\nu )$ are given by following expressions:

$$\begin{aligned} \Lambda _{1}(\nu )=\frac{-\left[ \left( t+\tilde{t}\right) \nu -\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}+\bar{\epsilon }_{\tilde{f}}t\right) \right] +\sqrt{\left[ \left( t-\tilde{t}\right) \nu +\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}-\bar{\epsilon }_{\tilde{f}}t\right) \right] ^{2}+4t\tilde{t}|\Delta |^{2}}}{4t\tilde{t}}, \end{aligned}$$

(62)

$$\begin{aligned} \Lambda _{2}(\nu )=\frac{-\left[ \left( t+\tilde{t}\right) \nu -\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}+\bar{\epsilon }_{\tilde{f}}t\right) \right] -\sqrt{\left[ \left( t-\tilde{t}\right) \nu +\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}-\bar{\epsilon }_{\tilde{f}}t\right) \right] ^{2}+4t\tilde{t}|\Delta |^{2}}}{4t\tilde{t}} \nonumber \\ \end{aligned}$$

(63)

and the functions $\chi _{i}\left[ \Lambda _{1}(\nu )\right]$ ($i=1,2$), in the denominators in the right-hand side in Eq. (61) are

$$\begin{aligned} \chi _{1}\left[ \Lambda _{1}(\nu )\right] =2\left( t+\tilde{t}\right) \nu +8t\tilde{t}\Lambda _{1}(\nu )-2\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}+\bar{\epsilon }_{\tilde{f}}t\right) \end{aligned}$$

(64)

and

$$\begin{aligned} \chi _{2}\left[ \Lambda _{2}(\nu )\right] =2\left( t+\tilde{t}\right) \nu +8t\tilde{t}\Lambda _{2}(\nu )-2\left( \bar{\epsilon }_{\tilde{c}}\tilde{t}+\bar{\epsilon }_{\tilde{f}}t\right) . \end{aligned}$$

(65)

Now, turning to the convolution forms for total fermionic and excitonic Green functions in Eqs. (48) and (49), we need an explicit expression for the phase-bosonic Green function $G_{z}\left( {\mathbf {k}},\omega _{n}\right)$. We will calculate it in the formalism of the effective phase action given in the quantum rotor model, discussed earlier in Ref. [47], where we have derived the effective phase-only action $S_\mathrm{eff}[\varphi ]$ by integrating the fermions in the partition function in Eq. (27). In the following, we cast $S_\mathrm{eff}[\varphi ]$ into the quantum rotor representation [47].

4.2 The Phase-Stiffness DOS

To proceed, we replace the phase degrees of freedom with the complex, unimodular field $z({\mathbf {r}}\tau )=e^{i\varphi ({\mathbf {r}}\tau )}$, which satisfies the time-periodic boundary condition $z({\mathbf {r}}\beta )=z({\mathbf {r}}0)$. The spherical constraint, imposed on a set of the unimodular variables is ${1/N}\sum _{\mathbf {k}}|z({\mathbf {r}}\tau )|^{2}=1$. Now, we introduce these new variables $z({\mathbf {r}}\tau )$ into the partition function in Eq. (27), in a way, consistent with the Faddeev–Popov ghost-field method [50]

$$\begin{aligned}&\int \mathcal{{D}}\bar{z}\mathcal{{D}}{z} \delta \left( \frac{}{}\sum _{{\mathbf {r}}}|z({\mathbf {r}}\tau )|^{2}-N\frac{}{}\right) \nonumber \\&\quad \times \,\delta \left( z-e^{i{\varphi }({\mathbf {r}}\tau )}\right) \delta \left( \bar{z}-e^{-i{\varphi }({\mathbf {r}}\tau )}\right) =1. \end{aligned}$$

(66)

The phase–phase propagator $G_{z}({\mathbf {r}}\tau ,{\mathbf {r}}'\tau ')$ will be rewritten in terms of $z({\mathbf {r}}\tau )$-variables as follows:

$$\begin{aligned} G_{z}({\mathbf {r}}\tau ,{\mathbf {r}}'\tau ')=\left\langle z({\mathbf {r}}\tau )\bar{z}({\mathbf {r}}'\tau ')\right\rangle . \end{aligned}$$

(67)

Furthermore, the variables $z({\mathbf {r}}\tau )$ play the role of the phase-flux attached to the fermions (see discussions in Sect. 2). In general case, the local expression of the phase–phase correlation function in Eq. (67) is equal to unity, but, at very low temperatures (especially at $T=0$), this law breaks down, because we have to consider the symmetry breaking related to the bosonic sector; thus, critically, we have fluctuation form $z({\mathbf {r}}\tau )=\left\langle e^{i\varphi }({\mathbf {r}}\tau )\right\rangle +\tilde{z}({\mathbf {r}}\tau )$, and the unimodularity constraint for $z$-field is violated.

Indeed, in the very low-temperature limit, considering the BEC of excitons, we have the spontaneous breaking of local U(1) gage symmetry, related to the phase field, leading to the non-vanishing expectation value of the $\left\langle e^{i\varphi }({\mathbf {r}}\tau )\right\rangle$. In order to demonstrate this, we separate the single-particle states ${\mathbf {k}}=0$ via the Bogoliubov displacement operation (see for details in Refs. [1, 47, 51]). This is so-called Bogoliubov phase coherence mechanism discussed in details in Ref. [1]. Then, we write for the complex variables $z({\mathbf {k}},\omega _{n})$

$$\begin{aligned} z({\mathbf {k}},\omega _{n})=\beta {N}\psi _{0}\delta _{{\mathbf {k}},0}\delta _{\omega _{n},0}+\tilde{z}({\mathbf {k}},\omega _{n})(1-\delta _{{\mathbf {k}},0})\times (1-\delta _{\omega _{n},0}), \end{aligned}$$

(68)

where $\psi _{0}$ is the BEC transition amplitude $\psi _{0}=\left\langle {z}({\mathbf {k}},\omega _{n})\right\rangle$. Next, $\tilde{z}({\mathbf {k}},\omega _{n})$ are the excitation part [54] (on-condensate) of effective bose-field. The Fourier transformation of the phase-phase propagator $G_{z}({\mathbf {r}}\tau ,{\mathbf {r}}'\tau ')$ in Eq. (47) is

$$\begin{aligned} G_{z}({\mathbf {k}}\omega _{n})=\frac{1}{\beta {N}}\Big \langle z\big ({\mathbf {k}},\omega _{n}\big )\bar{z}\big ({\mathbf {k}},\omega _{n}\big )\Big \rangle . \end{aligned}$$

(69)

We consider the expectation value $\left\langle z({\mathbf {k}},\omega _{n})\bar{z}({\mathbf {k}},\omega _{n})\right\rangle$ in the local limit, i.e., when ${\mathbf {d}}={\mathbf {r}}-{\mathbf {r}}'=0$ and $\tau -\tau '=0$ and we should draw the condensate part, by applying the transformation in Eq. (68). Hence, we have

$$\begin{aligned} G_{z}\big ({\mathbf {k}},\omega _{n}\big )&= \frac{1}{\beta {N}}\Big \langle z\big ({\mathbf {k}},\omega _{n}\big )\bar{z}({\mathbf {k}},\omega _{n})\Big \rangle \nonumber \\&= \beta {N}|\psi _{0}|^{2}\cdot \delta _{{\mathbf {k}},0}\delta _{\omega _{n},0}+\tilde{G}_{z}\big ({\mathbf {k}},\omega _{n}\big ). \end{aligned}$$

(70)

Thereby, in Eq. (70), we have defined the coherent macroscopic state for the excitonic system in the low-temperature limit, and an excitonic BEC is expected in the next. The Fourier Green function $G_{z}({\mathbf {k}},\omega _{n})$, defined in Eq. (47), could be calculated within the quantum rotor model, and we give here only the final result (for details, see in Ref. [47]).

$$\begin{aligned} G_{z}({\mathbf {k}},\omega _{n})=\frac{1}{\gamma ^{-1}(\omega _{n})-4J\epsilon ({\mathbf {k}})-\lambda }, \end{aligned}$$

(71)

where $\gamma ^{-1}(\omega _{n})$ is the inverse of the Fourier transformation of the two-point phase–phase correlation function [47]. We have

$$\begin{aligned} \gamma (\omega _{n})=\frac{8}{UZ_{0}}\sum ^{+\infty }_{{m}=-\infty }\frac{e^{-\frac{U\beta }{4}\left( {{m}}-\frac{2\bar{\mu }}{U}\right) ^{2}}}{1-16\left[ \frac{i\omega _{n}}{U}-\frac{1}{2}\left( {{m}}-\frac{2\bar{\mu }}{U}\right) \right] ^{2}}, \end{aligned}$$

(72)

where $Z_{0}$ is the partition function of the non-interacting Bose sector

$$\begin{aligned} Z_{0}=\sum ^{+\infty }_{{m}=-\infty }e^{-\frac{U\beta }{4}\left( {{m}}-\frac{2\bar{\mu }}{U}\right) ^{2}}. \end{aligned}$$

(73)

The summations, in Eqs. (72) and (73), run over topological winding numbers $m$ of the group U(1). Now, the local constraint on the $z$-variables will be rewritten as

$$\begin{aligned} 1-|\psi _{0}|^{2}=\frac{1}{\beta {N}}\sum _{\begin{array}{c} {\mathbf {k}}\ne 0 \\ \omega _{n}\ne 0 \end{array}}G_{z}({\mathbf {k}},\omega _{n}). \end{aligned}$$

(74)

Putting here the expression of the Fourier transform $G_{z}({\mathbf {k}},\omega _{n})$ from Eq. (71) and performing the Bose–Matsubara frequency summations in Eq. (74), we obtain the equation

$$\begin{aligned} 1-|\psi _{0}|^{2}=\frac{U}{4N}\sum _{{\mathbf {k}}} \frac{n_\mathrm{B}\left( \zeta _{1{\mathbf {k}}}\right) -n_\mathrm{B}\left( \zeta _{2{\mathbf {k}}}\right) }{\sqrt{\bar{\mu }^{2}+2Ug_\mathrm{B}{J}\left[ \epsilon ({{\mathbf {0}}})-\epsilon ({{\mathbf {k}}})\right] }}, \end{aligned}$$

(75)

where $n_\mathrm{B}\left( \epsilon \right)$ is the Bose–Einstein distribution function $n_\mathrm{B}\left( \epsilon \right) =1/\left( e^{\beta \epsilon }-1\right)$ and the variables $\zeta _{1{\mathbf {k}}}$ and $\zeta _{2{\mathbf {k}}}$ are defined as

$$\begin{aligned} \zeta _{\alpha {\mathbf {k}}}=-\bar{\mu }-(-1)^{\alpha }{\sqrt{\bar{\mu }^{2}+2U{J}g_\mathrm{B}\left[ \epsilon ({{\mathbf {0}}})-\epsilon ({{\mathbf {k}}})\right] }}. \end{aligned}$$

(76)

Here $\alpha =1,2$. We see also that at the fundamental state with ${\mathbf {k}}=0$, there is a residual gap

$$\begin{aligned} \Delta W=\zeta _{1}\left( {\mathbf {0}}\right) -\zeta _{2}\left( {\mathbf {0}}\right) =2|\bar{\mu }| \end{aligned}$$

(77)

developing, which is related to the condensate state. It equals exactly the binding energy of a molecule in the BEC limit $E_\mathrm{bind}\approx |2\bar{\mu }|$ [52, 53]. In fact, we have obtained in Eq. (75) a useful relation from which we can get the critical temperature $T_{c}$ of the excitonic condensate phase transition. Putting $|\psi _{0}|^{2}=0$ (thus passing to the critical line of the excitonic condensate phase transition) in Eq. (75), we obtain the equation for $T_{c}$

$$\begin{aligned} \frac{U}{4N}\sum _{{\mathbf {k}}} \frac{n_\mathrm{B}\left( \zeta _{1{\mathbf {k}}}\right) -n_\mathrm{B}\left( \zeta _{2{\mathbf {k}}}\right) }{\sqrt{\bar{\mu }^{2}+2Ug_\mathrm{B}{J}\left[ \epsilon ({{\mathbf {0}}})-\epsilon ({{\mathbf {k}}})\right] }}=1, \end{aligned}$$

(78)

The solution for $T_{c}$ for the hopping amplitude $\tilde{t}=-0.3$ is plotted in Fig. 2 in Sect. 3. For more details about the discussion of the phase diagram obtained, please see in Ref. [47].

The retarded bosonic Green function [54] is related to the Matsubara Green function, by the analytical continuation

$$\begin{aligned} G^{R}_{z}({\mathbf {k}},\omega )=G_{z}({\mathbf {k}},i\omega _{n})|_{i\omega _{n}\rightarrow \omega +i\eta }. \end{aligned}$$

(79)

And the ${\mathbf {k}}$-summed DOS for bosons reads as

$$\begin{aligned} \rho _{z}(\omega )=-\frac{1}{\pi }\sum _{{\mathbf {k}}}\mathrm{Im }G^{R}_{z}({\mathbf {k}},\omega ). \end{aligned}$$

(80)

After non -difficult algebraic manipulations and replacing the summation in Eq. (80) by integration with the help of 3D DOS $\rho _{3D}\left( x\right) =\frac{1}{N}\sum _{{\mathbf {k}}}\delta (x-\gamma _{{\mathbf {k}}})$, we get

$$\begin{aligned} \rho _{z}(\omega )&= -\int ^{+\infty }_{-\infty }\mathrm{d}x \frac{U\rho _{3D}\left( x\right) }{4}\cdot \left[ \frac{\delta \left[ \omega -\kappa _{1}(x)\right] }{\sqrt{\bar{\mu }^{2}+4UJ\left( 3-x\right) }}\right. \nonumber \\&\left. +\,\frac{\delta \left[ \omega -\kappa _{2}(x)\right] }{\sqrt{\bar{\mu }^{2}+4UJ\left( 3-x\right) }}\right] , \nonumber \\ \end{aligned}$$

(81)

where $\kappa _{i}(x)$ $i=1,2$ are given by the following relations $\kappa _{1,2}(x)=-\bar{\mu }\pm {\sqrt{\bar{\mu }^{2}+4U{J}\left( 3-x\right) }}$ and the stiffness parameter $J$ is given in Eq. (39) in Sect. 3. As it could be expected, the Bosonic DOS function Eq. (81) is negative $\rho _{z}(\omega )< 0$. This is consistent with the general considerations of the weakly non-ideal Bose gas [54].

4.3 Spectral Density Functions and Fermionic DOS

Furthermore, we separate the condensate modes $\left\{ {\mathbf {q}}={{0}},\omega _{n}=0 \right\}$ in Eqs. (48) and (49). We have

$$\begin{aligned} G_\mathrm{x,x}\big ({\mathbf {k}},\nu _{n}\big )&= |\psi _{0}|^{2} \cdot G_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n}) \nonumber \\&+\,\frac{1}{\beta {N}}\sum _{\begin{array}{c} {\mathbf {q}}\ne 0 \\ \omega _{n}\ne 0 \end{array}}\tilde{G}_{z}({\mathbf {q}},\omega _{n})\cdot G_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}}-{\mathbf {q}},\nu _{n}-\omega _{n}) \end{aligned}$$

(82)

and

$$\begin{aligned} G_\mathrm{c,f}({\mathbf {k}},\nu _{n})&= |\psi _{0}|^{2} \cdot G_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu _{n}) \nonumber \\&+\,\frac{1}{\beta {N}}\sum _{\begin{array}{c} {\mathbf {q}}\ne 0 \\ \omega _{n}\ne 0 \end{array}}\tilde{G}_{z}({\mathbf {q}},\omega _{n})\cdot G_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}}-{\mathbf {q}},\nu _{n}-\omega _{n}). \end{aligned}$$

(83)

As we see, the normal and excitonic propagators are composed of two parts, one responsible for the condensate state and the other on condensate excitation part (see discussion in Ref. [54] for the case of the pure Bose gas). Note also that first terms in the right-hand sides in Eqs. (82) and (83) consist of the condensate-transition probability function $|\psi _{0}|^{2}$, multiplied with the fermionic propagators $G_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n})$ and $G_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu _{n})$.

Now, we are ready to calculate the analytical forms of the normal excitonic spectral functions $A_\mathrm{x,x}({\mathbf {k}},\nu )$ ($x=f,c$) and anomalous excitonic spectral function $A_\mathrm{c,f}({\mathbf {k}},\nu )$ and, later on, the profiles of the respective DOS, including states of the condensate. We introduce here the spectral functions $A_\mathrm{x,x}({\mathbf {k}},\nu )$ and $A_\mathrm{c,f}({\mathbf {k}},\nu )$ that carries the same physical information as the correlation functions $G_\mathrm{x,x}({\mathbf {k}},\nu _{n})$ and $G_\mathrm{c,f}({\mathbf {k}},\nu _{n})$. We have

$$\begin{aligned} G_\mathrm{x,x}({\mathbf {k}},\nu _{n})=\int ^{+\infty }_{-\infty }\mathrm{d}\nu '\frac{A_\mathrm{x,x}({\mathbf {k}},\nu ')}{i\nu _{n}-\nu '}, \end{aligned}$$

(84)

$$\begin{aligned} G_\mathrm{c,f}({\mathbf {k}},\nu _{n})=\int ^{+\infty }_{-\infty }\mathrm{d}\nu '\frac{A_\mathrm{c,f}({\mathbf {k}},\nu ')}{i\nu _{n}-\nu '}. \end{aligned}$$

(85)

The integration here is over continuous frequencies. Note that $G_\mathrm{x,x}({\mathbf {k}},\nu _{n})$ and $G_\mathrm{c,f}({\mathbf {k}},\nu _{n})$ are total fermionic Green functions, including also the convolutions with bosonic parts. In the same way, we can introduce the spectral functions $A_{z}({\mathbf {k}},\nu )$, $A_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu )$ and $A_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu )$, associated with the charge and the pure fermionic parts (without bosonic sector). They correspond, respectively, to the correlation functions $G_{z}({\mathbf {k}},\omega _{n})$, $G_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n})$ and $G_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu _{n})$. We have the following equations for these counterparts:

$$\begin{aligned} G_{z}({\mathbf {k}},\omega _{n})=\int ^{+\infty }_{-\infty }\mathrm{d}\nu '\frac{A_{z}({\mathbf {k}},\nu ')}{i\omega _{n}-\nu '}, \end{aligned}$$

(86)

$$\begin{aligned} G_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu _{n})=\int ^{+\infty }_{-\infty }\mathrm{d}\nu '\frac{A_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu ')}{i\nu _{n}-\nu '} \end{aligned}$$

(87)

and

$$\begin{aligned} G_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu _{n})=\int ^{+\infty }_{-\infty }\mathrm{d}\nu '\frac{A_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu ')}{i\nu _{n}-\nu '}. \end{aligned}$$

(88)

Using these definitions, we get for the total spectral density functions $A_\mathrm{x,x}({\mathbf {k}},\nu )$ and $A_\mathrm{c,f}({\mathbf {k}},\nu )$ (see Appendix)

$$\begin{aligned} A_\mathrm{x,x}({\mathbf {k}},\nu )&= |\psi _{0}|^{2}\cdot A_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}},\nu ) \nonumber \\&-\,\frac{1}{N}\sum _{{\mathbf {q}}\ne 0}\int {\mathrm{d}\nu '}A_\mathrm{z}({\mathbf {q}},\nu ')A_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}({\mathbf {k}}-{\mathbf {q}},\nu -\nu ') \nonumber \\&\times \,\left[ n(\nu ')+f(\nu -\nu ')\right] \end{aligned}$$

(89)

and

$$\begin{aligned} A_\mathrm{c,f}({\mathbf {k}},\nu )&= |\psi _{0}|^{2}\cdot A_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}},\nu ) \nonumber \\&-\,\frac{1}{N}\sum _{{\mathbf {q}}\ne 0}\int {\mathrm{d}\nu '}A_\mathrm{z}({\mathbf {q}},\nu ')A_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}({\mathbf {k}}-{\mathbf {q}},\nu -\nu ') \nonumber \\&\times \,\left[ n(\nu ')+f(\nu -\nu ')\right] . \end{aligned}$$

(90)

The proof of these relations is given in Appendix. From the spectral functions, we can obtain the corresponding DOS, by summing over the reciprocal wave vectors ${\mathbf {k}}$; hence, the total DOS are $\rho _\mathrm{x,x}(\nu )=\frac{1}{N}\sum _\mathrm{{\mathbf {k}}}A_\mathrm{x,x}({\mathbf {k}},\nu )$ and $\rho _\mathrm{c,f}(\nu )=\frac{1}{N}\sum _\mathrm{{\mathbf {k}}}A_\mathrm{c,f}({\mathbf {k}},\nu )$.

Furthermore, using the expressions for $A_\mathrm{x,x}({\mathbf {k}},\nu )$ and $A_\mathrm{c,f}({\mathbf {k}},\nu )$ in Eqs. (89) and (90), we get for the total DOS functions

$$\begin{aligned} \rho _\mathrm{x,x}(\nu )=|\psi _{0}|^{2}\cdot \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu )+\tilde{\rho }_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu ) \end{aligned}$$

(91)

and

$$\begin{aligned} \rho _\mathrm{c,f}(\nu )=|\psi _{0}|^{2}\cdot \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )+\tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu ), \end{aligned}$$

(92)

where $\tilde{\rho }_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu )$ and $\tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )$ are DOS, corresponding to the excitation part of the system and are given, as convolutions, in terms of continuous frequency modes (see Appendix)

$$\begin{aligned} \tilde{\rho }_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu )&= -\int ^{+\infty }_{-\infty }\mathrm{d}\nu ' \rho _{z}(\nu ')\rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu -\nu ') \nonumber \\&\,\times \left[ n\left( \nu '\right) +f\left( \nu -\nu '\right) \right] \end{aligned}$$

(93)

and

$$\begin{aligned} \tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )&= -\int ^{+\infty }_{-\infty }\mathrm{d}\nu ' \rho _{z}(\nu ')\rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu -\nu ') \nonumber \\&\times \, \left[ n\left( \nu '\right) +f\left( \nu -\nu '\right) \right] . \end{aligned}$$

(94)

A key feature of the results in Eqs. (A.4) and (92) is that we have separated DOS contributions coming from the condensate and excitation parts. We define also the total DOS function as

$$\begin{aligned} \rho (\nu )=\sum _{x=f,c}\rho _\mathrm{x,x}(\nu ). \end{aligned}$$

(95)

In Sect. 4.4, we present the numerical evaluations of all discussed DOS functions.

4.4 Total DOS Functions

Employing Eqs. (56), (61), and (81), we can obtain the explicit analytical expressions for DOS functions in Eqs. (A.4) and (92). For the normal $f$- and $c$-band excitonic DOS functions, we obtain

$$\begin{aligned} \rho _\mathrm{x,x}(\nu )&= |\psi _{0}|^{2}\cdot \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu ) \nonumber \\&-\,U\int ^{+3}_{-3}\mathrm{d}x \frac{\rho _{3D}(x)}{4\sqrt{\bar{\mu }^{2}+4UJ\left( 3-x\right) }} \nonumber \\&\times \left\{ \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( \nu -\kappa _{1}\left( x\right) \right) \cdot \left[ n\left( \kappa _{1}(x)\right) +f\left( \nu -\kappa _{1}(x)\right) \right] \right. \nonumber \\&+\left. \rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}\left( \nu -\kappa _{2}\left( x\right) \right) \cdot \left[ n\left( \kappa _{2}(x)\right) +f\left( \nu -\kappa _{2}(x)\right) \right] \right\} . \end{aligned}$$

(96)

For the anomalous excitonic DOS function, we have

$$\begin{aligned} \rho _\mathrm{c,f}(\nu )&= |\psi _{0}|^{2}\cdot \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu ) \nonumber \\&-\,U\int ^{+3}_{-3}dx \frac{\rho _{3D}(x)}{4\sqrt{\bar{\mu }^{2}+4UJ\left( 3-x\right) }} \nonumber \\&\times \left\{ \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}\left( \nu -\kappa _{1}\left( x\right) \right) \cdot \left[ n\left( \kappa _{1}(x)\right) +f\left( \nu -\kappa _{1}(x)\right) \right] \right. \nonumber \\&+\,\left. \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}\left( \nu -\kappa _{2}\left( x\right) \right) \cdot \left[ n\left( \kappa _{2}(x)\right) +f\left( \nu -\kappa _{2}(x)\right) \right] \right\} . \end{aligned}$$

(97)

The numerical evaluations of calculated DOS functions at $T=0$ are given in Figs. 3, 4, 5, 6, and 7 for $\tilde{t}=-0.3$. The presence of singularities in the integration region causes that we used an adaptive 21-point integration routine combined with the Wynn $\epsilon$-algorithm [55] to calculate those integrals numerically. The accuracy for adaptive evaluations is achieved with an absolute error of order of $10^{-4}$ and with a relative error of order of $10^{-7}$.

Fig. 3

Single-particle normal DOS functions (incoherent) for different values of the Coulomb interaction parameter $U$ and for the case $T=0$. In panels I–III, $\tilde{x}=\tilde{f}$ and $\tilde{x}=\tilde{c}$ partial DOS structures are shown, and total DOS is plotted for $\tilde{t}=-0.3$. The chemical potential values are also shown in each panel (Color figure online)

Fig. 4

Single-particle normal DOS functions (incoherent) for different values of the Coulomb interaction parameter $U$ and for the case $T=0$. In panels I–III, $\tilde{x}=\tilde{f}$ and $\tilde{x}=\tilde{c}$ partial DOS structures are shown, and total DOS is plotted for $\tilde{t}=-0.3$. The chemical potential values are also shown in each panel (Color figure online)

Fig. 5

Condensate part $|\psi _{0}|^{2}\cdot \rho _{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}$ of the anomalous excitonic DOS function given in Eq. (97), for different values of the Coulomb interaction parameter $U$ and for the case $T = 0$. The case $\tilde{t}=-0.3$ is considered here, and the values of the functions $|\psi _{0}|^{2}$ were taken from the work in Ref. [47] (Color figure online)

Fig. 6

Total phase-coherent normal DOS function $\rho (\nu )$ given in Eq. (97) for different values of the Coulomb interaction parameter $U$ ($U=4$, $U=6$, $U=8$, $U=9$, $U=10.4$ in the figure) and for the case $T=0$. The case $\tilde{t}=-0.3$ is considered here, and the values of the functions $|\psi _{0}|^{2}$ were taken from the work in Ref. [47] (Color figure online)

Fig. 7

Incoherent (in green) and coherent (in black) normal DOS functions for the $c$-band and for different values of the parameter $U$ ($U=4$ in a, $U=6$ in b and $U=8$ in c). The case $\tilde{t}=-0.3$ is considered here (Color figure online)

Particularly, in Figs. 3 and 4, we have presented purely fermionic normal single-particle (incoherent) DOS $\rho _{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\omega )$. We examine their behavior over the entire BCS–BEC crossover region (i.e., for different values of the Coulomb interaction $U$). An artificial Lorentzian broadening $\eta =0.01$ is used in numerical evaluations for the incoherent partial and total DOS functions of $f$- and $c$ -orbitals. The chemical potential values are inserted along the upper-bound $\mu _\mathrm{max}$, where they are maximal (see in Ref. [47]). The principal reason of it is that the BEC transition amplitude $\psi _{0}$ has no physical solutions along the lower-bound $\mu _\mathrm{min}$ of the chemical potential. On the other hand, the values $\mu _\mathrm{max}$ are most convenient, because they are minimalizing the Hamiltonian of the system. Furthermore, as the reference for the coherent BEC transition amplitude $|\psi _{0}|^{2}$, we considered the self-consistent calculation results from the work in Ref. [47], where this function is calculated both analytically and numerically for different values of the Coulomb interaction parameter $U$ and for different temperatures, including the zero temperature limit. We use here the numerical data, which are evaluated there.

We see in Fig. 3 that in the small interaction limit, when $2 \leqslant U\leqslant 6$ (i.e., the BCS limit), the incoherent excitonic DOS exhibits a BCS-like double-peak structure (see the panels I–III in Fig. 3) and the peaks are separated with a well-defined hybridization gap. The principal reason of it is the non-vanishing Hartree-gap $\Delta _{H}\ne 0$ in the single-particle energy-spectrum discussed above, in Sect. 3. The hybridization gap is proportional to the parameter $U$, and it is increasing with the increase of parameter $U$. We observe also that the peaks become more separated when increasing of $U$. In the strong interaction limit, this displacement is stabilizing and we have practically constant value of the hybridization gap, when further increasing the interaction (see the panels I–III, in Fig. 4). The results in Figs. 3 and  4 are very similar with the previous theoretical results [15, 20, 52, 53]. Especially, they are close to those presented in Refs. [17] and [20], where the partial incoherent $f$- and $c$-band normal DOS functions, and total DOS is calculated using HF and SO(2)-invariant slave boson approaches.

In Fig. 5 (see the panels a–c), we have shown the coherent condensate part of the anomalous excitonic DOS given by the first term in the right-hand side in Eq. (A.4), corresponding to the fundamental state (${\mathbf {k}}=0$). Different values of the Coulomb interaction are considered, including the small and strong interaction cases. Here, we observe again the double-peak fermionic structure, but there is no the hybridization gap for the small and medium values of the Coulomb interaction parameter (see the panel-a and panel-b in Fig. 5), and we have an infinite number of states for all values of the frequency modes $\nu$. This is due to the coherence effects and the presence of the coherent excitonic condensate at the fundamental mode $\nu =0$. For higher values of $U$ in Fig. 5 (see the panel-c), this double-peak structure in DOS is smoothing, but, in contrast to the incoherent normal DOS behavior (see in Figs. 3, 4), here a small Mott-gap appears at the very high values of the Coulomb interaction parameter (see the DOS curves for $U=8$ and $U=9.6$ in the panel-c in Fig. 5). This is due to the fact that the very strong Coulomb interaction has a destructive role on the condensate state and, in the large-$U$ limit of interaction, we have the destruction of the excitonic condensate, and a very small Mott-type hybridization gap is enhanced. In this region of the interaction, we have the coherent exciton DOS separation into two separate parts, similar to the case of the incoherent DOS functions in Figs. 3 and 4.

In the positive (negative) frequency regions, the DOS spectrum in Figs. 3, 4, and 5 is slightly displacing and broadening into the direction of higher (smaller) frequencies, for both normal (incoherent) single-particle fermionic and coherent excitonic DOS functions and, for both, we observe also a gradual decrease in the DOS amplitudes, across the whole BCS–BEC crossover region, when increasing the Coulomb interaction parameter $U$. In Fig. 6, we have shown the total coherent fermionic DOS functions $\rho (\nu )$, given by Eq. (95), for different values of the Coulomb interaction parameter $U$ and in the limit of zero temperature. An artificial Lorentzian broadening $\eta$ = 0.01 is again used during the numerical calculations. Contrary the case of the incoherent DOS in Figs. 3 and 4, the coherent total fermionic DOS in Fig. 6 shows a different, gapless behavior. Here, again, as in the case of the incoherent DOS functions, for small and intermediate interactions $U$, we have typical double-peak fermionic structure in the DOS, which is smoothing and disappearing totally in the strong interaction limit. But, in difference with incoherent DOS functions, here we have not the presence of the hybridization gap in the spectra and we have always a finite number of states for all values of the frequency modes (see in Fig. 6). The reason of this gapless DOS behavior could be related to the strong coherence effects between two bands, which is due to the presence of the phase-stiffness mechanism considered here. It is worth to mention that another gapless-type behavior in the DOS spectrum has been found recently in Ref. [56], where this effect is associated with metallic charge-density-wave phase and is driven by strong electron correlations. In Fig. 6, all figures are combined together, in the way, to see the total DOS evolution with variation of the interaction parameter $U$. We see clearly how the tuning of the interaction parameter $U$ affects the general DOS behavior in the model, by reducing the DOS amplitudes with increasing of $U$ and reducing also the number of states at the Fermi level $\rho (\nu =\epsilon _{F})$ when increasing $U$. Thereby, in the case of normal fermionic DOS (both incoherent and coherent), presented in Figs. 3, 4, 5, and 6, the double-peak structure disappears for $U \gtrsim 6$, signaling the appearance of the SM (BEC) limit of the transition.

In Fig. 7, a comparison is given for the incoherent and coherent $c$-band normal DOS behaviors. Three different values of the Coulomb interaction $U$ are chosen (see the panels a, b, and c). We see clearly in Fig. 7 that how the coherence effects are reducing the incoherent DOS amplitudes and DOS spectra become broader, along frequency axis, and also, there is no gap near the Fermi level ($\rho (\nu \sim \epsilon _{F})\ne 0$).

In Figs. 8 and  9, we have presented the temperature dependence of the $c$-band normal fermionic (Fig. 8) and anomalous excitonic DOS functions (Fig. 9). They are given by the first terms in the right-hand side in Eqs. (96) and (97). In the left-panel in Fig. 8, we have presented the temperature dependence of the single-particle DOS $|\psi _{0}|^{2}\rho _\mathrm{c,c}(\nu )$, which corresponds to the fundamental state ${\mathbf {k}}=0$, and for $U=6$. The case $\tilde{t}=-0.3$ is considered. The corresponding values of the amplitude $\psi _{0}$ of BEC transition are taken again from the work in Ref. [47]. In the right panel in Fig. 9, the same function is plotted for the case $U=9$. As we see in Fig. 8, the hybridization gap is still open for all values of the temperature. So the system is always an insulator.

In Fig. 9, the temperature dependence of the phase-coherent anomalous excitonic condensate part of the DOS is presented (i.e., ${\mathbf {k}}=0$, and without excitation part) for $U=2$ and $U=6$. We observe in all Figs. 8 and  9 that the temperature has a destructive effect on the DOS amplitudes.

We realize also, from Figs. 8 and  9, that the anomalous excitonic DOS functions vanish at the temperatures that are far away from the region of the EI transition of about two orders of magnitude (compare the temperature scales in Figs. 8 and 9 with those given in Fig. 2 in Sect. 3). This result could be regarded as a good proof of the theory elaborated in Ref. [47], and we can theoretically clearly state that the excitonic BEC and EI states are not the same phases of matter. We see also, in Figs. 8 and 9, that the normal single-particle DOS persists for a rather large values of temperature than the anomalous condensate DOS, because of the presence of the hybridization gap in the low-energy spectra.

Fig. 8

Temperature dependence of the normal $c$-band DOS function in Eq. (96) for the case $\tilde{t}=-0.3$ and for different values of the Coulomb interaction parameter $U$ (Color figure online)

Fig. 9

Temperature dependence of the condensate part of the anomalous excitonic DOS function given in Eq. (97) for the case $\tilde{t}=-0.3$ and for different values of the Coulomb interaction parameter $U$ (Color figure online)

5 Final Remarks and Conclusions

We have studied 3D system of conduction band electrons and valence band holes in the frame of the extended Falicov–Kimball model. We have implemented the path-integral formalism, in which the Coulomb interaction term is expressed in terms of U(1) quantum phase variables $\varphi$ conjugated to the local particle number, providing a useful interpretation of the problem. In Sect. 3, we have shown that at low temperatures, the electron–hole system becomes unstable with respect to the formation of the excitons at $T=T_\mathrm{EI}$, and the local gap $\Delta$ is present in the excitation spectrum, controlled by the Coulomb interaction parameter $U/t$, which gives the relevant energy scales for the excitonic insulator state. Here, as a result of the spontaneous symmetry breaking, an expectation value of $\langle e^{i\varphi }\rangle \ne 0$ appears, which is signaling of the presence of the phase coherence in the system. Furthermore, pairing and condensation are not generally the same, as it was admitted in the literature, except the weak interaction limit, when we have a BCS-like condensation of excitonic pairs. However, in the excitonic system with the strong pairing, we have the situation where the pairs are strongly bound, but are uncorrelated one with each other, until they become phase coherent at temperatures $T\lesssim T_{c}$. This situation was studied in detail in Ref. [47].

We have evaluated the normal and anomalous excitonic spectral functions for the $f$- and $c$-band. We have determined, both analytically and numerically, DOS spectra, governed by the pure fermionic part (due to the condensate modes ${\mathbf {k}}=0$), and the excitation spectra. We have shown that there is a usual hybridization gap in the normal (incoherent) $f$- and $c$-band DOS structures. Contrary, in the case of the coherent normal DOS functions (presented in Fig. 6), this gap is lacking, and there is always a finite number of states at all frequency modes $\nu$. We associate this result to the strong coherence effects, which are present in the system at low temperatures.

In the anomalous excitonic DOS structure, we have found that the hybridization gap is absent for the weak and intermediate values of the Coulomb interaction parameter $U$ and this is due to the presence of the coherent excitonic condensate and strong coherence effects. A very small gap is opening in the spectra in the strong interaction limit, signaling the destruction of the coherence and condensate state.

The excitonic phase coherence may be evidenced by the coherence of their light emission, which can be studied by interferometry measurements [57]. Therefore, measurements of the intensity of the line-shape of the excitons decay (by emitting photons, upon electron–hole recombination) may be a powerful probe of the DOS spectra in the excitonic systems [58].

However, a final remark that should be featured is such that for the experimental determination of the excitonic Bose–Einstein condensate, the ARPES measurements should be provided at temperatures much lower than temperatures at which the excitonic insulator state is being determined [20]. This is important for the achievement of macroscopic phase coherence between excitonic pairs, and the strong coherence effects are manifesting at the very low temperatures.

As a continuation of our studies, we would like to consider the role of the charge-density-wave-like excitations in the excitonic systems, within the EFKM model and consider the temperature effects. It is especially interesting to find out how the coherent excitonic DOS will be affected in the case of the presence of such elementary excitations. This will give us a more complete picture on the excitonic phase transition scenario.

Appendix: Spectral Functions and DOS

1.1 A Convolution for DOS

Here, we give a short derivation of Eqs. (89) and (90) presented in Sect. 4.3. To this end, we use the denitions in Eqs. (86)–(88) and the convolution form of the Green functions given in Eqs. (82) and (83). Then, it is easy to see that the following identity holds for the integrals of the normal $x=f, c$-band fermionic spectral functions:

$$\begin{aligned} \int \mathrm{d}\nu \frac{A_{x,x}\left( {\mathbf {k}},\nu \right) }{i\nu _{n}-\nu }&= |\psi _{0}|^{2}\cdot \int \mathrm{d}\nu \frac{A_{\tilde{x},\tilde{x}}\left( {\mathbf {k}},\nu \right) }{i\nu _{n}-\nu } \nonumber \\&+\,\frac{1}{\beta {N}}\sum _{{\mathbf {q}}\ne 0, \omega _{n}\ne 0}\int \mathrm{d}\nu \int \mathrm{d}\nu '\left[ \frac{1}{i\omega _{n}-\nu '}\right. \nonumber \\&\left. -\,\frac{1}{i\omega _{n}-i\nu _{n}+\nu '-\nu }\right] \cdot A_{z}\left( {\mathbf {q},\nu '}\right) \cdot A_{\tilde{x}\tilde{x}}\left( {\mathbf {k}}-{\mathbf {q}},\nu -\nu '\right) .\nonumber \\ \end{aligned}$$

(A.1)

Furthermore, for calculating the Matsubara sums over bosonic frequencies $\omega _{n}$, we will use the property of the Bose–Einstein distribution function $n\left( \epsilon \right)$ (see Sect. 4.3) of complex argument

$$\begin{aligned} n\left( i\nu _{n}+\nu -\nu '\right) =-f\left( \nu -\nu '\right) , \end{aligned}$$

(A.2)

where $f\left( \epsilon \right)$ is the Fermi–Dirac distribution function. In Eq. (A.2), we used the fact that $\nu _{n}$ are even Fermion–Matsubara frequencies $\nu _{n}=\frac{\pi }{\beta }\left( 2n+1\right)$. Next, we sum the Bosonic Matsubara frequencies in Eq. (A.1) amd we rewrite the equality in Eq. (A.1) in the following form:

$$\begin{aligned} A_{x,x}\left( {\mathbf {k}},\nu \right)&=|\psi _{0}|^{2}\cdot A_{\tilde{x},\tilde{x}}\left( {\mathbf {k}},\nu \right) \nonumber \\&\quad -\frac{1}{N}\sum _{{\mathbf {q}}\ne 0}\int \mathrm{d}\nu 'A_{z}\left( {\mathbf {q}},\nu '\right) A_{\tilde{x},\tilde{x}}\left( {\mathbf {k}}-{\mathbf {q}},\nu -\nu '\right) \cdot \left[ n\left( \nu '\right) \!+\!f\left( \nu -\nu '\right) \right] . \end{aligned}$$

(A.3)

The derivation of the convolution form, given in Eq. (90), is exactly the same.

Summing over the wave vectors ${\mathbf {k}}$ in Eq. (A.3), we will have the forms of the normal $f$ and $c$-bands density of states $\rho _{x,x}\left( \nu \right)$ and also the excitonic DOS function $\rho _{c,f}\left( \nu \right)$ given in Eqs. (A.4) and (92) in Sect. 4.3.

1.2 Normalization Conditions for DOS

We have a composed nature of the interacting electron, given in Sect. 2; thus, we have various sum rules, corresponding to different counterparts of the total coherent DOS functions. As a general rule, the partial normal DOS and anomalous excitonic DOS functions $\rho _{x,x}\left( \nu \right)$ and $\rho _{c,f}\left( \nu \right)$ satisfy following normalization conditions:

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\nu \rho _\mathrm{f,f}\left( \nu \right) =n_\mathrm{f}, \nonumber \\ \int ^{+\infty }_{-\infty }\mathrm{d}\nu \rho _\mathrm{c,c}\left( \nu \right) =n_\mathrm{c}, \end{aligned}$$

(A.4)

and

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\nu \rho _\mathrm{c,f}\left( \nu \right) =1. \end{aligned}$$

(A.5)

For the total normal band DOS function $\rho (\nu )$, given in Eq. (95), we have

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\nu {\rho }\left( \nu \right) =1. \end{aligned}$$

(A.6)

The functions $\tilde{\rho }_{\tilde{\mathrm{x}},\tilde{\mathrm{x}}}(\nu )$ and $\tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )$ in Eqs. (93) and (94) are excitation parts [54] of the respective total DOS functions in Eqs. (A.4) and (92), i.e., without condensate modes ${\mathbf {k}}=0$; therefore, the following normalization conditions hold

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\nu \tilde{\rho }_{\tilde{\mathrm{f}},\tilde{\mathrm{f}}}(\nu )=n_{f}-|\psi _{0}|^{2}, \nonumber \\ \int ^{+\infty }_{-\infty }\mathrm{d}\nu \tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{c}}}(\nu )=n_{c}-|\psi _{0}|^{2} \end{aligned}$$

(A.7)

and

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\nu \tilde{\rho }_{\tilde{\mathrm{c}},\tilde{\mathrm{f}}}(\nu )=1-|\psi _{0}|^{2}. \end{aligned}$$

(A.8)

In contrary, the integral over the charge-DOS function $\rho _{z}(\omega )$, given in Eq. (76), vanishes, because it is antisymmetric: $\rho _{z}(\omega )=-\rho _{z}(-\omega )$. Then

$$\begin{aligned} \int ^{+\infty }_{-\infty }\mathrm{d}\omega \rho _{z}(\omega )=0. \end{aligned}$$

(A.9)

The norms of the full spectral density functions in Eqs. (89) and (90) depend on the constraints in the charge and spin bosonic sectors.