Density of states in the bilayer graphene with the excitonic pairing interaction

In the present paper, we consider the excitonic effects on the single particle normal density of states (DOS) in the bilayer graphene (BLG). The local interlayer Coulomb interaction is considered between the particles on the non-equivalent sublattice sites in different layers of the BLG. We show the presence of the excitonic shift of the neutrality point, even for the noninteracting layers. Furthermore, for the interacting layers, a very large asymmetry in the DOS structure is shown between the particle and hole channels. At the large values of the interlayer hopping amplitude, a large number of DOS at the Dirac's point indicates the existence of the strong excitonic coherence effects between the layers in the BLG and the enhancement of the excitonic condensation. We have found different competing orders in the interacting BLG. Particularly, a phase transition from the hybridized excitonic insulator phase to the coherent condensate state is shown at the small values of the local interlayer Coulomb interaction.


INTRODUCTION
The problem of excitonic pair formation in the graphene and bilayer graphene (BLG) systems is one of the longstanding and controversial problems in the modern solid state physics [1][2][3][4][5][6][7][8][9][10][11][12]. The chiral invariance of the free quasiparticle Hamiltonian, combined with electron Coulomb interaction term brought the idea about the possibility of spontaneous chiral symmetry breaking (CSB) in a single layer and bilayer graphene, reflecting in the form of the gapped states in the fermionic quasiparticle spectrum [1,[13][14][15][16][17][18]. The excitonic gap equation has been derived in the Refs. [5][6][7][8][9][10][11] at the Bardeen-Cooper-Schrieffer (BCS) limit of the excitonic transition scenario. The general idea used there is based on the supposition of the weak electronic correlations in the BLG, due to its large number of the fermionic flavors. It has been suggested in Refs. [1,16,17] that even a partial account of the wave vector or energy dependence of the gap function rules out the constant gap solution in the BCS limit. Meanwhile, it has been shown that even undoped graphene can provide a variety of electron-hole type pairing CSB orders especially for the strong Coulomb coupling case [16,17], which renders the treatments in [1,[13][14][15]18] to be obscure.
Concerning the excitonic condensation, the proper inclusion of the chemical potential fluctuation effects on the excitonic pairing gap and condensation in the BLG system, and also the role of the average chemical potential on the excitonic effects have been discussed in Refs. [10][11][12]. A very large excitonic gap, found in a sufficiently broad interval of the repulsive interlayer Coulomb interaction parameter, given in Ref. [12], suggests that the level of intralayer density and chemical potential fluctuations in the BLG system, discussed in Ref. [11], are sufficiently small and do not affect the robust excitonic insulator state in the BLG. The excitonic condensation in the single BLG is principally possible in the case of the static interlayer screening regime [10].
Moreover, it has been shown that there exists a critical value of the interlayer hopping amplitude γ 1 [15], which provides an interesting energy cutoff, below which the electron-hole correlations do not drive the system towards the CSB excitonic transition. More recently, it has been shown [19,20] that the single-particle coherent density of states (DOS) in the usual two-dimensional (2D) [19] and three-dimensional (3D) [20] semiconducting systems at the zero temperature limit is always finite, reflecting with the excitonic condensate regime in these systems. For the 3D semiconducting systems, the coherent DOS spectra survive also for the higher temperatures [20]. This situation is also typical for the double layer electronic structures at the half filling [21] when considering the electron-hole pair formation and condensation. In addition, it has been demonstrated that the excitonic insulator state and the excitonic condensation are two distinct phase transitions in the solid state [22,23], and the condensate states are due to the electronic phase stiffness [20,24], mechanism. Meanwhile, in the electronic bilayer systems, those mentioned phase transitions are indistinguishable as it was shown in Refs. [12,21]. On the other hand, the inclusion of the dynamic screening and the full band structure favors the pairing and condensation [25][26][27].
In this paper, we use the bilayer Hubbard model to calculate the single-particle DOS functions in the BLG and we examine the excitonic effects in the DOS. For the noninteracting layers, we show the principal modifications to the usual tight-binding single layer graphene's DOS behavior (with differences between the A and B DOS structures near the Dirac's neutrality point) and we estimate the excitonic blue-shift values in the normal DOS for the reasonable values of the interlayer hopping amplitude.
At finite values of the interlayer interaction parameter, the DOS shows always the remarkable four peaks structure and a very large interband hybridization gap opens for intermediate values of the interaction parameter. We have found the critical value of the interlayer interaction parameter at which the hybridization gap ∆ Hybr opens in the system. We estimate the values of ∆ Hybr for different strengths of the interlayer interaction parameter and for different values of the interlayer hopping amplitude γ 1 . We show the modifications of the van Hove singularity (vHs.) peaks in the DOS when augmenting the interlayer interaction parameter. An interesting interlayer hopping mediated crossover, from the insulating pairing state to the excitonic condensate state follows from our considerations. Moreover, we show that at any reasonable value of the interlayer interaction parameter W , there exist a critical value of the interlayer hopping parameter γ 1 at which the hybridization gap vanishes and the BLG pass to the coherent excitonic condensate state. We calculate the hybridization gap in the system, and we show that it becomes larger when increasing the interaction parameter W .
Furthermore, we estimate the excitonic shift energies, the intraband vHs. peaks separations and the values of the DOS at the Dirac's neutrality points for the zero temperature case. The full interaction bandwidth, considered here, mimics different limits of correlations in the BLG, which have been discussed only partially in the literature, where the discussions have been restricted only to the strong interlayer screening regime.
The paper is organized as follows: in the Section 2, we introduce the bilayer Hubbard model for our BLG system. In the Section 3 we give the form of the fermionic action in the Feynman's path integral formalism. Next, in the Section 4, we discuss the single-particle DOS in the BLG, for different interlayer interaction regimes. In the Section 5 we present the numerical results on the excitonic effects in the normal DOS and, in the Section 6, we give a short conclusion to our paper.

THE INTERLAYER HUBBARD MODEL
The Bernal Stacked (BS) BLG system is composed of two coupled honeycomb layers with sublattice sites A, B andÃ,B, in the bottom and top layers respectively, arranged in the z-direction in such a way that the atoms on the sitesÃ in the top layer lie just above the atoms on the sites B in the bottom layer graphene, and each layer is composed of two interpenetrating triangular lattices. We consider the electronic BLG structure with the equal chemical potentials and equal on-site quasienergies in each layer. When switching the local Coulomb potential between the layers, we keep the charge neutrality equilibrium across the BLG, by imposing the half-filling condition in each layer. The interlayer Hubbard model with the intralayer U and local interlayer W Coulomb interaction terms is subjected by the following Hamiltonian The first two terms describe the intralayer electron hopping with the hopping parameter γ 0 . The summation rr ′ , in these terms denotes the sum over the nearest neighbours lattice sites in the separated honeycomb layers in the BLG structure. We kept here the small letters a, b andã,b for the electron operators on the lattice sites A, B andÃ,B respectively. The third term corresponds to the hopping between two different monolayers in the BLG and the parameter γ 1 is the interlayer hopping amplitude. We neglect the nonlocal interlayer hopping terms because their role is essentially unimportant in the considered problem about the excitonic effects. When working with the grand canonical ensemble, we have added also the chemical potential terms for each layer and for each sublattice. We suppose here that the chemical potentials µ ℓ corresponding to different layers ℓ = 1, 2 and different sublattices A, B andÃ,B are the same (this is, of course, true, if we consider purely electronic layers and we neglect the effect of the disorder, caused by the additionally charged impurities). It is important to mention here that the chemical potentials of electrons on the nonequivalent sublattice sites, in the same layer, get different shifts in different layers due to the stacking order of the BLG structure. Next, n ℓσ (r) is the electron density operator for the fermions in the layer ℓ and with the spin σ. The four U -terms in the Hamiltonian in Eq.(1) describe the on-site intralayer Coulomb interactions in the BLG. The summation index η in the U -terms, in Eq.(1), refers to different sublattice fermions in a given layer, i.e., η = a, b, for the bottom layer with ℓ = 1, and η =ã,b, for the top layer with ℓ = 2. We will study the excitonic effects in the BLG at the half-filling condition in each layer, i.e., n ℓ = 1, for ℓ = 1, 2. The last term in Eq.(1), describes the local interlayer Coulomb repulsion, and the parameter W is the corresponding interlayer Coulomb interaction. We put γ 0 = 1, as the unit of energy, and we set k B = 1, = 1 through the paper.

THE FERMIONIC ACTION AND DIRAC REPRESENTATION
Here, we write the partition function of our electronic BLG system. For simplifying the notations we will introduce the two component fermionic fields, corresponding to different type of fermions in the two sublattices of the graphene monolayers: f c = (c, c), where c = a, b, for the layer with ℓ = 1 and c =ã,b for the layer with ℓ = 2. Next, the partition function will be written as The fermionic action in the exponential, in Eq.(2), is given by We have introduced here the imaginary-time variables τ [28,29], at each lattice site r in both layers of the BLG. The variables τ vary in the interval (0, β), where β = 1/T with T being the temperature. The first term in Eq.(3) describes the fermionic Berry-terms, corresponding to all fermionic flavors in the BLG system. They are given as Furthermore, in order to examine the excitonic effects in the BLG, we perform the real-space linearization of the four-fermionic terms, in Eq.(1). We do not present here the details of such a procedure and we refer to our recent work, in Ref. [12], where this procedure is presented in more details. For the next, we will consider the homogeneous BLG structure, when the pairing occurs between the particles with the same spin orientations, i.e., ∆ σσ ′ = ∆ σσ δ σσ ′ and we can assume that the pairing gap is real ∆ σ =∆ σ ≡ ∆. Then the excitonic pairing gap parameter is given by Next, we pass to the Fourier space representation, given by the transformation c σ (r, τ ) = 1 βN kΩn c σk (Ω n )e i(kr−Ωnτ ) , where N is the total number of sites on the η-type sublattice, in the layer ℓ, and we write the partition function of the system in the form Here, Ω n = π (2n + 1) /β with n = 0, ±1, ±2, . . . , are the fermionic Matsubara frequencies [29]. The four component Dirac spinors ψ σk (Ω n ), in Eq. (6), are introduced at each discrete state k in the reciprocal space and for each spin direction σ =↑, ↓. Being the generalized Weyl spinors [30,31], they are given as The matrix G −1 σk (Ω n ), in Eq. (6), is the inverse Green's function matrix, of size 4 × 4. It is defined as Indeed, the structure of the matrix does not changes when inverting the spin direction, i.e., The diagonal elements of the matrix, in Eq. (8), are the quasienergies where the effectve chemical potentials µ ℓ with ℓ = 1, 2 are defined as µ eff 1 = µ + U/4 and µ eff 2 = µ + U/4 + W . The parametersγ ℓk , in Eq. (8), are the renormalized (nearest neighbors) intralayer hopping amplitudes:γ ℓk = zγ ℓk γ 0 , where the k-dependent parameters γ ℓk are the energy dispersions in the BLG layers, Namely, we have The parameter z, is the number of the nearest neighbors lattice sites in the honeycomb lattice. The components of the nearest-neighbors vectors δ ℓ , for the bottom layer 1, are given by δ . For the layer 2, we have obviously δ Then, for the function γ 1k , we get where d is the carbon-carbon interatomic distance. It is not difficult to realize that γ 2k = γ * 1k ≡ γ * k , and we havẽ γ 2k =γ * 1k ≡γ * k , where we have omitted the layer index ℓ.
The form of the Green's function matrix, given in Eq. (8), has been used recently in the work of Ref. [12] in order to derive the self-consistent equations, which determine the excitonic gap parameter ∆ and the effective bare chemical potentialμ in the interacting BLG system. Particularly, this last one plays an important role in the BLG theory and redefines the charge neutrality point (CNP) in the context of the exciton formation in the interacting BLG, [12]. Quite interesting experimental results on that subject are given recently in Refs. [32][33][34].

THE SINGLE-PARTICLE DOS
A. The sublattice spectral functions In this section, we calculate the single-particle normal DOS for the BLG system, given by the Hamiltonian, in Eq.(1) and basing on the form of the fermionic action, given in Eq. (6). The normal single-particle DOS is straightforwardly defined as where the index c corresponds to the sublattice type in the BLG layers, i.e., c = A, B,Ã,B, and the spectral function S c (k, ω), in the right hand side in Eq. (12), is defined with the help of the retarded real-time Green's function Here, we have supposed that the spin variable σ is fixed in the direction spin-↑, being completely unimportant for the considered here problem. In turn, the retarded Green function G R c (k, ω), could be obtained after the analytical continuation into the real frequency axis, in the expression of the respective Matsubara Green's function G c (k, Ω n ), (see the similar procedures in Refs. [28,29]) The explicit calculation of the thermal normal Green's functions G c (k, Ω n ) follows from the definition of the normal Matsubara Green's function [29]. As usually, in the real space, they are defined as the statistical average of the product of an annihilation c and a creation typec operators, i.e., for our fermions, we have After transforming into the Fourier space the fermionic operators, entering in Eq.(15) the local expression of the Green's function (for the symmetry reasons of the action, in Eq. (6), we consider the time-and space-local expression of the single particle Green's function), in Eq. (15), is where the Fourier transforms G c (k, Ω n ) are given by In order to calculate the statistical average on the righthand side in Eq.(17), we will perform the Hubbard-Stratanovich transformation in the expression of the partition function. In the Dirac's spinor notations, the partition function in Eq.(2), in the Section 3, will be transformed as where the auxiliary source field vectors J kσ (Ω n ) are the subjects of the Dirac's spinors defined similar to the ψfields, in the Section 3 The matrixĜ σk (Ω n ) is defined as the inverse of the matrix given in Eq.(8) and we haveĜ σk . Furthermore, the calculation of the thermal Matsubara Green's function, in Eq. (17), is straightforward. For the sublattice-A, in the bottom layer of the BLG, we get Then it follows that After the inversion of the matrix given in Eq. (8), the Green's function matrix component G 11 k (Ω n ) takes the following form where the dimensionless coefficients α (1) ik are and P (3) (ε ik ) is a polynomial of third order in ε ik . Namely, we have with the coefficients ω ik , i = 1, ...3, given as and The quasiparticle energy parameters ε ik , in Eq.(23) are defined as follows where we have introduced a new bare chemical potential As we will show later on, the bare chemical potential has a fundamental impact on the DOS behavior in the BLG, and provide the excitonic shift on the frequency axis. The energy parameters in Eqs. (28) and (29) define the electronic band structure in the BLG system with the excitonic pairing interaction (see also in Ref. ([12])). It is not diffucult to verify that for the noninteracting BLG, i.e., when U = 0, W = 0 and ∆ = 0, the expressions in Eqs. (28) and Eq. (29) are reducing to the usual tight binding dispersion relations [35] in the context of the real-space Green's function study of the noninteractiong BLG. It is important to mention that the normal spectral functions in different layers of the BLG, coincide with each other when interchanging the sublattice notations in the monolayers. This follows from the symmetry of the action, given in Eq. (6). Particularly, we obtain S ℓ=2,c=ã (k, ω) = S ℓ=1,c=b (k, ω), (31) S ℓ=2,c=b (k, ω) = S ℓ=1,c=a (k, ω).
In Fig. 1, we have presented the variation of the effective chemical potentialμ/γ 0 as a function of the interlayer Coulomb interaction parameter W/γ 0 . The temperature dependence has been also shown. In the inset, in Fig. 1, we have shown the dependence ofμ/γ 0 on the local intralayer Coulomb interaction parameter U , given bȳ where µ(W, U, γ 0 , T ) is the exact solution of the chemical potential in the BLG. Different fixed values of the interlayer interaction parameters W are considered in the picture. The more detailed discussion on the role of the chemical potential is given in Ref. [12]. We see that the behavior of the effective chemical potential as a function of W/γ 0 shows a finite, very large jump at T /γ 0 = 0 and the shifted CNP corresponds well with the recent experimental observations of the behavior of bilayer's average chemical potential, given in Ref. [32], where a direct measurement of the chemical potential of BLG has been done, as a function of its carrier density n. Here, it is important to mention that the excitonic effects in the BLG lead to the significant shift of the double CNP in the BLG (see in Refs. [32][33][34]). Contrary, the intensity of the functionμ/γ 0 is much higher in our case, due to the single BLG considered here, while in the gated double BLG measurements, discussed in Ref. [32], the interlayer interaction is much weaker as compared with the double monolayer graphene (see about in Ref. [10]), and the reason for this is the effect of the finite amount of carrier density n T induced in the top BLG Ref. [32]. In Fig. 2, the solution for the chemical potential is shown as a function of the intralayer Coulomb interaction parameter U and for different values of the interlayer interaction parameter W . The linear slope of the chemical potential corresponds to the coefficient κ = 0.25, given in Eq. (33). In Fig. 3, the exact solution for µ is presented for different values of the intralayer interaction parameter U . The zero temperature case is considered in the picture and the interlayer hopping amplitude is fixed at the value γ 1 = 0.128γ 0 .

B. The sublattice density of states
With the help of the analytical continuation, given in Eq. (14) and by the use of the formula 1/(x + iδ) = P(1/x) − iπδ(x), we get for the single particle normal DOS for the sublattice A Here, we can transform the summation over the wave vectors, into the integration over the continuous variables, by introducing the 2D density of states corresponding to the noninteracting graphene layers The noninteracting DOS in the monolayer graphene beyond the Dirac's approximation [36,37], can be analytically expressed in terms of the elliptic integral of the first kind K(x) [38]. Namely, we have where, formally, we have enlarged the domain of variation of the argument, into ±∞. The function ϕ(x), in Eq.(36), is given by [36] ϕ After the definition in Eq.(37), we get for the normal A DOS function where the frequency dependent dimensionless parameters ǫ i (ω) with i = 1, . . . , 4, in Eq. (38), are given by the fol-lowing expressions Next, for the functions Λ ∓ (x), in the denominators, in Eq.(38), we have and it is clear that For the considered assumption of the halffilling in each layer, and as the theoretical and numerical calculations show, the normal single-particle DOS is not the same for different sublattices in the given layer with l = 1, 2 and near the shifted neutrality points. For the next, we will omit the layer indexes near the DOS functions notations (due to the relations in Eqs. (31) and (32), in the Section 3). Similarly, for the sublattice B, we have the following expression for the B DOS The coefficients β ik , in Eq.(41) with i = 1, ..4 are given by the following relations where P ′ (3) (ε ik ) in Eq.(42) is again a polynomial of third order in ε ik , namely with the coefficients ω ′ ik , i = 1, ...3, given as and We see that the coefficients ω ′ 1k , ω ′ 2k and ω ′ 3k could be obtained from the coefficients ω 1k , ω 2k and ω 3k , just by replacing the effective chemical potentials µ and by setting simultaneously ∆ = 0. Finally, for the single-particle B DOS we get We will examine numerically calculated DOS functions in the next Section of the present paper.

RESULTS AND DISCUSSION
In the panels a and b, in Fig. 4, we have presented the plots of the A and B DOS functions given in Eqs. (38) and (47). The zero interlayer interaction limit is considered W = 0. In the panel a, in Fig. 4, the plots of the A and B DOS functions are presented for the case of the zero intralayer Coulomb interaction U . It is clear in Fig. 4 that each band of the BLG, given in Eqs. (28) and Eq.(29), contributes with one vHs. peak inherited from the singlelayer graphene spectrum. We see that the DOS functions, given in Eqs. (38) and (47) reproduce correctly the tight binding graphene DOS behavior (see in Ref. [36]) with the difference that the neutrality Dirac's point is shifted toward the higher frequencies, ω 0 = 1.363γ 0 . For the realistic γ 0 = 3eV (see in Ref. [39]), we get for the shifted frequency ω 0 = 4.089 eV, which is sufficiently large as compared with the tight binding graphene's value. In the panel b, in Fig. 4, the same DOS functions are plotted for the case of the finite intralayer interaction parameter U = 2γ 0 . As we see, in this case, the separation between the vHs. peaks in the DOS is larger, and the value of the DOS at the neutrality point ω 0 , for the sublattice A, is higher than in the previous case, given in the panel a. The position of the neutrality point is unchanged, and we see also that the region, where the B DOS is drastically decreasing, is much larger in the case of the nonzero intralayer interaction parameter U = 0. The difference between the A and B DOS structures is clear in the panel b in Fig. 4. In Fig. 5, we have presented the DOS functions for different values of the intralayer Coulomb interaction parameter: U = 0, U = 1γ 0 and U = 2γ 0 . In all that cases the DOS functions, corresponding to different sublattices, are different near the neutrality point. It is important to mention that the DOS functions remain unchanged when U = 0, but they are different from the zero interaction case U = 0. The DOS behaviors presented in Figs. 4 and 5 are very similar to the DOS structures, discussed in Ref. [35], apart the shifted neutrality point. The shift effect of the neutrality point in the DOS structures is due to the strong excitonic effects in the BLG.
In Fig. 6, we have shown the the A and B DOS evolutions for different values of the interlayer hopping amplitude γ 1 . The intralayer Coulomb interaction parameter is fixed at the value U = 2γ 0 , and the zero interlayer Coulomb interaction case is considered in the picture. We observe in Fig. 6 (see in the panels b and c) that the increase of the interlayer hopping amplitude leads to a very large number of A DOS at the neutrality point. The vHs. peaks separations also become very large when increasing the parameter γ 1 . In Fig. 7, we have presented the evolution of the A DOS near the neutrality point for the same values of the interlayer hopping amplitude γ 1 . The A DOS behavior near the point ω 0 shows that the interlayer hopping amplitude could lead to the existence of the interlayer excitonic condensate states even at the zero value of the interlayer Coulomb interaction parameter. The very large A DOS value at the Dirac's point (see the panel c, in Fig. 6, and also in Fig. 7) is caused by the shift of higher situated energy states in the A DOS structure toward the neutrality point ω 0 , and mediated by the formation of coherent condensates states. This scenario of the excitonic condensation at the zero interlayer coupling is converging well with the general discussion about coherent excitonic density of states in the semiconducting systems (see in Ref. [20]), where it has been shown that a large amount of states in the DOS (without the hybridization gap) is the sign of the coherent excitonic condensates in these systems.

A. The hybridization gap in the BLG DOS
Here, we will examine the formation of the hybridization gap in the BLG system caused by the interlayer Coulomb interaction parameter W . For the convenience, we will fix the value of the interlayer hopping amplitude at the value γ 1 = 0.128γ 0 and the intralayer interaction parameter at U = 2γ 0 . In Fig. 8, we have presented the enlarged pictures of the A and B sublattice DOS functions near the neutrality point ω 0 . In Fig. 8, we consider two, very close, values of the interlayer interaction parameter W and we show how the hybridization gap ap-pear above the critical value W = 0.133γ 0 (for a given value of γ 1 ). We see in the upper panel a, in Fig. 8, that the behaviors of different sublattice DOS functions are drastically not the same near the neutrality Dirac's point ω 0 in the DOS. As the numerical calculations show, there is a critical value of the interlayer interaction parameter, above which the hybridization gap starts to open in the BLG (see in the lower panel b in Fig. 8) and the bilayer system passes to the regime of the excitonic pairing in the insulating state. A very small, but finite hybridization gap appears at W = 0.1331γ 0 , which is slightly higher than the critical value W = 0.133γ 0 , considered in the upper panel a, in Fig. 8. In Fig. 9, we have shown the A and B sublattice DOS functions for the large value of the interlayer Coulomb interaction parameter, corresponding to the maximum value of the excitonic pairing gap parameter (see in Ref. [12]). It is clear from the structure of the DOS functions, presented in Fig. 9, that the insulating state of the BLG is symmetric with respect to the hybridization gap formation, i.e., the A and B DOS functions goes to zero at the same values on frequency axes on the both sides of the hybridization gap ∆ Hybr . It is remarkable to note also that unlike the half-filling considered here, the vHs. in the DOS structures, corresponding to different particle channels in the band structure, are not symmetric with respect to the Dirac's point and the hybridization gap. The inter-peak separations in the particle or hole channel in the DOS become strongly asymmetric for the finite values of the interlayer interaction parameter W . The observed "blue" shift effect of the neutrality point and the strong asymmetries in the DOS structures are due to the strong excitonic effects in the BLG. At any finite value of the interlayer interaction parameter, the BLG system is in the excitonic pairing state. On the other hand, the excitonic condensation is impossible for the large values of W (even at large values of the parameter γ 1 ), because of the strongly hybridized states in the particle and hole channels and the very large hybridization gap.
In Fig. 10, we have shown the W dependence of the A and B sublattice DOS functions for the interlayer hopping amplitude γ 1 = 0.128γ 0 and for U = 2γ 0 . Five values of the interaction parameter are indicated in Fig. 10. Just for the interest, we have shown also the DOS functions at the critical value of the interaction parameter W c = 0.133γ 0 = 0.399 eV, above which the hybridization gap opens in the BLG. The principal observation in Fig. 10 is that the hybridization gap is increasing when increasing the interlayer interaction parameter. This fact is in good agreement with the principal results of the excitonic pairing scenario discussed in Ref. [12], where it has been shown that the interlayer coupling interaction favors the excitonic pairing state in the BLG.
On the other hand, in Fig. 11, we have shown how the insulating gap, in the A DOS spectrum, is closing when augmenting the interlayer hopping amplitude (the relatively small interlayer interaction parameter is considered in Fig. 11: W = 0.1331γ 0 ). In the inset, in Fig. 11, we have shown the A DOS spectrum with the hybridization gap of order ∆ Hybr = 0.00299γ 0 = 8.97 meV. The interlayer Coulomb interaction parameter is of order W = 0.1331γ 0 = 399 meV and the interlayer hopping amplitude is γ 1 = 0.128γ 0 = 384 meV. Then, in the picture in Fig. 11, we show the A DOS near the excitonic condensate states even at the non-zero value of the interlayer interaction parameter W . Thus at the large values of the interlayer hopping amplitude, the system BLG is passing from the insulating hybridized state into the possible excitonic condensate state. This improvement analog to this, about the excitonic condensate state in the BLG and mediated by the parameter γ 1 is also discussed in Ref. [12], where it has been shown how the excitonic condensation state is improved for the large interlayer hoppings. Let's mention also that the value of the excitonic shift-frequency ω 0 = 4.089 eV is very close to the absolute numerical value of the effective bare chemical potential solution in the BLG: |μ| = 1.37γ 0 = 4.11 eV, and which has been calculated in Ref. [12]. The important role of the non-zero chemical potential solution on the DOS behavior is discussed also in Ref. [40], concerning the single layer graphene, where a Drude peak arises in the longitudinal conductivity spectrum, and the DOS becomes finite at the Fermi level. We observe also in Fig. 11, that the Dirac's neutrality point ω 0 is shifting toward the lower frequency region (see the evolution from red to green lines in the picture). This red-shift effect and the excitonic shift observed in the previous pictures, presented here, are much more significant than the shift effects discussed in Refs. [36,40], which are due to the inclusion of the next nearest neighbors intra-and interlayer hoppings in the monolayer graphene and BLG, and also differ from the results on the single impurity problem, discussed in Ref. [41]. The increase of the interlayer interaction parameter above its critical value W c , with the appropriate highest value of the interlayer hopping, (see the green line, in Fig. 11 the excitonic absorption spectrum in the BLG both at zero interlayer coupling (zero applied bias) and nonzero coupling cases and should be verified ulteriorly (maybe the citations are needed). Particularly, we expect that in a large domain of the incident photon's energies, the BLG absorption spectrum, at the zero applied voltage, will show a sufficiently large absorption peak region in the case of the large interlayer hopping amplitude γ 1 . In contrast, for the finite interlayer Coulomb interaction, the A DOS has finite values at the neutrality points in the case of the relatively small Coulomb interactions W  [12]) and for the interlayer hopping amplitude fixed at γ1 = 0.128γ0 = 0.384 eV. The zero temperature case is considered in the picture. and relatively high values of the interlayer hopping γ 1 . This blue-shift effect, caused by the interlayer hopping amplitude, means that the excitonic condensate states survive for the higher values of the incident photon's energies, thus improving the excitonic insulator state at the large values interlayer hopping.

CONCLUSION
We have considered the density of states in the BLG system, by considering the bilayer Hubbard model at the half-filling condition in each layer, and by assuming the statistical equilibrium states for each value of the interlayer interaction parameter. The theoretical method considered here permits to obtain the important results for the effective chemical potential in the BLG, which shows the extraordinary close results with the recent experimental measurements of the chemical potential in the gated BLG and double BLG heterostructures. For the first time in the literature, we show theoretically, how the charge neutrality point is changing its position when considering the excitonic effects in the BLG system.
We have calculated the A and B sublattice DOS functions in the BLG for different interlayer interaction regimes and for different values of the interlayer hopping parameter. At the zero interlayer interaction case, we have obtained the results very similar to the usual tightbinding DOS in the BLG, and a very large "blue"-shift of the Dirac's neutrality point mediated by the strong excitonic effects in the BLG. At the zero interlayer cou-pling limit, we have shown the main modifications to the usual tight-binding DOS. We have shown that the excitonic condensation mechanism in the charge equilibrated BLG is possible even in the case of the noninteracting layers of the BLG. The principal tunable parameter, in this case, is the interlayer hopping amplitude γ 1 , the large values of which improves the excitonic condensate state. In addition, at any finite and realistic value of the interlayer interaction parameter, it is possible to find the critical value of the interlayer hopping amplitude that renders the BLG into the excitonic condensation state just by suppressing the hybridization gap, present in the system. For example, at the value, W = 0.1331γ 0 and at γ 1 = 0.128γ 0 a very small but finite hybridization gap is present in the BLG. When slightly augmenting the parameter γ 1 to γ 1 = 0.129γ 0 , the hybridization gap is suppressed and the system starts to pass into the excitonic condensate regime. The insulating excitonic state is also suppressed in this case. The principal consequence from this consideration is the following statement: Statement: at each fixed value of the parameter γ 1 , there is a critical value of the interlayer interaction parameter W c above which the hybridization gap opens in the BLG, and when the hybridization gap is present for a certain value of γ 1 then it is possible to find a realistic critical value of the parameter γ 1 itself, at which the hybridization gap closes, rendering the BLG into the possible excitonic condensate state. These statements are not valid only in the case of the very large interlayer interactions, for example at W = 2γ 0 , at which the very high, but approximatively realistic (for example γ 1 = 0.5γ 0 ), values of the interlayer hopping amplitude are not capable of suppressing the very large hybridization gap. Indeed, we think that the charge neutrality at very large W is rather not realistic and could not be achieved experimentally by anyway.
One of the principal achievements, which also ensues from our theoretical model, is the existence of the excitonic condensate states in the BS type bilayer graphene mediated by the interlayer hopping amplitude, even at the finite and relatively small values of the interlayer interaction parameter. The density of states calculations, effectuated in the present paper, show that the excitonic condensate and the excitonic pair formation are fully controlled by the interlayer Coulomb interaction and interlayer hopping. Moreover, in the limit W = 0, there exists an interesting inter-crossover from the hybridized insulating gapped state to the excitonic condensate states in the BLG, mediated by the interlayer hopping mechanism. Therewith, we have shown that the passage when W = 0, is not strictly speaking equivalent to the usual tight-binding description of the BLG.
The different interlayer interaction regimes have been considered in the paper, which correspond to different screening regimes in the bilayer graphene, and which have been discussed only partially in the known literature. From the experimental side of the problem, and taking into account the recent theoretical achievements on the bilayer graphene systems, [4,[9][10][11][12], the Coulomb drag measurements [28][29][30] are promising to observe the excitonic condensation in the pure BLG (without strong disorder) and double BLG heterostructures. For the future, the study the excitonic effects in the hBN intercalated multilayer graphene G/hBN/G could have a breakthrough impact in the technological applications and improvements of these materials as the solid state systems with the sufficiently large band gaps and also due to the recently growing interests in these materi-als for the potential interconnected circuit technologies with the improved high current capacities across these structures, approaching the pristine graphene's working performances.

AUTHOR CONTRIBUTION STATEMENT
All authors contributed equally to the paper.