Excitonic Tunneling in the AB-bilayer Graphene Josephson Junctions

We have considered the AB-stacked bilayer graphene tunnel junction construction. The bilayers are supposed to be in the charge equilibrium states and at the half-filling in each of the electronic layers of the construction and at each value of the external gate potential. By considering the interacting bilayers in both sides of the junction and by taking into account both intralayer and interlayer Coulomb interaction effects, we have calculated the normal and excitonic tunnel currents through the junction. The electronic band renormalizations have been taken into account, due to the excitonic pairing effects and condensation in the BLGs. The exact four-band energy dispersions, including the excitonic renormalizations, have been used for the bilayers without any low-energy approximation. The normal and excitonic tunneling currents have been calculated for different values of the gate potential and for different values of the interlayer interaction parameters in both sides of the tunnel junction. We demonstrate the existence of the excitonic Josephson current through the junction which persist even for the non-interacting bilayer graphene junction. For the non-interacting case, the mechanism of the excitonic condensates formation and tunneling between the condensates is attributed to the interlayer hopping between the layers. The role of the charge neutrality point has been discussed in details.

The excitonic gap formation and condensation has been examined also in the bilayer graphene structures [49][50][51][52][53][54][55][56][57]. Namely, the bilayer graphene is very promising for the optoelectronic applications due to its unique gate-controllable band structure properties [58]. The imposition of the external electrical field can tune the bilayer graphene from the semimetal to the semiconducting state. Nevertheless, the excitonic condensation in the bilayer graphene structures remains controversial in the modern solid-state physics because of the complicated nature of the single-particle correlations in these systems [49][50][51][52][53][54][55][56][57]. It has been shown recently [59] that the critical temperature, which describes the transition from the condensate state to the normal state in graphene double-layer structures, can be very high due to the extremely small effective mass of excitons. The coherence in exciton BEC condensates survives at the very high temperatures. An analog conclusion has been drawn in Ref. [56,57], concerning the bilayer graphene, where the condensate evolution has been analyzed as a function of the interlayer Coulomb interaction parameter in the BLG.
The excitonic condensation has been realized experimentally in the double bilayer graphene heterostructure, in the strong quantum Hall regime, with the help of a combination of the Coulomb drag and current counterflow measurements [60]. They have also found the evidence of strong interlayer coupling between graphene layers, thanks to the observation of quantized Hall "drag plateau". The zero-valued longitudinal resistance, measured there, confirms the dissipationless (friction-free) nature of the electron-hole condensate state. A quite simple experimental way to observe the excitonic condensate states in the bilayer structures is related to the possibility of engineering of a spatially confined excitonic condensates in the potential traps, and the investigation of the Josephson tunneling effects for excitons [61][62][63][64], related to the tunneling between two trapped Bose condensates [65] that possess the macroscopic phase coherence. The excitonic Josephson tunneling effects and thermal transport properties in the electron-hole type double-layer graphene junctions, separated by a dielectric layer, have been recently considered in Refs. [66,67].
In the present paper, we study the excitonic tunneling effects in the tunnel junction based on the AB-stacked bilayer graphene structures. We suppose the presence of the macroscopic phase coherent regime with the well-defined condensates phases and amplitudes and we consider only the local, on-site, interlayer excitonic pairing in each side of the junction. We treat the local intralayer and interlayer Coulomb interactions within the bilayer Hubbard model at the half-filling (in each layer of the BLG). Supposing the electronic bilayers, without the initial optical pumping mechanism, we study the normal and excitonic tunneling currents through the BLG/I/BLG junction for different values of the external gate voltage, applied to the heterostructure, including the dc limit (i.e., V = 0) when a dc excitonic Josephson current passes through the tunnel junction. We will assume the half-filling regime in each layer of the BLG structures, even in the presence of the applied gate potential, thus by supposing that not considerable changes occur in the electron density during the adiabatic switching of the external potential.
Here, by the adiabatic switching, we mean the very slow changes of the external switching potential applied to the junction, which causes the appearance of the ac Josephson current in the junction. The adiabatic changes of the external gate potential are important in many contexts concerning the tunnel junctions. First of all, here the adiabaticity of the external potential is important for the non-perturbative influence of the gate potential on electron densities in the individual bilayer structures, on both sides of the junction. Furthermore, this leads to the possibility of examination of the permitted values of interlayer Coulomb interaction parameters and excitonic gaps for a fixed value of the external gate potential. The half-filling condition in each layer of the BLGs is directly related to the adiabaticity of the external gate potential in the sense that the electronic reconfiguration time τ rec in the BLG system is much faster than the variation time of the external potential τ Field : τ rec << τ Field . Such instantaneous changes of the electron densities in the system lead to the well-defined Coulomb potentials, for each value of the external field configuration. On the other hand, we consider the explicit dependence of the condensates phase difference parameter on the external gate potential and we show that the excitonic dc Josephson current appears in the junction array at the zero value of the external gate potential. We show how a finite difference between the static phases of the coherent excitonic condensates in the BLG subsystems leads to the excitonic Josephson dc current through the tunnel junction even in the case of the non-interacting regime in the BLG systems. The adiabatic switching of the external gate potential could be very important in the context of the adiabatic graphene logic constructions similar to the superconductor logic constructions, which are actually under the strong theoretical and experimental investigations [68][69][70]. In such devices, the intrinsic switching speed of the junctions increases rapidly with the increase of the damping resistance, and the phase differences change more adiabatically, which results in a considerable reduction in the energy dissipation of the single logical gates.
Furthermore, we show also that any finite voltage leads to the ac Josephson current irrespective of the phases of coherent condensates in both sides of the junction and we calculate the ac Josephson current as a function of external gate potential, for different values of the interlayer Coulomb interactions in different sides of the junctions. We study the amplitude of the zero voltage Josephson dc current as a function of the interlayer Coulomb interaction parameters in the subsystems. The symmetric and asymmetric interaction cases have been considered straightforwardly. Also, we calculate the normal quasiparticle tunneling current and we show that normal tunneling in the BLG/I/BLG heterostructure is an interaction-protected process and the threshold frequency of the normal tunneling current strongly depends on the values of the Coulomb interaction parameters in the BLGs. We analyze the role of the charge neutrality point (CNP) on the behavior of the normal and excitonic tunneling currents.

Description of the Hubbard Interactions
We introduce here our model consisting of two Bernal stacking bilayer graphene (BLG) structures separated by a very thick dielectric layer (We suppose that the thickness of the insulating layer is such that we can neglect the quasiparticle scattering and recombination processes in the layer.) For the convenience, we will denote by a, b andã,b (and their conjugates a † , b † andã † ,b † ) the annihilation (creation) fermionic operators corresponding to different sublattice sites A, B in the left bottom, andÃ, B in the left top layer in the left BLG. Similarly, we denote by c, d andc,d (and their conjugates c † , d † andc † ,d † ) the annihilation (creation) fermionic operators corresponding to different sublattice sites C, D in the right bottom andC,D in the right top layer in the right BLG. In Fig. 1, we have presented the schematic setup of our BLG/I/BLG junction (here I represents the dielectric layer between the BLGs). We assume here that the layers, which have lattice sites A, B and C, D (bottom layer) are biased-V /2 (with = L, R), and the layers with the lattice sitesÃ,B andC,D (top layers) V /2, so that the potential difference between the two layers is V . For the first treatment of such a junction, we suppose the half-filling condition satisfied in both BLG systems, i.e., we suppose that n = 1, where n is the total particle number operator in each layer with = 1, 2 of each BLG with = L, R. Note, also that we attached the number = 1 to the bottom layers and = 2 to the top layers, in the heterostructure. The non-interacting tight-binding part of the total tunnel junction Hamiltonian could be written in the usual form with H 0 , given as where the fermionic operators X and Y refer to different sublattice fermions, i.e., for = L we have X = a,ã, and Y = b,b, while for = R we have X = c,c and Y = d,d. The intralayer hopping parameter γ 0 is supposed to be the same for both layers and for both BLGs. The fermionic operators P † and Q † in Eq. (2) are defined such that for = L we have P = b and Q =ã, whereas for = R we defined P = d and Q =c. The parameter γ 1 denotes the interlayer hopping amplitude between the layers of the BLGs, and it is supposed also the same for both sides of the junction. The variable σ , in all terms in Eq. (2), describes the fermionic spin variable, which takes two values σ =↑, ↓. The last term in Eq. (2) subjects the chemical potential terms with the chemical potentials coupled to the charge density operators in each layer of separated BLGs. Initially, we suppose that the BLG/I/BLG system is in the Grand canonical equilibrium state with the equal chemical potentials in both layers of each BLG structure. Thus, we have μ =1 = μ =2 for = L, R. The charge density operator n ,σ (r), as it was discussed above, describes the total electron density with the spin σ in the given layer = 1, 2 of a given BLG with = L, R. For example, for the layer = 1 in the BLG with = L, we have n (1) L,σ (r) = a † σ (r)a σ (r) + b † σ (r)b σ (r). It is important to notice here that we consider purely electronic graphene layers in both sides of the heterojunction and we do not suppose any initial optical pumping in the BLGs.
The fermionic operators, in the individual BLG, systems satisfy the following anticommutation rules The top indices, in Eq. (3), indicates the layers in the BLGs, and the bottom subscripts were introduced for the left-and right BLG, as in Eq. (2). Thus, we suppose the non-vanishing anticommutators only between the fermions on the same type of sublattices in the layers. The sign +, near the parenthesis [. . .] + , in the left-hand side in Eq. (3), indicates the anticommutation rule. We introduce here also the interaction terms by defining the generalized bilayer Hubbard model, in each side of the junction. Namely, we have for the interaction part where H i are the interaction Hamiltonians, for both sides of the junction where n (1) P σ and n (2) Q σ are the electron densities in the layers 1 and 2 in the BLGs ( = L, R) for the P and Q -type fermions, introduced above. The parameter U describes the intralayer Coulomb interactions in the layers of the BLG systems. We have attributed only one index notation to U because we suppose in the paper that U is the same for both layers in the individual BLG systems. Moreover, we put U L = U R = 2γ 0 throughout the paper, except the case of the non-interacting BLGs, when U = W = 0. As we have shown previously [56,57], this parameter does not play an important role in the formation of the excitonic condensates states in the individual BLG systems, and it leads only to the small modifications of the chemical potential behavior as a function of the interlayer Coulomb interaction parameter W . The interlayer Coulomb interaction W in the second term in Eq. (5) is considered as local, i.e., it is defined as the electrostatic Coulomb repulsion between the electrons at the same lattice position in different layers of the BLGs (with = L, R). Indeed, this type of local Coulomb interaction is possible to consider in the BLGs due to their stacking order. We have shown also in Ref. [56,57] that the interaction parameter W is responsible for the excitonic pair formation and excitonic condensates states in the BLG system. It is interesting to emphasize here that at the zero interlayer coupling limit W = 0 between the layers in the individual BLG system and at the zero value of the on-site intralayer Coulomb interaction U = 0 the BLG systems are in the condensates states. This becomes clear when calculating (for details see in Ref. [56,57]) the frequency-summed anomalous excitonic spectral functions g bã (k) (in the left BLG with = L) and g dc (k) (in the right BLG with = R) which give the momentum distribution of the spectral functions in the BLG systems. Indeed, the calculations show the existence of a very narrow but finite bandwidth in the momentum space where the excitonic condensate amplitude is not zero. Therefore, the fundamental excitonic state in the BLG system, at the zero interlayer interaction regime, is governed principally by the electron-hole pair condensed phase caused by the interlayer hopping γ 1 , and the individual free excitonic pairs are indeed absent in the BLG structure. This last effect will become clear in the present paper when calculating the excitonic Josephson current for the non-interaction BLGs, i.e., when U = 0 and W = 0. We should notice here that the interaction parameter W is not defined for example for the electrons on the lattice sites A (in the bottom layer = 1, in the left BLG) andB (in the upper layer = 2, in the right BLG), and, similarly, it is not defined for the electrons on the lattice sites C (in the bottom layer = 1, in the right BLG) andD (in the upper layer = 2, in the right BLG). Indeed, these lattice sites have not their local nearest neighbors sites on the opposite layer in the given BLG. Both interaction parameters U and W are positive because we consider, initially, the electronic layers in both BLG systems. Then, we will use the interaction representation for the fermions [71], in which the time dependence of the fermionic operators is given by the unperturbed Hamiltonian is the tight-binding Hamiltonian of the left side of the junction in the presence of the external gate voltage applied to the BLG/I/BLG (see). We assume that the gate voltage V (t) drops across the barrier and, in general, the Fermi levels in the left-and right-BLG structures will relatively shift by an amount proportional to the potential drop across the junction, i.e.,μ L −μ R = −eV (t). We should emphasize here that the Fermi levelsμ L andμ R are the functions of the shifted effective chemical potentials in the system which depend on the intralayer and interlayer Coulomb interaction parameters introduced in Eq. (5). For the reasons that will be clear in the following sections, we have denoted byμ l , l = L, R the exact Fermi levels in different sides of the tunnel junction. We suppose here the half-filling conditions for the total electron densities in each layer of the separate BLG structure and we assume that the influence of the external gate voltage on the charge densities in the layers is infinitesimally small, and the excitonic gap parameters L = W L ã † i b i (in the left bilayer) and R = W R c † i d i will not get modified by the external perturbation, i.e., δV δn δ , and consequently δn = 0, δ = 0. Such a non-perturbative effect on the electron densities and on the excitonic gap parameter permits to include properly the effect of the applied gate voltage on the excitonic properties in the system and do not affect the excitonic condensate state in the system. For a more sophisticated case that evolves the variation of the excitonic condensates states in the junction, one should include the influence of the gate voltage on the excitonic gap parameter and the half-filling assumption will be failed in this case. Such a treatment is out of the scope of the present paper. We will consider the hopping parameter γ 0 = 1 as the unit of energy measure in the system, and we set e = 1, k B = 1, = 1 in the paper.

The Tunneling Hamiltonian
The time evolution of the eigenstates is determined by the perturbation term, given via the quasiparticle tunneling Hamiltonian H T . In the following, we will transform the total Hamiltonian of the system H = H 0 + H int + H T into the Fourier space by introducing the Nambu fermionic spinors ψ L k,σ (t) and ψ R k,σ (t) in the left and right sides of the junction. We have The tunneling matrix Hamiltonian H T will include all possible quasiparticle tunnelings between the BLGs in the junction. Namely, it accounts all sublattice fermions in the hexagonal layers that participate to the tunneling process. Due to the AB stacking order in the BLGs and the form of the interlayer hopping and interlayer interaction terms in Eqs. (2) and (5), some of the sublattice fermions will not participate to the total excitonic tunneling (for example the fermions on the lattice sites a,b and c,d), and they will contribute only to the single-particle normal tunneling in the junction. Therefore, the tunneling Hamiltonian H T , which describe the tunneling of the particles between the similar lattice sites (here the pairs of similar lattice sites in the BLGs are (Ã,C), (B,D), (A, C) and (B, D), where the first letter in the pair (x, y) is the lattice site in the left BLG and the second letter is the lattice site which belongs to the right BLG) will be written in the form where t k,p is the tunneling matrix element, which describes the transition probability for the electron from the state p (in the right side of the barrier) to a state k (in the left side of the barrier). We will calculate the total tunneling current through the BLG/I/BLG junction by considering the rate of change of the total electron density in the top layer in the right BLG system.
where Ṅ =2 R (t) is the expectation value of the rate of change of the electron number operator N =2 R in the top right layer of the heterojunction, i.e., Indeed, we have considered here the change of the electron number operator in the top right layer of the heterojunction, because of the AB stacking order in the right BLG. The local excitons are created between the electrons, at the sitesÃ in the top right layer, and the holes, at the sites B in the bottom layer (or similarly between the sites C and D in the right BLG system). There are no excitons between the electrons at the sitesB (D) in the top layer and the holes at the sites A (D) in the bottom layer of the left (right) BLG. Thus, considering the zig-zag outermost sites as the sites which are close to the junction contacts with the insulating barrier, the tunneling of the electrons is coming only from the top layer of the right BLG system. Another equivalent geometry for the local excitons in the right BLG is also possible when the electrons at the sites D in the right bottom layer = 1 are coupled with the holes at the sitesC in the top layer = 2. In this case, we should consider the rate of change of the electron number operator N =1 R (t) in the right bottom layer. Those two descriptions are equivalent. We can consider also the tunneling current coming from the rate of change of the hole number operator in the right bottom layer, but, in fact, this contribution is exactly the same as that given by Eq. (8) except the sign in the right-hand side. The sign changing occurs because of the hole type of particles in the right bottom layer. Thus, the rate of total current between the layers will be twice as higher as the tunneling current given in Eq. (8) above. Let us mention also that the consideration of both the electron and hole particle number operator changes in the right BLG has nothing to do with the true excitonic tunneling process, and it defines only the amplitudes of the normal singleparticle and excitonic tunneling currents. In reality, the only electron tunneling from the right top layer is sufficient to consider the single-particle and excitonic tunnelings through the junction because the electron-hole recombination times is much large than the single-particle tunneling time τ rec >> τ T (Here τ T is the single-particle tunneling time.) The electron tunneling time is of order of τ T ∼ 10 −16 sec [72], while the time of the electron and hole recombination processes (Auger recombination time for Coulomb scattering) in graphene is of order of τ rec ∼ 1.1×10 −12 sec for electron-hole densities smaller than 10 12 cm −2 [73]. In order to not double counting the singleparticle and excitonic tunneling effects, we have not considered the rate of change of the left bottom electron number operator N =1 . Thus, as a result, we have the single exciton tunneling through the tunnel junction, from left to right. The expectation value Ṅ R is given by The notation [. . .] − refers to the usual Bose-type commutation rule. A very simple calculation shows that Furthermore, we apply the Rickayzen [74] transformation for the creation and destruction operators in the left-BLG system. Due to the shift of the Fermi levels μ L −μ R = −eV (the similar shit takes the place between the chemical potentials in the usual superconducting Josephson junctions and in the case of the one layer graphene junction arrays, which is not discussed here), the electrons in the left side of the junction gain the additional energy, driven by the external potential and the unperturbed Hamiltonian of the left BLG is given, in this case, as H L (V ) = H L (0) − eV N L . Here, N L is the total electron density operator in the left BLG: N L = =1,2 N L . The equations of motion for the electron in the left BLG at V = 0 will be written as while, for the case V = 0, we have We can suppose that the phase differences across the junction, between the upper and lower layers in the BLGs, evolve with voltage according to the relations known in the usual Josephson junctions theory [75], and a good reason for this is related to the fact that the total ground state energy of the separate non-interacting graphene layer equals exactly twice of the chemical potential in the layer E = 2|μ | with = L, R [76], which is the case of the usual superconductors according to the Gorkov's derivations of the energy spectrum of superconductors in the microscopic theory [77]. Therefore, for the phases in Eq. (14) we have the following relations where Then after integration in both sides in Eq. (15), from zero to t, we obtain Upper (t) and Lower (t). Furthermore, we will define the time-dependent phase difference parameter ϕ(t) as Next, we will consider the tunneling matrix term H T as the small perturbation turned on adiabatically from t = − ∞ and which determines the time evolution of the eigenstates in the Heisenberg picture. Therefore, in the first order in H T , we have The infinitesimal, positive constant η = 0 + has been introduced in the integral in Eq. (17) in order to assure the convergence of the integral. Here, the averages over the commutators are referred to the unperturbed Hamiltonian H 0 = H 0 + H int . The timedependent operators A † (t) and C † (t) in Eq. (17) have been introduced with respect to the reformulation of the tunneling Hamiltonian in Eq. (7) in terms of the composite tunnel-operators A(t), B(t), C(t) and D(t). Namely, we write where It is important to note that the phase factor in Eqs. (17) and (18) appears after the Rickayzen transformations, given in Eq. (14), above.

The Normal Quasiparticle and Excitonic Tunneling
Next, we calculate the commutators in Eq. (17) and we perform the statistical averaging of four-point correlation functions using the Wick's theorem at the finite temperatures [71]. After the whole averaging procedure, we keep only terms responsible for the normal and excitonic currents in the junction. We obtain for the total tunneling current the following expression where the first term, in the sum in the right-hand side in Eq. (20), gives the normal single-particle tunneling current and the second term in Eq. (20) is responsible for the excitonic tunneling. Concerning the excitonic part, we have neglected the term containing the average of type A † (t), B(τ ) − H 0 . This is because the AB-stacked geometry of the BLG structures and the convention that the current (included the excitonic one) is passing from left to right in the junction. We have neglected also all terms in which the following averages appear: Indeed, these mentioned averages are zero because of the symmetry of the Hamiltonian H 0 and the structure of the Nambu spinors in Eq. (6). The phase factor e i(ϕ L −ϕ R ) , in the last term in Eq. (20), is related to the broken symmetry states in the BLG systems in both sides of the junction and which is the consequence of the excitonic pair formation in the system. We have shown previously [56,57] that the individual BLG system undergoes the excitonic condensation phase transition when modulating the interlayer interaction parameter W . In this case, the macroscopic phase coherent state of the excitonic condensate appears in the BLG. Moreover, we showed that in the weak interaction limit, the pair formation in the BLG is BCS-like and a BCS-BEC crossover occurs when augmenting the interlayer interaction parameter. Therefore, as in the case of the usual superconducting junction array [78,79], here also we have the symmetry broken state here related to the local gauge invariance of the complex excitonic order parameter = | |e iϕ . The phase of the order parameter ϕ cannot be changed arbitrarily because of the presence of the macroscopic excitonic condensate state in each side of the junction. Indeed, the measurable physical quantities in the BLG/I/BLG junction arrays depend only on the difference ϕ 0 = ϕ L −ϕ R between the phases of the coherent condensate states in different sides of the junction. Furthermore, we can write Eq. (20) in the form Here, we have performed a change of variable of integration t = t − τ and we have introduced the time-dependent functions N (t) and F(t) for the normal and excitonic counterparts, which are defined with the help of the Kadanoff-Baym Green's functions [80]. For the normal function N (t), we have and for the excitonic function F(t), we get Th function n F (ω) = 1/(e βω +1) ( = L, R) in Eqs. (22) and (23) is the Fermi-Dirac distribution function. We have performed the Fourier transformation of the Kadanoff's Green's functions [80] in order to obtain the expressions of N (t) and F(t), given in Eqs. (22) and (23), and, also, we have introduced the single-particle spectral functions A X (k, ω) (X =ã,b for = L and X =c,d for = R) and excitonic spectral functions A L ab (k, ω), A R cd (k, ω). We have supposed in Eqs. (22) and (23) a simple form of the tunneling probability t k,p = δ k,p t. For a simple treatment, this approximation is sufficient to consider the excitonic effects.
The normal and anomalous (or excitonic) spectral functions in the separate BLGs with the presence of the excitonic pairing interaction have been discussed by us in Ref. [56,57], where the full four-band theory has been developed without low-energy assumption near the K -point in the Brillouin zone. From the form of the total fermionic action, derived there, it follows that the normal spectral functions in different layers in the left BLG are interconnected: The same is true also for the right-BLG system. Moreover, the explicit analytical expressions of the normal spectral functions in the bottom layer of the left BLG are [56,57] and the similar expressions can be written for the right-BLG structure = R. The spectral functions in the upper layers in both BLGs could be obtained after the relations in Eq. (24). The dimensionless coefficients α ik , in Eq. (25) with i = 1, . . . , 4 and = L, R, are given as where P (3) (κ ik ) is the polynomial of third order in κ ik , namely we have with the coefficients ω ik , i = 1, . . . , 3, given as The effective chemical potentials μ 1eff and μ 2eff are defined here as a function of the intralayer U and interlayer W (with = L, R) Coulomb interaction parameters as The energy dispersion parameters κ ik in Eqs. (25) and (26), define the band structures of the interacting BLGs in both sides of the tunnel junction [56,57]. We have The parameters β ik are given as where P (3) (κ ik ) is the polynomial of third order in κ ik , namely we have with the coefficients ω ik , i = 1, . . . , 3, given as As we have mentioned above, the intralayer and interlayer hopping parameters are supposed the same in both sides of the junction. It is particularly important to underline here the role of the bare chemical potentialsμ appearing in the expressions of the band structure parameters in Eq. (30). Indeed, they are playing the role of the exact Fermi energies in the BLG structures as it was pointed out in Ref. [56,57]. For the anomalous (or excitonic) spectral functions in Eq. (23), we have [56,57] and the parameters γ ik are given as Furthermore, the normal and excitonic tunneling currents can be simply expressed analytically after the Werthamer spectral decomposition [81] and by supposing simultaneously the case of the constant gate voltage, i.e., V (t) = V = const. We have and Here, we have used the definition of the phase difference function in Eq. (16). The total tunneling current through the tunnel junction will be where ω 0 is the field frequency: ω 0 = e V , and the normal tunneling is equivalent to the single-particle tunneling term I n (V , T ) = (N (iη + ω 0 )), while the coherent Josephson tunneling of excitons is given by the last two terms in Eq. (38). Thus, the total excitonic tunneling current is and we have I J 1 (V , T ) = (F(iη − ω 0 )) and I J 2 (V , T ) = (F(iη − ω 0 )).
Finally, after some calculations, we get for the normal tunneling current the following expression Here, we have introduced a product-function T (x, y, z), in the following way where each multiplier, in the total product, is given as a single θ -Heaviside step function. Furthermore, the interaction-dependent parameters a 1 , b 1 , c 1 and d 1 , in Eq. (40), are defined as Next, the frequency-dependent parameters ω 0 , ω 0 ,˜ ω 0 and˜ ω 0 , introduced in Eq. (40), form a detuning matrix and for each component we have where the detuning frequencies ω 0 ,ω 0 , ω 0 andω 0 depend explicitly on the difference between the Fermi energies in different sides of the tunnel junction. Namely, we get Next, eight parameters x i with i = 1, . . . , 8, have been introduced in Eq. (40), which are defined as a 1 , b 1 ), x 5,6 = ±X (˜ ω 0 , a 1 , d 1 ), whereas the function X (x, y, z) is defined as X (x, y, z) = (4γ 0 |x|) −1 x 2 − y 2 − z 2 2 − 4y 2 z 2 . The density of states (DOS) function ρ(x) in Eq. (40) appears after transforming the k-summation in Eq. (22) into the integration over the continuous variable, i.e., k . . . = dxρ(x) . . .. The DOS, in the non-interacting graphene layer, is defined as where γ k is the band dispersion in the non-interacting single graphene sheet, i.e., The parameter d, in Eq. (48), refers to the carbon-carbon distance in the graphene layers. Beyond the Dirac's approximation, the DOS can be analytically expressed [58,82,83] as where K(x) is the Elliptic integral of the first kind [84] K(x) = π/2 0 dt/ 1 − x 2 sin 2 t . The function (x), in Eq. (49), is given by [82,83] Furthermore, the functions f (x, y, z) in the denominators in Eq. (40) read as: Concerning the excitonic part of the tunneling current, the Josephson current I J 1 (V , T ), in Eq. (39), is given as: Furthermore, we will calculate numerically the principal value P.V., in Eq. (52), by using the explicit expressions of the excitonic spectral functions, given in Eq. (34), above. For the Josephson current I J 2 (V , T ), we get Again, transforming the summation over the wave vectors into the integration, we obtain for the Josephson current I J 2 (V , T ): In the following section, we will analyze numerically the obtained expressions for the normal and excitonic tunneling currents in the system.

The Normal Quasiparticle Tunneling
In Fig. 2, we have studied the evolution of the normal quasiparticle tunneling current [given in Eq. (40)] as a function of the applied gate voltage V and for two different limits of the right-BLG interlayer Coulomb interaction parameter W R : the weak and intermediate regime, starting from W R = 0 (solid black curve), W R = γ 0 (solid blue curve), W * R = 1.25γ 0 (bold-dashed darker green curve) and W R = 1.5γ 0 (dashed darker yellow curve) (the value W R = 1.5γ 0 is chosen very close to the CNP value W C R = 1.48999γ 0 ) and the high interaction limit, when W R = 1.8γ 0 (dashed darker blue curve), W R = 2γ 0 (solid green curve), W R = 3γ 0 (dot-dashed darker red curve) and W R = 5γ 0 (solid red curve). In all plots, in Fig. 2, we have set the intralayer Coulomb interaction parameters U η equal and U η = 2γ 0 , for both BLGs = L, R and for all layers in the individual BLGs, i.e., η = 1, 2. We denoted by W * R the value of the interlayer Coulomb interaction parameter at which the excitonic gap parameter is maximal: R = max R . The interlayer interaction parameter in the left BLG is fixed at the value W L = 0.5γ 0 , for all curves in Fig. 2. We see first of all, in Fig. 2, that the normal quasiparticle tunneling through the BLG/I/BLG heterojunction, accompanied with the excitonic pair formations in the system, is a threshold process, and the threshold frequency ω 0 = (e/ )V of the external field depends on the relative values of the Coulomb interaction parameters W L and W R at different sides of the construction. We also see that when augmenting the parameter W R in the interval W R ∈ [0, 1.5γ 0 ], the curves, corresponding to the positive part of the current function I n (V , T ), are shifting into left and the intensity of curves is increased. In turn, the threshold values of the normal tunneling current are also shifting left. Starting from the upper bound (UB) critical CNP value of W R (i.e., W R ≥ W C R (U B)), related to the upper bound solution of the chemical potential in the BLG (see in Ref. [56,57]), the normal tunneling current is shifting right, on the V -axis. Nevertheless, the same is not true for the negative part of the tunnel current. We see, in Fig. 2, that all curves of the negative part of the quasiparticle tunneling current are displacing to the right when increasing the interlayer Coulomb interaction parameter in the interval W R ∈ [0, 5γ 0 ]. Only a large jump of the threshold voltages occurs when passing across the lower bound (LB) CNP value W C R (L B), related to the lower bound chemical potential at the CNP. We observe also that for a very large disbalance between the values of the parameters and W R , the additional low-frequency peaks appear in the positive part of the current spectrum, and the threshold frequency values of V are gradually decreasing in these cases. We see also that the amplitudes of the low-frequency peaks are increasing with W R . It is interesting to note that for W R = 5γ 0 (see the solid red curve in Fig. 2), the threshold frequency ω 0 in the positive part of the normal current is of the order of ω 0 ∼ 2γ 1 = 0.256γ 0 = 0.76 eV. Another important observation in Fig. 2 is related to the formation of the fourpeak structures in the spectrum of the normal tunneling current (This is due to the strong excitonic excitations in the four-band structure of the BLGs.) Namely, for small values of W R , the spectrum of the normal current is stepwise, which furthermore transforms into the four-peak structure at the intermediate values of W R (W R = 1.2γ 0 , W * R = 1.25γ 0 , W R = 1.5γ 0 and W R = 1.8γ 0 ). This is more apparent in Fig. 3, where we have chosen very close values of W R , in order to demonstrate the gradual formation of the current four-peak structure.
We see in Fig. 3 that for W R = γ 0 the structure of the tunneling current is halfstepwise with the well-formed two peaks in the excitation spectrum and at the relatively high values of the gate potential. When slightly augmenting W R (see the curves at W R = 1.2γ 0 and W * R = 1.25γ 0 ), the current steps become more pronounced and, at the value W * R = 1.25γ 0 , the resonant tunneling peaks appear at the place of the current steps. Remember that at W * R = 1.25γ 0 the excitonic gap parameter attains its maximum value (see in Ref. [56,57]). Moreover, for higher values of W R , the fourpeak structure remains present and, additionally, the low-frequency resonant tunneling peaks appear in the electron tunneling spectrum. We relate the four-peak structure to the high energy strong resonant tunneling of single electrons, and the presence of the very large tunneling threshold is a direct consequence of the EI state in the bilayer graphenes. Thus, in order to do the tunneling, the electrons must break their contribution to the excitonic insulator state. Contrary, the additional low-frequency peaks at the large values of the parameter W R , are related to the excitonic condensate states in the BLGs. This is the manifestation of tunneling coming from the coherent excitonic condensates states in the system. Vis-a-vis the high interaction values of W R , accompanying the low-frequency peaks, coming from the excitonic condensates states, we can conclude that the excitonic insulator state and the excitonic condensates states in the system are not identical. These are two different states of matter, and the EI state does not survive at the high values of W R , considered in Fig. 2. This observation is in complete agreement with the ideas retrieved in Ref. [56,57]. This statement is also in agreement with the recent work in Ref. [66], where it has been shown that the excitonic Josephson current becomes extremely small before the EI state breaks down. A detailed analysis of the role of the CNP point W C = 1.48999γ 0 , for the case of interaction balanced BLGs, is given in Fig. 4. Particularly, in the upper panel in Fig. 4, we consider the equal values of W L and W R , below the upper bound critical value The latest was also considered in the figure.) We observe that the tunneling spectrum is perfectly symmetric in this case with respect to the origin. When augmenting the interaction parameter up to the upper bound solution at the CNP point, the positive part of the tunneling spectrum is shifted left (see the upper panel in Fig. 4), while for higher values of W R (W R = W L ≥ W C (U B)), the spectrum is shifted to right (This behavior is presented in the middle panel, in Fig. 4). In the lower panel, in Fig. 4, we have presented the smooth passage of the tunneling spectrum when crossing the CNP point. Both, lower bound and upper bound curves of the tunneling current have been considered there. It is worth to mention that the electronic band structure of the BLGs is doubly degenerated at W C because of the chemical potential solutions in the system (see in Ref. [56,57]). We see in the bottom panel in Fig. 4 that the left outermost curve is the UB normal tunneling current in the system. When passing from LB to UB, the tunneling current spectrum still shifted to the left, while for W > W C (U B) the spectrum is transferring to the right. In Fig. 5, we have presented the temperature dependence of the normal tunneling current for a special antisymmetric choice of the Coulomb interaction parameters W L and W R : W L = 0.5γ 0 and W R = γ 0 . As it is clear from the picture, the amplitude of the normal current decreases with increasing the temperature and also leads to the partial suppression of the threshold values of gate voltage, thus promoting a more flexible tunneling of the normal electrons.

The Excitonic Josephson Tunneling
The time dependence of the total excitonic Josephson current through the BLG/I/BLG heterostructure, given in Eq. (39), above, is evaluated numerically in Fig. 6, for different values of the applied gate potential, interaction parameters and condensates phase difference ϕ 0 = ϕ L − ϕ R . When calculating numerically the principal value in Eq. (52), three singular points (− 1; 0; 1) of the integrand have been straddled properly. For the left BLG, we have W L = 0.5γ 0 in all panels in Fig. 6, while for the right BLG we have chosen three different values W R = 0; 0.5γ 0 and W R = γ 0 , from top to bottom panels in Fig. 6. At V = 0 and ϕ 0 = 0, no tunneling current flows in the system (see the solid black lines on the time axis with the holly triangular plotmarkers), while a dc current appears at V = 0 if the phases of coherent excitonic states, in different sides of the junction, differ by ϕ 0 = π/2 or ϕ 0 = −π/2 (see the solid black lines in Fig. 6). Remarkably, the dc excitonic Josephson current changes the sign when changing the sign of ϕ 0 , i.e., ϕ 0 → − ϕ 0 ⇒ I Exc (t) → −I Exc (t). The change of sign in the total excitonic current is related to the fact that the second part of the total excitonic current is zero, i.e., I J 2 (V , T ) cos( ϕ 0 + 2ω 0 t) = 0 because ω 0 = 0 (This is the case V = 0) and ϕ 0 = π/2 (or ϕ 0 = −π/2 when reversing the sign of the condensates phase difference parameter ϕ 0 .) The first part in the total excitonic current at ω 0 = 0 is proportional simply to sin ϕ 0 which is odd function of its argument, and the factor I J 1 (V , T ) is given in Eq. (52), in Sect. 2.3. Such a finite dc current in the system suggests the degeneracy in the ground state of the system, i.e., the U(1) symmetry, and the presence of the excitonic condensates. Moreover, an ac excitonic current appears in the junction for any finite value of V (V = 0.5γ 0 and V = γ 0 , in the picture), even for the case ϕ 0 = 0 (see the blue dotted lines in Fig. 6 with the square plot-markers). The additional phase difference ϕ 0 = π/2 only amplifies the excitonic tunneling current amplitude and leads to a phase shift (see the black dashed lines in Fig. 6). We observe also that the amplitude of the tunnel current decreases when increasing the parameter W R (see the top panel with W R = 0 and the middle panel with W R = 0.5γ 0 , in Fig. 6). In the bottom panel in Fig. 6, we have chosen the larger value of the applied gate potential V = 1.5γ 0 , and W R = γ 0 . We see that the oscillations of the tunnel current are multiplied in this case within the same time interval. This effect is clearly seen in Fig. 7, where different values of the applied gate voltage are considered straightforwardly, and the interaction parameter W R is fixed at the value W R = 0, in correspondence with the plots of the excitonic Josephson current in the upper panel in Fig. 6. We see that when multiplying the value of V by an integer number V = nV , where n = 1, 2, 3, 4, we have for the current wavelength λ = λ/n, thus a relation of type λV = const, between the applied gate potential and the current wavelength, emerges naturally. In Fig. 8, the same function I Exc (t) is presented for the case of the fixed applied gate potential V = 0.5γ 0 , and for different values of the right-BLG interlayer Coulomb interaction parameter W R , below the critical CNP value W C R = 1.48999γ 0 . The condensates phase difference is fixed at the value ϕ 0 = π/2. We see, particularly, that when augmenting the interaction parameter W R (and keeping at the same time W L fixed at W L = 0.5γ 0 ) the amplitude of the excitonic tunneling current is decreasing considerably. The zeros of the current function do not shift in their positions, contrary to the case presented in Fig. 7, where the shift of the zeros is caused by the reduction of the current wavelength in the junction. Thus, we realize that the interlayer Coulomb interaction in the right BLG affects only the current amplitudes, while the changes of the applied gate potential modify principally the frequency of the excitonic current and have not a significant effect on the current amplitudes. Next, in Fig. 9, we have shown the time dependence of the excitonic tunneling current for the values of W R above the UB charge neutrality point W C R (U B). Contrary to the case, given in Fig. 8, the amplitude of the Josephson tunneling current is increasing with W R . This result is related again to the behavior of the chemical potential and Fermi energies in the BLGs (see in Ref. [56,57]). In order to compare the results in Figs. 8 and 9, we kept the curve for W R = 0, in both cases.
In Fig. 10, we have presented the excitonic Josephson current-voltage dependence on time at T = 0 and for the balanced values of the interaction parameters W L and W R . We put W L = W R = 0.5γ 0 , and we have considered the cases when ϕ 0 = 0 or ϕ 0 = 0. We observe, in the upper panel in Fig. 10, that for the case ϕ 0 = π/2 we have a finite excitonic Josephson tunneling current through the system at V = 0. At V = 0, we have principally the same qualitative behavior of the tunneling current for both cases ϕ 0 = 0 and ϕ 0 = π/2 apart the situation at t = 0, when there is a large threshold of V at the zero phase difference (see the solid black curve in the lower panel, in Fig. 10). Contrary, for ϕ 0 = π/2 and at t = 0 the excitonic current The value W R = 0 is also considered as a reference. The parameter W L is fixed at the value W L = 0.5γ 0 , and the condensates phase difference is fixed at ϕ 0 = π/2. Zero temperature case is considered here. The intralayer Coulomb interaction parameters U L and U R are set equal to U = 2γ 0 ( = L, R), for all plots in the figure (Color figure online) has been developed in the system (see the solid black curve, in the upper panel, in Fig. 10).
The I-V characteristics of the excitonic Josephson current, for the case of the equal values of the interlayer interaction parameters, i.e., W L = W R = W , is shown in Fig. 11. Different values of W are considered in the picture: W = 0.5γ 0 (solid black curve), W = 0.8γ 0 (solid blue curve), W = γ 0 (solid yellow curve), W = W * = 1.25γ 0 (solid red curve), W = 1.3γ 0 (dot-dashed darker green curve), W = 1.4γ 0 (large-dashed green curve), and W = W C (L B) = 1.48999γ 0 (bold-dashed darker yellow curve). We observe that up to the value W * = 1.25γ 0 , which corresponds to the maximum of the excitonic gap parameters in the BLGs, the amplitude of the excitonic Josephson current is increasing (including the dc values at V = 0). Furthermore, for W > W * , the amplitudes are continuously decreasing, for W up to the LB CNP value W = W C (L B) = 1.48999γ 0 . The further increase of W , above the LB CNP, leads to a drastic decrease (of about of one order of magnitude) of the excitonic tunnel current amplitude, and this is shown in Fig. 12. It is important to concentrate on another principal difference between the results presented in Figs. 11 and 12. This concerns the first deepest minima of the excitonic Josephson current. In Fig. 11, those minima appear for relatively small voltages, and the positions of the positive (negative) minima are shifting to right (left), for W ≤ W * = 1.25γ 0 . Then, with a further increase of W in the interval W * < W < W C (L B) = 1.48999γ 0 , they are slightly shifting to the left (right). On the other hand, the first minima for the curves in Fig. 12, for W C (U B) ≤ W ≤ 3γ 0 , appear for very large values of V , and the positions of positive (negative) minima are continuously shifting to right (left) in this case. Therefore, the very large values of V are very promising to observe the ac excitonic Josephson current in the system. Experimentally, this could be achieved with the appropriate choice of the left and right gate voltages, which will change the interlayer Coulomb interactions in the BLGs, until the expected effect takes place. In Fig. 13, we have presented the excitonic tunneling current for the asymmetric values of the interaction parameters W L and W R and we see how the increase of W R in the right BLG, above the symmetric value W R = 0.5γ 0 (the parameter W L is fixed at the value W L = 0.5γ 0 for all values of W R ), leads to the tunneling current transfer into right, on the positive axis of V . It is interesting and straightforward to consider the special case of the non-interacting BLGs and the excitonic Josephson current through the junction in this particular case. This result is shown in Fig. 14. We see that although the zero interaction limit, i.e., when U = 0 and W = 0, the excitonic Josephson current still present in the BLG/I/BLG junction, and the values of it are comparable to the case when W > W C (U B), shown in Fig. 12, above. The reason of such excitonic effect is related to the presence of the excitonic condensates states, and the principal mechanism for the formation of the excitonic condensates phase is due to the interlayer hopping between the layer in the individual BLG systems (see in Ref. [56,57]). In the context of the excitonic Josephson current, discussed here, the finite value of the excitonic current follows from the expressions of the functions I J 1 (V , T ) and I J 2 (V , T ), given in Eqs. (50) and (51), in Sect. 2.3. Indeed, the expression of the current amplitude I J 1 (V , T ), in Eq. (52) is proportional to the product of the anomalous spectral functions A L ab (k, ω) and A R cd , given in Eq. (34); in Sect. 2.3, thus, the function I J 1 (V , T ) is proportional to the product (γ 1 + L )(γ 1 + R ). The same is true for the function I J 2 (V , T ), given in Eq. (53). At the zero interlayer coupling limit L = R = 0 while the functions I J 1 (V , T ) and I J 2 (V , T ) remain finite, due to the finite interlayer hopping amplitude γ 1 . Thus, in the non-interacting regime, the fundamental state of BLG systems in both sides of the junction is represented as the phase coherent macroscopic state of the excitonic condensate and the principal mechanism for such a state is attributed to the finite interlayer hopping between the layers in the BLGs. Another spectacular effect, related to the non-interacting regime in the BLGs, is the presence of a finite solution of the chemical potential and the Fermi energy in the separate BLGs (see in Ref. [56,57], for more details). the threshold frequency have been obtained for the case of the interaction-symmetric junction W L = W R , below and above the charge neutrality point. Particularly, from the form of I-V spectrum, in the case of the interaction-asymmetric junction, it is clear that the excitonic insulator state and the excitonic condensate states are two different states of matter in the BLGs. We have shown that the low-frequency peaks appear from the excitonic tunneling in the condensates regime in BLGs, and at the large values of the Coulomb interaction parameter in the right BLG (when the excitonic insulator state breaks down). The temperature dependence of the threshold values of field frequency has been derived for a very large interval of temperature: from zero up to very high temperatures.
Furthermore, the time dependence of the excitonic Josephson tunnel current has been studied in details for different limits of interlayer Coulomb interaction parameter in the right BLG, and the excitonic Josephson dc effect has been found for the case ϕ 0 = 0 and V = 0. Ulteriorly, an ac excitonic Josephson effect has been shown in the junction for nonzero values of the applied gate potential. It has been shown that the phase difference between the excitonic condensates in the BLGs, only amplifies the tunneling current when keeping V constant. It has been shown that a relation of type λV = const is valid in the junction when changing the applied gate potential by the integer values. The time dependence of the excitonic Josephson tunneling current has been analyzed by considering different values of the right-BLG interaction parameter. The I − V characteristics of the excitonic Josephson effect has been found for both interaction-symmetric and interaction-asymmetric cases, and the effects of the charge neutrality point have been discussed in details. It has been shown that the very large values of V are promising to observe the ac excitonic Josephson current in the system. Finally, a particular case of the non-interacting BLGs junction has been considered separately, and the excitonic Josephson tunneling current has been calculated for this case. As a result, we have shown that the nonzero excitonic current is present in this case with the amplitudes comparable to that of the amplitudes of I-V spectrum at the high symmetric values of W . The principal mechanism for the excitonic condensate states and the nonzero excitonic tunneling current, in this case, is attributed to the interlayer hopping between the layers in the individual BLG systems.
We have presented a complete theory of the bilayer graphene-based tunnel junction, and the results, presented here, could represent a veritable framework on which the experimental setup could be made and the results could be compared straightforwardly. The theoretical results in the paper are especially important in the context of the long-standing problem about the excitonic condensation and pair formation in the BLGs and double-layer graphene systems. This is especially important for the theoretical understanding of the nature of the EI state and the coherent excitonic condensate states in such systems: their similarities and differences. The Josephson effect studied in the considered system could be furthermore an important construction in order to build the graphene-based quantum interference devices, which will bring a new idea about the ultra-sensitive magnetometers and voltmeters and will be more sensitive than the usual superconductors-based Quantum Interferometer Devices SQUIDs, due to the exceptional mobility of the electrons in graphene. The excitonic Josephson junctions would be also promising in the context of engineering a new type of ultrafast electronic circuits blocks, which could form a digital logic unit, in the modern ultrafast computers and fast electronics. The results, presented here, could be also important for the building of the new type of adiabatic graphene logic constructions based on the Graphene Quantum Interferometer Devices (Graphene-QUIDs) analog to the Josephson flux transfer logic devices such as the Adiabatic Quantum Flux Parametron (AQFP).
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.