Magnetic Field-Controlled Electrical Conductivity in AA Bilayer Graphene

Abstract

We consider the effect of the external magnetic field on the in-plane conductivity in the AA-stacked bilayer graphene system in the strong excitonic condensate regime. We include the effects of the applied inter-layer electric field and the Coulomb interactions. The on-site and inter-layer Coulomb interactions were treated via the bilayer Hubbard model. Using the solutions for the physical parameters in the system, we calculate the in-plane conductivity of the bilayer graphene. By employing the Green-Kubo formalism for the polarization function in the system, we show that the conductivity in the AA bilayer system is fully controlled by the applied magnetic field. For the partial filling in the layers, the electrical conductivity is different for different spin orientations, and, at the high values of the magnetic field, only one component remains with the given spin orientation. Meanwhile, for the half-filling limit, there is no spin-splitting observed in the conductivity function. The theory evaluated here shows the new possibility for spin-controlled electronic transport in the excitonic bilayer graphene system.

1. Introduction

Recent technological advances in two-dimensional electronic materials bring new insights into the methods and applications of modern nanoscale devices [1]. The AB-stacked bilayer graphene (BLG) is known for its extraordinary property of band gap formation when exposed to an external electric field potential [2,3,4]. In contrast with the AB-BLG systems, the single-particle excitation structure of the AA-stacked BLG has no band gap, and the energy spectrum is linear [4,5]. Nevertheless, a sufficiently large gap in AA bilayers was found in [6], induced by the mass terms via the dielectric medium. The authors in [6] analyzed the transmission and reflection probabilities in the AA bilayer with the layers encapsulated in the dielectric medium. The band gap was found to be of the order of 40 meV. The coexistence of antiferromagnetism with the excitonic ordering [7,8,9] and band gap opening in the doped AA bilayer [10,11] have been found recently, and a bilayer-based spin-valve device has been proposed in [12] when examining the charge-carrier transport dependence on the spin relaxation time. A recent theoretical treatment in [13] showed the possibility of metal–semiconductor transition and excitonic condensation in AA-BLG with an external magnetic field applied. The optical conductivity in AB- and AA-BLG systems has been analyzed in many studies [14,15,16,17,18,19,20,21].

In the present work, we calculate the electrical conductivity in the AA-stacked BLG beyond the Dirac approximation. We consider the system in the presence of the excitonic condensate regime, an external magnetic field, and electric field potential. In addition to the works performed on this subject, we show how the conductivity, for different spin directions, could be controlled by tuning the external magnetic field parameter. We use as a base our previous results in [13], where we showed the possibility of the formation of the strong excitonic condensate state with resulting large band gaps in AA-BLG at the zero temperature limit. We use the Green-Kubo imaginary-time formalism [22] to calculate the in-plane electric conductivity of the AA-BLG in the presence of excitons. Different spin orientations were taken into account. Both partial (with an average fractional number of particles at the given lattice site) and half-filling (with the average number of particles equal to one at each lattice site) regimes were considered. We show that, for sufficiently large values of the magnetic field and at the partial-filling regime, the x-component of the conductivity completely vanishes in one spin channel and remains sufficiently large in another spin channel. Additionally, we show the presence of optical gaps in the conductivity spectrum, which have a strong relationship with the excitonic condensate state in the system, and consider Coulomb interactions in the system. The results given here could open new possibilities for spin-controlled electronic transport applications of the AA-BLG system and allow its consideration as a new type of spin-valve device.

The paper is organized as follows: In Section 2, we introduce the bilayer Hubbard model and we discuss the Coulomb terms in it. In Section 3, we use the Kubo formalism for the electric conductivity concerning our system and we calculate it numerically for different values of the on-site filling coefficient. In Section 4, we give the results of our numerical calculations. Furthermore, in Section 5, we discuss the obtained results, and in Section 6, we give a short conclusion to our paper. Appendix A is devoted to the calculation of the polarization function in the system.

2. Generalization of Hubbard Model for AA Bilayer Graphene

We studied the electrical conductivity in the AA-BLG system with the help of the bilayer Hubbard model. The Hamiltonian of our system has the following form:

$$\begin{array}{ccc}& & {\widehat{\mathcal{H}}}_{\mathrm{AA}}={\widehat{\mathcal{H}}}_{\mathrm{t}}+{\widehat{\mathcal{H}}}_{\mathrm{int}},\hfill \end{array}$$

(1)

where ${\widehat{\mathcal{H}}}_{\mathrm{t}}$ is the tight binding part of the total Hamiltonian and ${\widehat{\mathcal{H}}}_{\mathrm{int}}$ is the interaction part. The Hamiltonian ${\widehat{\mathcal{H}}}_{\mathrm{t}}$ is given in terms of the intra-layer (${\gamma }_0$) and inter-layer (${\gamma }_1$) hopping amplitudes and

$$\begin{array}{ccc}& & {\widehat{\mathcal{H}}}_{\mathrm{t}}=-{\gamma }_0\sum _{〈\mathbf{r},{\mathbf{r}}^{\prime }〉}\sum _{\sigma }\left({\widehat{a}}_{\sigma }^{†}\left(\mathbf{r}\right){\widehat{b}}_{\sigma }\left({\mathbf{r}}^{\prime }\right)+h.c.\right)\hfill \\ & & -{\gamma }_0\sum _{〈\mathbf{r},{\mathbf{r}}^{\prime }〉}\sum _{\sigma }\left({\widehat{\tilde{a}}}_{\sigma }^{†}\left(\mathbf{r}\right){\widehat{\tilde{b}}}_{\sigma }\left({\mathbf{r}}^{\prime }\right)+h.c.\right)\hfill \\ & & -{\gamma }_1\sum _{\mathbf{r}\sigma }\left({\widehat{a}}_{\sigma }^{†}\left(\mathbf{r}\right){\widehat{\tilde{a}}}_{\sigma }\left(\mathbf{r}\right)+{\widehat{b}}_{\sigma }^{†}\left(\mathbf{r}\right){\widehat{\tilde{b}}}_{\sigma }\left(\mathbf{r}\right)+h.c.\right)-\mu \sum _{\mathbf{r}}\widehat{n}\left(\mathbf{r}\right).\hfill \end{array}$$

(2)

Here, the operators ${\widehat{x}}_{\sigma }\left(\mathbf{r}\right)$, ${\widehat{\tilde{x}}}_{\sigma }\left(\mathbf{r}\right)$ and ${\widehat{x}}_{\sigma }^{†}\left(\mathbf{r}\right),{\widehat{\tilde{x}}}_{\sigma }^{†}\left(\mathbf{r}\right)$ with $x=a,b$ are the electron annihilation and creation operators attached to different sublattice sites in the layers (see Figure 1). The spin $\sigma$ has two possible directions: $\sigma =\uparrow ,\downarrow$. The parameter ${\gamma }_0$ describes the hopping of electrons between the adjacent lattice sites, and ${\gamma }_1$ is the local hopping of electrons between the layers in the AA-BLG. Here, we put the hopping parameters equal to ${\gamma }_0=3$ eV and ${\gamma }_1=0.384$ eV, which are consistent with the values reported in [23]. The summation $〈\dots 〉$ in the first two terms in Equation (2) is over the nearest neighbor lattice sites $\mathbf{r}$, ${\mathbf{r}}^{\prime }$. The last term in Equation (2) describes the coupling of the chemical potential with the total electron density operator $\widehat{n}\left(\mathbf{r}\right)$, which is given as

$$\begin{array}{c}\hfill \widehat{n}\left(\mathbf{r}\right)={\widehat{n}}_1\left(\mathbf{r}\right)+{\widehat{n}}_2\left(\mathbf{r}\right),\end{array}$$

(3)

where the electron density operators ${\widehat{n}}_1\left(\mathbf{r}\right)$ and ${\widehat{n}}_2\left(\mathbf{r}\right)$ are defined as

$$\begin{array}{c}\hfill {\widehat{n}}_1\left(\mathbf{r}\right)=\sum _{x=a,b}{\widehat{n}}_x\left(\mathbf{r}\right),\\ \hfill {\widehat{n}}_2\left(\mathbf{r}\right)=\sum _{\tilde{x}=\tilde{a},\tilde{b}}{\widehat{n}}_{\tilde{x}}\left(\mathbf{r}\right).\end{array}$$

(4)

The interaction part of the Hamiltonian includes the on-site U-Hubbard interaction terms in the layers, the local inter-layer Coulomb coupling, and the interaction with the external electric field potential V and the magnetic field B. We have

$$\begin{array}{c}\hfill {\widehat{\mathcal{H}}}_{\mathrm{int}}=U\sum _{\mathbf{r}\eta }{\widehat{n}}_{\eta \uparrow }\left(\mathbf{r}\right){\widehat{n}}_{\eta \downarrow }\left(\mathbf{r}\right)+W\sum _{\mathbf{r}\sigma {\sigma }^{\prime }}{\widehat{n}}_{a\sigma }\left(\mathbf{r}\right){\widehat{n}}_{\tilde{a}{\sigma }^{\prime }}\left(\mathbf{r}\right)\\ \hfill +W\sum _{\mathbf{r}\sigma {\sigma }^{\prime }}{\widehat{n}}_{b\sigma }\left(\mathbf{r}\right){\widehat{n}}_{\tilde{b}{\sigma }^{\prime }}\left(\mathbf{r}\right)+\frac{V}{2}\sum _{\mathbf{r}}\left({\widehat{n}}_2\left(\mathbf{r}\right)-{\widehat{n}}_1\left(\mathbf{r}\right)\right)\\ \hfill -g{\mu }_B{B}_z\sum _{\mathbf{r}\eta }\left({\widehat{n}}_{\eta \uparrow }-{\widehat{n}}_{\eta \downarrow }\right).\end{array}$$

(5)

Here, the parameters U and W describe the intra-layer and inter-layer Coulomb interactions in the system. The parameter g, in Equation (5), is the Landé g-factor [24], and ${\mu }_B$ is the Bohr magneton (here, we use the conventions ${\mu }_B=1$ and $\hslash =1$). The magnetic field B is considered in the z-direction, perpendicular to the layers of the bilayer (see Figure 1). The summation index $\eta$, in the last term in Equation (5), is over the sublattice site variables a, b (for the layer $\ell =1$) and $\tilde{a},\tilde{b}$ (for the layer $\ell =2$).

The biquadratic fermionic terms in the Coulomb interaction parts of the interaction Hamiltonian in Equation (5) could be linearized via the Hubbard–Stratanovich decoupling procedure. This detail was reported in [13].

The effective Hamiltonian, after those decoupling procedures, is given as

$$\begin{array}{ccc}& & {\widehat{\mathcal{H}}}_{\mathrm{eff}}=-{\gamma }_0\sum _{\mathbf{r}\delta }\sum _{\sigma }\left(a_{\sigma }^{†}\left(\mathbf{r}\right)b_{\sigma }(\mathbf{r}+\delta )e^{\frac{ie}{\hslash c}{\int }_{\mathbf{r}}^{\mathbf{r}+\delta }A\left({\mathbf{r}}^{\prime }\right)d{\mathbf{r}}^{\prime }}+h.c.\right)\hfill \\ & & -{\gamma }_0\sum _{\mathbf{r}\delta }\sum _{\sigma }\left({\tilde{a}}_{\sigma }^{†}\left(\mathbf{r}\right){\tilde{b}}_{\sigma }(\mathbf{r}+\delta )e^{\frac{ie}{\hslash c}{\int }_{\mathbf{r}}^{\mathbf{r}+\delta }A\left({\mathbf{r}}^{\prime }\right)d{\mathbf{r}}^{\prime }}+h.c.\right)\hfill \\ & & -\sum _{\mathbf{r}\sigma }\left({\Delta }_{\sigma }+{\gamma }_1\right)\left(a_{\sigma }^{†}\left(\mathbf{r}\right){\tilde{a}}_{\sigma }\left(\mathbf{r}\right)+b_{\sigma }^{†}\left(\mathbf{r}\right){\tilde{b}}_{\sigma }\left(\mathbf{r}\right)+h.c.\right)\hfill \\ & & -\sum _{\mathbf{r}\sigma }\left(\mu -g{\mu }_B\mathbf{B}\sigma +\frac{V}{2}-\frac{U}{2}{\overline{n}}_a-2W\right)n_{a\sigma }\left(\mathbf{r}\right)\hfill \\ & & -\sum _{\mathbf{r}\sigma }\left(\mu -g{\mu }_B\mathbf{B}\sigma +\frac{V}{2}-\frac{U}{2}{\overline{n}}_b-2W\right)n_{b\sigma }\left(\mathbf{r}\right)\hfill \\ & & -\sum _{\mathbf{r}\sigma }\left(\mu -g{\mu }_B\mathbf{B}\sigma -\frac{V}{2}-\frac{U}{2}{\overline{n}}_{\tilde{a}}-2W\right)n_{\tilde{a}\sigma }\left(\mathbf{r}\right)\hfill \\ & & -\sum _{\mathbf{r}\sigma }\left(\mu -g{\mu }_B\mathbf{B}\sigma -\frac{V}{2}-\frac{U}{2}{\overline{n}}_{\tilde{b}}-2W\right)n_{\tilde{b}\sigma }\left(\mathbf{r}\right).\hfill \end{array}$$

(6)

We included in Equation (6) the vector potential $\mathbf{A}\left(\mathbf{r}\right)$ in order to consider the electric current response of the system. The vector potential is taken into account with the help of the Peierls–Onsager substitution [25,26]. We have

$$\begin{array}{c}\hfill \mathbf{B}\sigma =\left\{\begin{array}{cc}{\displaystyle }& B,\mathrm{if}\sigma =\uparrow ,\\ \\ \\ \\ {\displaystyle }& -B,\mathrm{if}\sigma =\downarrow .\end{array}\right.\end{array}$$

(7)

The parameter ${\Delta }_{\sigma }$, introduced in Equation (6), is the excitonic order parameter, defined as

$$\begin{array}{c}\hfill {\Delta }_{\sigma }=W〈a_{\sigma }^{†}{\tilde{a}}_{\sigma }〉=W〈b_{\sigma }^{†}{\tilde{b}}_{\sigma }〉,\end{array}$$

(8)

where the parenthesis $〈\dots 〉$ means the mean-field average of the product of two fermionic operators (see [13] for details). The averages ${\overline{n}}_{\eta }$, in Equation (6), (with $\eta =a,\tilde{a},b,\tilde{b}$) are over sublattice density operators, and come from the decoupling of the intra-layer on-site Coulomb interaction terms in the interaction part in Equation (5). We assume that the average electron concentrations on different sublattice sites in the layers are equal, i.e., ${\overline{n}}_a={\overline{n}}_b$ and ${\overline{n}}_{\tilde{a}}={\overline{n}}_{\tilde{b}}$. Additionally, we impose the condition for the occupation (particle filling) between adjacent lattice sites in different layers, i.e., ${\overline{n}}_a+{\overline{n}}_{\tilde{a}}=1/\kappa$. Thus, the half filling corresponds to $\kappa =0.5$ when a maximum of one particle occupies a given lattice site in the given layer.

3. The Green-Kubo Formalism and Electrical Conductivity in AA Bilayer Graphene

3.1. The Electric Current Operator beyond Dirac Approximation

Furthermore, by supposing the small variations of the vector potential (at the distances of the order of $\delta$) and expanding the exponential in the tight binding part of the Hamiltonian in Equation (6) up to the first order of the vector potential $\mathbf{A}\left(\mathbf{r}\right)$, we obtain, for the Hamiltonian,

$$\begin{array}{c}\hfill {\widehat{\mathcal{H}}}_{\mathrm{eff}}={\widehat{\mathcal{H}}}_{\mathrm{eff}}(\mathbf{A}\left(\mathbf{r}\right)=0)-\frac{1}{c}\sum _{\mathbf{r},i}\sum _{\ell =1,2}{A}_{\ell i}\left(\mathbf{r}\right)j_{\ell i}\left(\mathbf{r}\right),\end{array}$$

(9)

where $j_{\ell i}\left(\mathbf{r}\right)$ is the component of the electric current operator along the direction i. We can obtain the form of the in-plane operator $j_{\ell x}\left(\mathbf{r}\right)$ by functional differentiation of the Hamiltonian in Equation (9) with respect to the vector potential $A_{\ell x}\left(\mathbf{r}\right)$. We have, for the total current density operator in the AA-BLG, the following expression:

$$\begin{array}{c}\hfill j_x=-c\frac{\delta {\widehat{\mathcal{H}}}_{\mathrm{eff}}}{\delta A_{\ell x}\left(\mathbf{r}\right)}.\end{array}$$

(10)

After the Fourier transformation of the fermionic creation and annihilation operators into the reciprocal space, we obtain, for the current density operator,

$$\begin{array}{c}\hfill j_x=-\frac{e}{N}\sum _{\mathbf{k}\sigma }\left(v_{\mathbf{k}x}{a}_{\mathbf{k}\sigma }^{†}{b}_{\mathbf{k}\sigma }+v_{\mathbf{k}x}^{*}{b}_{\mathbf{k}\sigma }^{†}{a}_{\mathbf{k}\sigma }\right.\\ \hfill \left.+v_{\mathbf{k}x}{\tilde{a}}_{\mathbf{k}\sigma }^{†}{\tilde{b}}_{\mathbf{k}\sigma }+v_{\mathbf{k}x}^{*}{\tilde{b}}_{\mathbf{k}\sigma }^{†}{\tilde{a}}_{\mathbf{k}\sigma }\right),\end{array}$$

(11)

where we have introduced the velocity operator $v_{\mathbf{k}x}$

$$\begin{array}{c}\hfill v_{\mathbf{k}x}=-\frac{i{\gamma }_0}{\hslash }\sum _{\delta }{\delta }_x{e}^{i\mathbf{k}\delta }.\end{array}$$

(12)

It is easy to calculate explicitly the x-component of the velocity operator, given in Equation (12). Indeed, for the nearest neighbor vectors $\delta$ (which are the same for the AA-stacking configuration), we have the following expressions

$$\begin{array}{c}\hfill \delta =\left\{\begin{array}{cc}{\displaystyle }& \hfill {\delta }_1=\left(\frac{a}{2\sqrt{3}},\frac{a}{2}\right),\hfill \\ \\ \\ \\ {\displaystyle }& \hfill {\delta }_2=\left(\frac{a}{2\sqrt{3}},-\frac{a}{2}\right),\hfill \\ \\ \\ \\ {\displaystyle }& \hfill {\delta }_3=\left(-\frac{a}{\sqrt{3}},0\right),\hfill \end{array}\right.\end{array}$$

(13)

where $a=\sqrt{3}{a}_0$, in Equation (13), is the sublattice constant, while $a_0$ is the carbon–carbon length in the graphene layers (with $a_0=1.42$). Then, for the velocity operator, we obtain

$$\begin{array}{c}\hfill v_{\mathbf{k}x}=-\frac{i{\gamma }_{0}a}{\hslash }\left(e^{i\frac{k_{x}a}{2}}cos\frac{k_{y}a\sqrt{3}}{2}-e^{-i{k}_{x}a}\right).\end{array}$$

(14)

Hereafter, we use this expression for the velocity operator when calculating the electrical polarization function beyond the Dirac approximation.

3.2. The Polarization Function and Electrical Conductivity

The electronic polarization function can be calculated within the Matsubara imaginary time formalism [22] in the following way

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{ij}\left(i{\omega }_m\right)=-{\int }_0^{\beta }d\tau e^{i{\omega }_m\tau }〈T_{\tau }{j}_i\left(\tau \right)j_j\left(0\right)〉,\end{array}$$

(15)

where the indexes $i,j=x,y,z$ denote the components of the polarization operator. The frequencies ${\omega }_m$ are the Bosonic Matsubara frequencies, with ${\omega }_m=2m\pi /\beta$ (with $m=0$, $\pm 1,\pm 2,\dots$), and $\beta$ is $\beta =1/T$. In turn, the real part of the conductivity function is related to the imaginary part of the retarded polarization function, as

$$\begin{array}{c}\hfill \Re {\sigma }_{ij}\left(\mathrm{\Omega }\right)=\frac{\Im {\mathrm{\Pi }}_{ij}^R\left(\mathrm{\Omega }\right)}{\mathrm{\Omega }}.\end{array}$$

(16)

The retarded polarization function in Equation (16) could be obtained from that in Equation (15) via analytical continuation into the positive real frequency axis

$$\begin{array}{c}\hfill {\mathrm{\Pi }}_{ij}^R\left(\mathrm{\Omega }\right)={\mathrm{\Pi }}_{ij}(i{\omega }_m\to \mathrm{\Omega }+i{ϵ}^{+}).\end{array}$$

(17)

The average of the time-ordered product in Equation (15) could be evaluated using the Wick theorem [27]. Then, we obtain, for the x-component of the polarization function in Equation (15), the expression in terms of the fermionic and excitonic Green’s functions (the details of calculations and definitions of Green’s functions are given in Appendix A).

$$\begin{array}{ccc}& & {\mathrm{\Pi }}_{\mathrm{xx}}\left(i{\omega }_m\right)=-\frac{e^2}{\beta N}\sum _{\mathbf{k}{\nu }_n}\sum _{\sigma }{\left|v_{\mathbf{k}x}\right|}^2\times \hfill \\ & & \times \left({\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{aa}}(i{\omega }_m+i{\nu }_n){\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{bb}}\left(i{\nu }_n\right)\right.\hfill \\ & & +{\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{aa}}\left(i{\nu }_n\right){\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{bb}}(i{\omega }_m+i{\nu }_n)\hfill \\ & & +{\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{a}}\tilde{\mathrm{a}}}(i{\omega }_m+i{\nu }_n){\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{b}}\tilde{\mathrm{b}}}\left(i{\nu }_n\right)\hfill \\ & & +{\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{a}}\tilde{\mathrm{a}}}\left(i{\nu }_n\right){\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{b}}\tilde{\mathrm{b}}}(i{\omega }_m+i{\nu }_n)\hfill \\ & & +2{\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{a}\tilde{\mathrm{a}}}\left(i{\nu }_n\right){\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{b}}\mathrm{b}}(i{\omega }_m+i{\nu }_n)\hfill \\ & & \left.+2{\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{a}}\mathrm{a}}(i{\omega }_m+i{\nu }_n){\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{b}\tilde{\mathrm{b}}}\left(i{\nu }_n\right)\right),\hfill \end{array}$$

(18)

where ${\mathcal{G}}_{\mathbf{k}\sigma }^{\eta \eta }\left(i{\nu }_n\right)$ is the normal fermionic Green’s function and ${\mathcal{G}}_{\mathbf{k}\sigma }^{\eta \tilde{\eta }}\left(i{\nu }_n\right)$ is the excitonic Green’s function (see Appendix A). The frequencies ${\nu }_n$ in Equation (18) are the fermionic Matsubara freqeuncies ${\nu }_n=\pi (2n+1)/\beta$, with $n=0,\pm 1,\pm 2,\dots$. We present here their explicit forms following the work in [13]. We have

$$\begin{array}{c}\hfill {\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{aa}}\left(i{\nu }_n\right)={\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{bb}}\left(i{\nu }_n\right)=-\sum _{i=1}^4\frac{{\alpha }_{i\mathbf{k}}^{\left(\sigma \right)}}{i{\nu }_n+{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)}},\\ \hfill {\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{a}}\tilde{\mathrm{a}}}\left(i{\nu }_n\right)={\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{b}}\tilde{\mathrm{b}}}\left(i{\nu }_n\right)=-\sum _{i=1}^4\frac{{\beta }_{i\mathbf{k}}^{\left(\sigma \right)}}{i{\nu }_n+{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)}},\\ \hfill {\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{a}\tilde{\mathrm{a}}}\left(i{\nu }_n\right)={\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{b}}\mathrm{b}}\left(i{\nu }_n\right)={\mathcal{G}}_{\mathbf{k}\sigma }^{\tilde{\mathrm{a}}\mathrm{a}}\left(i{\nu }_n\right)={\mathcal{G}}_{\mathbf{k}\sigma }^{\mathrm{b}\tilde{\mathrm{b}}}\left(i{\nu }_n\right)=\\ \hfill =-\sum _{i=1}^4\frac{{\gamma }_{i\mathbf{k}}^{\left(\sigma \right)}}{i{\nu }_n+{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)}}.\end{array}$$

(19)

Here, the coefficients ${\alpha }_{i\mathbf{k}}^{\left(\sigma \right)}$, ${\beta }_{i\mathbf{k}}^{\left(\sigma \right)}$, and ${\gamma }_{i\mathbf{k}}^{\left(\sigma \right)}$ in the right-hand side in Equation (19) are given in Appendix A. Furthermore, we reinserted the expressions of Green’s functions in Equation (19) into Equation (18) and we performed the summation over the Fermionic Matsubara frequencies ${\nu }_n$. Then, we obtained, for the imaginary part of the polarization function, the following expression

$$\begin{array}{ccc}& & \Im {\mathrm{\Pi }}_{\mathrm{xx}}\left(\mathrm{\Omega }\right)=\frac{\pi e^2}{N}\sum _{\mathbf{k}\sigma }\sum _{i,j=1,\dots ,4}{\left|v_{\mathbf{k}x}\right|}^2\times \hfill \\ & & \times \left[2{\alpha }_{i\mathbf{k}}^{\left(\sigma \right)}{\alpha }_{j\mathbf{k}}^{\left(\sigma \right)}+2{\beta }_{i\mathbf{k}}^{\left(\sigma \right)}{\beta }_{j\mathbf{k}}^{\left(\sigma \right)}+4{\gamma }_{i\mathbf{k}}^{\left(\sigma \right)}{\gamma }_{j\mathbf{k}}^{\left(\sigma \right)}\right]\times \hfill \\ \hfill \times & & \delta (\mathrm{\Omega }-{\epsilon }_{j\mathbf{k}}^{\left(\sigma \right)}+{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)})\left(n_F(-{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)})-n_F(-{\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)}-\mathrm{\Omega })\right).\hfill \end{array}$$

(20)

Here, $n_F\left(x\right)$ is the Fermi–Dirac distribution function defined as $n_F\left(x\right)=1/(e^{\beta (x-\mu )}+1)$. The energy parameters ${\epsilon }_{i\mathbf{k}}^{\left(\sigma \right)}$ with the energy branches $i=1,\dots ,4$ define the electronic band structure in the AA-BLG with the excitonic pairing interaction. They are given and calculated in [13] for different values of the magnetic field parameter. Due to different physical solutions of the excitonic order parameter ${\Delta }_{\sigma }$, corresponding to different spin orientations, we expected different results for the imaginary part of the polarization function. Therefore, we separated

$$\begin{array}{c}\hfill \Im {\mathrm{\Pi }}_{\mathrm{xx}}\left(\mathrm{\Omega }\right)=\Im {\mathrm{\Pi }}_{\mathrm{xx}}^{\uparrow }\left(\mathrm{\Omega }\right)+\Im {\mathrm{\Pi }}_{\mathrm{xx}}^{\downarrow }\left(\mathrm{\Omega }\right)\end{array}$$

(21)

and we calculated both components separately. In the same manner, the real part of the electrical conductivity function $\Re {\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$ is

$$\begin{array}{c}\hfill \Re {\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)=\Re {\sigma }_{\mathrm{xx}}^{\uparrow }\left(\mathrm{\Omega }\right)+\Re {\sigma }_{\mathrm{xx}}^{\downarrow }\left(\mathrm{\Omega }\right).\end{array}$$

(22)

In the next section, we give the numerical results for the function in Equation (22), and we discuss different on-site filling regimes.

4. Results

In Figure 2 and Figure 3, we presented the numerical results for the energy dependence of the real part of the conductivity function normalized to the half of the conductance quantum ${\sigma }_0=e^2/\hslash$. On the axis of abscissa, we put the excitation energy, normalized to the intra-layer hopping energy ${\gamma }_0$, i.e., $E=\hslash \mathrm{\Omega }/{\gamma }_0$, where $\mathrm{\Omega }$ could be the frequency of the external radiating photons coming from the light source devices. The zero temperature limit was considered in the figures.

Both components of the sum, in the right-hand part of Equation (22), were calculated separately to show the difference related to the orientation of the spin. In Figure 2, the case of partial filling is considered with the inverse filling coefficient $\kappa =1$ (see [13]). When evaluating the $\mathbf{k}$-sum in Equation (20), the $\delta$-Dirac function was replaced by the Lorentzian function with the small broadening coefficient $\eta =0.05$. Different values of the normalized magnetic field parameter $\tilde{B}={\mu }_{B}B/{\gamma }_0$ are considered in panels (a)–(c) in Figure 2. The inter-layer Coulomb interaction parameter was set at the value $W=2{\gamma }_0$, and the electric field potential was fixed at the value $V=2{\gamma }_0$. In panel (a) in Figure 2, we consider the case of zero magnetic fields. The electrical conductivities coincide for both spin channels. The plot in black corresponds to the large intra-layer Coulomb interaction with $U=4{\gamma }_0$. The conductivity peak at the zero frequency limit $\mathrm{\Omega }\to 0$ corresponds to the Drude limit (when no interband transitions occur), and the width of the peak is equal to the inverse relaxation time of electrons. We see that the region where the conductivity function is non-zero is displaced to the right on the frequency axis for the small value of the intra-layer Coulomb interaction parameter $U={\gamma }_0$. Moreover, an optical gap appears in the conductivity spectrum for the case $U={\gamma }_0$ (see the low-frequency region in the red plot) with the value ${\Delta }_{\mathrm{opt}}={\gamma }_0$. The values of the conductivity function are the same for both spin orientations, i.e., ${\Re }^{\uparrow }{\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)={\Re }^{\downarrow }{\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$. The case of the non-zero magnetic field with $\tilde{B}=0.5$ is considered in panel (b) in Figure 2. Both large and small values of the parameter U are considered here. We see that the electrical conductivity is different in this case for different orientations of the spin. For $U=4{\gamma }_0$ the conductivity peaks, corresponding to the spin orientations $\sigma =\uparrow$ and $\sigma =\downarrow$, are separated on the $\mathrm{\Omega }$-axis and the conductivity ${\Re }^{\uparrow }{\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$ is placed in the high-frequency region, while the frequencies corresponding to the non-zero values of ${\Re }^{\downarrow }{\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$ are smaller. The optical gaps are present for both spin orientations, with ${\Delta }_{\mathrm{opt}}^{\uparrow }=2{\gamma }_0=6$ eV (corresponding to the case $\sigma =\uparrow$) and ${\Delta }_{\mathrm{opt}}^{\downarrow }=0.8{\gamma }_0=2.4$ eV (corresponding to the case $\sigma =\downarrow$). Meanwhile, we see that for $U={\gamma }_0$, the peaks of the conductivity function for $\sigma =\uparrow$ are displaced to the right and the optical gap is larger, in this case, ${\Delta }_{\mathrm{opt}}^{\uparrow }=2.8{\gamma }_0=8.4$ eV. The conductivity function for $\sigma =\downarrow$ is very different in this case (see the curve in green, in panel (b) in Figure 2). We can observe that the optical gap is absent in the case $\sigma =\downarrow$ and the conductivity maximum is displaced again to the right on the $\mathrm{\Omega }$-axis. Thus, by inducing the finite magnetic field, the conductivities corresponding to different spin orientations are completely different and get completely separated on the $\mathrm{\Omega }$-axis. The low-frequency region on $\mathrm{\Omega }$ corresponds to the conductivities with $\sigma =\downarrow$ and the high-frequency region corresponds to the conductivities with $\sigma =\uparrow$. Experimentally, this means that the electrical conductivity in the magnetic field $\tilde{B}=0.5$ and with $\sigma =\downarrow$ could be detected with the photons’ radiation with wavelengths $\lambda \in (137.76,6200)$ nm, i.e., covering the light spectrum from the near ultraviolet (UV) up to the far infrared regions (also covering the visible light spectrum). Meanwhile, the conductivity, corresponding to $\sigma =\uparrow$, could be observed in the presence of radiation with $\lambda \in (68,137.76)$ nm, covering the light spectrum from the soft X-ray up to the near UV regions.

A totally different picture appears for the higher values of the magnetic field parameter $\tilde{B}=1$. In this case, the conductivity becomes totally suppressed for $U={\gamma }_0$ and the results for the case of $U=4{\gamma }_0$ are given in panel (c) in Figure 2. We can observe in panel (c), in Figure 2, that in this case only the conductivity with $\sigma =\downarrow$ remains and $\Re {\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$ vanishes. The optical gap is also absent in this case. The results for the case of the half filling (with the inverse filling parameter $\kappa =0.5$) are presented in Figure 3 with the same values of the physical parameters in the system as in Figure 2. We can observe in panels (a)–(c), in Figure 3, that ${\Re }^{\uparrow }\left({\sigma }_{\mathrm{xx}}\right)\left(\mathrm{\Omega }\right)={\Re }^{\downarrow }{\sigma }_{\mathrm{xx}}\left(\mathrm{\Omega }\right)$ for all values of the intra-layer Coulomb interaction parameter U. We can see in panel (a), that for $\tilde{B}=0$ and $U=4{\gamma }_0$, the conductivity is suppressed apart from in the Drude part of the spectrum (see the black curve in the very low-frequency region). Meanwhile, the conductivity function for the case $U={\gamma }_0$ generates very large values on a very wide interval of frequencies. The conductivity functions at $\tilde{B}=0.5$ are shown in panel (b) in Figure 3, for $U=4{\gamma }_0$ and $U={\gamma }_0$. We can observe that the conductivity function for $U={\gamma }_0$ has practically the same shape as for the case of zero magnetic fields. Meanwhile, the large domain of conductivity appears for $U=4{\gamma }_0$ and $\tilde{B}=0.5$ (see the black curve in panel (b) in Figure 3). Furthermore, with a higher value of the magnetic field parameter equal to $\tilde{B}=1$, the amplitudes of the conductivity functions decrease drastically (see panel (c) in Figure 3), and also very large optical gaps appear in the conductivity spectrum, which transfer the spectrum to the large frequency region. For $U={\gamma }_0$, the optical gap is of the order of ${\Delta }_{\mathrm{opt}}^{\sigma }=1.75{\gamma }_0=5.25$ eV, and for $U=4{\gamma }_0$, we obtain ${\Delta }_{\mathrm{opt}}^{\sigma }=2.45{\gamma }_0=7.35$ eV. We can relate the obtained high values of the optical gap parameters to the presence of the strong excitonic condensate in the system (which has been proven to be possible in [13]) and the considered limits of the Coulomb interaction parameters in the AA bilayer graphene.

As an example, we also considered in Figure 4 a very high value for the magnetic field, equal to $\tilde{B}=2$, and a case of partial filling in the layers with $\kappa =1$. The obtained result is very similar to that obtained for the case of $\tilde{B}=1$ (see panel (c) in Figure 2).

5. Discussion

Indeed, the results of our work are based, fundamentally, on the initial calculations that started in [13], where the important physical parameters, such as the excitonic gap parameter, average charge density imbalance, and chemical potential, were obtained self-consistently within the framework of the same Hamiltonian that was used in the present work. In this sense, in general, it could be slightly hard to capture the complete background on which the results here are based. For this reason, it is perhaps recommended to read the above-mentioned paper.

Concerning the numerical values of the interaction parameters, given in the Hamiltonian in Equation (1) and used in the present work, it is useful to notice that the most realistic values of the intra-layer (U) Coulomb interaction parameters are given in the intervals 5–7 eV [4]. In practice, the on-site Coulomb interaction parameter U is the strongest interaction constant. The local inter-layer Coulomb interaction W is, in general, smaller than U. In order to measure the interaction constant W in units of U, we can choose the approximate formula $W/U=a_0/d_0$, where $a_0$ is the lattice constant in the layers and $d_0$ is the inter-layer separation distance. By the way, all mentioned estimations are especially true for the pristine AB-staked bilayer graphene. For the artificially constructed AA-BLG systems, we can change the inter-layer separation distance $d_0$ experimentally and obtain different interaction regimes (which are desired by the authors). In this sense, the artificially built AA bilayers (which consist of two separate graphene monolayers and are not naturally obtained bilayers with the given fixed value of the inter-layer distance) are more purposeful to verify the results presented in the frames of our theory. Moreover, our BLG system consists of such artificial constructions (for this reason we sometimes call them “double-layer” systems and not “bilayer”) with variable inter-layer distances. The parameters V and B could be changed experimentally or arbitrarily, within the reasonable bandwidths (see, in [13], about the relation of B and V).

To our knowledge, the experimental results for U and W were only known (until now) for the AB-stacked BLG structure (see [28] or in [29]). The authors introduced a fundamental emphasis on the control of the charge carrier concentration (by changing the doping level in the layers), which permits the measurements of parameters, such as U and ${\gamma }_1$. They synthesized the BLG structure on the silicon carbide (SiC) substrate and measured the electronic properties using angle-resolved photoemission spectroscopy. The most important parameter, which merits measurement, is the local Hubbard interaction U for different carrier concentrations. For the AA-BLG, we calculated, numerically [13], the exact average carrier concentration difference $\delta \overline{n}$ (charge-imbalance) between layers. which gives fundamental information about the electron or hole doping in the individual layers. The nature of the electron or hole type of the layers was determined by the sign of the parameter $\delta \overline{n}$. If $\delta \overline{n}>0$, then the upper layer in BLG was populated by the electrons, while the lower layer was populated by the holes. For $\delta \overline{n}<0$, the holes are dominating in the upper layer and the electrons in the lower layer. The other parameter W (especially in the case of AA-stacked BLG) could be measured just by varying the inter-layer distance in the bilayer construction with the simultaneous ARPES measurements of the electron concentrations in the layers. Furthermore, the results for W could be compared with those obtained for U, and then the energy scales of W could be extracted subsequently. On the other hand, the charge density modulations due to the electric field and in the example of the AB-BLG structures have been measured recently via low-temperature Raman spectroscopy techniques [30]. Particularly, the authors analyzed and measured experimentally the charge density non-uniformity of the BLG caused by phonon anharmonicity.

Indeed, very good estimates for the magnetic fields, needed to employ the spin-valve regime on several spin transport devices, are given in a series of works in [31,32,33] and in the references therein. Particularly, the order of magnitudes of the magnetic fields at which the suitable spin-valve effect takes place are very small and of the order of mT. Nevertheless, our estimations for the magnetic field values at which the AA-BLG behaves like a spin-valve device are presented in panels (a)–(c) in Figure 2, and in Figure 4 for the case of the partial-filling regime. We should remark that the spin-valve effect takes place in our BLG at the partial-filling regime and at sufficiently high values of the external magnetic field parameter, at the order of 1.5–6 T. It is notable that the spin-valve effect observed here occurs at the partial-filling regime, which is also mentioned in each panel in Figure 2 and Figure 3. The partial-filling regime deviates considerably from the half-filling regime. In the case of partial filling, the average number of particles (electrons or holes) at the single atomic lattice site position is not one (as in the case of half filling) but could be fractional and less than one. The most prominent effect of the spin-split conductivity was attained at the value of the dimensionless magnetic field parameter $\tilde{B}=0.5$ and is presented in panel (b) in Figure 2. This regime corresponds to the real magnetic field with the value $B={\gamma }_0\tilde{B}=1.5$ T, which is completely achievable in experimental situations. Moreover, such realistic values of the magnetic field have been discussed in many works, such as in [16,34,35].

It is remarkable here to note that the spin-valve effect does not take place in the case of the half filling. Indeed, as shown in Figure 3, there is no spin-splitting of the electrical conductivity function; nevertheless, the behavior of the conductivity function was different at the small and large U limits. This does not occur at the high value of the magnetic field, and at the half filling we found the opening of the large optical gap in the system [36], which is due to the existence of the excitonic condensate state in the AA-BLG [13].

6. Conclusions

In this paper, we considered the bilayer Hubbard model to study the electronic conductivity of AA-stacked bilayer graphene in the strong excitonic condensate regime. The applied electric field potential and the magnetic field were considered. The effect of the in-plane electric field was taken into account in the form of the Peierls phase variables in the hopping terms of the initial Hamiltonian of the system. The Green-Kubo approach was used to calculate the electrical conductivity function in the system with the presence of excitons at the zero temperature limit. Both partial and half-filling cases were considered, and we used the numerical results obtained in [13] for the physical parameters in the system, such as the chemical potential, the average charge density difference between the layers, and the excitonic order parameters for different spin orientations $\sigma =\uparrow$ and $\sigma =\downarrow$. We numerically calculated the conductivity functions for both spin channels and different values of the magnetic field parameters and interaction parameters. For the partial-filling case and when changing the magnetic field, the conductivity function obtained different values for different spin orientations. Particularly, for $\sigma =\uparrow$, the conductivity function vanished at the high values of the magnetic field, while for $\sigma =\downarrow$ it remained finite for all values of the applied magnetic field. Contrarily, for the half-filling case, the conductivity functions nearly coincided for different spin orientations and at small and intermediate values of the applied magnetic field. Moreover, in this case, the obtained large values of the optical gap parameters were strongly related to the presence of an excitonic condensate state in the system.

The results obtained in this paper, show how one can control the electrical conductivity function in AA-stacked bilayer graphene by changing the magnetic field. After an appropriate tuning of the interaction parameters, one can obtain the electrical conductivity in the system with the required spin orientations. The results obtained here could be important for spintronic and optoelectronic applications of the AA bilayer graphene and give new insight into the possible applications of the AA bilayer graphene as a new type of spin-valve device.