Collective resonances near zero energy induced by a point defect in bilayer graphene

Intrinsic defects give rise to scattering processes governing the transport properties of mesoscopic systems. We investigate analytically and numerically the local density of states in Bernal stacking bilayer graphene with a point defect. With Bernal stacking structure, there are two types of lattice sites. One corresponds to connected sites, where carbon atoms from each layer stack on top of each other, and the other corresponds to disconnected sites. From our theoretical study, a picture emerges in which the pronounced zero-energy peak in the local density of states does not attribute to zero-energy impurity states associated to two different types of defects but to a collective phenomenon of the low-energy resonant states induced by the defect. To corroborate this description, we numerically show that at small system size N, where N is the number of unit cells, the zero-energy peak near the defect scales as 1/lnN for the quasi-localized zero-energy state and as 1/N for the delocalized zero-energy state. As the system size approaches to the thermodynamic limit, the former zero-energy peak becomes a power-law singularity 1/|E| in low energies, while the latter is broadened into a Lorentzian shape. A striking point is that both types of zero-energy peaks decay as 1/r2 away from the defect, manifesting the quasi-localized character. Based on our results, we propose a general formula for the local density of states in low-energy and in real space. Our study sheds light on this fundamental problem of defects.

. A nanotorus of BLG with a point defect on the upper layer. The top layer is sketched in a blue shadow, as compared to the bottom one. The black line represents an unit cell, which encloses four sites in sublattices A, B and layers 1, 2. The intra-layer hopping is denoted as t, and two different inter-layer hoppings are respectively represented as γ 1 between B 1 and A 2 (connected sites) and γ 3 between A 1 and B 2 (disconnected sites), in the unit of t. The defect potential is denoted as V 0 . , with an impurity potential V 0 . We note that a vacancy, which corresponds to the elimination of lattice sites without lattice relaxation, is equivalent to the unitary limit of impurity potential, V 0 /t → ∞.
Before we consider the model with a point defect, it is instructive to understand the electronic structure of a pristine BLG in low energy. In the basis of Ψ = † † † † † , the Hamiltonian in the momentum space is represented as , where the 4 × 4 matrix  0 has the following form, at Dirac points k D = (k x , k y ) = (±4π/3, 0), the eigenstates at zero energy are here q± states correspond to the gapless continuum with quadratic dispersion E(k) = k 2 /(2m * ) and m * = γ 1 /2, and g± states correspond to the bands with finite gap ±γ 1 . It is evident that the q± states have large amplitudes at the disconnected sites A 1 /B 2 and the g± states have large amplitudes at the connected sites B 1 /A 2 . After introducing a point defect to a BLG, we can distinguish two different types of LDOS, due to contribution from the gapless continuum or the gapped one, associated with the position of the defect. These will then be investigated analytically and numerically in the following sections.
Finite-size Calculation. For a defect placed at a connected site B 1 , we study the LDOS at the first-nearestneighbor A 1 site in a reference frame centered at the defect position. The N dependent spectral weight of the LDOS represents the intriguing localization character of defect-induced states. First let us consider γ 3 = 0. As illustrated in Fig. 2 our numerical results reveal that the height of the peak at zero energy scales as 1/lnN when N is smaller than 10 6 . The 1/lnN behavior is a consequence of the quasilocalized zero-energy state existing around the point defect, which was discussed in the previous study 56 . However, as N increases to ∼ N 10 7 , the spectral weight would eventually saturate at a finite value. For γ 1 = γ 3 = 0.1, as shown in Fig. 2, the weight of the zero-energy peak decreases with strong oscillations. Nevertheless, the spectral weight still saturates at a finite value near ∼ N 10 7 . To understand the saturation, one might observe the low-energy behavior of the LDOS with respect to different system sizes. Similar to previous discovery in a MLG 55 , we find that a point defect also generates a lots of resonant peaks with large spectral weights near zero energy. When the system size approaches to infinity, the defect-induced resonant peaks crowd to zero energy, and eventually the spectral weight at zero energy saturates. The collection of these resonant states constitutes the zero-bias anomaly in the LDOS. In Sec. II C, we will analytically compute the LDOS in low-energy and in the thermodynamic limit and show that the peak is a power-law singularity.
For a defect placed at a connected site A 1 , Fig. 2 shows the N dependent zero-energy peak of LDOS at the B 1 site nearest to the defect. For γ 3 = 0 and < N 10 6 , the height of zero-energy peak scales as 1/N, which is attributed to a delocalized zero-energy state induced by a defect at a A 1 /B 2 site 56 . When N is larger than 10 6 , the spectral weight at zero energy saturates at a finite value. By investigating the low-energy behavior of the LDOS with respect to different system sizes, we find that the defect-induced resonant states crowd into the zero-energy regime when N increases to infinity. When γ 1 = γ 3 = 0.1, the height of the zero-energy peak decreases with oscillations, as shown in Fig. 2. Nevertheless, the spectral weight still saturates at a finite value.
We emphasize that the spectral weights of the induced resonant states are much smaller than those around a defect at B 1 /A 2 site. Therefore the height of the zero-energy peak is saturated at a relatively small value, as shown in Fig. 3. In the following, we will analytically compute the LDOS in the low-energy and the thermodynamic limit.
Analytical Computation in the Thermodynamic Limit. In previous section, we focused on an exact numerical evaluation of the LDOS in finite-size systems. For small N the finite-size-scaling of the zero-energy peak follows the localization character of zero-energy states induced by a point defect. When N becomes large, collective phenomena from the induced resonant states near zero energy are expected to become of particular relevance. Eventually the zero-energy peak does not vanish but saturates at a finite value. Now we will analytically derive the change of the LDOS in the low-energy and thermodynamic limit.
The Green's function techniques allow us to obtain the change of the LDOS at i site nearby the defect at j site, where i, j = A 1 , B 1 , A 2 , B 2 are within the same unit cell. In the last step, we used time-reversal symmetry and took the unitary limit V 0 /t → ∞, which makes Δρ i independent on the strength of the impurity potential. We are interested in the low-energy regime near the Dirac points k D = (±4π/3, 0), where the Hamiltonian, Eq. (2), is written as , where q± correspond to the quadratic bands ε q± = ±q 2 /γ 1 , and g± correspond to the gapped bands ε g± = ±γ 1 . Evaluating the Green's functions analytically (for details see Method), we obtain the energy dependance of LDOS for two different types of defects.

A Point Defect at a Connected Site B 1 /A 2 .
To the leading order, the retarded Green's functions are approximated as , where Λ is a high-momentum cut-off. If a point defect is placed at a connected B 1 site, the change of the LDOS at the nearest-neighbor A 1 site is approximated as To confirm our analytical result, we compute the LDOS numerically in the thermodynamic limit. In the left panel of Fig. 4, we show the LDOS of a A 1 site before and after placing a point defect at the nearest B 1 site. As one can discern from Fig. 4, there exists large spectral weight transferred from the high-energy regime to the low-energy one in the presence of a point defect. Because the logarithmic correction is very weak, our numerical results, shown in the right panel of Fig. 4, exhibit a 1/E power-law singularity for the change of the LDOS. We find that the power-law singularity are robust for both γ 1 = γ 3 = 0.1 and γ 1 = 0.1, γ 3 = 0. These results on the power-law singularity allow us to draw connections to the previous study in MLG 55 . We will develop a simple interpretation in terms of Harper equations in the subsection 2.3.3.
A Point Defect at a Disconnected Site A 1 /B 2 . When a point defect is located at a A 1 site, we approximate the retarded Green's function to the leading order, , as shown in Method. Thus, the change of the LDOS at the nearest-neighbor B 1 site is expressed as It is remarkable that the change of the LDOS takes the form of Lorentzian function with the broadening factor γ 1 , the inter-layer hopping. To confirm the analytical result, we show the numerical study of the LDOS in the thermodynamic limit in it is reasonable to ignore the inter-layer hopping γ 1 , γ 3 momentarily. Therefore, the problem is reduced to a point defect problem in MLG. As the system size grows to infinity, the zero-energy singularity 54,55 induced by a point defect at For a point defect at a B 1 site, the Harper equations involving ϕ A 1 do not change at all if we turn on the inter-layer hopping γ 1 but keep γ 3 = 0. Instead, turning on γ 3 will cause the wave function ϕ A 1 spreading to A 2 sites. The effect is described by the Harper equation at B 2 sites: where δ i and δ ′ i are the three displacement vectors pointed from a B 2 site to the nearest A 1 sites on the top layer and to the nearest A 2 sites on the bottom layer, respectively. Because of γ  1 3 , the spatial wave function ϕ r ( ) A 1 remains more or less robust, and ϕ r ( ) is of the order of γ 3 . Thus, ϕ r ( ) accounts for a rather small spread of the wave function from A 1 sites to A 2 sites. Since in the low-energy limit excitations living on A 2 sites are gapped, there is no significant effect for the zero-energy impurity states being coupled to a gapped continuum. Therefore, the LDOS in low energies can be understood in a simple monolayer picture, where the gapless continuum living on A 1 sites is disturbed by a point defect. Based on the previous study of MLG 55 , significant defect-induced resonant states from the continuum lead to a power-law singularity in the LDOS. Our numerical and analytical results indeed confirm our argument in the thermodynamic limit for BLG. Now we consider the case of a point defect at a A 1 site. When we gradually turn on the inter-layer hopping γ 1 , this gives rise to the Harper equation of A 2 site as where δ i represents the set of the three displacement vectors pointed from a A 2 site to the nearest B 2 sites on the top layer. Following the same argument, ϕ r ( ) is of the order of γ 1 . The hopping amplitude γ 1 can be viewed as a coupling between the zero-energy state living on B 1 sites and the gapless continuum on B 2 sites of the bottom layer. It explains that a sharp delta-function peak from a single state in the spectral function is broadened into a Lorentzian shape when coupling to a continuum. Meanwhile, the half-width of the Lorentzian is proportional to the coupling, the interlayer hopping amplitude γ 1 . This is indeed what happens in the our numerical and analytical results.
The spatial profiles of the zero-energy peak. While the numerical and analytic results in previous sections focused on LDOS in the energy domain, now we study the zero-energy peak in the spatial domain. In Fig. 6, we identify that the spatial dependance of the spectral weight at zero energy is proportional to 1/r 2 for two different types of point defects, located at A 1 or B 1 site. This 1/r 2 behavior manifests the quasi-localized character of the zero-energy peak. In previous sections we have excluded a single impurity state, either quasi-localized or extended, at zero energy as the cause for those peak in the thermodynamic limit. Figure 6, however, suggests that the induced resonant states which crowd to zero energy share the same spatial profile with the quasi-localized state.
Although the spatial dependence of the resonant states is inaccessible by analytic approach, our numerical support allows us to propose an asymptotic formula for the LDOS in low energies The factor F(r) = 1/r 2 shows the spatial profile with the quasi-localized character. The first term represents the quasi-localized zero-energy state with the normalization C N ln A 1 = for a defect at a connected site B 1 , or the delocalized one with the normalization C N B 1 = for a defect at a disconnected site A 1 56 . The second term Δρ i (E), defined by Eqs (6) or (7), is the contribution from the resonant states in the low-energy regime induced by the two different types of defects. We note that here the LDOS from a pristine BLG is neglected, because it is much smaller than above two terms. The LDOS in Eq. (10) shows that as the system size approaches to the thermodynamic limit, the contribution from zero modes fades away and the LDOS is dominated by infinite resonant peaks. We emphasize that the formula, Eq. (10), is proposed based on our numerical observation. To fully understand the LDOS in low energies, analytic solutions for the resonant states are still necessary.

Discussions
We have studied the zero-energy peak induced by a point defect in a Bernal stacking BLG. We numerically computed the LDOS at the first-nearest-neighbor of a point defect and investigated the system-size N dependence of the LDOS. For small size, the zero-energy peak of the LDOS scales as 1/lnN for the induced quasi-localized state or as 1/N for the delocalized state at zero energy. When N approaches infinity, the defect-induced resonant states crowd into zero energy and lead to non-vanishing zero-energy peaks for both cases. To further support our numerical findings, we analytically evaluated the change of the LDOS in the thermodynamic limit. We found that the zero-energy peak for a defect at a B 1 /A 2 site becomes a power-law singularity, while the peak for a defect at a A 1 /B 2 site is broadened into a Lorentzian shape. By studying the spatial dependence of the zero-energy peak, we showed that the zero-energy peaks for both cases decay as 1/r 2 away from the point defect. Combing all above discoveries, we proposed a formula for the LDOS in low energies and real space.
The previous theoretical study in monolayer graphene shows that a finite density of vacancies leads to a sharp peak of LDOS exactly at the Fermi level, superimposed upon the flat portion of the DOS 49 . This impurity band can be observed from the numerical calculations. The quasi-localization nature of defect states enables an impurity band form near zero energy, and it is shown that defects induce an impurity band with density of state characterized by the Wigner semi-circle law 55 . The band width of the impurity band is proportional to the square root of the density of defects. Since the impurity band is a flat-band, which supports ferromagnetism when electrons correlation effect is included 59 . In BLG, we expect that an impurity band will appear with finite defect concentration, and the peak of the impurity band will still pin at the Fermi level. The quasi-localization nature of LDOS in BLG discussed in our study may, although further study is needed, exhibit similar physics which a monolayer graphene has.
In the present work we do not investigate correlation effects in the problem of BLG with a monovacancy. Here we would like to discuss the possibility of defect-induced magnetism due to the Coulomb interaction and the large enhancement of the local density of states near the vacancy. By including inter-electronic interactions on a pure π-band model of BLG, the enhanced LDOS around impurity sites implies that defects may lead to the formation of local magnetic moments 60 . Eventually the spin-split DOS could be observed and the peak of the enhance LDOS will not exactly locate at the Fermi level. However, in real graphene we should consider the sp 2 σ-orbital electrons and the lattice reconstruction near vacancies. To understand a reconstructed single-atom vacancy in graphene, we shall rely on the results from first-principle calculations 61,62 .
In principle, near a single atom vacancy there are three unsatisfied sp 2 σ-orbital electrons and one π-orbital electron. To maximize its spin configurations, a net magnetic moment is 4 μB. If the dangling σ orbitals are not passivated, the vacancy reconstructs. If we consider a planar structure, π and σ bonds do not mix. The ab initio calculations have showed that these three impurity levels from the three dangling σ bonds split due to the crystal field and a Jahn-Teller distortion 61,62 . The π bond state, however, remains at the Fermi energy, being introduced in the midgap of the π bands, which refers to a zero-energy impurity state 62 . If the relaxed structure for the vacancy allows non-coplanar structures with out-of-plane displacements, the σ orbitals near the Fermi energy will be able to hybridize with the π band. Since the zero-energy state decays as 1/r, its overlap with σ orbital states is large. This produces the dominant exchange interaction from Coulomb interaction between local zero-energy state from the π orbital electron and the σ orbital electron near the Fermi energy. This exchange, due to overlap between the π impurity state with the σ orbital state, is responsible to the spin-splitting of the vacancy-induced zero-energy π states. However, in the literature, the predicted magnetic moment varies widely for the defect [61][62][63][64][65][66][67][68] . The formation of the enhanced LDOS at Fermi level is crucial for determining the hybridization function between the impurity and itinerate states in the picture of the Anderson-Kondo model 69 . For BLG, further investigations will be on exploring how the two types of zero-energy peaks evolve when electronic correlations being considered. Since the correlation effects may be significant near defects, our study on LDOS here becomes important on the idea of fabricating spin qubits by defects 71 . The inclusion of realistic interactions to defect-induced resonant states is expected to become of particular relevance for maintaining quantum coherence for a qubit. Moreover, recent experiments 72 showed that in graphene localized states induced by disorders can significantly enhance electron-phonon coupling and become a local drain of hot carriers, which is crucial for electronic transport at low temperatures 73 . For BLG, two types of defects could play different roles on dissipation. The thermal imaging of dissipation shall be revealed by a superconducting quantum interference nano-thermometer 72 .

Method
Green's function. To calculate the LDOS, we solve the Dyson equation for BLG in the presence of a point defect with a potential V 0 . The LDOS can be expressed in terms of the non-interacting retarded Green's function. Assuming the defect placed at j site of the unit cell at the origin r = 0, we use standard Green's function techniques to formulate the LDOS at i site of the unit cell at r as Here the LDOS and the change of the LDOS are represented respectively as where we sum over all discrete k x and k y points, k x = 4πn x /N x , π = k n N 2 / 3 y y y and n x/y = 0, 1, 2, 3, ..., N x/y − 1. We investigate how the LDOS changes as the size N = N x × N y evolves from 10 2 to 10 7 .
Retarded Green's functions around the Dirac points. Here we will elaborate the calculation of the Green's functions in the thermodynamic limit. Near the Dirac points, the Hamiltonian, Eq. (5), be diagonalized by a unitary transformation , where the unitary matrix U connects the eigenbasis q The retarded Green's functions in the eigenbasis can be computed straightforwardly, i.e.
By employing time-reversal and structure symmetries of a bilayer, we obtain Integrating all states near the Dirac points within a momentum cutoff Λ, we can obtain the Green's functions in real space, 2  cos  1  8  cos  3  8  sin  3  2  sin  3  4 sin cos , where θ is the angle centered at k = (4π/3, 0 where δ i and δ i′ are the three displacement vectors pointed from a B 2 site to the nearest neighbors A 1 sites on the top layer and A 2 sites on the bottom layer respectively.