Conductivity scaling and the effects of symmetry-breaking terms in bilayer graphene Hamiltonian

We study the ballistic conductivity of bilayer graphene in the presence of symmetry-breaking terms in effective Hamiltonian for low-energy excitations, such as the trigonal-warping term ($\gamma_3$), the electron-hole symmetry breaking interlayer hopping ($\gamma_4$), and the staggered potential ($\delta_{AB}$). Earlier, it was shown that for $\gamma_3\neq{}0$, in the absence of remaining symmetry-breaking terms (i.e., $\gamma_4=\delta_{AB}=0$), the conductivity ($\sigma$) approaches the value of $3\sigma_0$ for the system size $L\rightarrow{}\infty$ (with $\sigma_0=8e^2/(\pi{}h)$ being the result in the absence of trigonal warping, $\gamma_3=0$). We demonstrate that $\gamma_4\neq{}0$ leads to the divergent conductivity if $\gamma_3\neq{}0$, or to the vanishing conductivity if $\gamma_3=0$. For realistic values of the tight-binding model parameters, $\gamma_3=0.3\,$eV, $\gamma_4=0.15\,$eV (and $\delta_{AB}=0$), the conductivity values are in the range of $\sigma/\sigma_0\approx{}4-5$ for $100\,$nm$\<L<1\,\mu$m, in agreement with existing experimental results. The staggered potential ($\delta_{AB}\neq{}0$) suppresses zero-temperature transport, leading to $\sigma\rightarrow{}0$ for $L\rightarrow{}\infty$. Although $\sigma=\sigma(L)$ is no longer universal, the Fano factor approaches the pseudodiffusive value ($F\rightarrow{}1/3$ for $L\rightarrow{}\infty$) in any case with non-vanishing $\sigma$ (otherwise, $F\rightarrow{}1$) signaling the transport is ruled by evanescent waves. Temperature effects are briefly discussed in terms of a phenomenological model for staggered potential $\delta_{AB}=\delta_{AB}(T)$ showing that, for $0<T\leqslant{}T_c\approx{}12\,$K and $\delta_{AB}(0)=1.5\,$meV, $\sigma(L)$ is noticeably affected by $T$ for $L\gtrsim{}100\,$nm.


I. INTRODUCTION
The universal conductivity of monolayer graphene (MLG), σ MLG = 4e 2 /(π h) (with the elementary charge e and the Planck constant h), accompanied by the pseudodiffusive shot noise (quantified by the Fano factor F = 1/3), is one of the most recognizable landmarks of the Dirac nature of electrons that dwell in this material [1][2][3][4][5]. These unique characteristics are linked to the dominant role of transport via evanescent waves in graphene near the charge-neutrality point [6]. What is more, the effective Hamiltonian for low-energy excitations, where v F = √ 3 t 0 a/(2h) ≈ 10 6 m/s is the energyindependent Fermi velocity (with t 0 ≈ 3 eV the nearestneighbor hopping integral and a = 0.246 nm the lattice spacing), p j = −ih∂ j are in-plane momentum operators, and σ j are the Pauli matrices acting on the sublattice degree of freedom (with j = x, y), possesses several symmetries which are crucial for the simplicity of transport properties. (We further note that the absence of valley-coupling factors is supposed throughout the paper, and the discussion is limited to the K valley.) These include the rotational invariance (RI), the electron-hole symmetry (EHS), and the sublattice equivalence (SE), which is embedded in the so-called symplectic symmetry (or time-reversal symmetry in a single valley) [7,8].
In bilayer graphene (BLG) the situation is more complex due to the couplings between the layers [9][10][11][12][13]. Historically, the effective Hamiltonians for BLG were constructed by taking only the leading tight-binding parameters of the Slonczewski-Weiss-McClure model [14,15], which are indicated in Fig. 1(a).
Even in the simplest possible approach [9,10], including the nearest-neighbor interlayer hopping γ 0 (being numerically different than t 0 ) and the direct interlayer hopping γ 1 , SE is already eliminated due to the inequivalence of sites connected by γ 1 (dimer sites) and the remaining ones (nondimer sites), giving an opportunity to open the band gap by perpendicular electric field introducing the layer inequivalence [16]. (The second-nearest-neighbor interlayer hopping, formally breaking EHS, is usually omitted as-in the low-energy limit-it only shifts the charge-neutrality point by a constant value; see Ref. [17].) Quite surprisingly, the approach of Refs. [9,10] leads to the conductivity σ 0 = 2σ MLG = 8e 2 /(π h) and F = 1/3 (in the absence of a gap), as one could expect for two decoupled layers. The results are also size independent, provided that W L l ⊥ , with the sample width W , the length L marked in Fig. 1(b), and l ⊥ = √ 3 aγ 0 /(2γ 1 ) ≈ 1.77 nm being a new length scale due the coupling between the layers.
Next, skew-interlayer hopping (or the trigonal-warping term) γ 3 [18] breaks RI, leading to the appearance of three additional Dirac cones at each valley [16]. The effect of γ 3 on quantum transport is also significant [11][12][13]; namely, the conductivity σ (L) is no longer universal but length dependent, approaching the value of 3σ 0 for large L [19]. In contrast, the Fano factor is unaffected (i.e., F = 1/3), showing that the pseudodiffusive nature of charge transport in BLG cannot be attributed any particular value of σ . In this paper, we complement the previous studies of ballistic charge transport in BLG by examining numerically the effect of EHS-breaking interlayer hopping γ 4 on the σ (and F ) dependence on L. The results show that for γ 4 = 0, σ (L) may be either divergent (for γ 3 = 0) or vanishing (for γ 3 = 0) with L → ∞ (with F ≈ 1/3 in the first case or F → 1 in the second case), as marked schematically in Fig. 1(c). These findings extend the collection of possible behaviors associated with transport via evanescent waves in graphene-based systems. The role of an intrinsic (i.e., not related to the external electric field but rather interaction induced) band gap reported by some experimental works [20][21][22][23] (and parametrized here by the staggered potential δ AB ) is also discussed.
The remaining part of the paper is organized as follows. In Sec. II we present the model Hamiltonian and discuss how each of the symmetry-breaking terms (γ 3 , γ 4 , or δ AB ) affects the low-energy dispersion relation. Then, in Sec. III we demonstrate, by means of numerical mode matching for the Dirac equation, the behavior of σ and F with growing L separately in the presence and in the absence of each symmetry breaking. The concluding remarks are given in Sec. IV.

II. THE MODEL
We start from the minimal version of the four-band Hamiltonian [16], in which all the symmetry breakings mentioned in Sec. I are quantified by independent parameters: where π = e −iθ (p x + ip y ), π † = e iθ (p x − ip y ), with the angle θ (between an armchair direction and the x axis) defining the crystallographic orientation of the sample, v 0 = √ 3aγ 0 /(2h), v 3 = v 0 γ 3 /γ 0 , and v 4 = v 0 γ 4 /γ 0 . In the forthcoming numerical discussion, we set θ = π/4, γ 0 = 3.16 eV, and γ 1 = 0.381 eV [24]; for each of the remaining parameters the cases of zero and nonzero value are studied independently to demonstrate the impact of a particular symmetry breaking on ballistic transport. Namely, we took γ 3 = 0 or 0.3 eV, γ 4 = 0 or 0.15 eV, and δ AB = 0 or 1.5 meV.
Our specific choice of the staggered potential δ AB in the Hamiltonian H BLG (2) follows from the demand that it opens a band gap without breaking EHS, which is solely controlled by γ 4 . (In the parametrization of Ref. [16] the energy difference between dimer and nondimer sites V also breaks EHS; here we set V = 0). Physically, δ AB represents the irreducible part of a gap (i.e., one that cannot be closed by external electric fields) and can be attributed to charge or spin order which may appear in the BLG ground state when electron-electron repulsive interactions are taken into account [25,26].
In Fig. 2 we present the low-energy band structure following from the Hamiltonian H BLG (2) by displaying the cross sections, for p y =hk y = 0, of dispersion relations for eight different combinations of symmetry-breaking parameters γ 3 , γ 4 , and δ AB . An apparent feature visible in Fig. 2(a) is the energy shift of a secondary Dirac cone (same for all three secondary cones) due to EHS breaking for γ 4 = 0 [see Fig. 2(c) for a comparison], making it impossible (for γ 3 = 0 and δ AB = 0) to achieve the exact zero-doping case, in which the transport is fully carried by evanescent waves. The consequences of these features for BLG transport properties are discussed next.
being the k x position of a secondary Dirac cone calculated for the parameters as listed in (c). Each panel displays the cross section taken at k y = 0.

A. Zero-temperature charge transport
We employ the Landauer-Büttiker expressions for zerotemperature conductivity and the Fano factor in the linearresponse regime [27], namely, with the conductance quantum g 0 = 4e 2 /h accounting for spin and valley degeneracies. In order to determine the transmission matrix at a given Fermi energy, t = t(E ), for a rectangular sample attached to the two heavily doped regions, we employ the computational scheme similar to that presented in Ref. [13], with a numerical stabilization introduced in Ref. [28]. In brief, at finite-precision arithmetics, the modematching equations may become ill defined for sufficiently large L, as they contain both exponentially growing and exponentially decaying coefficients. This difficulty is overcome by dividing the sample area into N div consecutive, equally long parts and matching the wave functions for all (i.e., N div + 1) interfaces. Typically, using the double-precision arithmetic, we put N div = L/(40 l ⊥ ) + 1, with x denoting the floor of x.
Our numerical results for E = 0 are presented in Figs. 3 and 4. As the debate on the ground-state nature of BLG is currently ongoing [26]  The behavior of transport properties is relatively simple for δ AB = 1.5 meV (coinciding with the gap reported in Refs. [22,23]); we observe a fast decay of σ (L) with growing L, accompanied by F → 1 (see the blue lines in Figs. 3 and 4), indicating the insulating (or semiconducting) behavior. The remaining parameters (γ 3 and γ 4 ) are essentially meaningless in such a case; a slightly elevated conductivity (namely, σ > σ 0 ) is visible for γ 3 = 0.3 eV and L < 100 l ⊥ , due to the finite-size effects.
In a gapless case (δ AB = 0) we identify four apparently different behaviors of σ (L), depending on whether the remaining parameters (γ 3 and γ 4 ) take zero-or nonzero values.
where m = γ 1 /2v 2 0 and k l = γ 1 v 3 /hv 2 0 , leading to a nonzero number of propagating modes (open channels) at zero energy, which can be approximated as [28] Subsequently, the excess conductivity from secondary Dirac cones can roughly be bounded as . The transmission reduction for propagating modes, approximately by a factor of 2, can be attributed to the additional backscattering appearing in the double-barrier geometry, which is usually much weaker for a single barrier [31]. A secondary feature of σ (L) is a quasiperiodic oscillation due to the Fabry-Pérot resonances appearing for L = n L, where n = 1, 2, . . . , and (up to the order of magnitude) L ∼ πhv 3 / E l ≈ 340 l ⊥ [28]. None of these effects is present for γ 4 = 0, for which the conductivity follows the scenario earlier described in Refs. [12,13]. (In Appendix A, we present the analytical derivation explaining why σ (L) → 3σ 0 for L → ∞ and arbitrarily small γ 3 = 0.) We further notice that the available experimental value in Ref. [29], reporting σ ≈ 2.5σ 0 for L ≈ 400 nm = 226l ⊥ [red circle in Fig. 3(b)], is equally close to both the results for γ 4 = 0 and 0.15 eV, and the determination of γ 4 via conductivity measurements requires a sample length exceeding L 2 μm.
For δ AB = γ 3 = 0, the conductivity behavior with growing L is a bit more peculiar.
If In contrast, if γ 4 = 0.15 eV (solid black lines) we observe a slow power-law decay of σ (L) with growing L, which can be approximated as σ (L) ∝ L −2.0 for L 1000 l ⊥ , accompanied by F → 1. Notice that the Fano factor is F ≈ 1/3 in the range of L 300 l ⊥ shown in Fig. 4(d); the convergence to 1 becomes visible for L 1000 l ⊥ [see Fig. 4(c)]. Unlike for nonrelativistic electrons [32], we still obtain a finite σ (L) in the limit of infinite doping in the leads at fixed W and L, signaling the relativistic nature of charge carriers. The vanishing conductivity for L → ∞ at a fixed W/L, in the absence of a gap, clearly represents a remarkable feature of the results, providing an opportunity to verify the γ 3 = 0 model as put forward in Ref. [21] within ballistic transport experiments. A further reasoning that such a behavior appears generically for γ 3 = 0 and γ 4 = 0 is given in Appendix B.

B. Finite-temperature effects
For T > 0 and in the linear-response regime the electronic noise is dominated by the Nyquist-Johnson term of S(0) ≈ 4k B T σW/L [27], and the Fano factor becomes irrelevant. Therefore, we limit our discussion to the temperaturedependent conductivity, which is given by being sufficiently high to reach a convergence up to the machine roundoff errors. Additionally, when calculating the transmission matrix t(E ), we parametrize the staggered potential in the effective Hamiltonian H BLG (2) as follows: with T C = 12 K and δ AB (0) = 1.5 meV reproducing the temperature dependence of a gap reported in Refs. [22,23]. (The gapless case δ AB (0) = 0 is considered separately.) Our numerical results, for T = 0 and the selected temperatures 0 < T T C , are presented in Fig. 5. Similarly as in the previous section, the data sets for γ 3  An apparent feature of the results presented in Fig. 5(a) is that the curves for different temperatures closely follow each other up to L 120 l ⊥ , above which σ (L) grows noticeably faster with L for higher T . [Notice that the T = 0 curve also shows approximately linear growths with L, which manifests itself for L 10 3 l ⊥ ; see Fig. 3(a).] The position of a coalescence point, L ≈ 120 l ⊥ , can be attributed to fact that above such a length the quantum-size effects are less significant, allowing the finite-temperature effects to dominate transport properties. This can be rationalized taking into account the time-energy uncertainty relation limiting the energy resolution, with τ flight ≈ L/v 3 being the ballistic time of flight [28], together with the fact that the energy of thermal excitations is k B T E l for T 4 K, with E l given by Eq. (5). Subsequently, one can expect that the propagating modes in secondary Dirac cones are employed (by thermal excitations) provided that δE E l k B T , leading to For γ 3 = 0 the above reasoning no longer applies; however, a relatively flat σ dependence on L for T > 0 [see Fig. 5(c)] coincides with the divergent lower bound for L in Eq. (11). In such a case, one should rather estimate the time of flight (up to the order of magnitude) as τ flight ∼ L/v 0 . In turn, the condition k B T δE allowing the conductivity enhancement by thermal excitations is equivalent to giving, for instance, L 440 l ⊥ for T = 5 K. The lower bound for L in Eq. (12) allows one to understand why temperature effects on σ (L) are noticeably weakened for γ 3 = 0, comparing to the γ = 0 case. In the presence of a staggered potential, δ AB (0) = 1.5 meV, the primary temperature effects on σ (L) visible in Figs. 5(b) and 5(d) can be attributed to the gap closing for T approaching T C [see Eq. (9)]. The characteristic system length, above which the value of δ AB becomes significant, derived from the condition for the energy uncertainty δE δ AB (0), reads where we have estimated . This time, our numerical results show that the temperature effects are visible for significantly shorter systems in the γ 3 = 0.3 eV case, comparing to the γ 3 = 0 case, in a qualitative agreement with the estimation given in Eq. (13).

IV. CONCLUDING REMARKS
We have investigated, by calculating ballistic transport characteristics within the Landauer-Büttiker formalism, the role of symmetry-breaking terms in the effective Hamiltonian for bilayer graphene. Three such terms, the trigonal warping (γ 3 ), the electron-hole symmetry-breaking interlayer hopping (γ 4 ), and the staggered potential (δ AB ) quantifying a spontaneous band gap, are independently switched on and off, resulting in different behaviors of the conductivity (σ ) and the Fano factor (F ) with the increasing system length (L) at a fixed width-to-length ratio (W/L).
In the presence of a staggered potential at T = 0 (δ AB (0) > 0), a semiconducting behavior is observed regardless of the remaining parameters (γ 3 and γ 4 ); i.e., σ (L) → 0 (showing the exponential decay) and F → 1 for L → ∞. For T > 0, a zero-gap behavior is gradually restored (for any combination of γ 3 and γ 4 ) when the energy of thermal excitations k B T δ AB (0).
We hope that our numerical results will help verify the bilayer graphene models proposed in the literature, as soon as ballistic samples of the length L 1 μm become available. So far, conductivity measurements for shorter samples [29] suggest that the models neglecting the trigonal warping (γ 3 = 0) cannot correctly reproduce transport properties in the mesoscopic range, but conclusive information concerning the value of γ 4 is missing.
Apart from the material-science aspects outlined above, the asymptotic conductivity behavior suggests that bilayer graphene represents a model case when discussing the generality of spontaneous symmetry breaking in quantum systems [33,34]. When σ is considered as an order parameter, our findings can be summarized by putting forward the noncommuting order of limits, as L → ∞ and the relevant symmetry breakings vanish, namely, where we have introduced the distance between the layers d (with d → ∞ corresponding to simultaneous limits of γ 3 → 0 and γ 4 → 0 [35]), and the dots [· · · ] in Eq. (16) mark that the order of the two remaining limits is arbitrary in such a case. From this perspective, it becomes clear that both the sublattice and the combined rotationalelectron-hole symmetry breakings may appear spontaneously, as consequences of the layer stacking in graphene (d = const < ∞).
The peculiar cases of γ 3 = 0 or γ 4 = 0 in the absence of other symmetry breakings (i.e., δ AB = γ 4 = 0 or δ AB = γ 3 = 0) do not seem to have as clear a physical interpretation. However, in heterostructures containing graphene, a variety of spontaneous symmetry breakings may appear due to the couplings to surrounding layers, encouraging one to consider also anomalous parameter configurations.
Our considerations are limited to Bernal-stacked bilayer graphene, since other stacking usually opens the band gap leading to the vanishing conductivity. One notable exception is the magic-angle twisted bilayer graphene [36][37][38][39][40] showing a gapless band structure. However, the role of electron-electron interaction in such a system is more significant; in particular, the ground state near the chargeneutrality point is a correlated insulator [38][39][40]. Therefore, we expect the vanishing zero-energy conductivity to appear generically in bilayer graphene for non-Bernal stackings. Here, we present the analytical derivation of the total transmission (i.e., transmission probability summed over normal modes), coinciding with the Landauer-Büttiker conductivity [see Eq. (3) in the main text] σ (L) → 3σ 0 in the limit of L, W → ∞, at W/L = const 1. Some partial results were earlier reported in Ref. [12], but the full derivation, to our best knowledge, is missing in the literature.

ACKNOWLEDGMENTS
The dispersion relation for the Hamiltonian given by Eq. (2) in the main text, with δ AB = γ 4 = 0, takes the form where p = √ p 2 x + p 2 y and we have set θ = 0 for simplicity (later, we show that the physical results are independent of the lattice orientation in the L → ∞ limit).
In the vicinity of zero energy (|E | → 0), there are four solutions of the above equation corresponding to four Dirac cones: the central cone, located at p = (p x , p y ) = (0, 0), and three satellite cones, located (in polar coordinates) at p = γ 1 v 3 /v 2 0 , ϕ = 0, 2π/3, 4π/3. Below, we calculate the transmission of the system assuming that the states corresponding to different Dirac cones do not interfere among themselves. Physically, such a supposition corresponds to the conditions for the energy and system sizes where the Lifshitz energy E L = 1 4 γ 1 (v 3 /v 0 ) 2 . For γ 0 = 3.16 eV, γ 1 = 0.381 eV, and γ 3 = 0.3 eV [24], we have E L ≈ 1 meV, and the last two conditions in Eq. (A3) are equivalent to L, W 4ł ⊥ γ 3 /γ 0 = 75 nm. Expanding the dispersion relation given by Eqs. (A3) and (A4) up to second order around p = (0, 0), we obtain Thus, the central Dirac cone has an isotropic dispersion relation, closely resembling the dispersion relation following from the monolayer graphene Hamiltonian [see Eq. (1) in the main text]; in fact, the only difference is the proportionality coefficient v 3 instead of v F . Now, we write the effective single-cone Hamiltonian, corresponding to the dispersion relation given by Eq. (A4): Solving the scattering problem for a rectangular sample described by the above Hamiltonian with heavily (infinitely) doped leads one gets the formula for the transmission coefficient as a function of the transverse momentum (k y = p y /h): For the periodic boundary conditions, the transverse momentum gets quantized values, k ( j) y = 2π j/W , with j = 0, ±1, ±2, . . . . For W L, one can approximate the sum over k ( j) y by an integral, obtaining the total transmission Next, we expand (up to second order) the dispersion relation around p = (γ 1 v 3 /v 2 0 , 0) (i.e., the satellite Dirac cone at 125425-7 ϕ = 0), arriving at The corresponding single-cone Hamiltonian reads This time, solving the scattering problem for a rectangular sample we get the transmission coefficient and the integration over k y leads to The calculations for remaining Dirac cones at ϕ = 0 are more involved, yet straightforward. Generalizing the above reasoning for p = γ 1 (v 3 /v 2 0 )(cos ϕ, sin ϕ), we get and where we have defined ζ = k y L.
Changing the variables according to T (k y ) ≡ T (ζ , L) we find that the Landauer-Büttiker conductivity, for a fixed W/L 1, is bounded by vanishing in the L → ∞ limit. For μ 4 = 0, the conductivity σ (L) ≈ (π/4)σ 0 , and the Fano factor F ≈ 1 − 2/π for L l ⊥ , being numerically close to the results by Snyman and Beenakker [10].
The approximate upper bound for μ 4 = 0, given in Eq. (B4), is further supported with the numerical results presented in Fig. 6, where we have set the doping in the leads such that (k F l ⊥ ) 2 = 0.2 (after Ref. [10]), and W/L = 20. Numerical calculations for the full four-band model given the Hamiltonian H BLG (B1) leads to a noticeably faster but also a power-law decay of the conductivity, which can be approximated as σ (L) ∝ L −2.0 for L 1000 l ⊥ .