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 a e-mail: v.apinyan@int.pan.wroc.pl range of the spectrum, the band gap was calculated in reference [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 reference [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 electronhole (e-h) pairs are of particular importance because of their well-defined binding energies, which decides the efficiency of photovoltaic solar cells [43,44].
A new type of high-frequency excitonic effects have been obtained in references [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 interaction 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 references [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 Section 2, we introduce the Hubbard model for the BLG. In Section 3, we obtain the general expression for the longitudinal optical polarization function and the ac optical conductivity. In Section 4, we discuss the numerical results. Finally, in 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 Figure 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 equation (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 reference [31]. The chemical potential term in equation (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 equation (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 equation (1). Namely, we have Parameters U and V in equation (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 equation (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.
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, τ ), The Hamiltonian H (τ ) of the BLG system, in the last term in equation (6), is described in equation (1), and here we will write it in the more convenient form See equation (9) above.
We have introduced in equation (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 equation (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 equation (9), where χ σσ (r, τ ) is defined as The Hamiltonian, in the form given in equation (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 real-space 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 equation (9) we have The integral in the right hand side (r.h.s.), in equation (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 equation (11), can be evaluated by the steepest descent method. We get Here, in order to obtain the r.h.s. in equation (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 equation (6). Thus, we havē n 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, +∆ η c,l (r,τ )(n η l,↑ (r,τ )−n η l,↓ (r,τ )) (13) in equation (9), is also straightforward. Namely, for the l-layer and η-type sublattice variables we have See equation (13) above.
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 equation (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 equation (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 equation (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 reference [11]. The parametersγ lk , in equation (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 l = 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 Fig. 1). The vectors δ and δ are the nearest neighbor vectors in different layers. The components of δ, for the bottom layer-1, are given by is the sublattice constant (with a, being the carbon-carbon length in the graphene sheets). In the layer-2, we have 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: 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 equation (20) with i = 1, . . . , 4, are given as with and The coefficients β ik in equation (21), with i = 1, . . . , 4 are given by the relations The function n F (x), in equations (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 equations (20) and (21), and the changes in the electronic band structure of BLG system, in the presence of the excitonic pairing, are discussed in references [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]).

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 equation (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 equation (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 equation (29), into the Hamiltonian in equation (1) and then expanding 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 equation (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 equation (1). This term is responsible for the paramagnetic part of the current operator, and we neglected the second order diamagnetic contribution in equation (31).
3 The UV optical conductivity

The polarization function
We will consider the retarded polarization function (see 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 B-sublattice 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 The explicit analytical forms of the Fourier transformed single-particle Green's functions in equation (33) where the energy the parameters ε ik are defined in equations (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 equation (22). Next, for the B-sublattice Green function G bb (k, ν n ), we have [52] The coefficients γ ik , in the nominators in equation (35) are defined as with the coefficients a ik , i = 1, . . . , 3, given as We already showed, in reference [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 reference [52]. Next, for the product of the Fourier transformed single-particle Green's functions in equation (33), we get .
Furthermore, we perform the summation over the fermionic Matsubara frequencies ν n , in all terms in equation (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 σ 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 equation (A.2), in 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 equation (40) is the Dirac's delta function. The summation over the reciprocal lattice vectors k, in equation (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 non-interacting 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 equation (42), is given as [54] ϕ For the real part of the optical conductivity function, we obtain finally The function P ij (x), in equation (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 equations (22) and (36).

The optical conductivity
In order to perform explicitly the integration in equation (44), we will use the following rule for the composite Dirac's function: xx ≡ 0. This effect leads furthermore (see Sect. 4) to the total suppression of the Drude region in the spectrum of the optical conductivity 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 σ M G = e 2 /4 [20,56]). After performing the explicit analytical integration in equation (44), we get See equation (46) next page.
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 The frequency dependent functions F mn (Ω) in equation (27) with m, n = 1, 4 are defined as:

See equation (49) next page
where the index permutation function P mn (x), in equation (49) is given previously in equation (45), in the precedent section. The arguments ξ mn 1,2 , in equation (49), are the solutions of the equation Ω + ε n − ε m = 0. These solutions, together with the functions Λ mn (ξ mn 1 ) in equation (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 equation (46) with m, n = 1, 4 are expressed with the help of the parameters a 1 and b 1 , given in equation (28). We obtain We see here that all functions Λ mn (ξ mn 1 ) entering in the equation of the optical conductivity in equation (46), depend on frequency Ω via the solutions ξ mn 1,2 , given in equation (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][60]. 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 [61], the electronic self-energy 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 Figure 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 Figure 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 1, 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 Figure 2. From the results, given in Figure 2, and Table 1, 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 reference [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 Figure 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,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 reference [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 reference [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 σ xx (Ω)/σ Bi = 0. Furthermore, starting from the CNP value of V , the optical gap is increasing (see 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,12] and also in Tab. 1), 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 Ω ∈  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 Figure 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 Fig. 3. The imaginary part of the normalized longitudinal conductivity function in the bilayer graphene, given in equation (53). Different values of the interlayer Coulomb interaction parameter are considered (from zero up to Vc).
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 equation (53) are given in Figure 3. The function σ 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 Fig. 3) are displacing to the lower-Ω regions when augmenting the interaction parameter V up to critical value V c = 1.49γ 0 which is related to the critical value of the lower bound solution of the chemical potential at the CNP: µ c = −0.499γ 0 = −1.498 eV (see in Ref. [11]).
The amplitudes of peaks are increasing when increasing the interaction parameter from 0 up to V c . The multiple peak structure in the imaginary conductivity spectrum (see the peaks structures in the curves corresponding to the negative values of σ xx (Ω)/σ Bi ) is the artifact of the strong excitonic effects in the bilayer graphene. The important observation that could be gained from the results in Figures 2 and 3 and from the values given in Table 1 is that the optical gap fits well with the formula Fig. 4. The real part of the optical conductivity function, given in equation (46) for different values of normalized temperature T /γ0. The interlayer Coulomb interaction parameter is set at V = γ0 = 3 eV.
The shift by γ 1 in equation (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 reference [23], where the Keldysh transport formalism has been used to consider the influence of the band renormalization and excitonic electron-electron 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 Figure 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 Figure 2 and the curve in red in Figure 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, 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 Figure 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 Figure 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.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Appendix A: The current-current operator
The paramagnetic current operator j i (r), given in equation (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] σ ij (Ω) = Π R ij (Ω) Ω , 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 current-current 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 equation (A.3) corresponds to the particlehole bubble, in terms of Matsubara Green's functions. The current operator j i (τ ), in equation (A.3), could be obtained after summing over the lattice sites in the operator given in equation (31). We have j i (τ ) = r j i (r, τ ). (A.4) Here, we have attributed the imaginary time variables to the creation and annihilation operators, in equation (31). After transforming the r.h.s. in equation (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 charge-fluctuations 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 equation (1). 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 equation (A.5) in the expression of the polarization function, given in equation (A.3). 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 correlation function kx v k x G aa (k 0, kτ )G bb (kτ, k 0) +v kx v * k x G bb (k 0, kτ )G aa (kτ, k 0) +v * kx v k x Gãã(k 0, kτ )Gbb(kτ, k 0) +v kx v * k x Gbb(k 0, kτ )Gãã(kτ, k 0) .(A.7) 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 equation (6), in Section 2.1 of the present paper.