Gate-induced magneto-oscillation phase anomalies in graphene bilayers

The magneto-oscillations in graphene bilayers are studied in the vicinity of the K and K' points of the Brillouin zone within the four-band continuum model ased on the simplest tight-binding approximation involving only the nearest neighbor interactions. The model is employed to construct Landau plots for a variety of carrier concentrations and bias strengths between the graphene planes. The quantum-mechanical and quasiclassical approaches are compared. We found that the quantum magneto-oscillations are only asymptotically periodic and reach the frequencies predicted quasiclassically for high indices of Landau levels. In unbiased bilayers the phase of oscillations is equal to the phase of massive fermions. Anomalous behavior of oscillation phases was found in biased bilayers with broken inversion symmetry. The oscillation frequencies again tend to quasiclassically predicted ones, which are the same for $K$ and $K'$, but the quantum approach yields the gate-tunable corrections to oscillation phases, which differ in sign for K and K'. These valley-dependent phase corrections give rise, instead of a single quasiclassical series of oscillations, to two series with the same frequency but shifted in phase.


I. INTRODUCTION
In solids subject to a magnetic field B, the energy spectrum of charge carriers is quantized into Landau levels (LLs).The magneto-oscillations (MOs) observed in the Shubnikov-de Haas and de Haas-van Alphen effects reflect the oscillations of the density of states (DOS) with the field intensity.The DOS reaches maxima at magnetic fields, B n , for which the LLs with the index, n, cross the Fermi energy, E F .
The Landau plot is a plot of inverse magnetic fields, 1/B n , versus the LL index, n.It is a standard tool used to determine the frequency and phase of MOs, and the related important characteristics of the investigated systems.
The construction of the Landau plot is based on the Onsager-Lifshitz quasiclassical quantization rule, 1,2 A(E F ) = 2π|e|B (n + γ) , where A(E F ) is an area of the extremal cross-section of the Fermi surface (FS) cut by the plane perpendicular to the magnetic field direction, e is the electron charge and γ is a constant which describes the phase of MOs.It follows from Eq. ( 1) that MOs of DOS are periodic in 1/B and their frequency F is related to A(E F ) by The Onsager-Lifshitz quantization rule has been originally designed for three dimensional metals, where the validity of the quasiclassical approximation is guaranteed by a large number of LLs bellow E F in accessible magnetic fields.However, the method is also widely used when two-dimensional (2D) systems are investigated.Here, the importance of F is stressed out by the fact that the carrier concentration is proportional to the area surrounded by the Fermi contour.
In general, the rule should not be applicable to 2D systems.Subject to strong magnetic fields, the quantum limit with only one LL below E F can be easily reached.But in the majority of such systems, the periodicity of MOs is preserved due to the simple parabolic (Schrödinger-like) energy spectra of the 2D electron layers in the semiconductor structures, which yields the LL energies proportional to B.
In 2004 a single sheet of graphene was separated from bulk graphite by micromechanical cleavage. 3It was confirmed experimentally that electrons in graphene obey a linear energy dependence on the wave-vector k, as predicted many years ago by the band structure calculation. 4oth electron and hole charge carriers behave like massless relativistic particles -Dirac fermions (DFs), and there is no gap between the valence and conduction bands.The electron and hole Dirac cones touch at a neutrality point.
Subject to a magnetic field B, the DFs form LLs with energies proportional to √ B. In the seminal papers 3,5,6 the Shubnikov-de Haas MOs in graphene were found periodic in 1/B, similarly as in the 2D gas of Schrödinger fermions (SFs) with the parabolic energy spectra, but with the phase shifted by π.The shift, which was clearly demonstrated by the Landau plot of magneto-resistance oscillations, is due to the existence of the zero-energy LL in the linear Dirac spectrum, shared by electrons and holes.Note that γ = 1/2 for SFs, and γ = 0 for DFs.
In addition to a single layer graphene, also a few layer graphene samples can be prepared.Among them a bilayer graphene (BLG), in which two carbon layers are placed on top of each other with a standard Bernal stacking, is of particular interest.8][9] .][12][13] Note also that the application of the gate voltage is a necessary condition for the experimental observation of MOs in BLG.Without a gate voltage, the sample is neutral, the Fermi energy is located in the neutrality point, and no free charge carriers should be present in perfect samples.
There are two ways of how to apply the gate voltage.If the external voltage is applied symmetrically from both sides of a sample, just E F and the concentration of carriers are varied, and no gap is opened.The tunable gap appears in the presence of external electric field resulting from the asymmetrically applied gate voltage.
Let us point out that the charge carriers in BLG are neither SFs nor DFs, and therefore it is of interest to construct the corresponding Landau plots to see how far the bilayer energy spectra from these two simplest possibilities are.
This task is simplified by the fact that the electrochemical potential (i.e., also E F ) is kept constant during magnetic field sweeps in gated samples.According to Ref. 14 and 15 carrier density oscillations are compensated by gate current oscillations in the case of fixed E F .Note that in bulk samples, where the charge neutrality must be preserved, the carrier concentration is considered to be fixed.
To construct the Landau plot, we will first calculate the quasiclassical frequencies of MOs in BLG, based on their zero-magnetic-field electronic structure.
Later on we will compare these quasiclassical frequencies with results of the quantum-mechanical calculation of the electronic structure of BLG subject to a perpendicular magnetic field.

II. ZERO-FIELD ELECTRONIC STRUCTURE
][18][19][20][21][22] A single layer honeycomb lattice, with two atoms per unit cell, results from two superimposed triangular lattices labeled A and B. The unit cell is defined by the lattice vectors a 1 and a 2 , making the angle 60 • , the lattice constant a is equal to 2.46 Å.The bilayer is formed by two graphene sheets, 1 and 2, arranged in the Bernal stacking.The distance between layers is 3.37 Å.Thus the unit cell of a bilayer has four atoms, its lattice structure is sketched in Fig. 1.
In addition to the intralayer parameter γ 0 and the interlayer parameter t, the corresponding Hamiltonian depends on the potential energy difference between the two layers, which we denote 2u.The parameter γ 0 ≈ 3.1 eV yields the Fermi velocity v F ≈ 1.0 × 10 6 m/s, defined by v F = γ 0 √ 3a/2.We further consider that t ≈ 0.39 eV, and the energy 2u varies between 0 and 250 meV. 23hile γ 0 and t are fixed by nature, we assume that u and E F are the adjustable parameters.
If we employ the continuum approximation, 4 the Hamiltonian H in the vicinity of the K point can be written as where the matrix elements of the first layer are given by Similarly, the matrix elements corresponding to the second layer read There are only two nonzero interlayer matrix elements The Hamiltonian H ′ in the vicinity of the K ′ point has a similar structure, the matrix elements of H ′ are complex conjugates of the matrix elements of H.
The above Hamiltonians can be diagonalized analytically. 16,17,20,22,24The zero-field energy branches of the conduction band, E c1 (k) and E c2 (k), and the valence band, E v1 (k) and E v2 (k), of a bilayer result from hybridization of Fermi cones of layers 1 and 2, mediated by the interlayer matrix element t.Note that and that the valley degeneracy is preserved, i.e., we get the same bands in valleys K and K ′ .
Fermi contours (Color online) The ,,mexican hat" shape of the valence and conduction bands of a biased bilayer.The blue and red colors correspond to higher probability of finding charge carriers in the layers 1 and 2, respectively.Three groups of the Fermi contour are possible depending on the value of EF : the double circles (A1), the circles (A2), and the Fermi rings (A3).
For u = 0 two Fermi cones are replaced by four bonding and antibonding hyperbolic bands.The bonding valence and conduction bands, E v1 (k) and E c1 (k), touch at k = 0, the separation between bands of a bondingantibonding pair is equal to t on the energy scale.
When the interlayer voltage is applied, the Fermi cones of two layers are shifted along the energy axis, and the separation of the neutrality points becomes equal to 2u.The hybridization due to the interlayer parameter t is strongest near the cone cross-points.The resulting four bands are shown in Fig. 2. It turns out that for any finite u a gap is open between the topmost valence band E v1 (k) and the bottom conduction band E c1 (k).The conduction band acquires a ,,mexican hat" shape with energy minima at nonzero k and a local maximum at k = 0. We can write Note that for large k the band E c1 (k) describes electrons localized mostly in the layer 1. Near the local maximum at k = 0 the holes in the layer 2 prevail.Similar conclusions can be drawn for the topmost valence band E v1 (k).
As mentioned in Introduction, the quasiclassical frequencies of the bilayer, F 1 and F 2 , are proportional to areas surrounded by the Fermi circles, which depend, for a given u, on the Fermi energy value E F .Three different possibilities are depicted in Fig. 2 for the case of conduction/valence bands.(For the valence bands E F should be replaced by −E F .) The analytic expressions for the quasiclassical frequencies F 1 and F 2 read the frequencies F 1(2) are even functions of variables E F and u.The frequency F 2 is equal to zero at the local maximum E max c1 (0), and at the minimum E min c2 (0).For a finite u, the frequency F 2 approaches F 1 at E min c1 (k).Three forms of the Fermi contour are possible depending on the value of E F .First, the large E F cuts both conduction bands and F 1 > F 2 > 0. The frequency F 1 corresponds to electron orbits localized mainly in the layer 1, the frequency F 2 corresponds to hole orbits localized mainly in the layer 2. Second, only the band E c1 is cut by E F .Then F 1 > 0 and F 2 < 0. In this case F 2 is just a parameter and does not have the meaning of a true frequency.At last, the E F cuts the bottom conduction band E c1 (k) twice, if it is less than a local energy maximum, E max c1 (0).Then again F 1 > F 2 > 0. In that case F 1 is the frequency of an electron orbit in the layer 1 while F 2 belong to a hole orbit in the layer 2. Close to the local minima the difference between electrons and holes is smeared and charge carriers are present in both layers as indicated by the change of line colors in Fig. 2.
For the special case of u = 0, Eq. ( 5) reduces to Then the gap between the valence and conduction bands as well as the local maximum E max c1 (0) all disappear.The quasiclassical phases of MOs are not accesible via the Onsager-Lifshitz quantization rule, Eqs. ( 1) and (2).To find the energy spectra beyond the quasiclassical approximation, we need to diagonalize the magnetic Hamiltonians H and H ′ .

III. MAGNETIC FIELD EFFECTS
The magnetic Hamiltonians can be obtained from the zero-field Hamiltonians by modification of matrix elements H A1B1 , H B1A1 , H A2B2 and H B2A2 . 16,25,26The matrix elements of the magnetic Hamiltonian in the vicinity of the K point are ).The other matrix elements remain the same as in the zero-field Hamiltonian.Near the K ′ point, We need not diagonalize these Hamiltonians to construct the Landau plot.If we look for magnetic fields B n at which the LLs cross E F , it is enough to find the poles of the resolvent G(z) = (z − H) −1 , as it defines the density of states g(E F ) through The easiest way to find the poles is to solve the corresponding secular equation for B n assuming the fixed E F .We start with the simplest case of the unbiased BLG (u = 0).Then the secular equations can be given a very convenient form, utilizing the quasiclassical frequencies of MOs, presented in the previous paragraph, Eq. ( 6), The Hamiltonians H and H ′ yield identical equations for valleys K and K ′ .While the secular polynomial is quartic in energy it is only quadratic in B. Therefore, to construct the Landau plot it is enough to solve the quadratic equation to find The quasiclassical phase γ can be easily obtained from Eq. ( 8).For a large number of LLs below E F one may assume that n(n + 1) → (n + 1/2) 2 , and then Eq. ( 8) can be written in the form From here we obtain the asymptotic quasiclassical Landau plots i.e., we found that the phases of MOs correspond to SFs with γ = 1/2, in agreement with quasiclassical treatments of systems with inversion symmetry.Note that F 2 is positive only in the rather unrealistic case |E F | > t.
To get Landau plots for an arbitrary n we can express the solution of Eq. ( 8) as we can write (see also Ref. 27) Here the negative sign in the numerator/denominator corresponds to the frequency F 1 in the quasiclassical limit, and the positive sign to the quasiclassical frequency F 2 .It is obvious that for δ = 0 the MOs are not periodic in 1/B.The case of the biased BLG (u = 0) must be treated separately, as the presence of the electric field perpendicular to layer planes breaks the inversion symmetry and lifts the valley degeneracy. 11,12e secular equations can be given again a form quadratic in B, but the coefficients do not depend exclusively on the quasiclassical frequencies as in Eq. ( 8).We can write for the Hamiltonian H in the vicinity of K.In comparison with Eq. ( 8) there is an extra term In the vicinity of K ′ we obtain a very similar equation from the Hamiltonian H ′ , the only difference is that F 0 is replaced by −F 0 .The extra term, F 0 , is the reason of the valley asymmetry.It is obvious that Eq. ( 14) gives two different series of solutions, B n , for positive and negative F 0 .The quasiclassical frequencies F 1 and F 2 are even functions of E F and u.It means that there are the same frequencies not only for K and K ′ , but also for the electrons and holes with energies E F and −E F , respectively.Note also that F 1 and F 2 do not depend on the sign of u.On the other hand, F 0 is an odd function of E F and u.Thus F 0 breaks the K -K ′ symmetry, and also the symmetry between the electron and hole oscillations with the same quasiclassical frequencies.The change of sign of u also reverts the roles of K and K ′ valleys, i.e., what is valid for K with u > 0 is valid for K ′ with u < 0. Also Eq. ( 14) can be rewritten to an equation similar to Eq. ( 13), but with an additional dimensionless parameter λ, which depends on F 0 , Then the analytic solution reads (17) This equation reduces to Eq. ( 13) for λ = 0. Now it is a more difficult task to find an asymptotic expression for the oscillation phase than in the previous case u = 0.If we solve Eq. ( 14) for n + 1/2 we get which for B approaching zero yields Here ξ is a gate-tunable correction to the oscillation phase, given by This correction differs in sign for K and K ′ and also differs for electrons and holes from the same valley with the same absolute value of energy.

IV. RESULTS AND DISCUSSION
In the unbiased BLG the energy u is equal to zero and the parameter δ, which appears in Eq. ( 13), has a particularly simple form For small Fermi energies only the bottom branch E c1 (k) of the conduction subband is cut by E F and only the frequency F 1 is defined.For E F approaching zero, the parameter δ diverges.This implies that for energies close to the band bottom Eq. ( 13) can be written as  the form found for the extremal electron and hole orbits in graphite, 27 which clearly indicates the aperiodicity of oscillations.
The Fermi energies greater than t are rather unrealistic.Nevertheless we can consider this hypothetical case in our theoretical treatment.We can write, for E F = t and δ = 1, Eq. ( 13) as follows The Landau plots calculated for two selected values of E F , E F < t and E F > t, which cross the dispersion curves are presented in Fig. 3.The Landau plots are the same for E F and −E F due to the inversion symmetry conservation.One can observe that in the unbiased bilayer the phases of MOs, corresponding to both frequencies F 1 and F 2 , approach the phase of massive fermions,γ = 1/2, for higher quantum numbers of LLs.
In BLG, an applied electric field leads to asymmetry between K and K ′ valleys that gives rise to nontrivial oscillation phenomena in magnetic fields.To illuminate an anomalous behavior of oscillations, we plotted in Fig. 4 the field dependence of LLs in BLG.
In a single layer graphene the LL fans of electrons and holes start at the zero-field neutrality point.The neutrality points of two independent layers are shifted by 2u and the LLs of holes from the layer 1 cross the LLs of electrons from the layer 2, as shown in Fig. 4 by thin brown lines.
In BLG the shape of LL spectrum results from hybridization of LL spectra of layers 1 and 2. Due to the interlayer interaction, represented by the matrix element t, we have four fans of LLs which start at zero-field energies E v2 (0), E v1 (0), E c1 (0) and E c2 (0).The electron and hole Landau levels (in the K ′ valley) of two layers are mixed by the interlayer interaction, t.For the energy range corresponding to Fermi rings in zero magnetic filed (see A3 in Fig. 2), EF cuts the Landau levels twice.This is the reason for the anomalous phase in the quasiclassical limit B → 0. The hole levels from the layer 1 and the electron levels from the layer 2 avoid to cross, and the low-field hole LLs smoothly turn into the electron LLs as B increases.This is indicated in Fig. 4 by the change of LL color from red to blue.The LLs from a fan starting at zero-field energy E c1 (0) have minima in their field dependence and, therefore, can be cut twice by a single E F .Moreover, the minima are not the same for all levels and, consequently, not all levels are cut by a single E F .This is reflected in the quasiclassical approach as the gate-dependent correction to the MO phase, ξ, which is related to the energy difference 2u between two layers.Note that in the region of energies corresponding to the Fermi rings the expression (20) diverges at E min c1 (k) and is equal to 1/2 for E F = E max c1 (0).As the many body effects can play a role in this low concentration range, the above one-electron picture is probably oversimplified.We start our discussion with the simplest case, E c1/v1 (0) < |E F | < E c2/v2 (0), when only one conduction/valence energy band is cut by E F (see Fig. 2).The single quasiclassical frequency F 1 corresponds to a single Fermi area, which is the same for K and K ′ .
According to Eq. ( 19), the electron peaks in DOS are shifted by ξ to the left in the K valley, whereas the peaks in the K ′ valley are shifted by ξ to the right.The shift magnitude ranges from zero to 1/2 depending on the energy difference between two layers.In Fig. 5 the ,,phase" γ 1 is plotted as the function of n for three different cases with u equal to 0, 0.05 and 0.125 eV.The Fermi energies are chosen to keep the same Fermi area (and the same fixed F 1 ) in all three systems.Only for u = 0 the curves are identical for K and K ′ , for u = 0 the curves are substantially different.
In Fig. 6 the shifts of peaks in the above three cases are shown explicitly.There is a single series of oscillations for the unbiased bilayer, as the valley degeneracy is preserved in a system with the inversion symmetry.
The effect of the gate-tunable valley splitting originates in two series of oscillations which differ for the different u.Let us emphasize that all series of oscillations have the   same quasiclassical frequency F 1 , but the quasiclassical phases depend on the choice of u and on the choice of the valley index.We complete our discussion with cases when the Fermi energy cuts the conduction/valence bands twice.The Landau plots of the biased bilayer with u = 0.125 eV, which is probably the highest experimentally accessible value, 23 calculated for four selected Fermi energies, two in the conduction band, and two in the valence band, are presented in Fig. 7.In accordance with types of the Fermi contours in Fig. 2, the first and second cases are depicted.
The situation is more complicated in the region of energies, ∆ < |E F | < u, for which E F cuts the lowest subband E c1/v1 (k) twice, which is characteristic for the third type of the Fermi contours, as shown in Fig. 2. The bottom of E c1 (k) is at ≈ 0.105 eV.The parameter ξ is far from the values expected for the phase of quasiclassical oscillations, it reduces/grows heavily when E F approaches the bottom of E c1 (k)/E v1 (k).For E F = u, ξ becomes closer to −1/2 for E F in the conduction band and 1/2 for E F in the valence band.

V. CONCLUSIONS
Using a four-band continuum model, we calculated analytically the Landau plots in biased and unbiased BLG subject to external perpendicular magnetic fields.
It turns out that the magneto-oscillations are only asymptotically periodic, and that in the unbiased bilayers their phase is equal to the phase of massive fermions.The convergence to the quasiclassical limit is slow, and depends strongly on the the value of E F .The convergence is slower for higher values of E F .
Anomalous behavior of oscillation phases was found in biased bilayers with broken inversion symmetry.The oscillation frequencies again tend to quasiclassically predicted ones, which are the same for K and K ′ , but the quantum approach yields the gate-tunable corrections to oscillation phases, which differ in sign for K and K ′ .These valley-dependent phase corrections give rise, instead of a single quasiclassical series of oscillations, to two series with the same frequency but shifted in phase.
We also found that for E F in the region of energies corresponding to the Fermi rings in the quasiclassical approach, only a limited number of LLs can cut the Fermi energy and thus a limited number of magneto-oscillations can be achieved.Moreover, their quasiclassical phases reach very large values.As the many body effects can play a role in the corresponding concentration range, the above one-electron picture is probably oversimplified.

FIG. 1 .
FIG. 1. (Color online) Lattice structure of a graphene bilayer.The unit cell is a green parallelepiped.

FIG. 4 .
FIG. 4. (Color online)The electron and hole Landau levels (in the K ′ valley) of two layers are mixed by the interlayer interaction, t.For the energy range corresponding to Fermi rings in zero magnetic filed (see A3 in Fig.2), EF cuts the Landau levels twice.This is the reason for the anomalous phase in the quasiclassical limit B → 0.
FIG. 5. (Color online) The ,,phase" γ1 = F1/B − n calculated for the fixed quasiclassical frequency F1 = 70 T and various u, for the electron K and K ′ valleys as a function of the LL index, n.

FIG. 6 .
FIG. 6. (Color online)The DOS of the unbiased (a) and biased (b, c) BLG versus dimensionless value of the Landau plot, F1/B, with the fixed quasiclassical frequency F1 = 70 T.The frequency F1 corresponds to the situation when F1 > 0, F2 < 0, i.e., only the lowest conduction/valence energy band is cut by EF .In (b, c) the blue peaks show DOS calculated for the K valley, whereas the red ones are related to the K ′ valley.

FIG. 7 .
FIG. 7. (Color online) (a) The electron/hole bands of the biased graphene bilayer with the gap 2u = 0.25 eV at k = 0.The horizontal lines denote the Fermi energies which cross the electron (solid lines) and hole (dashed lines) dispersion curves.(b) The ,,phases" γ 1(2) (EF ) = F 1(2) /B − n calculated for EF depicted in (a) plotted as functions of the LL index, n.