Evolution of quasi-bound states in the circular n–p junction of bilayer graphene under magnetic field

Electron in gapless bilayer graphene can form quasi-bound states when a circular symmetric potential is created in bilayer graphene. These quasi-bound states can be adjusted by tuning the radius and strength of the potential barrier. We investigate the evolution of quasi-bound states spectra in the circular n–p junction of bilayer graphene under the magnetic field numerically. The energy levels of opposite angular momentum split and the splitting increases with the magnetic field. Moreover, weak magnetic fields can slightly shift the energy levels of quasi-bound states. While strong magnetic fields induce additional resonances in the local density states, which originates from Landau levels. We demonstrate that these numerical results are consistent with the semiclassical analysis based on Wentzel–Kramers–Brillouin approximation. Our results can be verified experimentally via scanning tunneling microscopy measurements.


Model and quasi-bound states spectrum
In this section, we consider the scattering of a plane wave electron on a CNPJ in BLG and calculate the local density of states (LDOS) of QBSs based on the two-band continuum Hamiltonian. Then we use the LDOS map to analyze how the QBSs' properties change with different potential barrier radius and strength.
The Bernal ( A − B ′ ) stacked BLG is shown in Fig. 1a. Taking into account in-plane hopping parameter γ AB = γ A ′ B ′ ≡ t and inter-layer coupling parameter γ A ′ B ≡ t ⊥ for undoped BLG, four bands model can be obtained by considering one 2p z orbital on each of the four atomic sites in the unit cell ( A, B, A ′ , B ′ ) 2,23 . Near the Dirac point K and K ′ , the two low-energy bands are touched and can be approximated as E ± (k) = ± 2 k 2 2m , where m = |t ⊥ | 2v 2 is the effective mass, v = 10 6 m/s is fermion velocity of electron, and a ≈ 1.42Å is the nearestneighbour distance.
We consider a CNPJ of gapless BLG and model a circular potential barrier with the step-like potential V = V (r)σ 0 = V 0 �(R − r)σ 0 , as shown in Fig. 1b,c. Focusing on the dynamics near a single Dirac point at K, the full two-band Hamiltonian is given by 38 : where p ± = p x ± ip y . The validity of Eq. (1) are discussed in Section "Discussion".
To solve this equation, we start by writing the canonical momentum operators as p ± = i e ±iφ ∂r ± i r ∂φ . The Hamiltonian commutes with the pseudo angular momentum operator J z = L z + σ z due to the radially symmetric potential. Here, L z = (r × p) z , and σ z is the third Pauli matrix. Then, we need look for eigenfunctions of J z = L z + σ z with eigenvalues j = l + 1 = 1, 2, 3,…, where l = 0, 1, 2,…. Assuming wavefunction is where the phase factor e ±iφ is derived from the BLG Hamiltonian, the coupled eigen equations are obtained: (3)   (4) by considering a scattering processes 25,39 , that an incident electron with energy E in BLG is scattered by a CNPJ created by a gate-induced circular potential barrier V(r). The incident plane wave in the n region can be written as a certain combination of cylindrical waves, then we utilize the scatter theory and the properties of and Bessel functions to get the wavefunctions outside the barrier (n region) as 25 : Within the barrier (p region), the regular eigenfunctions of the Hamiltonian H with energy E are 25 These functions are simultaneously eigenfunctions of H and J z , with eigenvalues E and j , respectively. Here, H (1) j , H (2) j are the Hankel functions of first and second kind, K j , I j are the modified Bessel functions, wavevectors in the n, p region, respectively. And we use α ′ = sgn(E) and α = sgn(E − V 0 ) to ensure the proper signs for electrons and holes.
In the n region, H j , H j , and K j are effective eigenfunctions which are bounded for large arguments, and we disregard other eigenfunctions that diverge for large arguments. Similarly, in the p region, we consider J j and I j but ignore the other eigenfunctions which are divergent at the origin. Thus, the complete wavefunction can be written as 25 in the n and p region, respectively. The coefficients S j , A j , B j and C j can be obtained from the boundary conditions at the interface of the CNPJ: the wavefunctions and their derivatives at r = R are continuous Therefore, we can calculate the local density of states by LDOS(j, r, E) ∝ |Ψ (r, E)| 2 . In the n region, where In the p region, where (4) j+1 (k n r)e iφ e ijφ , i j (r, φ) = I j−1 k p r e −iφ αI j+1 k p r e iφ e ijφ .
(10) ψ (n) j =h (2) j + S j h (1) j + A j k j ,    for resonant states at different angular momenta. As j increases, the resonant modes shift from the center of the quantum dot and gradually move outward, which acts similarly to those in monolayer graphene 10 . However, the QBSs in BLG for l = 0 are narrower compared to the l = 1 , l = 2 modes. This feature is different from monolayer graphene due to their different band structures 10,26 . The QBSs spectra can be measured experimentally via STM. Besides, we notice that the QBSs are formed in the p region, which can be regarded as BLG quantum dots. Figure 2b and c show the QBSs energy levels change with different potential barrier radius R and strength V 0 for the j = 1 mode at position r 0 = 10 nm. Here, the peaks of LDOS in the p region represent the energy levels of QBSs. In general, the higher barrier potential trap more QBSs along with the wider energy spacings. Likewise, larger bilayer graphene quantum dot can trap more QBSs with the narrower energy spacing. Additionally, the trapping time can be obtained through half-width of energy levels by τ = △E . For larger R and higher V 0 , QBSs can be trapped longer. These results suggest that we can confine specific energies and angular momentum modes by adjusting the potential size and depth. Note that the above calculations neglect valley mixing, for reasons explained in Section "Discussion".

Energy spectra of quasi-bound states in magnetic fields
In this section, we numerically solve the radial equation for a CNPJ of BLG in the presence of an external perpendicular magnetic field. Following, we focus on the case that the magnetic field is not sufficiently strong to make system fully evolve into Landau levels. In order to provide a simpler and more intuitive physical picture, we also give a semiclassical analysis of QBSs based on the WKB approximation at the end of this section.
When a magnetic field is applied perpendicularly on the graphene surface, the orbital motion of electrons in two-dimension is quantized and the spectrum becomes discrete, called Landau levels. These Landau levels F + k p r =B j J j+1 k p r + C j I j+1 k p r Numerical solution under magnetic fields. To obtain the QBSs under the magnetic field, we start by solving the Eqs. (22) and (23) via the two sides finite difference method discretized in 600 sites in the interval 0 + < r < L ( L > R is truncation position). The initial wavefunction of finite difference method is ψ 0 at r ≃ 0 + , while ψ L at r = L . According to the finite difference method, ψ 0 and ψ L evolve with the formula from two sides to n-p junction boundary r = R . Then, analogy to the analytic method at B = 0 in Section "Model and quasibound states spectrum", we apply the boundary condition of ψ 0 and ψ L at r = R to obtain the new coefficients under the magnetic field 11,26,35 . To be specific, at r ≃ 0 + side, owing to eBr ≪ 2j−1 r and e 2 B 2 r 2 4 2 + eB(j−1) ≪ j 2 −1 r 2 , Eqs. (22) and (23) can be reduced to Eqs. (3) and (4) by neglecting the magnetic terms. Consequently, we can directly use the analytical solution of zero magnetic field case, a set of Bessel function [Eqs. (5)(6)(7)(8)(9)], as ψ 0 . For other side r = L , the magnetic field become dominant to give rise to the Landau level spectrum. Thus, we can not directly utilize the Eqs. (5)(6)(7)(8)(9) as initial wavefunctions at r = L . Here, we consider the low magnetic field case, and under the influence of disorder the wave function acquires a Lorentzian weight n △ 2 (ε−εn(B)) 2 +△ 2 on the zero magnetic field ψ 0 to produce ψ L . This treatment continuously returns to the zero magnetic field case, and reflects the Landau levels effect on large distance induced by magnetic field. Note that this approximation works well when the magnetic length l B = eB is comparable to barrier radii R. Because small l B will make the whole system evolves into Landau levels, which beyond our research. Here, △ = 0.5 meV is a broadening parameter from disorder scattering, ε n (B) = ω c √ n(n − 1) is spectrum of BLG under magnetic field 2 , and ω c = eB mc is the cyclotron frequency of non-relativistic electrons with effective mass m. Moreover, the results are proved to be insensitive to the details of the cutoff, as an example, we take a cutoff at L = 1.5R below.
At relatively weak magnetic field B < 0.5 T, magnetic field only slightly shifts QBSs, as shown in Fig. 3a,b. Here, we plot the LDOS in logarithmic scale to clearly display the subpeaks induced by weak magnetic field. The energy shift is about 0.06 meV for 0.1 T, which evaluated from Fig. 3c. Furthermore, we notice the system has time-reversal symmetry at B = 0 T, which guarantees the degeneracy pair E K (j) = E K (−j) . However, the finite magnetic field break time-reversal symmetry of the system. The degenerate states of opposite angular momentum separate and the energy splitting enlarges as B increases. Specifically, with B increasing, the energies of QBSs for j = +n mode decrease while for j = −n increase.
For relatively strong magnetic field, we plot the evolution of QBSs spectra under B = 1.1 T, 1.3 T and 1.7 T, as shown in Fig. 3d,e. Comparing with the weak magnetic case, a stronger magnetic field have a more obvious effect on the QBSs due to the appearance of Landau levels and this effect becomes larger as B increases. The Landau levels appear next to the QBSs energy levels, and it seems to be no arresting interplay between them from numerical results. The QBSs and Landau levels do not merge or simply repel but coexist in this transition region. The LDOS peaks are the superposition between the confined state and the Landau levels. If we continue to increase the magnetic field B, when it exceeds the critical magnetic field B c , the QBSs will disappear and the whole system will evolves into Landau levels. Here, B c can be evaluated from the magnitude of Landau levels and QBSs energy, and it is about 3 T at R = 40 nm for j = 1 mode. Therefore, the QBSs in a CNPJ of BLG can be tuned by adjusting the magnetic field strength as well.
WKB approximation for zero and weak magnetic fields. Aforementioned numerical results cannot directly show the impact of particular parameters on the QBSs. To obtain a better understanding of numerical results, we can analyze the QBSs with semiclassical method based on WKB approximation. Reconsidering Eqs. (22) and (23), we firstly rewrite these equations in the matrix form 37,40,41 : www.nature.com/scientificreports/ (29) χ(r) = ψ(r)exp i y(r)dr .  www.nature.com/scientificreports/ we can obtain the relation between energy level E n,j of QBSs and quantum number n. Here, we take phase factor δ = 0 , because the step-like circular potential can be regarded as two vertical barrier potential in radial direction, and this geometrical shape of boundary corresponds to δ = 0 44 . In Fig. 4, the blue curves in rightmost panels depict the relation of E j and n, and the red dots represent the position of QBSs. Setting B = 0 T in the above formula, the energy levels of QBSs are consistent with the rigorous results shown in Section "Model and quasi-bound states spectrum". Likewise, at relatively weak magnetic field B < 0.5 T, the WKB solution are in accordance with the numerical results in Section "Numerical solution under magnetic fields". These results verify the availability of the WKB approximation. Thus, WKB approximation provides an easier way for predicting the QBSs energies.

Discussion
Throughout the above analysis of CNPJ, we have modeled the electrostatic potential as a step-like function of position. This assumption is justified if a ≪ R ensures the absence of inter-valley scattering at the interface, where a is the lattice constant and R is the characteristic length representing the width of the transition region between the junction's n and p sides. The inter-valley scattering is inevitable in experiments, but as investigated in Fig. 3 of Ref. 15 , the authors demonstrated that the inter-valley scattering caused by the step potential is very weak and QBSs are nearly insensitive to the smoothness of boundary. These results are in good agreement with those of experiments 10,13,14 . These previous investigations justify the validity of our approximation. Regarding the feasibility of the two-band Hamiltonian in Eq. (1), there are two points needs to address. Firstly, we require E F < t ⊥ /2 ≈ 200 meV to ensure the quadratic dispersion relation holds. Secondly, we neglect the trigonal warping term. Thus, our model is reliable for quasi-particles in the energy regime under 200 meV 25 . Furthermore, we calculate the LDOS as a function of energy derived from the two-band Hamiltonian and fourband Hamiltonian to verify the validity of the two-band Hamiltonian in Eq. (1), as shown in Fig. 5.
Besides, the QBSs of CNPJ in our study differ from those of Coulomb potential. The QBSs of long range Coulomb potential exhibit the dramatic property of discrete scale invariance 45,46 . In contrast, the circular potential added on our system is confined potential and the QBSs of them don't show the discrete scale invariance. Thus, the results of these two kinds of potential are different in the qualitative and quantitative studies. Moreover, the construction of circular np junction in graphene has already been achieved experimentally, and our theoretical results can be useful for qualitative analysis to them.

Summary
In this paper, we have studied the quasi-bound states in a circular n-p junction of bilayer graphene and their evolution under the magnetic field numerically. We have shown that the quasi-bound states spectra can be controlled by adjusting the potential barrier radius and strength. These energy spectra are quantitatively different www.nature.com/scientificreports/ from those for monolayer graphene due to different band structures. We also have demonstrated that the energy level degeneracy of opposite angular momentum states breaks under magnetic field and the energy splitting enlarges as magnetic field strength increases. Moreover, applying weak magnetic fields on system leads to slight shift of quasi-bound states. While the strong magnetic fields induce additional resonances beside the quasi-bound states. These additional resonances originate from the Landau levels. The evolution of quasi-bound state spectra under magnetic field is also supplemented with semiclassical analysis based on the WKB approximation. Our results are highly relevant to recent experiments and can be verified in STM measurement. www.nature.com/scientificreports/