Full electric-field tuning of the nonreciprocal transport effect in massive chiral fermions with trigonal warping


 In a noncentrosymmetric system, an intrinsic electric
polarization is allowed and may lead to
unusual nonreciprocal charge transport phenomena. As a result, a current-dependent resistance,
arising from the magnetoelectric anisotropy term of k · E × B, appears and acts as a current
rectifier with the amount of rectification being linearly proportional to the magnitude of both
current and applied magnetic field. In this work, a different type of nonreciprocal transport effect
was demonstrated in a graphene-based device, which requires no external magnetic field. Owing to
the unique pseudospin (valley) degree of freedom in chiral fermions with trigonal warping, a large
nonreciprocal transport effect was uncovered in a gapped bilayer graphene, where electric-field
tunabilities of the band gap and valley polarization play an important role. The exact cancellation
of nonreciprocal effect between two different valleys is effectively removed by breaking the inversion
symmetry via electric gatings. The magnitude of the current rectification appears to be at a
maximum when the Fermi surface undergoes a Lifshitz transition near the band edges, which is
proportional to the current and the displacement field strength. The full electric-field tuning of the
nonreciprocal transport effect without a magnetic field opens up a new direction for valleytronics
in two-dimensional based devices.


INTRODUCTION
From a symmetry argument, a magnetoelectric anisotropy term in the form of k · E × B can generally exist in all systems [1], which gives rise to additional contributions to various magnetoelectric properties in matters [2][3][4][5]. Among them, the nonlinear and nonreciprocal transport effect (NRTE) is one of the intriguing subjects recently, where an unusual current-direction dependent on the charge transport is a direct consequence of the NRTE. Phenomenologically, it can be described by a current dependent resistance of R = R 0 (1+γ I · ( P × B)), where γ is the NRTE coefficient and P is the charge polarization of the system. As is evident from the equation, the NRTE appears only when current ( I), magnetic field ( B) and P are all nonzero and also perpendicular to each other, setting stringent criteria for the experimental realizations. So far, a number of discoveries for NRTE * These authors contributed equally to the work. † wlee@phys.sinica.edu.tw had been made in various systems, such as bulk-Rashba polar semiconductors [6], chiral mangets [7], noncentrosymmetric superconductors [8][9][10][11][12], and magnetic bilayer heterostructures [13][14][15]. In theories, the microscopic mechanism of NRTE is typically associated with the velocity asymmetry and the scattering lifetime asymmetry in the Boltzmann transport equation, giving rise to a nonzero transport coefficient in the E 2 term. Therefore, such a NRTE can be readily probed by measuring the second harmonic responses to an alternating driving current with a frequency of ω. In general, the nonlinear (E 2 term) resistance (R 2ω ) is negligibly small as compared to E-linear resistance (R ω ), but the ratio of R 2ω /R ω can be as large as 0.01 or higher in systems with NRTE.
As it turns out, the Dirac-like systems with chiral fermions are good platforms for the study of NRTE. The trigonal warping effect [16,17] and the chiral property with (pseudo)spin-momentum locking give rise to the asymmetries of energy dispersion and velocity on the Fermi surface. A large current dependent resistance was reported in magnetic topological insulator (TI) / TI heterostructure, which was attributed to the Dirac surface states and also magnon-driven asymmetric scattering [14]. On the other hand, in a thin film of nonmagnetic Bi 2 Se 3 , a sizable R 2ω also appears when the charge transport is dominated by the Dirac surface states [18]. The Dirac Fermi surface with trigonal warping is further distorted by an applied magnetic field parallel to the film surface, giving rise to the observed R 2ω signals. In other words, the applied magnetic field acts effectively as a spin-to-charge converter via the nonlinear resistance channel due to the spin-momentum locking. The aforementioned systems all require an external magnetic field for the appearance of NRTE, which put a severe limitation for the practical device applications. Here, we report a new type of NRTE that is rooted in the valley-momentum locking effect of massive and trigonal-warped chiral fermions with Berry phase 2π in a dual-gated bilayer graphene (BLG) device. The valley degree of freedom in a BLG is a counterpart of the spin degree of freedom, but the manipulation of the valley polarization, in principle, can be made possible by electric gating only. Large R 2ω signals were revealed in a dual gate BLG device whenever the Fermi energy (ε F ) was tuned to be near the band edges. Its magnitude grows larger with increasing displacement field strength |D| and also the bias current. We note that the |D| not only opens up a band gap in a BLG but also effectively polarizes the valley degree of freedom, and thus the exact cancellation of the NRTE from two different valleys of K+ and K− is removed to result in a large R 2ω response. This effect is fully adjustable by electric gating alone and requires no external magnetic field. Those findings can be understood qualitatively from the velocity asymmetry in the Boltzmann transport equation for massive chiral fermions with trigonal warping. Figure 1 (a) shows an optical image of a dual-gated BLG device we fabricated, where the measurement circuits were also illustrated (See Methods for more information). A BLG was sandwiched between two hexagonal boron nitride (hBN) microcrystals with thicknesses of about few tens of nanometers as illustrated in Fig. 1(b). The device geometry was designed to give more than 80 % of the BLG area between the voltage leads to be topgate effective. Therefore, the measured voltage signals are mostly dominated by the BLG region that can be simultaneously tuned by the top gate voltage (V tg ) and bottom gate voltage (V bg ), enabling the introduction of a large displacement field D on a BLG. By neglecting the screening effect, the strength of the displacement field can be calculated via

EXPERIMENTAL RESULTS
for bottom(top)-gate, and d b (d t ) is the dielectric layer thickness of bottom(top)-gate. The operation of a dual-gated device not only capacitively varies the density of the system but also changes the D values that acts on the BLG [19,20]. As demonstrated in Fig. 1(b), a band gap ε g appears due to a finite potential difference (∆) between the upper and lower layer of BLG due to D, giving rise to a rapid increase of the zero-density sheet resistance (R ω0 ) with increasing D as show in Fig.1(c) and (e). The color codes in Fig. 1(c) represent the values of sheet resistance R ω in logarithmic scale at T = 5K, and the corresponding R ω -V bg curves at different V tg values ranging from +3 V to -3 V are shown in Fig. 1 where V bg0 and V tg0 were determined to be ≈ -1.86 V and -0.15 V, respectively. The R ω0 equals about 2.2 kΩ at zero V tg , and it increases by more than two orders of magnitude of R ω ≈ 285.6 kΩ at V tg = -3 V that corresponds to a D ≈ +0.58 V/nm. The large variation in R ω indicates a high device quality, which is also supported by the high Hall mobility of ≈ 10,000 cm 2 /Vs as determined from the Hall measurements (See supplementary information section 1). Surprisingly, large second harmonic signals R 2ω were uncovered and shown in R 2ω values appear near the zero-density regime with magnitudes rapidly growing for |V tg | ≥ 2 V. As V bg sweeps from + 50 V to -50V, R 2ω first reaches a local maximum and then rapidly falls off to a local minimum with an opposite sign. At V tg = -3V, the magnitude of R 2ω is as large as 5.0 kΩ. The bias current I equals 100 nA for all the data shown in Fig.   1, and the curves at different V tg values in Fig.1 Both |R 2ωp(v) | and R ω0 are exponentially growing with D in low D regime, and it slows down for |D| ≥ 0.4 V/nm. By tuning D up to about a magnitude of ≈ 0.58 V/nm, we remark that the |R 2ωp(v) | spans more than four orders of magnitude, which is more than two-fold larger than that for R ω0 . The inset cartoon in the lower-panel of Fig. 2(c) illustrates the sign of D in relation to the device geometry.
In a dual-gated BLG device, the unscreened D values as function of V bg and V tg can be calculated as shown in Fig. 3(a). The sweeping of gate voltages along a D-constant contour line ensures the shifting of only the ε F in BLG while keeping the degree of inversion asymmetry unchanged. As an example, a sweeping along a black line, connecting point 1 and point 2 in Fig. 3 (a) along a D-constant contour line, will result in a shifting of ε F from the conduction band to the valence band as illustrated in the lower panel of Fig. 3(a). We carried out such sweepings along D-constant contour lines at 10 different D values ranging from +0.5 V/nm to -0.5 V/nm. The upper-panel and lower-panel of Fig. 3(b) show the measurement results for R ω and |R 2ω |, respectively, as a function of the corresponding sheet density n 2D . When ε F = 0 (n 2D = 0), R ω attains a maximum value of R ω0 , but R 2ω is nearly zero. As ε F moves toward conduction (valence) band with n 2D ≤ 0 ( n 2D ≥ 0), |R 2ω | first increases rapidly to a local maximum, and then it falls off to a vanishingly small value for |n 2D | ≥ 3 × 10 11 cm −2 . ∆n 2D is the sheet density at which the |R 2ω | reached a maximum as illustrated in the lower panel of Fig. 3 (b). An enlarged view of the R 2ω -n 2D data in a linear scale is shown in Fig. 3(c), which shares the same color codes as that in Fig. 3 As noted before, by reversing the sign of D (dashed lines in Fig. 3 (b) and (c)), consistent behaviors of R ω and |R 2ω | are observed with a minor difference in magnitude, but no sign changes in R 2ω were found due to the sign reversal of D. The extracted ∆n 2D is plotted as a function of D as shown in Fig. 3(d). For |D| ≤ 0.4 V/nm, ∆n 2D tends to drop slightly with an increasing |D| to a value of ≈ 3.7 × 10 10 cm −2 at D = 0.4 V/nm.
The T -dependent measurements were carried out to extracted the band gap ε g using a thermal-activation model for the conductance. The normalized zero-density conductance of R ω0 (D=0)/R ω0 versus 1/T curves for different D values up to ≈ 0.58 V/nm are plotted in  where it decays from a value of about 2.8 % at T = 5K to 0.28 % at T = 80 K. The ∆V bg , defined as the V bg difference between the local extremes in R 2ω -V bg curves, is shown in Fig.   4(d), which shows similar behavior as ∆n 2D in Fig. 3 (d). At T = 5 K, ∆V bg first decreases with increasing |D| up to a critical value of D c ≈ 0.4 V/nm. For D ≥ D c , ∆V bg shows a small increase with |D|. We note that D c turns out to move to a higher value as T increases.
For all the T -dependent data, the R 2ωp(v) versus the corresponding R ω collapsed on a single curve as demonstrated in Fig. 4(e), giving a relation of R 2ωp(v) ∝ R 2.02 ω .

NUMERICAL CALCULATIONS OF NRTE
As illustrated in the lower panel of Fig. 5(a), the sublattices A1 and B2 on different layers form an effective honeycomb lattice in the reduced 2-component pseudospin representation with two different valleys of K+ and K−. The reduced 2×2 Hamiltonian of a BLG [21] can be described by , m * is the effective mass, and σ is the Pauli spin matrices. The resulting 2-component pseudospin states, referring to the electron probabilities on sublattices A1 and B2, are thus chiral fermions with the valley-momentum locking effect. By further considering trigonal warping and an asymmetry parameter ∆ [22], the energy dispersion can be calculated and expressed as where (p =hk, φ) is the momentum in polar coordinate, γ 1 is the interlayer nearest neighbor hopping amplitude, ξ = +1(−1) refers to the K+(K−) valley, and υ is the Fermi velocity.
The trigonal warping derives from the additional velocity term of a is the lattice parameter for the sublattice, and γ 3 is the next nearest interlayer hopping amplitude as illustrated in the inset cartoon of Fig.1 x with respect to k x , but it cancels exactly while summing up contributions from both valleys of K+ and K−. Therefore, an additional breaking of the valley symmetry (i.e., valley polarization) is needed in order to have a non-zero contribution due to the asymmetry terms of ∂ε/∂k x and ∂ 2 ε/∂k 2 x , which is required for the appearance of NRTE to be shown as follows.
Based on the Boltzmann transport theorem [23], the nonequilibrium distribution function (g k ) up to the second order of driving electric-field (E 2 ) can be described by , where we used E = E xî . The resulting second harmonic conductance σ 2ω due to NRTE at T = 0 K for the low energy bands (α = 1) can be expressed as (See supplementary information section 3) where τ is the transport lifetime, p =hk, and υ k,ξ = ∂ε 1,ξ /∂p. For simplicity, we neglect the complex asymmetric scattering effects for chiral fermions in a BLG, and Eq.2 reduces to a simple form of σ 2ω,ξ = − e 3 τ 2 2πh 2 I 2ω,ξ , where ε ξ -dependent I 2ω,ξ term can be numerically calculated from Eq. 1 and shown in Fig. 6(a)  increases, reaching a maximum value near the band edge at ε ≈ 25 meV (dashed-line in Fig.   6(a)). For ε ≤ 25 meV, a Lifshitz transition occurs [21,24,25], and the single roundedtriangle Fermi surface transformed into multiple separated ellipses as illustrated on the right panel of Fig. 6(a) (See supplementary information section 3). The magnitude of I 2ω,+1 drops rapidly to zero for ε ≤ 24.3 meV, where the total area of ellipses and thus the density of states are negligibly small. As expected, the I 2ω,+1 turns out to be canceled exactly by I 2ω,−1 for K− at the same ε (red line in Fig. 6(a)). We remark that such a Lifshitz transition is a direct consequence of the trigonal warping effect that becomes significant at lower energy, and it happens for all ∆ values near the band edges [25]. The corresponding energy ε L for the occurrence of Lifshitz transition equals ≈ 1 meV for D = 0, and ε L ≈ ∆/2 for a finite D.
For more rigorous calculations, the asymmetric scatterings in τ and correlation effect should be considered [14,26,27], which is beyond the scope of current work. However, we argue that the inclusion of the asymmetric scatterings in τ should further amplify the asymmetry terms in Eq. 2 and thus may not affect the qualitative scenario that we described here.

DISCUSSIONS
The breaking of the inversion symmetry by D not only results in a finite ε g but also gives rise to a valley polarization in a gapped BLG [28]. A large nonlocal resistance has been reported in a gapped BLG and is attributed to the inverse valley Hall effect with Denhanced valley polarization [29][30][31]. As illustrated in Fig. 6(b), the D-induced potential difference between the upper and lower layer of a BLG causes a small energy difference for the electrons in different sublattices (valleys). The energy difference can be described by ∆ε ≈ eDdδn/n c , where d is the distance between layers, n c ≡ γ 2 1 πh 2 υ 2 , and δn is the D-induced density difference between layers. The resulting Fermi energies in K+ (black solid line) and K− (red solid line) are then shifted in the opposite direction as illustrated in the right panel of Fig. 6(b). The total second harmonic conductance is then described by As a result, a non-zero σ 2ω appears when ∆ε and thus D is finite. The measured second harmonic resistance can be linked to the calculated σ 2ω via R 2ω =− I 0 . This is an outcome from the band calculations at zero temperature and, in principle, can be smoothed out when taking into account a finite temperature effect. In general, ∂I 2ω,+1 /∂ε appears to be the largest when ε F is close to the band edge with multiple separated Fermi surfaces.
It then rapidly falls off as ε F ≥ ε L . The decrease of ∂I 2ω,+1 /∂ε with increasing ε becomes faster as D increases, which qualitatively agrees with the data shown in the lower-panel of Fig. 3 (b). We also note that the density for the maximum R 2ω shown in Fig. 3(d) gives ∆n 2D ≈ 3 -8 × 10 10 cm −2 , which is within the same order of magnitude for the density at E L ≈ ∆/2 (See supplementary information section 3). On the other hand, by observing Eq. 3, the sign of total σ 2ω is dictated by ∆ε ≈ Dδn. When ε F moves from the conduction band to the valence band with a constant D, δn changes sign, and thus R 2ω undergoes a sign reversal. On the contrary, reversing the sign for D also causes a sign reversal in δn, and it thus gives no sign changes for the resulting R 2ω . Those behaviors are well consistent with our observations shown in Fig. 1(f) and Fig. 3(c). For an order of magnitude estimation of R 2ω , we consider the case with a D = 0.6 V/nm and a total density of n 2D ≈ 6 × 10 10 cm −2 based on our experimental data, and it gives a ∆ε ≈ 0.8 meV, which takes into account the self-consistent screening effect [21]. A rough estimation of the transport lifetime, using an effective mass of 0.024 m e [35] and a sheet resistance of 200 kΩ, gives τ ≈ 7 fs. From Eq.3 , the calculated R 2ω for a bias current of 100 nA equals ≈ 2 kΩ with a ratio of R 2ω /R ω ≈ 1%, which is in the same order of magnitude as the measured R 2ω .
, the measured R 2ω is expected to scale with R ω with exponent larger than 3. This is, however, not the case as shown in Fig. 4(e), giving R 2ω ∝ R λ ω with an λ ≈ 2.02. As it turns out, the exponent λ is device dependent (See supplementary information section 4). The charge-puddle effect is known to give rise to potential fluctuations and a finite minimum conductance at charge neutral point in graphene-based systems [32][33][34].
It may enhance the conductance near the band edges, resulting in a much lower R ω and spoiling the scaling of R 2ω = R λ ω . A rough estimate from V bg0 value in our device gives an impurity-charge density n imp ≈ 1.3 × 10 11 cm −2 , which is larger than the ∆n 2D shown in Fig. 3(d). This indicates the importance of the charge-puddle effects in our device, and more rigorous theoretical descriptions that include charge puddle effects are needed. On the technical side, we also noticed a ω-dependent on the measured R 2ω (See supplementary information section 5). This is attributed to the presence of certain capacitance effect from the contact and gate electrodes, which can become larger as ε g increases in high D regime.
Such a finite capacitance in the measurement circuits can shunt the higher frequency signals and thus cause a much reduced amplitude of R 2ω as compared to R ω . Even though the issues of the magnitude of R 2ω and its relation with R ω require further investigations, the observed large R 2ω /R ω is well reproducible in dual-gated BLG devices, and, in principle, a larger R 2ω /R ω can be achieved with a further increase in D. The upper panel of Fig. 6(d) illustrates a scheme for a practical electric-field tuning of the current rectification effect. By sourcing a DC voltage (V dc ), the measured DC current (I dc ) is found to give a high-current and low-current states, corresponding to the setting values of V bg for extremal R 2ωp and R 2ωv , respectively. As demonstrated in the lower panel of Fig. 6(d), the difference in the high-current and low-current states scale linearly with the magnitude of V dc and thus I dc .
Upon reversing the sign of V 2ω , the high-current and low-current states are also reversed as expected. This is consistent with a current-dependent resistance in the form of R = R 0 (1+γ (D)I), where γ (D) is a D-dependent factor.

CONCLUSION
In conclusion, we uncovered an unusual large R 2ω signal in dual-gated BLG devices, which turns out to grow with increasing bias current and also displacement field D. Such a NRTE in a zero magnetic field sets a clear distinction from other NRTE systems that require the breaking of the time reversal symmetry by magnetic field. The largest ratio of R 2ω /R ω for I = 100 nA was found to be about 2.8% at T = 5 K and D = + 0.58 V/nm. The extremal R 2ω occurs whenever the Fermi energy is shifted to the band edges, and its magnitude quickly declines when Fermi energy moves away from the Lifshitz transition energy.
Those behaviors qualitatively agree well with the calculated velocity asymmetry contribution in the Boltzmann transport equation, giving a total second harmonic conductance ∂ε ] · ∆ε with ∆ε ≈ eDdδn/n c . The breaking of the inversion symmetry removes the exact cancellation between I 2ω,+1 and I 2ω,−1 by D, which is crucial for the appearance of NRTE in a BLG. The scaling relation of R 2ω ∝ R λ 2ω gives a device dependent values for the exponent λ ranging from 1.3 to 2, which infers somewhat reduced increase of R 2ω with D in our dual-gated BLG devices. Nevertheless, the revelation of a large R 2ω for the massive chiral fermions with trigonal warping enables a full electric-field tuning of the NRTE and opens up a new possibility for the valleytronics in graphene-based 2D devices.

METHODS
High quality natural graphites and hexagonal boron nitrides were mechanically exfoliated by tapes and then transferred onto 300 nm SiO 2 /Si wafers. From optical contrast analyses and atomic force microscope scannings, high-uniformity flakes with proper thicknesses were selected and sequentially stacked using the dry transfer technique with a polycarbon-ate/PDMS stamp. For the device shown in Fig. 1(a), the thicknesses for the top hBN-layer and bottom hBN-layer are about 22 nm and 6 nm, respectively. Contact electrodes of Ti(5 nm)/Au(30 nm) were then fabricated by electron beam lithography facility. The device measurements were carried out in a closed-cycle cryostat that provides stable sample temperature above T = 5 K. A lock-in amplifier was used for the measurements of R ω and R 2ω simultaneously with a bias current at a low frequency of 18.111 Hz.

DATA AVAILABILITY
All the supporting data are included in the main text and also in supplementary information. The raw data and other related data for this paper can be requested from W.L.L.
(wlee@phys.sinica.edu.tw).     to the Fermi surfaces at 6 different energies for K+ (K−). A transformation into three separated ellipses was found when ε F is close to the band edge. Upon introducing D on a BLG, a small energy difference ∆ε ≈ eDdδn/n c between ε K+ and ε K− appears as illustrated in (b). It thus breaks the exact cancellation of I 2ω terms for K+ and K−, giving rise to a nonzero σ 2ω . The corresponding term of ∂I 2ω /∂ε as a function of shifted energy (ε − ∆/2) is shown in (c) for three different ∆ values of 50 meV, 30 meV and 10 meV. In general, the magnitude of σ 2ω is at maximum when ε F is near the band edge, and it falls off rapidly as ε F moves away from the Lifshitz transition energy (shaded region in (c)). A demonstration of the current rectification in a dual-gated BLG device is illustrated in the upper-panel of (d), where a DC current (I dc ) drawn by a biased voltage (V dc ) was measured (lower-panel of (d)) with alternating V bg values at which the R 2ω attains local extremal values.