Excitonic Phase Transition in the Extended Three-Dimensional Falicov–Kimball Model

We study the excitonic phase transition in a system of the conduction band electrons and valence band holes described by the three-dimensional (3D) extended Falicov–Kimball (EFKM) model with the tunable Coulomb interaction U\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U$$\end{document} between both species. By lowering the temperature, the electron–hole system may become unstable with respect to the formation of the excitons, i.e, electron–hole pairs at temperature T=TΔ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T=T_{\Delta }$$\end{document}, exhibiting a gap Δ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta $$\end{document} in the particle excitation spectrum. To this end we implement the functional integral formulation of the EFKM, where the Coulomb interaction term is expressed in terms of U(1) phase variables conjugate to the local particle number, providing a useful representation of strongly correlated system. The effective action formalism allows us to formulate a problem in the phase-only action in the form of the quantum rotor model and to obtain analytical formulas for the critical lines and other quantities of physical interest like charge gap, chemical potential and the correlation length.

excitonic condensate are not exactly the same [32][33][34]. The author in Refs. [32] and [33] shows from general considerations that in the low density limit of the excitonic pairs, the critical temperature T c of excitonic BEC should be much smaller than the temperature T of the EP formation. This is in contrast with previous treatments [18,23,24,28], where the EI state is associated with the BEC state of excitons. Similarly, in Ref. [34] it is shown that the EI state is an excitonium state, where the 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 preformed excitonic pairs. Obviously, in the low density limit, the gas of free excitons undergoes the BEC phase transition at very low temperatures, and the BEC temperature transition line is not coinciding with that of EP formation. The Bose condensation of the excitonic pairs is possible only when the macroscopic phase coherence is present in the system [32].
Contrary, at high e-h density, where the mean distance between the particles is shorter than the excitonic Bohr radius, the weakly bound e-h pairs behave like the Cooper pairs in the conventional superconductors and at sufficiently low temperatures, i.e., the BCS state of e-h pairs [8,11,35]. Therefore, an expected BCS-BEC crossover represents actually a fascinating problem typical to the excitonic systems. Especially, it is interesting from the viewpoint of the difference from similar crossover in superconductors, or the trapped atomic Fermi gases [35][36][37][38]. The transition to e-h pair condensed phase, in the weak-coupling limit, is related to the relative motion between electrons and holes [34], 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. 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 EP formation temperature [34].
In the present paper we explore the quantum collective behavior of the excitons in 3D system going beyond the simple HF method. To this end, we study the excitonic phase transition in a system composed of the conduction band electrons and valence band holes, described by the 3D extended Falicov-Kimball [39,40] model with tunable Coulomb interaction U between both species of particles. We implement the quantum rotor approach, where the Coulomb interaction of the EFKM model is expressed in terms of U(1) quantum phase variables conjugate to the local particle number, providing a useful representation of strongly correlated systems. This allows us to obtain the analytical formulas for the critical lines and other quantities of physical interest like the charge gap, chemical potential and the correlation length. We present also the numerical evaluations of all physical quantities discussed in the paper.
The plan of the paper is as follows: in the Sect. 2 we provide the Hamiltonian of the model EFKM, then in the Sect. 3 we introduce the new decoupling potentials and we handle with four fermion interaction term in the initial Hamiltonian. Furthermore, in the Sect. 4, we obtain the transition temperature of the excitonic pair formation, excitonic gap parameter, the charge gap and other important physical quantities. At the end of that sect. we discuss the numerical results. In Sect. 5 we obtain the effective phase action in the context of the quantum rotor approach and we derive the equation for the excitonic BEC transition critical temperature. Numerical results are also discussed there. The momentum distribution functions and the excitonic coherence length are calculated in Sect. 6. Conclusions are given in Sect. 7. A number of technical details is given in the "Appendices".

The Hamiltonian
We consider the Hamiltonian of the extended FKM wherec(r) (c(r)) are the creation (annihilation) operators of the electron of the corbitals at the site with the position r and rr runs over pairs of the nearest neighbor (n.n.) sites on a 3D cubic lattice. Furthermore t is the hopping integral for the celectrons and c is the on-site energy level. Similarly,f (r) ( f (r)) are the creation (annihilation) operators of the f -orbital electrons andt is the hopping integral for the f -electrons. 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, thus replacing the fermion Hilbert-space by the pseudo-fermionic one, and by linearizing the interaction term via bosonic states [23]. Furthermore, f is the on-site energy level of the f -orbital and μ is the chemical potential. The equilibrium value of chemical potential μ will be determined from the half-filling condition, i.e., n c = 1 − n f , where n x ≡ n x (r) is the average particle density with x = c, f for the c and f -orbital electrons respectively. Furthermore, we suppose that the chemical potentials of both orbitals are the same, as in the work in Ref. [25]. The parameter U , which enters in the last term of the Hamiltonian, is the Coulomb repulsion between two types of electrons. Furthermore n c (r) and n f (r) are the c-and f -electron density operators and they are defined as usual by the relation n x (r) =x(r)x(r). We consider also the following values for the band parameters c = 0 and f = −1. With this consideration the c-electrons are itinerant and the f -electrons are quasilocalized on the atomic sites. Throughout the paper we set k B = 1 andh = 1 and lattice constant a = 1.

The Method
In the first step, we transform the fermionic interaction term in the Hamiltonian by rewriting the density product in the last term in Eq. (1) in the equivalent form where we introduced the short-hand notations With the new notations we can rewrite the Hamiltonian of the system in Eq. (1) as We have put hereμ = μ −¯ and¯ = c + f /2 is the average energy level parameter. The Hamiltonian in Eq. (5) is now suitable for decoupling quadratic density terms using the Gaussian path integral method [41].

Functional Integral Formalism: Decoupling of Interactions
Dealing with fermions within the path integral method, requires introduction of the Grassmann variables c(rτ ) and f (rτ ) at each site r and at each imaginary time τ . The latest is varying in the interval 0 ≤ τ ≤ β, where β = 1/T (with T being the thermodynamic temperature). The variables c(rτ ) and f (rτ ) satisfy the anti-periodic boundary conditions x(rτ ) = −x(rτ + β). The partition function of system of the fermions, written as a functional integral over the Grassmann field, is where the action in exponential is given in the path integral formulation as Here S B [x, x] is the fermionic Berry term for the c and f -electrons. It is defined as Next, we decouple quadratic density terms in Eq. (5) using the Hubbard-Stratonovich (HS) transformation [41] and by introducing the new variables V (rτ ) and (rτ ) conjugated to the density terms n(rτ ) andñ(rτ ) respectively. For the quadratic term proportional to n 2 (rτ ), in the exponential of the partition function in Eq. (6), we have After combining the exponential in Eq. (9) with the term linear in total electron density n(r) in Eq. (5), we can decompose the variables V (rτ ) into the static and periodic parts where β 0 dτṼ (rτ ) = 0. As a result, the integration over V (rτ )-variables becomes now the integration over the scalar static variables V 0 (r) and the integration over the periodic fieldṼ (rτ ): For the periodic part in Eq. (10), we introduce the U(1) phase field variables φ(rτ ) using a Faraday-type relation [42] V (rτ ) = ∂ϕ(rτ ) ∂τ ≡φ(rτ ).
Thus, for the dynamic part, we transform the integration over the gauge variables V (rτ ) into the integration over the generic phase variables ϕ(rτ ) The periodicity ofṼ (rτ ) implies that φ (rβ) = φ (r0). The integration measure in Eq. (13) over φ variables is defined as where the notations φ i and φ f mean the initial and final paths. The path integral in Eq. (14) could be transformed into path integration over the compact U(1) group manifold, since the electromagnetic U(1) group governing the phase field is compact, i.e. φ(rτ ) has the topology of a circle (S 1 ), thus we have a non-homotopic mapping of the configuration space onto the U(1) gauge group S 1 → 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 m (r). 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. Accordingly, the path integral in Eq. (14) is transformed as The integration measure in Eq. (15) is now In performing the integration over the phase field one should take into account that the field configurations satisfy the boundary conditions [43,44] ϕ(rβ) − ϕ(r0) = 2π m(r).
Thus, integration over all phases φ(rτ ) amounts the integration over the β-periodic field ϕ(rτ ) and the summation over a set of U(1) winding numbers m(r). For the scalar static part V 0 (r), we have the following functional integral The saddle-point value of V s. p.
where n is total average particle density n = n c + n f . And we have the contribution in the partition function in Eq.
where the effective phase-only action S[ϕ] is given as (20) and the effective chemical potential μ n attached to the total density operator is introduced as μ n = Un 2 −μ. The decoupling of the quadratic term proportional toñ 2 (rτ ) in the exponential of the partition function in Eq. (6) is also straightforward. We obtain By combining the expression in the exponential in Eq.
The saddle-point evaluation gives for whereñ = ñ(rτ ) . As the result of decoupling, we obtain following "Zeeman-like" contribution in the partition function with attached effective chemical potential μñ = c − f 2 − Uñ 2 . To summarize, the partition function of the system after the decoupling procedures will be where the action S[c, c,f , f, ϕ] in the exponential is The action in the form given in Eq. (26) is suitable for derivation of the effective phase action and the fermionic action.

The U(1) Gauge Transformation
In the perspective of treating the local and non-local correlations in the excitonic system it is important to separate the U(1) gauge degrees of freedom related to the phase sector. To this end, we perform the local gauge transformation to new fermionic Grassmann variables a(rτ ) and b(rτ ). This procedure will automatically eliminates also the last imaginary term appearing in the expression of the phase action in Eq. (20). For the electrons of c and f -orbitals the U(1) transformation is whereÛ(ϕ) is the U(1) transformation matrixÛ(ϕ) =Î · cos ϕ(rτ ) + iσ z · sin ϕ(rτ ) with the unit matrixÎ andσ z being the Pauli matrix. The variablesx = a, b. We used the bosonic phase variables ϕ (rτ ) introduced in Eq. (15). In fact, the electron factorization in terms of two variables has an unprecedented impact on the whole theory. Especially, the emergent bosonic gauge sector, related to the phase variables, leads to a Bose-type of band bandwidth-renormalization factor (see in the Sect. 4). The action in the Eq. (26) after transformation procedure takes the following form with the new phase action S 0 [ϕ] Then the partition function of the system in new variables is This form of the partition function will be the starting point for deriving the effective actions for the fermions and for the phase sector.

Excitonic Gap
The EI low-temperature phase is characterized by local excitonic order parameter (excitonic gap). The nonvanishing of the expectation value signals the appearance of the electron-hole bound pairs, which manifests as a gap in the excitation spectrum and signals the presence of the EI state. The EI state develops from the local on-site electron-hole correlations.
We start with derivation of the EI gap equation. First, we apply the tranformation given in Eq. (27) to fermionic variables in the initial Hamiltonian of the system in Eq.
(1). Then, we decouple four fermionic interaction term within the HF approach [31] by applying Bogoliubov MF approximation. We have Here n a (rτ ) and n b (rτ ) are the electron densities after the U(1) gauge transformation.
The Fourier transformation of fermionic variables a(rτ ) and b(rτ ) is given by with x = a, b for the a and b type electrons. N is the number of lattice sites and ν n = π(2n + 1)/β are the Fermi-Matsubara frequencies with n = 0, ±1, ±2, . . .. Furthermore, we will integrate out the phase variables in the expression of the partition function given in Eq. (30), we obtain where the effective phase-averaged fermionic action in the exponential is given by Now, using Eq. (33) we can write the action where the effective chemical potentials μ a eff and μ b eff have been introduced as The factors n a and n b in Eqs. (37) and (38) are the average fermion densities n x = n x (rτ ) . Next, t k andt k are band-renormalized hopping amplitudes t k = 2tg B (k) andt k = 2t g B (k), where g B is the bandwidth renormalization factor The explicite expression of this important factor will be given in the Sect. 5, within the quantum rotor representation. (k) is the 3D lattice dispersion relation with d α (α = x, y, z), being the components of lattice spacing vector d = r − r with r and r n.n. positions For the simple cubic geometry they are all equal: d α ≡ a. Employing the vector-space notations, we can rewrite the action in Eq. (36) in more compact form withĜ −1 k (ν n ) inverse of the Green function matrix where the single-particle quasienergies E a k (ν n ) and E b k (ν n ) are given after Eq. (36) The general form of the normal fermionic propagator Gxx (rτ, r τ ), defined in terms of the transformed fermionic variablesx = a, b is and the anomalous or, the excitonic propagator, is given by The averages in Eqs. (45) and (46) are defined with the help of the effective fermionic action in Eq. (41) As a consequence, using Eqs. (41) and (47) we have A similar expression for G bb (rτ, r τ ) could be obtained with the simple replacement Furthermore, for the anomalous propagator we obtain while G ba (rτ, r τ ) is obtained by the substitution¯ → .

Self-consistent Solution for , g and c
Using the local expressions of the Green functions in Eqs. (48) and (49) obtained above, we have the equations for average electron densities n a and n b corresponding to the a and b-orbitals respectively, and also a self-consistent equation for the excitonic order parameter . We have Then, summing over the fermionic Matsubara frequencies, we can rewrite an equivalent system of equations It is worth to mention that the only difference between the obtained MF-like equations Eqs. (53)(54)(55) and the usual HF theory results given in Ref. [25] lies in the presence of the bandwidth renormalization factor g B attached to the c and f -band's hopping amplitudes t andt. In the low-temperature limit this factor goes to 1 and for T = 0 g B = 1 for all values of the Coulomb interaction parameter (see also the discussion at the end of the Sect. 5.2). Here, we assumed the half-filled band case n = n a + n b = 1 and we defined the fermion density differenceñ = n a − n b . Without any loss of generality, we have supposed the case of the EI state with the uniform real gap parameter . Furthermore, n F denotes the Fermi-Dirac distribution function n F ( ) = 1/ e β + 1 . Next, we have the band-energy parameters E + k and E − k defined as with the quasiparticle dispersion ξ k The defines the charge-transfer gap, which we will discuss later on in this Section.

Numerical Results and Discussion
The quantities n a , n b , , μ and c can be determined by solving numerically of Eqs. (53)-(55) in a self-consistent way. We start with the discussion of the stability region for the EI phase on the T − U plane, when approaching EI gap to zero: → 0. The temperature T of the excitonic pair formation, the functionñ and the chemical potential μ are considered here. The summations over the wave vectors in Eqs. (53)(54)(55) can be simplified by introducing the appropriate density of states (DOS) ρ 3D (x) for the 3D lattice. Using Eq. (40) we have For the simple cubic lattice the density of states is given as where (x) is the Heaviside step function and k(x) is the elliptic function of the first kind [45]. In Fig. 1 we have presented the solution for the EI stability region in 3D EFKM by solving the equation (T, U ) = 0, which determines the temperature T for which the pairing gap vanishes. The lowest curve in the Fig. 1 corresponds to the case of the vanishing narrow-band hoppingt = 0 [46]. In this case the critical temperature T still finite. Above this temperature, i.e. when T T we are in the normal Band-Insluator (BI) regime, and = 0. Just below the temperature T , i.e. when T T , the pair formation began and the system is passing into the EI regime. Our calculations, regarding the temperature T of the pair formation, agree very well with the analogous results in previous works (see Refs. [23][24][25][26][27][28]). For the completeness, the density difference between the conduction band and valence band, and the solutions of chemical potential at the EI transition boundary are presented in Figs. 2 and 3.
In Fig. 4 the solution for the excitonic pairing gap is plotted as a function of U/t for different values of the f -band hopping amplitudet and for T = 0. The excitonic gap is non-zero for a rather large domain of the Coulomb interaction in agreement with  [23] and in contrast with the results for the 2D square lattice in Ref. [25]. The obtained values for the lower and upper bounds of the Coulomb interaction in Ref. [25] are about (U c1 , U c2 ) = (0.66, 6.95) and, as it could be expected, they differ considerably from our results, especially for the large hopping. The solution forñ is plotted in Fig. 5 as a function of the dimensionless Coulomb interaction parameter U/t. It is clear in Fig. 5 that in the strong coupling limit U/t 1 the system is in the BI regime, because at the upper bound of the Coulomb interaction the f -band is fully occupied (n b = 1) and the c-band is totally empty (n a = 0). In the inset in Fig. 5 the plot of the function ρ 3D (x) is presented.
The exact numerical solutions for the chemical potential at T = 0 in the intermediate and strong interaction limits (for example 1.8 ≤ U/t ≤ 12 fort = −0.4t) form  Fig. 6) for all values oft, and a single particle excitation gap g = μ max − μ min is opening, where μ max and μ min are the upper and lower bounds of the chemical potential. The evolution of the upper bound of the chemical potential, as a function of Coulomb interaction parameter U/t, is presented in Fig. 7.
By moving from weak into intermediate coupling regime, the single-particle gap g and the pairing gap parameter , both are increasing, while in the strong coupling limit (U/t > 8 fort = −0.3t as an example) decreases rapidly with increasing U/t while g remains open (the Hartree-like gap structure). In the case of vanishing of the pairing gap = 0, the single particle gap g collapses g → 0 and the solution for the chemical potential is single valued (see Fig. 3) (this case corresponds to the case of the boundary of the EI state and is discussed above in Figs. 1, 2 and 3). In other words, we can conclude that in case of intermediate and strong Coulomb interaction parameter, the pairing interaction (when = 0) removes in some sense the degeneracy related to the chemical potential μ. Indeed, the difference between Figs. 3 and 6 is due to the pairing interaction . The charge-transfer gap c defined in Eq. (58) is calculated as a function of the Coulomb interaction parameter U/t. The results are presented in Fig. 8a, b. We see in Fig. 8a that for the small values of the Coulomb interaction, the charge-transfer gap is nearly zero. The small value of it is the manifestation of the semimetallic limit or the BCS limit. By augmenting the interaction parameter U , the gap c is gradually opening. In Fig. 8b we presented the charge-transfer gap for a smaller value of the hopping amplitudet = −0.1t. With decreasing the hopping amplitude we are decreasing also the charge-transfer gap. This is consistent with the results for the excitonic gap parameter presented in Fig. 4 and with the behavior of the single particle excitation gap ( g ) given in Fig. 6.

Effective Phase Action
We are interested now in purely phase action and thus, we will integrate out the fermions in Eq. (30) to obtain the effective phase action in the model. The partition function of the phase-only model is where and the action S 0 [ϕ] is given in Eq. (29). The detailed calculation of the average of second order term in Eq. (62) is given in the "Appendix 1". As result, we have for the phase-only action where, for the parameter J , we have the expresison (see "Appendix 1") with the parameters 1 (x, y) and 2 (x, y) Here E + (x) and E − (x) are the contineous versions (we have just replaced here the index k in Eq. (56 by the contineous variables x) of the similar parameters given in Eq. (56). From the expression of the parameter J in Eq. (64) it follows that the non-zero value of this quantity is directly linked with the pairing gap since J ( = 0) = 0. In Fig. 9 we presented the parameter J as a function of Coulomb interaction parameter U/t and for two different values of the f -band hoppingt. As figures show, the values of J are strictly positive for all regions of the normalized Coulomb interaction parameter. Indeed, the parameter J , in units of the hopping parameter t is very small J/t 1, but is persistent in the whole interaction region with non-vanishing values of the pairing gap . It is also important to emphasize on the form of the phase-stiffness parameter J in Eq. (64). Especially, it follows from Eq. (64) that the macroscopic phase coherence in the system is characterized by an energy scale J ex ∼ ( t e t h )/(t e + t h ) for all values of the Coulomb interaction parameter U , which is related to the motion of the center of mass of e-h composed particle, because (t e t h )/(t e + t h ) ≈ (m e + m h ) −1 . For the strong interaction limit we are converging with the hard core Boson model, with the kinetic energy proportional to t e t h /U ( being the local excitonic order parameter). Thereby, we have shown that non-local correlations between the electrons and holes of different n.n. excitonic pairs, are related with the excitonic BEC condensation.
In the discussion above, we have derived the effective phase-only action S eff [ϕ] = S 0 [ϕ] + S J [ϕ]. In the following, we cast the S eff [ϕ] into the quantum rotor representation [43]. To proceed, we replace the phase degrees of freedom with complex unimodular field z(rτ ) = e iϕ(rτ ) which satisfies the periodic boundary condition z(rβ) = z(r0). The spherical constraint, imposed on a set of the unimodular variables z(rτ ) is The spherical constraint in the Eqs. (67) and (68) can be resolved by introducing the Lagrange multiplier λ resulting from the Laplace transform of functional delta representation [43] δ r |z(rτ ) This adds a quadratic term (in the z-field) to the phase action. Next, the phase action in Eq. (62) can be rewritten in more convenient form using the trigonometric half-angle transformation formula Then, in terms of complex variables z(rτ ), the transformation in Eq. (70) leads to a biquadratic term in the phase action in Eq. (62). We have We can rewrite now the partition function in the form Furthermore, we linearize the action in Eq. (71) (for details see in "Appendix 2") and after we integrate out the phase variables in Eq. (60). Then, after the Fourier transformation of z-variables z(rτ ) = 1 β N k,ω n z(kω n )e i(kr−ω n τ ) with ω n , being the Bose-Matsubara frequencies ω n = 2π n β with (n = 0, ±1, ±2, . . .), the partition function assumes the form with the action S λ [z, z] where Furthermore, g b stands for the bandwidth-renormalization factor, γ −1 (ω n ) in Eq.
(75) is the inverse of the Fourier transformed two-point phase correlation function γ (rτ, r τ ) where Z 0 is the statistical sum of the noninteracting set of quantum rotators The calculation of the Fourier transformation γ (ω n ) of the function in Eq. (76) is straightforward [43] γ where The summations in Eqs. (78) and (79) are over the winding numbers m of the U(1) group (see the Sect. 3.1).

Exciton Condensate at T ∼ T c
In the thermodynamic limit N → ∞ the integration over λ-field in Eq. (72) can be performed exactly using the saddle-point method As a result, one can write the constraint for the saddle-point value of the Lagrange multiplier λ 0 where the average in Eq. (81) is defined as Then, with the help of the Eq. (75) we can write After Eq. (76), the equation Eq. (83) takes now the following explicit form The explicite value of the parameter λ 0 could be determined with the help of the Thouless criterion [48]. It states that the uniform static order parameter susceptibility diverges at the phase transition. Thus G −1 z (k = 0, ω n = 0) = 0 from which we can derive the critical value of the Lagrange multiplier Furthermore, we find After performing the Bose-Matsubara frequency summations in Eq. (83), we obtain the equation for the excitonic BEC transition critical temperature T c where n B ( ) is the Bose-Einstein distribution function n B ( ) = 1/ e β − 1 and the variables ζ 1k and ζ 2k are given by where α = 1, 2. We see also, that at the fundamental state with k = 0 there is a residual gap W = ζ 1 (0) − ζ 2 (0) = −2μ related to the condensate, which equals the binding energy of a molecule in the BEC limit E bind ≈ |2μ| [49][50][51]. At the end of this section we give also the analytical expression of the bandwidthrenormalization factor g B r − r In fact, the calculation of the factor g B r − r could be done alternatively within the self-consistent-harmonic-approximation (SCHA) [52,53]. In this approximation the quantum rotor description is reduced to classical Hamiltonian one. We do not present here the SCHA results for g B r − r . This could be a subject of a future investigation. The calculation of the factor g B r − r shows that, at T = 0, it is equal identically to 1. The numerical solution of the Eq. (87) is presented in Fig. 10. We see that the excitonic BEC transition critical temperature T c is much smaller than the EP formation critical temperature T discussed in the Sect. 4. This conclusion is in the good agreement with the previous theoretical investigations [32][33][34][35].

BEC Transition Amplitude at T T c
In general case, the local constraint in Eq. (67) for bosonic unimodular variables z(rτ ) breaks down at very low temperatures, (especially at T = 0) because we have to consider the symmetry breaking related to the Bosonic sector, thus, critically, we have the fluctuation form z(rτ ) = e iϕ (rτ ) +z(rτ ) and the unimodularity constraint is broken. In the limite of very low temperatures, considering the BEC of excitons, we have the spontaneous breaking of local U(1) gauge-symmetry related to the phase field, leading to the nonvanishing expectation value of z(rτ ). In order to demonstrate this, we separate the single particle state k = 0 by using Bogoliubov displacement operation (see, for details in Refs. [1] and [54]). Then, we write for the complex variables z(k, ω n ) where ψ 0 is the condensate transition amplitude ψ 0 = z(k, ω n ) of the bosonic field. Nextz(k, ω n ) is the excitation part of effective Bose-field. The general form of the bosonic charge propagator is given by where The average in Eq. (92) is defined in Eq. (82). We consider the expectation value z(k, ω n )z(k, ω n ) and we draw the condensate part by applying the transformation in Eq. (90). Hence, we have HereG z (k, ω n ) is related to the on-condensate exctitation part of the bosonic sector Finally, the local constraint in Eq. (67) will be rewritten as then, we get The obtained values for |ψ 0 | 2 are plotted in Fig. 11. With increasing the temperature, BEC transition probability decreases and disappears for the high temperature limit.

Momentum Distribution Functions and Exciton Coherence Length
To proceed we define frequency-summed normal and anomalous momentum dependent functions n a (k) = 1 β ν n G aa (kν n ), where G aa (kν n ) and G ab (kν n ) are the Fourier transformations of the local, normal and anomalous, propagators. Using Eqs. (48) and (49) we obtain Summing over the fermionic Matsubara frequencies ν n , we get while the function n b (k) for the b-orbital is simply n b (k) = 1−n a (k). The Bogoliubov coefficients appearing in Eqs. (100) and (101) are given by The plots of the normal and anomalous functions n a (k) and F(k) are given in Figs. 11, 12c. The k -summations in the analytical expression of ξ c were done with the (100 × 100 × 100) k-points in the FBZ. In the weak coupling regime the normal distribution function n a (k) drops at k F (see the top plots in Figs. 11, 12c) and anomalous momentum function (the bottom plots in Figs. 11, 12c) is picked at the Fermi level. With increasing the Coulomb interaction n a (k) spread out in the k-space and also k F becomes broad with the Fermi level k F displaced to the value (0, 0, 0) in the momentum space. Across the crossover regime, the anomalous momentum function decreases for all momenta of the reciprocal space and this is consistent with the behavior of the excitonic gap parameter in the strong coupling regime presented in Fig. 2. Subsequently, in Fig. 12b, we have presented the temperature dependence of the anomalous distribution function.
The spatial coherence of a fermionic system is encoded in its one-body density matrix, therefore, the anomalous momentum function is directly related to the excitonic coherence length. We can associate a characteristic decay of F(k) with the coherence length ξ c , defined by the relation [25] The quantity ξ c provides the quantitative information about the properties of the system. By calculating the coherence length given by Eq. (105) for different values of the Coulomb interaction parameter U/t, we can see directly the spatial extension of a single exciton. The results are given in Fig. 12c, d, where a rapid growth of the coherence length, for the small values of the Coulomb interaction parameter, is anticipated with the excitons cooled down below the temperature of their quantum degeneracy and the system is in the macroscopic phase-coherent regime. On the other hand, opposite to this behavior, the coherence length decreases rapidly with increasing U/t. A very similar decrease of the coherence length of the Cooper pairs is also proved in exact-diagonalization study on the attractive Hubbard model [55]. The coherence length has a minimum in the intermediate coupling regime and this is due to the denominator k |F(k)| 2 in Eq. (105), which is largest in this case. In the strong interaction region, ξ c slightly increases with increasing U/t. Our results are in good agreement with the HF results discussed earlier [16,25,28] ( Fig. 13).

Summary and Outlook
Now it is interesting to relate results of our calculations on the 3D excitonic system to the experimental results, e.g., for the compound TmSe 0.45 Te 0.55 , which is an intermediate valent semiconductor [4][5][6][7]. The hopping parametert is estimated for |t| = 0.3|t| = 5 meV (see Ref. [28]). By using these values, we find for the maximum of the excitonic pair transition temperature T max = 186.6 K at U = 8|t|, while the maximum of the exciton BEC transition temperature is found to be smaller of about two orders of magnitude at T max c = 0.44 K for U = 4.8|t|. Furthermore, the charge-gap bandwidth was found as to be W = | min c | = 0.0682 eV and the single particle excitation gap is of order g = 0.057 eV at U = 10.6|t|. The In conclusion, we have studied the excitonic phase transition in a system of conduction band electrons with transfer parameter t, and the valence band holes, described by 3D extended Falicov-Kimball model with the tunable Coulomb interaction U between both species. To this end we implement the functional integral formulation of our model, where the Coulomb interaction term is expressed in terms of U(1) quantum phase variables ϕ conjugated to the local particle number, providing a useful representation of strongly correlated systems. At low temperatures, the electron-hole system may become unstable with respect to the formation of the excitons at T = T , exhibiting a gap in the particle excitation spectrum controlled by the parameter U/t, which gives the relevant energy scale for the excitonic insulator state. In the weak coupling limit, U/t 1, the binding energy of excitonic pairs is small, thus pair breaking effect controls the excitonic phase transition in analogy to that what happens in a standard BCS superconductor. In the excitonic system with the strong pairing U/t 1, we have the situation, where the pairs are strongly bound and localized which diminish, T for large U/t.
We have shown that the excitonic BEC transition temperature T c is much smaller than the critical temperature T of the excitonic pair formation in good agreement with the discussions in Refs. [32][33][34]. Our results are in good agreement with the previous theoretical and experimental results. A possible direction for future work will be the determination of the single-particle excitation spectra and the excitonic density of states, which would be instrumental for interpretation of the coherent light emission measurements in the excitonic system.
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Fermionic Action
We would like now to derive the effective action for fermions. Our starting point is the partition function given in the Eq. (29) derived with the help of the U(1) gauge transformation as it is introduced in the Sect. 3 and which can be written as where the effective fermionic action S eff [ā, a,b, b] in the exponential is defined as Furthermore, we expand the logarithm keeping only the terms up to second order in S. As a result we obtain Here, the averages with respect to the phase variables are defined as

Phase Action
In a similar way, the integration over the fermions in Eq. (29) gives the effective action for the phase sector. The partition function in this case is where the effective phase action in the exponential is Again, by expanding the logarithm in the Eq. (110), we will have up to second order in S whereS 0 is an unimportant constant. Here the fermionic average · · · S eff is given by The Eq. (108) is important for deriving the excitonic phase-stiffness parameter. We present here derivation of the terms in the effective phase action, which are proportional to tt, in Eq. (111). The derivation of the other term, proportional tott, is very similar. Using Eqs. (111) and (27) we get for the mixed term tt − 1 2 As an example, we give the Wick averaging result of the first four-fermion term in the expression of Eq. (114) We kept in Eq. (114) only the terms proportional to excitonic gap. We neglected other terms like ā(rτ )b(r τ ) or a(rτ )b(r τ ) , which vanish due to the symmetry of the action in Eq. (63). Contributions, proportional to fermionic densities ā(rτ )a(r τ ) and b (rτ )b(r τ ) could be also omitted, since they are not contributing directly to the excitonic pair formation. After calculating all averages in Eq. (114) and recombining them with the similar terms coming from the component proportional tott, we obtain the relevant portion of the phase action in the form which contains the phase-stiffness parameter J Furthermore, in order to simplify the non-local (in time variables) effective phase action in Eq. (116), we resort to the gradient expansion of the phase field in the form As a result we can deduce the phase-stiffness parameter J in the form while the phase action in Eq. (116) simplifies to that given in Eq. (62), which now is local in time variable τ . For the product of the anomalous propagators in Eq. (119) we have

.(120)
Here z = 6 is the number of the n.n. sites on the 3D cubic lattice. After integrating over the imaginary time τ in Eq. (119), we perform the Matsubara frequency summations in Eq. (120) and obtain the phase-stiffness parameter J in Eq. (119) in the final form The parameters 1 (k, k ) and 2 (k, k ) entering in Eq. (121) are given by The summations over the k wave vectors in Eq. (121) could be transformed into the integrations with the help of the density of states given in Eq. (60) (see the Sect. 4.2). As we see, Eq. (121) relates the parameter J with the local pairing gap . The numerical evaluations of the expression in Eq. (121) for T = 0 K are presented in Fig. 9 and discussed in the Sect. 5 of the present paper.

Appendix 2: The Action S λ [z, z]
The action in Eq. (71) is quartic in unimodular z-field and could be decoupled with the help of the MF-like decoupling procedure Then we get We can write also an analogue expression for the exponential e i β 0 dτ r η z−e −iϕ(rτ ) . Thereby, we have where we introduced the phase-correlation function γ rτ, r τ = e i[ϕ(rτ )−ϕ(r τ )] . Now, we integrate out the bosonic η-field by employing the HS complex transformation for bosons where G −1 z (rτ, r τ ) is the inverse of the real-space bosonic Green-function matrix.