Electric field control of spins in bilayer graphene: Local moment formation and local moment interactions

We study local moment formation for adatoms on bilayer graphene (BLG) within a mean-field theory of the Anderson impurity model. The wavefunctions of the BLG electrons induce strong particle-hole asymmetry and band dependence of the hybridization, which is shown to result in unusual features in the impurity model phase diagram. We also study the effect of varying the chemical potential, as well as varying an electric field perpendicular to the bilayer; the latter modifies the density of states of electrons in BLG and, more significantly, shifts the impurity energy. We show that this leads to regimes in the impurity phase diagram where local moments can be turned on or off by applying modest external electric fields. Finally, we show that the RKKY interaction between local moments can be varied by tuning the chemical potential (as has also been suggested in monolayer graphene) or, more interestingly, by tuning the electric field so that it induces changes in the band structure of BLG.


Introduction
Single layer graphene hosts a plethora of phenomena that arise from the Diraclike band dispersion and chirality of its low-energy quasiparticle excitations [1,2].It is interesting to explore how these unusual single particle properties impact the physics of adatoms on graphene.The combination of adatom-graphene hybridization and Hubbard-like interactions on the adatom has been studied in the context of local moment formation [3,4], Kondo physics [5,6,7,8,9], RKKY interactions [10,11,12,13,14,15,16], and adatom positional ordering [17,18,19,20].The study of adatoms on monolayer graphene is also of interest to the nanoscience and quantum computation communities given the possibility to control local moment physics, and adatom-adatom spin and density interactions, by varying the carrier concentration via gating [21].
In contrast to monolayer graphene, bilayer graphene (BLG), which has Bernal stacking of single layers, has an extra tuning parameter.Using a dual-gate geometry, shown in figure 1a, enables one to separately tune the chemical potential and an electric field perpendicular to the layers, which is equivalent to separately tuning the potential on each of the two layers of BLG.While tuning the chemical potential modifies the carrier concentration, applying an electric field normal to the layers generates a gap in the band structure of BLG [22,23,24,25,26,27] (Figure 1c).(We will refer to the potential difference between the two layers, induced by this electric field, as the 'bias'.)Such a tunable gap system enables one to envision device applications and the ability to dynamically control various states in BLG [28,29] This tunability also allows for the study of interesting fundamental physics -for instance, it has been shown that engineering the electric field to flip direction (from pointing up to pointing down) as a function of position leads to localized one-dimensional modes at the kink in the bias [30,31,32].We have shown in recent work that incorporating interaction effects converts this 'nanowire' into a 2-band Tomonaga-Luttinger liquid whose properties, such as Luttinger parameters and mode velocities, can be controlled by the bias strength [33].
In this paper, we study adatoms in BLG.We examine local moment formation on the adatoms, RKKY interaction between such local moments, and how these effects can be controlled by tuning the chemical potential and a applying perpendicular electric field.
Our work goes beyond Ref.[34], which studied local moment formation for sitecentered adatoms on BLG, in several important respects.(i) We consider adatoms that are positioned at the center of a hexagonal plaquette on one of the layers.The study of this configuration is motivated by a recent ab initio study of adatoms in monolayer graphene that indicates plaquette centered impurities are generally more energetically favourable than on-site impurities [35].We expect a similar situation to hold in BLG.(ii) An applied electric field is shown to directly tune the impurity energy.This is because an impurity position will, in general, be located closer to the top layer of BLG.Accounting for this impurity energy shift allows us to identify regions of the phase diagram where local moment formation can be turned on and/or off by the application of a perpendicular electric field.(iii) For a particular impurity level chosen so that its renormalized (with self-energy corrections) energy level lies in the middle of gap in presence of the bias, we construct phase diagrams at zero, positive and negative bias by sweeping the chemical potential.The resulting phase diagram exhibits the onset of a Coulomb-blockade phase where any arbitrarily small U results in the formation of local moments.(iv) As a consequence of the chiral wavefunctions of BLG and the fact that the plaquette centered impurity adatom couples to many sites, the coupling between the impurity and the quasiparticles of BLG has strong momentum and band dependence.This affects many of the details of the phase diagram.For instance, the self-energy develops a large real part that has nontrivial frequency dependence, and substantially renormalizes the position of the impurity spectral peak in a manner that depends on the chemical potential and the applied bias.We provide detailed a physical explanation for how this affects the resulting phase diagrams, which were not provided in reference [34].Furthermore, to better illustrate the effect of the wavefunctions and chirality of BLG on the phase diagrams, the BLG system is contrasted with a fictitious system of non-chiral fermions with the same DOS and dispersion relation.(v) We go beyond the issue of local moment formation to address the tunable RKKY interactions between such local moments on BLG.We begin, in Section II, by introducing the Anderson impurity model specific to BLG.Section III summarizes the Anderson mean field theory formalism.Armed with this background, in Section IV we construct the impurity model phase diagrams for plaquette-centered adatoms (shown schematically in figure 1b).To highlight some of the unusual features of these phase diagrams, we contrast it with an impurity model of a fictitious system of electrons that have an identical dispersion but a band-independent coupling to the adatom.Finally, in Section V, we discuss the RKKY interaction, and its tunability, for local moments on BLG.

Adatom model in bilayer graphene
Consider an adatom on BLG, described by the Anderson impurity model [36], Here ks is the BLG electron dispersion for electrons labelled by momentum k and band index s.We assume a minimal model for the BLG dispersion that includes a nearest-neighbor hopping amplitude, t, to sites on the same layer, and an interlayer hopping amplitude, t ⊥ , between the two sites that sit one on top of the other.Henceforth, we set t = 1 and note that t ≈ 3 eV and t ⊥ /t ≈ 0.15 in BLG.In H imp , we denote the impurity energy by d , while U denotes the electron-electron repulsion on the impurity site.BLG electrons at sites r can hop on or off the adatom impurity with an amplitude χ r .We assume a common equilibrium chemical potential µ for the impurity and BLG electrons.The complete Hamiltonian for unbiased BLG is then given by H = H BLG + H imp + H mix .Electronic structure studies of transition metal adatoms on monolayer graphene suggest that the low-energy configuration of many types of impurities corresponds to the adatom residing at the center of a hexagonal plaquette [35].We therefore fix the adatom position to be at the plaquette center on the top layer (labelled = 1) of BLG, as shown in the schematic diagram on the right in figure 1.For simplicity, we assume that χ r =χ for the set of sites {r n }, which includes the six nearest neighbor plaquette sites in layer-1 and the site on layer-2 that lies directly below the adatom, and χ r =0 for all other sites.This simplifying assumption about the impurity model allows us to focus on unconventional features of local moment formation intrinsic to bilayer graphene.Future density functional studies would be useful in incorporating details of the impurity atomic orbitals.Turning to the mixing Hamiltonian H mix which allows the impurity electrons to hybridize with the BLG electrons, let us set where φ ks (r) denotes the wave function at site r for electrons in band-s and momentum k.We then obtain While the impurity model Hamiltonian in BLG looks similar to that in conventional systems or monolayer graphene, there are two important new ingredients in the impurity physics of BLG with plaquette centered impurities.First, for BLG (or multilayer graphene), as opposed to monolayer graphene, one can tune the density of states by applying an electric field perpendicular to the layers.Let ∆ denote the potential imbalance between the top and bottom layer induced by the electric field.Assuming that the adatom is at the same height as the top layer, this leads to an extra term in the Anderson Hamiltonian where r denotes the sites in the top ( = 1) and bottom ( = 0) layers.In writing this modification to the Hamiltonian, we have assumed that χ and t remain unchanged in the presence of an electric field.If intercalation of the impurity occurs, this will reduce the shift in the impurity energy, but will always be nonzero on grounds of the crystal symmetry.Incorporating the bias in this way thus has three effects: (i) a renormalization of the BLG dispersion; (ii) a modification of the hybridization V ks through a change in the BLG quasiparticle wavefunctions; and (iii) a shift the impurity energy to d + ∆/2.We will refer to the renormalized BLG dispersion and the hybridization as ks (∆) and V ks (∆) respectively.It is well-known that such a bias in BLG can open a band gap and significantly change the low-energy density of states; what is perhaps not appreciated is that this also effectively tunes the impurity energy in multilayer graphene.The last term in equation 6 describing this effect was not present in reference [34] and it will be shown to have a remarkable effect on local moment formation in presence of a bias.
A second important difference arises from the tunneling matrix elements, V ks , for the four bands of the bilayer.As shown in figure 2, these matrix elements display strong band-and momentum-dependence, which does not appear for the site-centered impurities discussed in reference [34].The rich structure of the coupling between the chiral BLG quasiparticles and the impurity site leads to a number of differences in the impurity model phase diagram when compared with conventional non-chiral fermions with a similar density of states, where we simply replace φ ks (r) ∼ exp(ik • r) in equation 4.

Mean field theory
A mean field treatment of the adatom impurity model is obtained, following Anderson [36], by setting where ρ d = σ n dσ , and With this mean field approximation, the entire Hamiltonian splits into two single particle impurity Hamiltonians, one for each spin, with These are coupled together by the self-consistency conditions that fix ξ dσ via m d and ρ d .The single particle Green function for the impurity is given by where the impurity self-energy is given by We can analytically continue this to the real frequency axis by setting iω n → ω + i0 + to obtain the real and imaginary parts of the self-energy Σ d (ω).We can then compute at T = 0 Within this mean field approach, the presence of a local moment on the impurity is signalled by a self-consistent solution with a nonzero m d .Alternatively, it is possible to self-consistently solve the mean field Hamiltonian using exact diagonalization for small system sizes.All of the phase diagrams in the next section were checked for consistency using this method.

Local moment formation
Using the above mean field theory enables us to study local moment formation on an impurity atom residing on BLG.Since the BLG band structure can be tuned by the electric field, we choose to define Γ 0 ≡ πχ 2 /t as a rough scale for the impurity level broadening in the absence of interactions.Thus, Γ 0 remains fixed for a given χ even as the electric field and chemical potential are varied.In this section, we begin by discussing the case when ∆ = 0 (i.e.without an applied electric field perpendicular to the layers).Phase diagrams are constructed by varying d and U for fixed χ = 0.3t (which implies χ ∼ 1 eV in conventional units) with various choices of the chemical potential.Next, we consider how varying ∆ can be used to tune the phase diagrams.After which, we discuss an alternative phase diagram for an impurity with a fixed bare energy level (although the actual energy level will be modified in the presence of a bias) with various choices of ∆.To construct these phase diagrams, µ and U are varied, while d and χ (= 0.3t) are kept fixed.
We have checked that varying χ modestly makes no qualitative changes to various features in the phase diagram, although it does shift the phase boundaries as expected.We ascribe the complexities of the impurity model phase diagram in BLG to the effective momentum-and band-dependent mixing V ks .As we discuss below, the strong variation of this coupling between different bands results in particle-hole asymmetry of the impurity model phase diagram via the impurity self-energy.This is despite the fact that in the simplest tight-binding parameterization, which we have considered, the BLG band dispersion itself is particle-hole symmetric for µ = 0.

Phase diagram in the unbiased case: ∆ = 0
The T = 0 mean field phase diagram for a plaquette-centered impurity embedded in 'intrinsic' (µ = 0) bilayer graphene with ∆ = 0 is shown in figure 3(a).The phase diagram shares some qualitative features with that of local moment formation in a typical host metal.Namely, there exists a critical ratio of Γ 0 /U before the onset of mean field magnetization and a clear Coulomb staircase in the small Γ 0 /U limit.Despite these similarities, there are two unusual aspects to this phase diagram.We next start by highlighting these novel features and then clarify their physical origin.
(i) As seen from figure 3(a), there is an extreme skewing of the magnetic regime from being centered at (µ − d )/U ∼ 0.5 for small Γ 0 /U to being centered around large positive values of (µ − d )/U with increasing Γ 0 /U .This strong particle-hole asymmetry arises from the fact that the impurity couples asymmetrically to the two layers of BLG, leading to a significant real part of the impurity self-energy Σ d (ω).The effect of which is to strongly renormalize d , which causes the observed skewing.In order to eliminate this large skewing in later plots, we split the real part of the impurity self-energy as and absorb Σ d (0) into the impurity energy, defining a renormalized impurity energy ) then vanishes at ω = 0, and remains small but nonzero away from ω = 0. Plotting the impurity model phase diagram in terms of the renormalized impurity energy ¯ d , to a large degree but not completely, removes the strong particle-hole asymmetry for µ = 0; this can be seen in figure 3(b).Of course, strong particle-hole asymmetry continues to exist away from µ = 0 even after accounting for the impurity energy renormalization, as shown in figure 3(c),(d); this can be ascribed to the particle-hole asymmetry of the BLG dispersion at nonzero µ.
(ii) As seen from figure 3(a), there is a dramatic elongation of the magnetic region to large values of Γ 0 /U ∼ 10, which one can partially attribute the small density of states at µ = 0.However, the phase diagram is also influenced by the wavefunctions of the BLG quasiparticles.A close inspection of the phase boundaries reveals that they are not symmetric about µ = 0 even after accounting for the selfenergy correction discussed above.We understand that this residual particle-hole symmetry breaking arises from the asymmetric broadening of impurity level caused by the disparate effective hybridizations with the different bands.This effect is also seen in the phase diagrams for systems by comparing the µ = 0.05t and µ = −0.05tphase diagrams.While one might naïvely expect that the symmetry between the valence and conduction dispersions would lead to symmetric phase diagrams for positive and negative chemical potential, subtle di fferences between the two regions again reflect the influence of the wavefunctions of the electrons that hybridize with the impurity level.
It is, in fact, extremely instructive to compare the complete impurity phase diagram of bilayer graphene with a fictitious system of electrons obtained by setting φ ks (r) = exp(ik • r)/ √ N s in equation 4, where N s is the total number of sites in the bilayer.These fictitious electrons are chosen to have the same dispersion as the BLG quasiparticles, but their coupling to the impurity does not account for the chirality or the band dependence of the quasiparticle wavefunctions.We find that some of the unusual features of the BLG impurity phase diagram, discussed above, are eliminated upon making this change.
Most noticeably, the phase diagram of the fictitious fermions is not skewed when µ = 0 even when plotted in terms of the unrenormalized impurity energy, indicating that ¯ d = d , so that Σ d (0) = 0.This stems from a symmetry Σ d (−ω) = −Σ * d (ω) in the expression for the self-energy in equation 14 upon assuming a band-independent V ks .Moreover, it also follows that the phase diagram for the fictitious fermions is exactly particle-hole symmetric, in contrast to the case of BLG.
Similar arguments also exactly relate the non-chiral phase diagrams of systems with corresponding chemical potential µ and −µ by noting that the self-energy at finite chemical potential can be obtained by Σ d (ω + µ) of the self-energy at µ = 0. Consequently, the phase diagram of the −µ system is obtained by reflecting the phase diagram of the µ system.This is again in contrast to the phase diagram of BLG where there is no such relation between systems with positive and negative µ.In BLG, particle-hole excitations in a system with positive chemical potential favour different bands than those of a system with negative chemical potential.Since each band has a unique effective coupling to the impurity in BLG, the hybridization of the impurity states will depend on the sign of the chemical potential and so the phase diagrams will be different.Finally, the other major distinction between the finite chemical potential phase diagrams of the two systems is that, once again, the BLG phase diagram is more strongly skewed, even when plotted in terms of ¯ d .This confirms that the band-and momentum-dependence of the hybridization to the BLG quasiparticles is responsible for sizeable shift in the impurity energy via a sizeable real self-energy.

Phase diagram in the biased case: ∆ = 0
We now turn our attention towards a BLG system in a dual-gate configuration.This setup allows one to continuously tune the layer bias and the average chemical potential independently by applying an external electric field perpendicular to the layers.In the presence of a symmetric interlayer bias, the chemical potential remains fixed while a band gap opens in the bulk electronic spectrum of BLG.In the context of local moment formation, this modification to the density of states is expected to substantially change the extent to which an impurity state hybridizes with the BLG electrons.In addition to this, the impurity energy levels also shift up or down depending on the potential of the layer in which it resides.This remarkable ability to alter the energy of an impurity level with respect to the chemical potential through the application of an external electric field is unique to multilayer systems, and has no analog in monolayer graphene.
In the first part of this section, we explore how biasing the layers affects local moment formation by reconstucting phase diagrams similar to those above, but for gated systems with different layer bias and fixed µ = 0. Doing so allows us to identify regions of impurity parameters where local moment formation can be turned on and/or off by the electric field.In the subsequent part of this section, we consider the ability to tune both the chemical potential and bias by constructing alternative phase diagrams where µ and U are varied and it is the bare impurity energy which is fixed.This is again done for a selection of values for the bias.5 is the phase diagram of the impurity model for experimentally accessible values of ∆ plotted in terms of the redefined impurity energy ¯ d introduced above (χ = 0.3t and µ = 0).The bias has two effects on the impurity model: (i) it opens a band gap ∼ ∆ in the BLG dispersion, and (ii) it shifts the impurity energy by ∆/2.Let us discuss, in turn, the impact of these two effects on the phase diagram.

Impurity energy variation Figure
(i) First consider restricting the effect of turning on a bias to opening a gap in the BLG spectrum so that the impurity energy level remains unaltered.Then, the dominant effect of a large ∆ is the elongation of the phase boundary to large Γ 0 /U , regardless of the parity of ∆.This occurs because the bias induces a large band gap and, when the impurity spectral peak lies in this gap, the coupling between the impurity and the extended states becomes negligible because the density of states vanishes.(We have to be careful that the renormalized impurity energy, taking selfenergy corrections into account, should lie in the gap; this renormalization is small if the band gap is large compared to Γ 0 .)Hence, the impurity spectral functions become simple delta functions and if we vary d for fixed Γ 0 /U the local moment phase boundary resembles that of a simple Coulomb staircase in the atomic limit.
(ii) The effect of shifting the energy of the impurity level is similar to the effect of the real part of the self-energy in the phase diagram; it dramatically skews the local moment phase about d /U = 0.5.The direction of the skewing depends on the parity of the bias, as this determines the direction of the impurity energy shift.
If the impurity energy shift and opening of a band gap are taken together, both skewing and elongation of the local moment phase boundary occur.As the electric field is increased from zero to large field strengths, the local moment phase continuously elongates and 'peels' away from the zero bias boundary.Although slight, it is important to note that the phase diagrams with opposite bias parity are not symmetric but have slight differences that arise from the breaking of layer symmetry by the impurity.One of the key new results is the identification of regions in the impurity parameter space where local moments can be turned either on and/or off by adjusting the electric field.The region where local moments survive both in the presence and absence of the electric field are simply where the phases overlap.

Chemical potential variation
Now we explore the possibility of tuning the chemical potential of the system to control local moment formation both in the unbiased and biased cases.To do this, we construct phase diagrams for a given d and ∆, and we now vary µ and U .We do this for ∆ = 0, ±0.2t, for a choice of the bare impurity energy such that the noninteracting impurity spectral peak appears in the midgap when ∆ = −0.2t,which we do by choosing d + ∆/2 = −Σ(ω = 0, ∆).
In figure 6, the phase diagram is plotted in terms of a redefined impurity energy ¯ d = d + Σ(ω = 0, ∆ = 0).It is important to emphasize that the location of the spectral peaks mostly do not correspond to ¯ d .The real part of the self-energy has significant frequency dependence that shifts the location of the spectral peak, whose effect must also be accounted for in order to fully understand the phase diagrams.
When the system is unbiased (i.e.∆ = 0) the impurity energy level lies within the conduction bands (see reference [1] or reference [23] for details on the band structure).In this case, the phase diagram is qualitatively similar to that of a single site impurity (see reference [34]).When a positive bias is in place, ∆ = 0.2t, a band gap opens and the impurity energy shifts deeper into the conduction band.Consequently, the phase boundaries for local moment formation are significantly reduced because of the enhanced broadening due to the increase in the density of states at higher energy in the conduction bands.
The more interesting case is when ∆ = −0.2tand the impurity spectral peak shifts down into the middle of the gap.Then, if U is small enough so that the doubly occupied state also lies at sub-gap energies, both the singly and doubly occupied states can no longer hybridize with the BLG states and the impurity spectral function reduces to delta functions.Hence, we again recover local moment formation very similar to the atomic limit, but now in the large Γ 0 /U limit.However, in this limit Thus, this phase diagram is unusual in the sense that it has two regimes resembling the atomic limit at large and small Γ 0 /U .Separating these regimes is the part of the phase diagram where the doubly occupied state's energy lies beyond the band edge and hybridizes with the conduction states.This occurs at about Γ 0 /U ∼ 2.5 when U ∼ 0.1t, precisely where the unusual 'hump'-like feature is seen in the upper phase boundary.The cause of the feature can again be attributed to level repulsion, as it becomes very strong for states close to the the gap edge and Σ (ω) exhibits a large peak.

RKKY interaction between local moments
In this Section, we the RKKY coupling between local moments [37,38,39] and study how it can be tuned by varying the band gap and chemical potential using a dual-gate configuration.We have seen in the previous section that such variations will, in general, modify the local moment.Here, we focus on changes to the RKKY coupling induced purely by changes in the bulk band structure and filling.
We consider two classical local moments that couple to the set of sites {r} and {r }, respectively, where is the strength of the exchange coupling of an electron's spin, s r , at site r with the magnetic impurity S a .Upon integrating out the itinerant electrons and retaining only those terms that are second order in J (a) r , one obtains a reduced Hamiltonian for the local moments, The coupling J RKKY is given by where m/n are band indices, i/j are the combined sublattice and layer label, n F is the Fermi distribution, and M ij (q) is a matrix describing the Fourier transform between different sites weighted by J (1) r J (2) r .The explicit form of M ij for the case of interest is provided below.
For monolayer graphene, it has been shown that a perturbative treatment in the continuum low-energy theory [13] produces approximate results that match closely with exact diagonalization [14] and lattice Green's functions methods [16], as long as an appropriate high-energy cutoff scheme is applied.In the above perturbative treatment, the entire band structure is used in the calculation so as to avoid any cutoff dependence and the RKKY coupling is accurately reproduced for monolayer graphene.We therefore expect this perturbative calculation to also be a reasonable approach to study the RKKY coupling in BLG in the dual-gate configuration.
We analyzed various moment configurations such as single site (AA, BB, AB) and plaquette coupled moments both along the zigzag and armchair directions (see figure for labelling conventions).The effects of varying the chemical potential and layer bias were seen to be qualitatively similar for each case, so we have chosen to present only the results for plaquette centered moments that lie along the zigzag direction.The impurity atom is taken to lie above an A 2 in the center of a hexagonal plaquette in layer 1.For simplicity, we assume the coupling to each of the seven sites is equal so that J (1) r = J (2) r ≡ J, although we have checked that the results are qualitatively unaffected if there is an unequal coupling to the site below the impurity on the other layer (A 2 ).In this case, the components of M ij are M A1A1 = J 2 4 + 2 (cos(q a ) + cos(q b ) + cos(q a − q b )) M A1B1 = J 2 2 + 2 cos(q a ) + cos(q b ) + e i(qa−q b ) + e −i(2qa−q b ) M A1A2 = J 2 1 + e iq b + e i(qa−q b ) M B1A2 = J 2 e iq b + e i(qa−q b ) + e i(qa−2q b ) where q a = q • a, q b = q • b and a = x and b = x/2 + √ 3ŷ/2.To demonstrate the ability to tune the RKKY interaction using the dual-gate configuration, the J RKKY coupling, normalized to its value at µ = 0 and ∆ = 0, is plotted in figure 8 as a function of the interlayer bias ∆ for two moments separated by 10 lattice spacings.This is done for µ = 0 and µ = 0.05t.For experimental considerations, one must keep in mind that the RKKY coupling is quite small in BLG.As an example, a bare exchange term equal to J (1) = J (2) = 0.2t produces an effective coupling J RKKY = 1.3 × 10 −4 t (∼ 4.4 K) at 4 lattice spacings, and just J RKKY = 7.5 × 10 −6 t (∼ 0.3 K) at 10 lattice spacings.However, similar to monolayer graphene, electron interactions are expected to make the coupling strength more long ranged [15].At shorter distances, the RKKY interaction is enhanced, but the tunability is reduced.
Before describing the tunable features of the RKKY coupling, it is important to first understand that the wavefunctions of a given band are sensitive to the parity of the bias between the layers, even though the dispersion is not.Their dependancy on the parity can significantly influence how J RKKY changes with bias, as explained below.When a positive bias is present, states in the upper two bands are more heavily weighted to layer 1 sites, while states in the lower band are more heavily weighted to the layer 2 sites.This weighting is reversed when the parity of the bias is negative.In contrast, when there is no bias the weighting of the wavefunction is the same for each layer.
With this background, it is possible to explain the symmetry/asymmetry between the two curves.When µ = 0, the chemical potential lies between the valence and conduction bands, and so particle-hole excitations can only occur between them.This corresponds to one of the states being localized to layer 1 and the other localized to layer 2. It follows that the coupling strength J RKKY is parity invariant and so it is symmetric for positive and negative biases.
In contrast, when µ = 0, the coupling is sensitive to the parity of the bias.If µ = 0.05t, ∆ > 0 and µ is less than the band gap, the chemical potential lies in the third band where the states tend to localize to layer 1, the layer in which the moments reside.The finite chemical potential causes some of the particle-hole excitations between the lower and upper two bands to be suppressed by Pauli-blocking, and also introduces low-energy excitations between the two upper bands where the wavefunctions are weighted to layer 1.If however, µ = 0.05t, ∆ < 0 and µ is less than the band gap, the chemical potential lies in the third band, but now these states tend to localize to layer 2. Although the energetics of the scattering processes remain the same, the matrix elements do not.The excitations between the upper two bands now have matrix elements whose weighting on layer 1 is much less.Thus, the coupling is dependent on the relative parity o f µ to ∆.Hence, if we consider a system with µ = −0.05t, the J RKKY curve will be reflected about ∆ = 0.
In addition to the effects described above, the density of states about the chemical potential tends to increase for small biases, as the band edge flattens and is pushed closer to the chemical potential.At large bias strengths, the dispersion close to the band-edge resembles that of a 'mexican-hat', leading to further complexity in the density of states.Furthermore, at finite chemical potential, the Fermi points extend out to form a Fermi-surface symmetric about the K-points.The combination of all these effects lead to the non-trivial changes seen the RKKY coupling in figure 8. Interestingly, at a distance of 10 lattice spaces, the coupling remains antiferromagnetic when µ = 0.05t and ∆ > 0 and tends to increase with bias strength.However, when ∆ becomes increasingly negative, the antiferromagnetic coupling strength is reduced to zero then switches to a ferromagnetic coupling about ∆ ∼ −0.12.For the case of µ = 0.05t considered here, the ability to fully turn off the coupling and switch the sign of J RKKY with an applied electric field sets in when the moments are separated by at least 9 lattice spacings and persists to about 20 lattice spacings.This window of tunability may be augmented by carefully adjusting µ.Regardless of the sign of the bias, once the band gap exceeds the chemical potential, the µ = 0.05t curve begins to merge with the µ = 0 curve, as expected at low temperature.
This ability to dynamically tune the strength of the RKKY interaction, as well as its sign from being antiferromagnetic to being ferromagnetic, for two fixed moments using just an electric field (rather than doping) is perhaps the most interesting feature of this system.

Discussion
In this paper, we have discussed the physics of local moments on bilayer graphene.We have, in our discussion, ignored the effects of electron-electron interactions among the BLG electrons.These are known to be important for quadratic band touching [40,41,42,43,44], leading to symmetry breaking and many-body gaps for zero doping and zero electric field.However, so long as there is nonzero doping or the presence of an electric field that gaps out the low-energy BLG states, the perturbative effects of electron-electron interactions are benign and will not lead to qualitatively new many-body effects.Using a clean substrate may also be a viable route to mitigating the effects of rippling and disorder.Moreover, the substrate may partially screen the BLG electron-electron interactions and reduce the many-body gap recently observed in suspended BLG [45].The competition between impurity physics and many-body interactions in BLG deserves a careful separate investigation.The ability to turn on/off local moments, and the ability to tune the sign and magnitude of the RKKY coupling between local moments using electric fields perpendicular to the bilayer, which we have studied, constitutes physics beyond what has been discussed for monolayer graphene.The experimental realization of such tunable local moments in BLG is a compelling prospect.It would be interesting to study such a system using scanning tunnelling spectroscopy and to probe the quantum dynamics of interacting local moments in experiments.We expect that thermal fluctuations will only slightly alter the phase boundaries as long as the temperature is below the Hubbard gap.Quantum fluctuations are expected to lead to Kondo screening or valence fluctuations in points of the phase diagram -this is an interesting direction for future research.

Figure 1 .
Figure 1.(a) Bilayer graphene in a dual-gate configuration.(b) Schematic diagram of a plaquette-centered (large, red) adatom impurity on the top layer of bilayer graphene.(c) Cross-section of the dispersion relation for unbiased (Left) and biased (Right) graphene close along ky = 0 through the K-point (∆ = 0 and ∆ = 0.025t, respectively).Inset: The two unique K-points and the cross-sectional cut are indicated in the Brillioun zone.

Figure 4 .
Figure 4.The phase diagram of local moment formation for fictitious fermions with the same dispersion as bilayer graphene for (above) µ = 0, and (below) µ = −0.05t,plotted in terms of ¯ d .(Note, ¯ d = d when µ = 0.) The phase diagram for µ = 0.05t is related to that of µ = −0.05by a reflection about µ − ¯ d /U = 0.5.

Figure 6 .
Figure 6.Local moment phase diagram for biased bilayer graphene.The impurity level at ∆ = −0.2tbias was chosen so that its spectral peak lies in the middle of the gap.

Figure 7 .
Figure 7. Crystal structure of bilayer graphene and site labelling convention.The primitive lattice vectors are a and b and the armchair and zigzag directions are indicated by the arrows.The local moments considered here are plaquette centered and reside on the same layer above an A 2 atom.

Figure 8 .
Figure 8. Normalized RKKY coupling strength between to classical plaquette centered moments at a distance of 10 lattice spacings along the zigzag direction.Both moments are located on the same layer and the chemical potential is chosen to be µ = 0 and µ = 0.05t at a temperature T = 0.002t (∼ 70 K).