High energy shift in the optical conductivity spectrum of the bilayer graphene

We calculate theoretically the optical conductivity in the bilayer graphene by considering Kubo-Green-Matsubara formalism. Different regimes of the interlayer coupling parameter have been considered in the paper. We show that the excitonic effects substantially affect the optical conductivity spectrum at the high-frequency regime when considering the full interaction bandwidth, leading to a total suppression of the usual Drude intraband optical transition channels and by creating a new type of optical gap. We discuss the role of the interlayer coupling parameter and the Fermi level on the conductivity spectrum, going far beyond the usual tight-binding approximation scheme for the extrinsic bilayer graphene.


INTRODUCTION
The introduction of two-dimensional materials and the possibility to control their optical properties brings the new and novel high valued technological applications in nanophotonics, optoelectronics and solar cells [1]. The optical properties of the monolayer (MG) and bilayer graphene (BLG) structures are of great importance in the context of the modern technological applications in the infrared, visible and terahertz range of the frequency spectrum. By applying the external gate voltage one can modify the density of charge carriers and the position of the Fermi level in these systems [2]. The imposition of external electric field tunes the bilayer graphene from the semimetallic to the semiconducting state [3], by allowing for novel terahertz devices [4] and transistors [5]. On the other hand, doped, or electrically tuned graphene and bilayer graphene allow for the spontaneous chiral symmetry breaking states, reflecting in the form of the gapped states in the fermionic quasiparticle spectrum [6][7][8][9][10][11][12]. The spectacular optical properties of bilayer graphene make it as a promising material for infrared optoelectronics. The optical transitions in the system can be alternated, after electrical gating of graphene and BLG [13], and the effects are very similar to the case of the charge transport in the field-effect transistor constructions [14,15]. The optical and charge transport properties in the BLG structures have been widely studied both theoretically [16][17][18][19][20][21][22][23][24][25][26][27][28][29][30] and experimentally [31][32][33][34]. The optical conductivity has been crucial for the optical determination of the band gap formation in the BLG [35][36][37][38]. In the visible range of the spectrum, the band gap was calculated in Ref. 39. A very comprehensive comparison of different microscopic models for the BLG and also the analysis of the optical transitions in the BLG heterostructures is done in Ref. 40. In particular, the effects of short ranged scatterers and screened Coulomb-impurities has been discussed and both intrinsic and extrinsic case of bilayer graphene has been considered in details. The existence of the resonant excitons and bound electron-hole (e-h) pairs in graphene and bilayer graphene structures has been confirmed after recent first principle calculations and experimental studies [41,42]. In the main part of those studies, the excitonic effects enhance from broadly resonant excitonic states, consisting of π and π * bands corresponding the low-frequency regime (up to 10 eV) and with the extremely short lifetimes. The bound electron-hole (e-h) pairs are of particular importance because of their welldefined binding energies, which decides the efficiency of photovoltaic solar cells [43,44].
A new type of high-frequency excitonic effects have been obtained in Refs. 45, 46, in the high-frequency regime (9∼ 20 eV), in the optical spectra of BLG. For the intrinsic graphene, those excitonic effects are the consequence of the unique parallel σ and π * bands, which results in a giant joint density of states. In contrast to ordinary semiconductors or insulators, the excitonic effects in graphene are stronger in the high-energy range of the spectrum [46]. Nevertheless, the most of the theoretical studies on the excitonic effects, discussed above, are limited only to the intrinsic case of the BLG, i.e., when the Fermi level is set to be zero. Therefore, it is of a fundamental interest to consider the excitonic effects in the high-frequency regime at which the BLG shows significant optical activities.
In the present paper, we consider the excitonic effects on the optical conductivity in the BLG and we show how the solution of the Fermi energy in the system affects the optical conductivity spectrum when passing through the balanced charge neutrality point (CNP) in the interacting BLG. We considered the optical conductivity in the BLG by taking into account the intralayer Coulomb inter- action and a wide range of the local interlayer Coulomb repulsion between the electrons on different sublattices of BLG. We show how the excitonic effects shift the spectrum of the optical conductivity and absorption spectrum towards the UV region of photonic frequencies. We suppose the half-filling condition in each of the BLG sheets and for all different values of the interlayer coupling parameter. When evaluating the ac-conductivity in the bilayer graphene, we use the numerical values for the excitonic gap parameter and chemical potential, obtained in Refs.11, 12. We calculate the optical conductivity using the Kubo-Green-Matsubara formalism [47] and we include the local excitonic pairing between the layers of the BLG. We consider different regimes of the interlayer Coulomb interaction parameter, corresponding to different states in the BLG: from semimetallic to semiconducting. The structure of the paper is following: in the Section 2, we introduce the Hubbard model for the BLG. In the Section 3, we obtain the general expression for the longitudinal optical polarization function and the ac optical conductivity. In the Section 4, we discuss the numerical results. Finally, in the Section 5, we give a conclusion to our manuscript. The Appendix, at the end of the paper, is devoted to the calculation of the current-current correlation function.

THE MODEL HAMILTONIAN
The vertical from top-view of the AB-stacked bilayer graphene is presented in Fig. 1. Different lattice site positions in the layer-1 are shown by letters A and B and in the layer-2 by tilde lettersÃ andB. The bilayer Hubbard Hamiltonian for the presented AB-stacked bilayer graphene structure without the external electric field is given by First two terms, in Eq.(1), form the well-known tight binding model. The first term in the Hamiltonian describes the usual hopping of the electrons between the nearest neighbor lattice sites, in a given layer, i.e., Here, a † σ (r) (a σ (r)) and b † σ (r) (b σ (r)) are the electron creation (annihilation) operators and the same operators with the tilde notations refer to the layer-2 in BLG. The parameter γ 0 describes the intralayer hopping in the graphene sheets (the most realistic value of it is γ 0 = 3 eV, given in Ref. 31. The chemical potential term in Eq.(2), has been added in order to deal with the Grand canonical ensemble. Initially, we postulate the equilibrium state in the BLG, i.e., µ 1 = µ 2 ≡ µ. The density operator n ℓσ (r), in the last term in Eq.(2), is n ℓσ (r) = η ℓ n η ℓ σ (r), where we putted η ℓ = a, b for ℓ = 1 and η ℓ =ã,b for ℓ = 2 and n η ℓ σ (r) = η † ℓσ (r)η ℓσ (r) is the usual density operator for a given sublattice, in the layer ℓ. The second term H ⊥ is responsible for the interlayer hopping in the BLG and is given as The parameter γ 1 describes the interlayer hopping between different layers in BLG. The interaction part in the system is given by the last term in the Hamiltonian in Eq.(1). Namely, we have Parameters U and V in Eq.(4) signify local intralayer and interlayer Coulomb interactions in the BLG structure. Indeed, as we will see later on, the consideration of the local interlayer coupling simplifies the problem substantially. This becomes clear after transforming the Hamiltonian in Eq.(1) into the Fourier space representation with appropriate linearisation of the fermionic action of the BLG. We consider the parameter γ 0 as the unit of the energy scale in the considered problem.
A. The interaction term U − V Here, we will show how the interaction terms will be handled in the fermionic-field path integral formalism [50]. For this, we pass into the Grassmann representation for the fermionic variables, and we write the partition function of the system in the imaginary time fermion path integral method. We introduce the imaginary-time variables τ , at each lattice site r and the variables τ vary in the interval (0, β), where β = 1/T with T being the temperature. The grand canonical partition function of the system is and the fermionic action S X , X,Ȳ , Y is given as follows Here, the first two terms are the Berry terms for the layers with the indices ℓ = 1, 2 where we have introduced the following notation for the fermionic operators: X 1,σ (r, τ ) = a 1,σ (r, τ ), X 2,σ (r, τ ) = a 2,σ (r, τ ), Y 1,σ (r, τ ) = b 1,σ (r, τ ) and Y 2,σ (r, τ ) = b 2,σ (r, τ ). The Hamiltonian H (τ ) of the BLG system, in the last term in Eq.(6), is described in Eq.(1), and here we will write it in the more convenient form We have introduced in Eq.(9) the z-component of the generalized spin operator S η l (r, τ ) = 1/2 α,β=↑,↓η l,α (rτ )σ αβ η l,β (r, τ ), for different sublattices in the layers of the BLG structure. It is defined as S η l,z (r, τ ) = 1/2 (η l,↑ (r, τ ) − η l,↓ (r, τ )). The chemical potentials µ 1 and µ 2 in Eq.(9) are the shifted chemical potentials, defined as: µ 1 = µ + U/2, µ 2 = µ + U/2 + V . Indeed, the chemical potentials of electrons on the nonequivalent sublattice sites get different shifts in different layers due to the stacking ordering of the BLG structure. We have introduced the new complex variables χ σσ ′ (r, τ ) and their complex conjugatesχ σσ ′ (r, τ ) in the last term in Eq.9), where χ σσ ′ (r, τ ) is defined as The Hamiltonian, in the form given in Eq. (9), is more suitable for further decouplings of four-fermionic terms, within the Hubbard-Stratanovich saddle-point linearisation procedure. We give here the procedure of the realspace linearization of the U − V terms for the case of the a-sublattice in the layer-1 in the BLG. Namely, for the U -term in Eq. (9) we have The integral in the right hand side (r.h.s.), in Eq. (11), is over the decoupling field variables V a 1 (r, τ ) (which are introduced at each sublattice site position r and at each time τ ) coupled to the density term n a 1 (r, τ ). The field integral, in r.h.s., in Eq.(11), can be evaluated by the steepest descent method. We get Here, in order to obtain the r.h.s. in Eq.(12), we have replaced the field integration over V a 1 (r, τ ) by the value of the function in the exponential at the saddle-point value of the decoupling potential υ a 1 = iU/2 (n a 1 − 2µ 1 /U ), where the average densityn a 1 is defined with the help of the total action of the system, given in Eq. (6). Thus, we haven a 1 = n a 1,↑ (r, τ ) + n a 1,↓ (r, τ ) . The same procedure could be repeated also for the nonlinear density terms on the other sublattices in the BLG structure, containing nã 2 (r), n b 1 (r) and nb 2 (r) density terms. The decoupling of the nonlinear density-difference term, in Eq. (9), is also straightforward. Namely, for the l-layer and η-type sublattice variables we have +∆ η c,l (r,τ )(n η l,↑ (r,τ )−n η l,↓ (r,τ )) .
The saddle-point values of the variables ∆ η c,l (r, τ ) are given by δ η c,l = U/2 n η l,↑ (r, τ ) − n η l,↓ (r, τ ) . Thus, it is proportional to the difference between the electron densities with the opposite spin polarizations. For simplicity, we suppose the case of the spin balanced BLG layers, with equal density numbers for each spin direction, i.e. n η l,↑ (r, τ ) = n η l,↓ (r, τ ) and the quantities δ η c,l vanish in the case. For the case of the half-filling occupation, considered here, we put n η l,σ (r, τ ) = 1/2, for each spin direction σ =↑, ↓.
Next, for decoupling the last term interaction V -term in Eq.(9), we apply the complex form of the Hubbard-Stratanovich transformation [50] for the one-component fermion-field In fact, the saddle-point value of the decoupling field Γ σσ ′ (r, τ ), introduced in Eq. (14), is directly related to the excitonic gap parameter. Indeed, we have We consider here the homogeneous BLG structure with the pairing between the particles with the same orientation of spin variables, i.e. ∆ σσ ′ = ∆ σ δ σσ ′ . Furthermore, we can write the total action of the system in the Fourier representation, given by the transformations η l,σ (r, τ ) = 1 βN k,νn η k,σ (ν n )e i(kri−νnτ ) , where ν n = π (2n + 1) /β, with n = 0, ±1, ±2, ..., are the fermionic Matsubara frequencies [51], and N is the total number of sites on the η-type sublattice, in the layer l. We introduce the four component Nambu-spinors at each discrete state k in the reciprocal space and for each spin direction σ =↑, ↓, Then the action of the system reads as Here, G −1 k,σ (ν n ), is the inverse Green's function matrix, of size 4 × 4. It is defined as The diagonal elements of the matrix in Eq.(17) are the energy parameters where, the effective chemical potentials µ eff 1 and µ eff 2 , are defined with the help of the intralayer and interlayer interaction parameters U and V as Thus the effect of the interaction parameters is such that the chemical potential µ gets shifted in BLG system and also the position of the CNP point gets shifted, as it was discussed in Ref. 11. The parametersγ lk , in Eq. (17), l = 1, 2, are the renormalized (nearest neighbors) hopping amplitudesγ lk = zγ lk t, where the k-dependent parameters γ 1k and γ 2k are the energy dispersions in the BLG layers with l = 1 and l = 2, respectively. We have γ 1k = 1/z δ e −ikδ for ℓ = 1 (and γ 2k = 1/z δ ′ e −ikδ ′ for the layer with l = 2). The parameter z is the number of the nearest neighbors lattice sites on the honeycomb lattice for a given sublattice variable and z = 3 for each monolayer (see in Fig. 1). The vectors δ and δ ′ are the nearest neighbor vectors in different layers. The components of δ, for the bottom layer-1, are given by , and a 0 = √ 3a is the sublattice constant (with a, being the carbon-carbon length in the graphene sheets). In the layer-2, we havebmδ where a is the carbon-carbon interatomic distance. By the convention, we put a ≡ 1, for both layers. For a given geometry of the AB-stacked BLG it is not difficult to realize that γ 2k = γ * 1k ≡ γ * k and it follows that γ 2k =γ * 1k ≡γ * k , where we have omitted the layer index l.
We assume here that the pairing gap is real and is not spin-dependent (∆ σ ≡ ∆ =∆). Therefore, the structure of the Green's function matrix does not changes for the opposite spin direction:Ĝ −1 k,↑ (ν n ) ≡Ĝ −1 k,↓ (ν n ). Next, we use the half-filling condition in each layer of the BLG system which determine the chemical potential in the BLG.
For the layer-ℓ, this condition holds thatn a ℓ +n b ℓ = 1. Here, we present the resulting system of coupled nonlinear self-consistent equations for the chemical potential µ and excitonic pairing gap parameter ∆. We get where the dimensionless coefficients α ik , in Eq. (20) with i = 1, ..4, are given as with and The coefficients β ik in Eq. (21), with i = 1, ..4 are given by the relations The function n F (x), in Eqs. (20) and (21), is the Fermi-Dirac distribution function n F (x) = 1/ e β(x−µ) + 1 .
The energy parameters ε ik define the interacting band structure in the BLG in our problem of the excitonic effects in the BLG. They are given by the following relations The exact numerical solution of Eqs.(20)- (21), and the changes in the electronic band structure of BLG system, in the presence of the excitonic pairing, are discussed in Refs.11, 52. The parameterμ plays the crucial role in the whole physics related to the excitonic BLG. It plays the role of the exact Fermi energy in the BLG and it is defined with the help of the effective chemical potentials asμ = 1/2 µ eff 1 + µ eff 2 . As the results show, the Fermi energy in the bilayer graphene as a function of the interaction parameter V and at T = 0 has also a very large jump at the CNP, similar to the exact chemical potential in BLG, and it happens nearly at the same value of the interlayer interaction parameter V = 1.49γ 0 . Moreover, at the zero interlayer interaction limit, the bare chemical potentialμ coincides with the Dirac's crossing energy level ε D [11]. The Coulomb interaction parameter U redefines the Fermi level in the BLG by the way that ε F =μ = µ + κU + 0.5V (with κ = 0.25) [11,12] and the chemical potential µ could be calculated self-consistently (see in Ref. 11).

B. The effect of the electric field
We suppose that the bilayer graphene gets excited in the external electromagnetic field, with the electric field component polarized along the x-axis in the graphene's sheet. Then, in order to consider the electric current response of the system, we include the vector potential A(r) in the tight binding part of the Hamiltonian in Eq.(2). This could be done via the Peierls-Onsager substitution [48,49] with C 1 = a, b for ℓ = 1, and C 2 =ã,b for ℓ = 2. The electron operators in the interaction terms do not get modified in that case because of the local nature of those terms. The shift of the operators product in Eq.(29) by a phase factor is related to the electronic Wannier states modifications in that case (for a detailed description, see in Ref. 49). After inserting the transformations, given in Eq. (29), into the Hamiltonian in Eq.(1) and then expand-ing it up to first order in vector potential A(r), we get where c is the speed of light and j ℓi (r) is the component of the current operator along the direction i, in the given sheet ℓ. It could be obtained from the expression in Eq. (30), after the functional differentiation of the Hamiltonian with respect to A ℓi (r) (here i indicates the component of the vector potential A, along the considered direction i). For the total current density operator in the BLG, we get straightforwardly We consider here only the linear contribution of the external vector potential into the Hamiltonian in Eq.(1). This term is responsible for the paramagnetic part of the current operator, and we neglected the second order diamagnetic contribution in Eq.(31).

THE UV OPTICAL CONDUCTIVITY
A. The polarization function We will consider the retarded polarization function (see in the Appendix A) in order to calculate the optical conductivity in the system. In the continuum approximation at the K point, we have for the velocity operator |v kx | 2 = v 2 F (see also in Ref. 20), where v F is the Fermi velocity, which relates to the intralayer hopping parameter γ 0 , i.e., v F = √ 3a 0 γ 0 /2 . For the A-sublattice in the bottom layer-1, we have defined the normal Green's functions as follows The similar expression could be written also for the Bsublattice Green's function G bb (kτ, kτ ′ ). For the top layer sublattice Green's functions Gãã and Gbb we have the relations Gãã = G bb and Gbb = G aa [52]. Then, for the component of the polarization function, along the current direction, we get +Gãã(k, ν n )Gbb(k, ν n − ω m ) +G aa (k, ν n )G bb (k, ν n + ω m ) +Gãã(k, ν n )Gbb(k, ν n + ω m ) . (33) The explicit analytical forms of the Fourier transformed single-particle Green's functions in Eq.(33) could be obtained after the functional derivation techniques [52] where the energy the parameters ε ik are defined in Eqs. (27) and (28). Here, ε 1k and ε 4k are the split valence and conduction bands, and ε 2k and ε 3k are the low energy conduction and valence bands, according to the usual tight-binding definitions [12,40]. The coefficients α ik , figuring in the nominator in the sum, in Eq. (22). Next, for the B-sublattice Green function G bb (k, ν n ), we have [52] G bb (k, The coefficients γ ik , in the nominators in Eq.(35) are defined as with the coefficients a ′ ik , i = 1, ...3, given as We already showed, in Ref. 52, that the semiconducting (or insulating) state could be reached from the semimetallic limit in the BLG, thus leading to the enhancement of the Bardeen-Cooper-Schrieffer (BCS)-Bose-Einstein-Condensate (BEC) type crossover mechanism in the interacting BLG system. A detailed description of such a transition is given in Ref. 52. Next, for the product of the Fourier transformed single-particle Green's functions in Eq.(33), we get Furthermore, we perform the summation over the fermionic Matsubara frequencies ν n , in all terms in Eq. (33). Then, we obtain for the polarization operator The retarded polarization function enters into the expression of the real part of the longitudinal conductivity function Re{σ} xx and could be obtained after the standard analytical continuation techniques into the real frequency axis of the upper-half complex semi-plane. This procedure is given in Eq.(A2), in the Appendix A. Then, for calculating the imaginary part of the polarization function, we use the real line version of the Sokhotskii-Plemelj identity, i.e., where P denotes the Cauchy principal value, and the parameter ζ takes two values ±1. The function δ(x) in Eq. (40) is the Dirac's delta function. The summation over the reciprocal lattice vectors k, in Eq.(39), could be replaced by the integration over the continuous variables via the introduction of the density of states (DOS) function ρ(x) for the non-interacting bilayer graphene sheets, i.e., k . . . = dxρ(x).... The DOS in the noninteracting graphene layer, is defined as Beyond the Dirac's approximation it could be analytically expressed [3,54] as where K(x) is the Elliptic integral of the first kind [55] K(x) = π/2 0 dt/ 1 − x 2 sin 2 t . The function ϕ(x), in Eq. (42), is given as [54] ϕ For the real part of the optical conductivity function, we obtain finally The function P ij (x), in Eq. (44), is the index-permutation function, given as where the coefficients α i (x) and γ i (x) are the continuous versions of the coefficients, given in Eqs. (22) and (36), above.

B. The optical conductivity
In order to perform explicitly the integration in Eq.(44), we will use the following rule for the composite Dirac's function: Then, after some calculations, we get the analytical expression for the real part of the conductivity function, in which only the interband transitions contribute (here, we normalize conductivity function in units of σ Bi = e 2 /2 , i.e., which is twice of the dc conductivity in the monolayer graphene σ MG = e 2 /4 [20,56]). After performing the explicit analytical integration in Eq.(44), we get Here, Θ(x) is the Heaviside's unit step function, and with α = a, b. The frequency dependent parameters a(Ω), b(Ω) and the interaction induced parameters a 1 , b 1 are defined in the following way a(Ω) = 2 (Ω + ∆ + γ 1 ) , b(Ω) = 2 (Ω − ∆ − γ 1 ) , The frequency dependent functions F mn (Ω) in Eq. (27) with m, n = 1, 4 are defined as: where the index permutation function P mn (x), in Eq. (49) is given previously in Eq. (45), in the precedent Section. The arguments ξ mn 1,2 , in Eq. (49), are the solutions of the equation Ω + ε n − ε m = 0. These solutions, together with the functions Λ mn (ξ mn 1 ) in Eq. (48), are different for different optical transitions in the system, and for each transition m → n they should be separately specified. We have Next, the functions Λ mn (x) in Eq. (46) with m, n = 1, 4 are expressed with the help of the parameters a 1 and b 1 , given in Eq. (28). We obtain We see here that all functions Λ mn (ξ mn 1 ) entering in the equation of the optical conductivity in Eq. (46), depend on frequency Ω via the solutions ξ mn 1,2 , given in Eq. (50).
We have considered here the pure BLG system. The inclusion of the effects of the scattering by dilute charged impurities could be done in the self-consistent Born approximation [57][58][59]. In this case, the effect of the finite external bias voltage should be included directly and the dependence of optical conductivity on the impurity concentration and bias voltage should be considered. The universal, disorder independent, electrical conductivity at low temperatures will be relevant in this case, leading to the localized states near the Fermi level in the electronic density of states. According Born approximation in the scattering theory [60], the electronic selfenergy matrix Σ(iν n ) of disordered system in the presence of finite but small density impurity atoms is expressed in terms of the local Green function of clean system G(k, iνn). Then the perturbative expansion for the Green's functionḠ(k, iν n ) of disordered system could be obtained via the Dyson equation given as The scattering rate generated by the effect of the disorder has the energy of the order of the Fermi energy. Therefore the inclusion of the effects of disorder will have a strong impact on the optical properties in the BLG, particularly on the optical conductivity. The study of the effects of the disorder is extremely important for the technological applications of the BLG because it is almost sensitive to the unavoidable disorder created by the substrate on which it is deposited, adatoms, ionized impurities, dopping, etc. The inclusion of the effect of disorder in our problem of the optical conductivity in the presence of the excitonic pairing could be done within the bilayer Hubbard model considered here. Especially, the large bandwidth of the Coulomb interaction potentials U and V is possible to consider within the present approach which will correspond to different screening regimes of the external potential by the impurity atoms in the BLG. Unfortunately, this subject it is out of the scope of the present paper.

THE NUMERICAL RESULTS
In Fig. 2, we have presented the numerical results for the real part of the longitudinal optical conductivity function, normalized to the dc conductivity in the bilayer graphene σ Bi = e 2 /2 . Different values of the interlayer Coulomb interaction parameter are considered (from zero up to intermediate V c = 1.49γ 0 = 4.47 eV and higher values) in Fig. 2 and the zero temperature limit is set for the presented plots. The intralayer Coulomb interaction parameter is fixed at the value U = 2γ 0 = 6 eV, and the interlayer hopping amplitude is γ 1 = 0.128γ 0 = 0.384 eV. In table I, we give the obtained values of the chemical potential [52], the threshold frequency values (the optical gaps) for the real part of the conductivity function, and the Fermi energy in the BLG corresponding to the interlayer interaction parameters, considered in Fig. 2. From the results, given in Fig. 2, and table I, we can see the general behavior of the real part of the optical conductivity function, as a function of the chemical potential µ and Fermi level ε F in the bilayer graphene. First of all, let's mention that the BLG system is automatically extrinsic in our case (this is the case when ε F = 0 because of the finite number of the electron-hole pairs), in difference with the discussion in Ref. 40, where both intrinsic and extrinsic cases have been considered equally, and the optical gap has been attributed to the first large peak in the intrinsic optical conductivity spectrum. Particularly, for all values of the interaction parameter V , we observe the presence of a very large optical gap in the real part of the conductivity spectrum in Fig. 2. Considering the case of the noninteracting layers, i.e., when V = 0, we see that the highest peak in optical conductivity spectrum (see the deepest black plot in Fig. 2) is situated at the value Ω = 2.856γ 0 , or at the high-energy resonant position at Ω = 8.57 eV, which is very close to the excitonic resonance at 8.3 eV obtained in recent abinitio many-body calculations of the high energy optical absorption in graphene [46]. Moreover, when augmenting the interaction parameter from zero up to the critical value V c = 1.49γ 0 = 4.47 eV the optical gap is decreasing, while the optical conductivity is largely increased, and the peaks positions become more apparent. The value V c = 1.49γ 0 , at which the optical gap attains its minimum ∆ ω0 = 1.337γ 0 = 4.011 eV, (corresponding to the photon wavelength λ = 309.15 nm, in the near UV range of the spectrum), is the threshold value of V , above which the chemical potential and the Fermi level jump to their upper bound solutions (for details, see in Refs.11, and 52) and the critical value V c = 1.49γ 0 plays the role of the new CNP in the interacting BLG, as it was explained in Ref. 11. Let's mention also that the value of ∆ ω0 at V c coincides exactly with the value of the hybridization gap ∆ H in the BLG, obtained in Ref. 52). Thus, we observe that for V < V c the interaction effect acts to redshift the conductivity peaks in the BLG.   For Ω < ∆ ω0 , we have Re σ xx (Ω)/σ Bi = 0. Furthermore, starting from the CNP value of V , the optical gap is increasing (see in Fig. 2) and becomes very large at V = 4γ 0 = 12 eV, of the order of ∆ ω = 12.48 eV.
On the other hand, there are no remarkable changes in the conductivity peaks amplitudes, in this case, corresponding to the different optical transitions in the bilayer graphene. The other important feature in the strong interaction regime is related to the blue-shift effect of the conductivity peaks above the CNP value of V , i.e. when V > V c . Thus, the general conclusion is that when the chemical potential and the Fermi energy in the BLG pass to their upper bound solutions as a function of the interlayer interaction parameter (see about in Refs.11 and 12 and also in Table I), the excitonic pairing causes to change the red-shift effect of the conductivity peaks into the blue-shifted ones.
From zero up to critical value of the interlayer interaction parameter V ∈ [0; V c ], the optically active photon frequency region is given by the frequency interval Ω ∈ [4.011; 19] eV, and corresponds to the photon's wavelengths λ ∈ [65.26; 309.15] nm, thus situating in the range, which starts from the extreme ultraviolet part of the light spectrum and increases up to near UV region. We observe also that inside this region of the interaction parameter, there is no conductivity spectrum displacement in the far UV side of the spectrum, and the observable changes occur only in the near UV part of the photon energies, governing the changes of the optical gap ∆ ω . Contrary, for the strong interlayer interaction regime, the considerable changes occur on both sides of the spectrum. For the large interaction values, considered in Fig. 2, the conductivity spectrum is squeezing when increasing the interaction parameter V , and the near UV part of the spectrum gets displaced into the smallest wavelengths sides, i.e., along with the UV-c part of the spectrum. For the very strong interaction V = 4γ 0 = 12 eV, the corresponding active optical frequencies are given by the energy interval Ω ∈ [12.47; 22] eV and the corresponding wavelengths are given in λ ∈ [56.17; 99.43] nm. Thus, in this case, the strong excitonic effects are present in the system and, experimentally, it would be quite difficult to observe them in this limit because of the very narrow region of the permitted photon's wavelengths.
The imaginary part of the optical conductivity function σ xx (Ω) can be easily calculated using the Kramers-Kronig formula, which relates the real and imaginary parts of the complex optical conductivity function. It is given as where a special attention should be paid to the singularity points Ω ′ = ±Ω when performing the numerical integration. The plots of the imaginary part of the conductivity function, normalized to the dc conductivity of the BLG σ Bi and calculated with help of the relation in Eq.(53) are given in Fig. 3. The function Im σ xx (Ω)/σ Bi has been evaluated for different values of the interlayer interaction parameter (from zero up to the critical value V c = 1.49γ 0 ) and is presented as a function of the normalized frequency Ω/γ 0 . We see that the positions of the negative conductivity peaks (see in The shift by γ 1 in Eq.(54) is due to the interlayer hopping. This result corresponds approximatively to the interband absorption threshold of the optical conductivity in the statically screened Thomas-Fermi regime, which has been well discussed in Ref. 23, where the Keldysh transport formalism has been used to consider the influence of the band renormalization and excitonic electronelectron interaction effects on the optical conductivity in doped bilayer graphene. The Mahan's exciton bound state vanishes in our case due to the strong interaction effects, present in the BLG. We considered here the suspended bilayer graphene where the interaction effects have a spectacular impact on the conductivity spectrum because the dielectric environment portion of the screening is absent.
In Fig. 4, we have presented the optical conductivity function for different values of temperature. The interlayer interaction parameter is fixed at the value V = γ 0 , corresponding to the plot in orange, in Fig. 2 and the curve in red in Fig. 3. We observe that when increasing the temperature, the far-UV part of the spectrum is unchanged, while the near UV-c side of the spectrum (simultaneously with the optical gap) is considerably modified and spread over the visible region of the photon's spectrum, i.e., for example, at the very high temperature T = 0.5γ 0 the photon wavelength corresponding to the optical gap ∆ ω = 0.822γ 0 = 2.46 eV, is of order λ = 1240/∆ ω = 502.8 nm, situating in the visible range of the light spectrum. We see that the optical gap and the amplitudes of the conductivity peaks are decreasing when increasing the temperature. In Fig. 5, we give the explicit temperature dependence of the real part of the optical conductivity function for the fixed interaction parameter V = γ 0 = 3 eV and for the incident photon wavelength of the order of λ = 165.651 nm. For that given wavelength, the optical conductivity function attains its dc limit σ Bi only at the very high temperatures, of the order of T ∼ 0.7γ 0 . The temperature dependence of the optical gap ∆ ω (within the same temperature range) is shown in the inset, in Fig. 5, for the case V = 3 eV. We see that for the temperatures T > 0.07γ 0 , the optical gap varies very slowly with temperature.

CONCLUDING REMARKS
We have calculated the optical conductivity function in the interacting bilayer graphene system, exposed in the external longitudinal electric field and in the presence of the excitonic pairing interaction in the system. We have used the exact four-band model of the BLG system, by avoiding low-energy approximation in the band dispersion spectrum. We have shown that for a fixed value of the interlayer hopping amplitude γ 1 , the excitonic pairing interaction affects considerably the optical conductivity spectrum in the bilayer graphene, which becomes more pronounced when increasing the interlayer interaction parameter V . We have shown that the inclusion of the excitonic effects in the system leads to the total suppression of the low-frequency Drude-contribution, related to the intraband optical transitions in the system. We have calculated the optical conductivity function for different values of the interlayer interaction parameter, and we have shown the existence of the optical gap in the extrinsic case of the bilayer graphene. The spectacular features in the conductivity spectrum are related to the behavior of the chemical potential and Fermi energy in the BLG when modifying the interlayer coupling parameter V . At the CNP, corresponding to the critical interlayer coupling strength V c , the Fermi energy passes into its upper bound solution through CNP, and the optical conductivity spectra get switched inversely, turning from the redshifted lines into the blueshifted ones. This effect dominates along the whole optical spectrum, and we suggest that it should affect considerably also on the other optical properties in the bilayer graphene such as the refractive index, the optical absorption spectra, and the optical reflectance. The general rule, obtained in the present paper, is that the optical conductivity spectrum and the optical gap parameter are strongly governed by the position of the Fermi energy in the interacting bilayer graphene.
The obtained results could represent a considerable interest when studying experimentally the excitonic effects in the optical conductivity spectrum in the strongly interacting bilayer graphene. Furthermore, the results obtained here would represent as the reference for the calculation of the other optical properties in the system such as the optical absorption spectra, reflectivity, and the electron energy loss spectra. The results obtained in the present study could spread a new insight into the possible optoelectronic and nano-optical applications of the bilayer graphene.

AUTHOR CONTRIBUTION STATEMENT
All authors contributed equally to the paper.
Appendix A: The current-current operator The paramagnetic current operator j i (r), given in Eq. (31), describes the in-plane electric current density, driven in the BLG layers by the external electric field component. We present here the calculation of the real part of the ac conductivity σ ij (Ω) in the BLG with the help of the retarded polarization function Π R ij (Ω) and by using the standard Kubo-Green-Matsubara formalism [47] Re{σ ij (Ω)} = Im Π R ij (Ω) Ω , where Π R ij (Ω) = Π R ij (iω m → Ω + iη + ) (A2) is the real frequency retarded function, which is obtained after the analytical continuation of the bosonic Matsubara frequencies iω m in to the upper half of the real axis in the complex plane, i.e., iω m → Ω + iη + with η + being the infinitesimal positive constant, and the bosonic Matsubara frequencies ω m are give as usual by ω m = 2πm/β, where m = 0, ±1, ±2, ..., and β = 1/k B T . The currentcurrent correlation function, in turn, is given in the Kubo-Green-Matsubara formalism [47] as where τ is the imaginary time, and T τ is the chronological time ordering operator [47,50]. The response function Π R ij (iω m ) in Eq.(A3) corresponds to the particle-hole bubble, in terms of Matsubara Green's functions. The current operator j i (τ ), in Eq.(A3), could be obtained after summing over the lattice sites in the operator given in Eq. (31). We have j i (τ ) = r j i (r, τ ). (A4) Here, we have attributed the imaginary time variables to the creation and annihilation operators, in Eq. (31). After transforming the r.h.s. in Eq.(4) into the Fourier k-space, and after summing over the lattice sites positions, we get the following expression for the total current density operator is the electron velocity operator in the individual graphene sheet. There exist several theoretical approaches for treating the BLG in the presence of the external fields. Usually one considers the applied external bias voltage [20,24], which makes the considerable potential difference between the layers in the BLG. Moreover, it leads to the charge imbalance between the layers. The electron Coulomb interaction between the layers was not considered is such approaches [20,24] and the tight-binding description is discussed only. The authors found a very large asymmetry gap, induced in the quasiparticle excitation spectrum, which leads to the chargefluctuations in the system. In another treatment [11,12], the intralayer and interlayer Coulomb interactions have been included properly in the Hamiltonian of the bilayer graphene, given in Eq.(1), here. We follow here the second routine, described above, in order to calculate the ac conductivity in the interacting bilayer graphene with the presence of the excitonic pairing interaction. We perform the full bandwidth 4-band calculation scheme, which goes beyond the usually known Dirac-cone approximation [3]. We will concentrate here on the calculation of the four-point fermionic correlation functions that appear after putting Eq.(A5) in the expression of the polarization function, given in Eq.(A3). We consider here only the x component of the vector potential, thus the current density operator is given in the same direction j x (τ ) as the electric field component E x , which oscillates parallel to the graphene's layers. Next, after performing the Wick averaging procedure [51], we get a very compact expression for the current-current cor-relation function The single-particle Green's functions G aa and G bb are non zero only if k = k ′ , therefore, G aa (k ′ 0, kτ ) = δ kk ′ G aa (k0, kτ ). This follows from the symmetry of the action, given in Eq. (6), in the Section 2 A of the present paper.