Focusing RKKY interaction by graphene P-N junction

The carrier-mediated RKKY interaction between local spins plays an important role for the application of magnetically doped graphene in spintronics and quantum computation. Previous studies largely concentrate on the influence of electronic states of uniform systems on the RKKY interaction. Here we reveal a very different way to manipulate the RKKY interaction by showing that the anomalous focusing - a well-known electron optics phenomenon in graphene P-N junctions - can be utilized to refocus the massless Dirac electrons emanating from one local spin to the other local spin. This gives rise to rich spatial interference patterns and symmetry-protected non-oscillatory RKKY interaction with a strongly enhanced magnitude. It may provide a new way to engineer the long-range spin-spin interaction in graphene.


I. INTRODUCTION
Graphene [1], a single atomic layer of graphite, is featured by ultra-high mobility and electrical tunability of carrier density and hence provides an attractive platform for studying the unique electron optics of Dirac fermions owing to its gapless and linear dispersion. Cheianov et al. [2] proposed the interesting idea that an interface between electron (N)-doped and hole (P)-doped regions in graphene can focus an electron beam, which may lead to the realization of an electronic analog of the Veselago lens in optics [3][4][5][6]. This anomalous focusing effect has motivated many new ideas and device concepts [7][8][9][10][11]. Very recently, this effect was observed experimentally [12,13], which paves the way for realizing electron optics based on graphene P-N junctions (PNJs). A common feature of these works is that they concentrate on focusing the electrons themselves, leaving its potential applications to other fields unexplored.
In this work, we explore a very different direction by showing that the anomalous focusing effect could be utilized to manipulate the carrier-mediated Rudermann-Kittel-Kasuya-Yosida (RKKY) interaction [14][15][16] between magnetic moments (spins), with potential applications in spintronics [17][18][19], scalable quantum computation [20,21], and majorana fermion physics [22] as the RKKY interaction enables longrange correlation between spatially separated local spins [23][24][25][26][27], a crucial ingredient in these developments. In recent years, a lot of efforts have been devoted to characterizing the RKKY interaction in different 2D uniform systems such as two-dimensional electron gases [28,29], graphene [30][31][32][33][34][35][36], and the surface of topological insulators [37][38][39]. There are also many interesting schemes to manipulate the RKKY interaction [40][41][42][43][44][45][46][47][48]. Recent experimental advances further enable * wenyang@csrc.ac.cn † kchang@semi.ac.cn the RKKY interaction to be mapped out with atomic-scale resolution from single-atom magnetometry using scanning tunnelling spectroscopy [49][50][51]. However, in a d-dimensional uniform system, the RKKY interaction usually decays at least as fast as 1/R d . This rapid decay -a common feature of the previous studies mentioned above -makes the RKKY interaction very short-ranged and may hinder its applications. This motivates growing interest in modifying the long-range behavior of the RKKY interaction, e.g., the 1/R 3 long-range decay in undoped graphene can be changed by thermal excitation [52] and even be slowed down by electron-electron interactions [53]. Here we show that the graphene PNJ allows the diverging electron beams emanating from one local spin to be refocused onto the other local spin, thus the electronmediated RKKY interaction between these two local spins can be strongly enhanced and tuned beyond the 1/R d limit of noninteracting uniform systems. The graphene PNJ also gives rise to symmetry-protected non-oscillatory RKKY interaction as a function of the distance, in sharp contrast to the "universal" oscillation of the RKKY interaction in uniform systems with a finite carrier concentration. This may provide a new way for engineering the correlation between spatially separated local spins for their applications in spintronics and quantum computation.
This paper is organized as follows. In Sec. II, we explain intuitively how to utilize the graphene PNJ to manipulate the RKKY interaction and highlight the symmetry-protected non-oscillatory RKKY interaction. Then in Sec. III, we perform numerical simulation based on the tight-binding model to demonstrate this all-electrical manipulation and discuss the experimental feasibility. Finally, we present a brief summary in Sec. IV.

II. RKKY INTERACTION IN GRAPHENE P-N JUNCTION: PHYSICAL PICTURE
Let us consider two local spinsŜ 1 (located at R 1 ) andŜ 2 (located at R 2 ) coupled to the spin densityŝ(x) ≡ŝδ(r − x) of itinerant carriers via the exchange interactionV ex = −J 0Ŝ1 · s(R 1 ) − J 0Ŝ2 ·ŝ(R 2 ), whereŝ (r) is the carrier spin (position) operator. The total Hamiltonian of the coupled system is the sum ofV ex and the carrier HamiltonianĤ. The carriermediated RKKY interaction originates from the local excitation of carrier spin density fluctuation by one local spin and its subsequent propagation to the other spin. At zero temperature, the effective RKKY interaction between the local spins is obtained by eliminating the carrier degree of freedom through second-order perturbation theory as [14][15][16]54] (1) E F is the Fermi energy of the carriers, Tr traces over the carrier spin, andĜ(r, r 0 ; E) ≡ r|(E + i0 + −Ĥ)|r 0 is the unperturbed (i.e., in the absence of the local spins) propagator of the carriers in real space. In general,Ĝ(r, r 0 ; E) is still an operator acting on the carrier spin degree of freedom. In a ddimensional uniform system, the carriers excited by the first local spin at R 1 propagate towards the second local spin at R 2 in the form of an outgoing waveĜ( k is a characteristic wave vector of the carriers with energy E, and the denominator R (d−1)/2 ensures the conservation of probability current. The integration over the energy in Eq. (1) yields another factor 1/R from the oscillating phase factor e ikR , so J αβ ∝ 1/R d . This provides a rough explanation for the "universal" 1/R d decay of the RKKY interaction, as discovered in a great diversity of materials by previous studies. It also reveals a very different way -tailoring the carrier propagation and interference -to manipulate the RKKY interaction beyond this constraint, as opposed to previous studies that exploit the electronic states and energy band structures of different uniform materials. The anomalous focusing effect in graphene PNJs [2] provides a paradigmatic example for this manipulation.
As shown in Fig. 1(a), the honeycomb lattice of graphene consists of two sublattices (denoted by A and B) and each unit cell contains two carbon atoms (or π z -orbitals), one on each sublattice. Let us use R to denote the location of each carbon atom, |R for the corresponding orbital, and s R (= A or B) for the sublattice on which R locates. For carriers in the graphene PNJ, the tight-binding Hamiltonian is the sum ofĤ 0 for uniform graphene andV J for the on-site junction potential: where R, R ′ denotes nearest neighbors, t ≈ 3 eV is the nearest-neighbor hopping [1], and V R is equal to −V 0 (+V 0 ) when R locates in the left (right) of the junction [shaded stripe in Fig. 1(a)] with V 0 ≥ 0. As shown in Fig. 1(b), the zero point of energy is chosen such that the Dirac point in the left (right) of the junction lies at −V 0 (+V 0 ). Uniform graphene corresponds to V 0 = 0, while nonzero V 0 corresponds to a junction, e.g., N-N (P-P) junction corresponds to E F > V 0 (E F < V 0 ). Here we consider the P-N junction (PNJ), corresponding to E F ∈ [−V 0 , +V 0 ]. In uniform graphene (V 0 = 0), the doping is determined by E F . In graphene PNJ, the electron doping in the N region (left) is V 0 + E F , while the hole doping in the P region (right) is V 0 − E F , e.g., E F = 0 corresponds to the electron doping in the N region being equal to the hole doping in the P region. Due to the absence of spin-orbit coupling in the carrier HamiltonianĤ, the carrier-mediated RKKY interaction at zero temperature between one local spinŜ 1 at R 1 (sublattice s R 1 ) in the N region and another local spinŜ 2 at R 2 (sublattice s R 2 ) in the P region [see Fig. 1(a)] assumes the isotropic Heisenberg form [33,34]:Ĥ RKKY = JŜ 1 ·Ŝ 2 , where the range function is determined by the unperturbed propagator (i.e., in the absence of the local spins) of the carriers from R 1 to R 2 : In arriving at Eq. (3), we have used G(R 2 , R 1 , E) = G(R 1 , R 2 , E) due to the time-reversal invariance of the graphene HamiltonianĤ. Note that the propagator G(R 2 , R 1 , E) and hence the RKKY range function J are very sensitive to the sublattices on which R 1 and R 2 locate (i.e., s R 1 and s R 2 ), e.g., when R 2 moves from an atom on the A sublattice (s R 2 = A) to a neighboring atom on the B sublattice (s R 2 = B) [see Fig. 1(a)], the propagator and hence the RKKY interaction may change significantly. Now we discuss how the anomalous focusing effect in graphene PNJs [2] can be utilized to manipulate the carrier propagator and hence the RKKY interaction beyond the "universal" 1/R d long-range decay as encountered in previous studies. Ever since the pioneering work of Cheianov et al. [2], there have been many studies on the anomalous focusing effect, either based on the classical analogy to light propagation in geometric optics or based on the scattering of the electron wave functions. Below, we provide a physically intuitive analysis on how the graphene PNJ focuses the carrier propagator and hence the RKKY interaction based on the continuum model of graphene [1]. The purpose is to provide a qualitative picture for the focusing of the RKKY interaction and establish its effectiveness for an arbitrary direction of the P-N interface.
A. Anomalous focusing of carrier propagator: continuum model The low-energy physics of graphene is described by two Dirac cones located at K and K ′ = −K, which form a Kramer pair. For clarity, we first analyze the focusing effect based on the K-valley continuum model, leaving the discussion including both valleys to the end of this subsection. Using the band-edge Bloch functions |Φ K,A (from π z -orbitals on the A sublattice) and |Φ K,B (from π z -orbitals on the B sublattice) of the K valley as the basis, the continuum model for the K valley reads [1],ĥ wherep is the momentum relative to the K valley and v F is the Fermi velocity. Note that the continuum model regards the two atoms of the same unit cell to locate at the same spatial point, so each spatial point contains two sublattices/orbitals and the Hamiltonianĥ is a 2 × 2 matrix. Correspondingly, the carrier propagator from R 1 to R 2 is also a 2×2 matrix: For uniform graphene, the K-valley continuum model [Eq. (4) with V 0 = 0] leads to a massless Dirac spectrum E ± (q) ≡ ±v F |q| and chiral eigenstates |±, q = e iq·r |u ± (q) , where |u ± (q) is the two-component spinor for the sublattice degrees of freedom. The conduction (valence) band state The RKKY interaction is usually dominated by the contributions from carriers near the Fermi surface [the energy integral in Eq. (3) merely produces a multiplicative factor ∝ 1/R], so we focus on the carrier propagator on the Fermi level E F . For E F > 0, the Fermi momentum is q F ≡ E F /v F , and the right-going eigenstates |+, q on the Fermi contour are characterized by the momentum q ≡ ((q 2 F − q 2 y ) 1/2 , q y ). The 2×2 propagator from R 1 to R 2 (on the right of R 1 ) in uniform graphene can be expressed in terms of these eigenstates as The RKKY interaction in uniform graphene is given by Eq.
For the graphene PNJ described by the K-valley continuum model in Eq. (4), the first local spin S 1 in the N region excites a series of outgoing plane wave eigenstates on the Fermi contour, but only the right-going eigenstates, i.e., |+, q N with momentum q N ≡ ((q 2 N − q 2 y ) 1/2 , q y ), can transmit across the P-N interface, becomes a right-going eigenstate |−, q P with momentum q P ≡ (−(q 2 P − q 2 y ) 1/2 , q y ) on the Fermi contour of the P region, and finally arrive at Fermi momenta in the N and P regions, respectively. In terms of these local, right-going eigenstates on the Fermi contours, the carrier propagator from R 1 to R 2 is (see Appendix A): where w(q y ) is the transmission coefficient of the incident state |+, q N across the PNJ. The RKKY interaction in the graphene PNJ is given by Eq.
The key difference between the carrier propagator in the graphene PNJ [Eq. (6)] and that in uniform graphene [Eq. (5)] is the change of the propagation phase factor from e iq·(R 2 −R 1 ) ≡ e iφ 0 (q y ) for uniform graphene to e i(q P ·R 2 −q N ·R 1 ) ≡ e iφ NP (q y ) for the graphene PNJ. This change is responsible for the anomalous focusing of the diverging carrier wave into a converging one. For an intuitive analysis of this behavior, we discretize the q y axis into grids m∆ (m ∈ Z), where the spacing ∆ ≪ size of the graphene Brillouin zone. Then Eq. (6) gives g(R 2 , R 1 , E) = m g (m) (R 2 , R 1 , E) and g (m) is the contribution from the q y integral over the mth segment [(m − 1/2)∆, (m + 1/2)∆], corresponding to a wave packet characterized by the center momentum q y = m∆. In other words, the entire propagator is the sum of contributions from all these wave packets characterized by different center momenta q y 's on the Fermi contour. The same analysis is applicable to the propagator g uniform in uniform graphene [Eq. (5)]. Since the integrand is the product of a slowly-varying part and a rapidly oscillating propagation phase factor, the contribution from a given wave packet characterized by the center momentum q y is appreciable only when the propagation phase is stationary: ∂ q y φ 0 (q y ) = 0 (for uniform graphene) or ∂ q y φ NP (q y ) = 0 (for graphene PNJ). This first-order stationary phase condition determines the most probable (or classical) trajectory of a wave packet emanating from R 1 .
In uniform graphene, the classical trajectory of a given wave packet characterized by the center momentum q = ((q 2 F − q 2 y ) 1/2 , q y ) on the Fermi contour [Eq. (5)] is a beam emanating from R 1 and going along the wave vector q. The classical trajectories of different wave packets on the Fermi contour form many outgoing beams emanating from R 1 , which where R 2 is on the cusp of the caustics. The PNJ potential V 0 = t/2 for (a) and t/5 for (b) and (c).
manifests the diverging propagation of the carriers in uniform graphene and leads to g uniform ∝ 1/R 1/2 . By contrast, in the graphene PNJ, the classical trajectory of a given wave packet characterized by the center momentum q N = ((q 2 N − q 2 y ) 1/2 , q y ) with incident angle θ N = tan −1 (q y /q N,x ) consists of the incident beam along q N , the reflection beam with a reflection angle θ N , and the refraction beam with a refraction angle θ P ≡ tan −1 (q y /q P,x ), as sketched in Fig. 1(a) and further visualized in Fig. 2(a). Here the refraction angle θ P is determined by the Snell law sin θ N = n sin θ P with a negative effective Let us use (X 1 , Y 1 ) and (X 2 , Y 2 ) to denote the Cartesian coordinates of R 1 and R 2 , respectively. In the graphene PNJ, when R 2 locates on the caustics [2] the wave packet going from R 1 to R 2 obeys not only ∂ q y φ 0 (q y ) = 0, but also ∂ 2 q y φ 0 (q y ) = 0, so that its contribution to the propagator g(R 2 , R 1 , E F ) is enhanced. The most interesting case occurs at E F = 0 or equivalently n = −1.
In this case, we have q N,x = −q P,x , thus for R 2 at the mirror image of R 1 about the PNJ, i.e., R 2 = R m 1 ≡ (|X 1 |, Y 1 ), the phase φ NP (q y ) vanishes for all q y , so that the integrand in Eq. (6) no longer suffers from the rapidly oscillating phase factor e iφ NP (q y ) . This corresponds to constructive interference of all the transmission waves at R m 1 or equivalently perfect focusing of the diverging electron beams emanating from R 1 onto R m 1 [2]. This not only lead to strong local enhancement of the propagator g(R 2 , R 1 , E F = 0) when R 2 locates in the vicinity of R m 1 [see Fig. 2(b)], but also makes g(R m 1 , R 1 , E F = 0) independent of the distance R [black, solid line in Fig. 2(c)], in sharp contrast to the 1/R 1/2 decay in uniform graphene. This distance independent propagator can be attributed to the existence of a hidden symmetry on the Fermi contours of the host material [55].
When E F 0, the anomalous focusing locally enhances the propagator from R 1 to its caustics and slows down its decay with the distance to a slower rate ∼ 1/R ξ (ξ < 1/2) compared with the 1/R 1/2 decay in uniform graphene, as a consequence of imperfect focusing away from n = −1, e.g., for R 2 on the cusp (|nX 1 |, 0), we have ξ ≈ 0.24 nearly independent of E F , as shown in Fig. 2(c). By contrast, for R 2 far away from the caustics, the propagator recovers the 1/R 1/2 decay of uniform graphene. In addition to locally enhancing the propagator, the PNJ also slightly decreases the propagation amplitude via the finite transmission w(q y ). However, this effect is of minor importance because Klein tunneling [56] allows carriers with a small incident angle θ N to go through the PNJ almost completely, as demonstrated by the weak reflection in Fig. 2(a) when the incident angle is small.
From the above analysis, it is clear that the PNJ qualitatively changes the diverging spherical propagation of carriers in graphene into a converging one. This greatly enhances the propagator near the caustics, so that its decay with interspin distance R slows down (for E F 0) and even ceases (for E F = 0). At large distances, the non-decaying propagator could lead to 1/R decay of the RKKY interaction between two mirror symmetric spins about the PNJ, as we demonstrate shortly (see Sec. IIIA).
Before concluding this subsection, we emphasize that the above analysis is based on the K-valley continuum model, featured by a single Dirac cone and a circular Fermi contour. This simplified model ignores two important effects: the trigonal warping at high energies and the presence of another Dirac cone at K ′ = −K. At high Fermi energies |E F | ∼ t, the trigonal warping leads to non-circular Fermi contours, while the inter-valley scattering may decrease the transmission probability across a sharp PNJ [57]. The former makes the focusing effect no longer perfect even when E F = 0, while the latter decreases the carrier propagator and hence the RKKY interaction. Even in the linear regime |E F | ≪ t, the presence of two inequivalent valleys still gives rise to inter-valley interference that significantly affect the RKKY interaction, as discussed by Sherafati and Satpathy [33,34] for uniform graphene. Here we briefly discuss this issue for the graphene PNJ. First, we assume that the P-N interface does not induce inter-valley scattering. Then, when both K and K ′ valleys are included, the 2×2 carrier propagator from R 1 to R 2 would be where g (g ′ ) is the propagator in the presence of the K (K ′ ) valley alone [see Eq. (6) for the expression of g]. The K ′ -valley continuum modelĥ ′ = −v Fσ * ·p + sgn(x)V 0 is the time reversal of the K-valley model [Eq. (4)], so that g ′ s 2 s 1 (R 2 , R 1 , E) = g s 1 s 2 (R 1 , R 2 , E), i.e., the propagator of a K ′ -valley electron from the sublattice s 1 at R 1 to the sublattice s 2 at R 2 is equal to the propagator of a K-valley electron from the sublattice s 2 at R 2 back to the sublattice s 1 at R 1 . This also ensures the time-reversal invariance of the total propagator: G s 2 s 1 (R 2 , R 1 , E) = G s 1 s 2 (R 1 , R 2 , E). The RKKY interaction is given by Eq. (3) with G(R 2 , R 1 , E) replaced by the (s R 2 , s R 1 ) matrix element of G(R 2 , R 1 , E). Consequently, the RKKY interaction consists of the intra-valley contributions and the inter-valley interference term. The former oscillates slowly on the length scale of the Fermi wave length, while the latter oscillates rapidly as e i(K−K ′ )·(R 2 −R 1 ) on the atomic scale, similar to the case of uniform graphene [33,34]. In the presence of inter-valley scattering by the P-N interface, a quantitative description is very difficult within the continuum model, but we still expect the contribution from the inter-valley interference to be rapidly oscillating on the atomic scale. Therefore, the slowly-varying envelope of the RKKY interaction is always determined by the intra-valley contributions, which are independent of the direction of the P-N interface with respect to the crystalline axis of graphene. In other words, the continuum model suggests that the focusing of the RKKY interaction should occur for an arbitrary direction of the P-N interface, as confirmed by our subsequent numerical calculation based on the tight-binding model.

B. Symmetry-protected non-oscillatory RKKY interaction
For uniform graphene, the tight-binding HamiltonianĤ 0 [see Eq. (2a)] possesses electron-hole symmetryPĤ 0P −1 = −Ĥ 0 . [31,36,58], whereP inverts all the π z -orbitals on sublattice B but keeps all π z -orbitals on sublattice A invariant, i.e., P|R = ±|R , with the upper (lower) sign for s R = A (s R = B). For undoped graphene, this makes the RKKY interaction between local spins on the same (opposite) sublattice always ferromagnetic (antiferromagnetic), irrespective of their distance. However, once the graphene is doped, the RKKY interaction recovers its "universal" oscillation with a characteristic wavelength λ F /2 (λ F is the Fermi wavelength) between ferromagnetic and anti-ferromagnetic couplings, as also found in many other materials.
For the graphene PNJ, the presence of the junction potential breaks the electron-hole symmetry of uniform graphene. Moreover, since both the N region and the P region are doped, the RKKY interaction is also expected to oscillate between ferromagnetic and antiferromagnetic couplings with the distance. Interestingly, we find that the electron-hole symmetry can be restored under certain conditions. Let us consider a general graphene PNJ described by the tight-binding Hamiltonian Eq. (2) with a general on-site junction potential V R and define the mirror reflection operatorM that maps the π zorbital |R to another π z -orbital |R m at the mirror image location R m about the P-N interface, i.e.,M|R = |R m . The key observation is that as long as the mirror reflectionM about the P-N interface keeps the graphene lattice invariant but inverts the junction potential (i.e., V R = −V R m ), the PNJ HamiltonianĤ possesses a generalized electron-hole symmetry: (PM)Ĥ(PM) −1 = −Ĥ. This ensures the eigen-energies of the PNJ to appear in pairs (ε, −ε) and the corresponding eigenstates |φ ε and |φ −ε obey |φ −ε =PM|φ ε , similar to the electron-hole symmetry in uniform graphene [31,36,58]. As a consequence of this symmetry, when E F = 0, the Matsubara Green's function G(R, R ′ , τ) [59] obeys G(R 1 , R m 1 , −τ) = ±G(R m 1 , R 1 , τ), with the upper (lower) sign for R 1 and R m 1 on the opposite (same) sublattices. The time-reversal symmetry further dictates the Matsubara Green's functions to be real. According to the imaginary-time formalism for the RKKY interaction [31,36,58] [equivalent to the real-time formalism in Eq. (3)], the sign of the RKKY range function is determined by G (R 1 , R 2 R 1 , τ). Therefore, when the two local spins are mirror symmetric about the PNJ, i.e., R 2 = R m 1 , their RKKY interaction is always ferromagnetic (antiferromagnetic) on the same (opposite) sublattices, irrespective of their distance. If the P-N interface is along the zigzag direction, then R 1 and R m 1 are always on opposite sublattices, so the RKKY interaction is antiferromagnetic. If the P-N interface is along the armchair direction, then R 1 and R m 1 are always on the same sublattice, so the RKKY interaction is ferromagnetic.

III. NUMERICAL RESULTS
Here we calculate the RKKY interaction in the graphene PNJ numerically based on the tight-binding model [Eq. (2)], with the on-site junction potential V R = −V 0 (V R = +V 0 ) in the N (P) region. For convenience, we introduce the dimensionless RKKY range function J ≡ tJ/J 2 0 . Due to the transformationĤ → −Ĥ upon (V 0 , t) → (−V 0 , −t) and the time-reversal symmetry, J is invariant upon (E F , V 0 ) → (−E F , −V 0 ) (see Appendix B), so we need only consider V 0 > 0. Our results show that for different orientations of the PNJ (e.g., along the zigzag direction, the armchair direction, and a slightly misaligned direction) and different sublattice locations of R 1 and R 2 , the RKKY interaction exhibits similar anomalous focusing behaviors, consistent with our previous analysis based on the continuum model in Sec. IIA. For specificity, we present our results for a PNJ along the zigzag direction and always take R 1 on the A sublattice and R 2 on the B sublattice.

A. Anomalous focusing of RKKY interaction
For the first local spinŜ 1 fixed at R 1 = (−91a, 0) [a is the C-C bond length, see Fig. 1(a)], the spatial map of the scaled range function JR 2 /a 2 as a function of the location R 2 = (X 2 , Y 2 ) of the second local spinŜ 2 is shown in Fig.  3(a) for uniform graphene and Fig. 3(b)-(d) for graphene PNJ. Here we follow Ref. 60 and use the multiplication factor R 2 /a 2 to remove the intrinsic decay ∝ 1/R 2 of the RKKY interaction in uniform graphene [30,34]. This helps us to present an overall view of the spatial texture of the RKKY interaction over both the N region and the P region in a single contour plot and highlights the focusing of the RKKY interaction by the P-N interface. For example, in uniform graphene with V 0 = 0 and E F = 0.2t [ Fig. 3(a)], the scaled range function does not decay, manifesting the intrinsic 1/R 2 decay of the RKKY interaction. By contrast, in the graphene PNJ [ Fig.  3(b)-(d)], the P-N interface induces strong local enhancement of the RKKY interaction in the P region, but it has a negligible influence in the N region. This is an obvious consequence of anomalous focusing: in the N region, the carrier propagation remains diverging, similar to uniform graphene shown by Fig.  3(a); while in the P region, the carrier wave is refocused by the PNJ [2]. For E F = 0 in Fig. 3(b), corresponding to an effective refraction index n = −1, the RKKY interaction is significantly enhanced whenŜ 2 locates near the mirror image ofŜ 1 about the PNJ. For E F = 0.02t in Fig. 3(c) [E F = −0.02t in Fig.  3(d)], corresponding to n ≈ −0.82 (n ≈ −1.2), the maximum of the RKKY interaction shifts towards (away from) the PNJ, consistent with the shift of the caustics [see Eq. (7)].
Now we discuss two fine features in Fig. 3. First, for uniform graphene, Fig. 3(a) reproduces the C 3v spatial symmetry at short distances [60], the slow oscillations with a characteristic wavelength λ F /2 [30,34,60], and the rapid oscillations on the atomic scale due to the inter-valley interference [31, 34, 60], as described by the dimensionless range function at large distances (q F R ≫ 1) [34]: where q F is the Fermi momentum and θ R is the angle between R ≡ R 2 − R 1 and K − K ′ . For the graphene PNJ in Fig. 3(b)-(d), the RKKY interaction also consists of a slowlyvarying envelope and a rapidly-varying part that oscillates on the atomic scale. As discussed at the end of Sec. IIA, the former comes from the intra-valley contributions, while the latter comes from the inter-valley interference and hence oscillates with a momentum K − K ′ , similar to Eq. (8). Notice that in Figs. 3(a)-(d), the most rapid atomic-scale oscillation occurs along the y axis [i.e., the zigzag direction, see Fig. 1(a)], consistent with previous studies in uniform graphene [33,34].
Second, in Fig. 3(b)-(d), there is no local enhancement of the RKKY interaction near the P-N interface. According to Eq. (3), this manifests the fact that there is no local charge accumulation near the P-N interface, since the incident electron wave emanating from R 1 either reflects back or transmits through the P-N interface. The interference between the incident wave and the reflection wave in the N region (near the P-N interface) can be clearly seen by comparing Fig. 3(b)-(d) to Fig. 3(a).
Let us consider the RKKY interaction between two mirror symmetric spins in graphene PNJ at E F = 0, i.e., electron doping V 0 in the N region and hole doping V 0 in the P region. According to the symmetry analysis in Sec. IIB, for the PNJ with its interface along the zigzag direction, the RKKY interaction between two mirror symmetric spins is always antiferromagnetic. This feature is demonstrated by Fig. 4. As shown in Fig. 4(a), in the linear regime (V 0 ≪ t), the RKKY interaction J ∝ V 2 0 increases quadratically with V 0 . As shown in Fig. 4(b), the scaled RKKY interaction strength JR/a is Rindependent at large distances. This indicates that J follows 1/R asymptotic decay due to the perfect refocusing, in sharp contrast to the 1/R d asymptotic decay in a great diversity of doped d-dimensional materials, as well as the 1/R 3 asymptotic decay in undoped graphene. Therefore, the RKKY interaction at V 0 ≪ t and R ≫ λ F can be well approximated by the analytical expression J ≈ 0.012(V 0 /t) 2 /(R/a), where the constant 0.012 is obtained by fitting the data in Fig. 4(a)-(b). For com-parison, in uniform graphene with the same doping level as the PNJ, the envelope of the RKKY interaction [see Eq. (8)] has the asymptotic form J uniform ≈ 0.037(V 0 /t)/(R 2 /a 2 ). The RKKY interaction in the graphene PNJ differs qualitatively from that in uniform graphene in the scaling with both the distance (J uniform ∝ 1/R 2 vs. J ∝ 1/R) and the junction potential V 0 or equivalently carrier concentration (J uniform ∝ V 0 vs. J ∝ V 2 0 ). Since the localized spins are usually fixed, dynamic tuning of the PNJ by electric gating potentially allows for selective control of localized spins, an important ingredient for spin-based quantum computation.

B. Experimental feasibility and generalization to other materials
Since the focusing of the RKKY interaction is dominated by the contribution from the electron states near the Fermi level, and these states cannot not "feel" any potential variation on the length scale ≪ Fermi wavelength λ F , the finite width of the PNJ has a small influence as long as it is much smaller than λ F . Taking V 0 = 0.1t and the first local spin at R 1 = (−151a, 0) as an example, using the experimentally fabricated linear PNJ [61] of width 29a ≈ 4.1 nm, instead of a sharp PNJ, only reduces the magnitude of the RKKY interaction by ∼ 25% without changing the 1/R asymptotic scaling.
In the presence of a finite gap ∆ (e.g., due to substrate mismatch [62]) in the Dirac spectrum of graphene, as long as ∆ ≪ |V 0 |, the gap does not significantly influence the states near the Fermi level, which dominates the anomalous focusing effect. This has been confirmed by our numerical calculation using E F = 0, V 0 = 0.2t, and a typical gap ∆ = 0.03t: no appreciable change of the focusing behavior occurs. As a matter of fact, from our analytical analysis following Eq. (6), it is clear that perfect focusing of the PNJ with E F = 0 essentially arises from the opposite momenta q N,x = −q P,x = [(V 0 /v F ) 2 − q 2 y ] 1/2 in the N region and P region of the PNJ, which suppresses the rapidly oscillating phase q P · R 2 − q N · R 1 = 0 [see Eq. (6)] as long as R 1 and R 2 are mirror symmetric about the PNJ. The key ingredients of this effect are the circular Fermi contours and the opposite group velocities in the N region and P region, although the linear dispersion and gapless feature of graphene allows high transmission of electron waves (i.e., Klein tunneling) and hence quantitatively stronger focusing effect. Consequently, similar principle should lead to similar effect in other materials with a nonlinear dispersion and/or a finite gap (e.g., silicene PNJ [63]). For example, we have numerically verified that the local enhancement of the RKKY interaction remains effective even when the gap of the graphene PNJ increases to ∆ = 0.1 t.
The anomalous focusing is pronounced at low temperature and persists up to nitrogen temperature [2]. To observe the 1/R long-range RKKY interaction across the graphene PNJ, experiments should be carried out at a temperature higher than the Kondo temperature to avoid the screening of the local spins by the carriers [35]. The anomalous focusing in graphene PNJ has been demonstrated by two recent experi-ments [12,13]. Thus we expect our theoretical results to be experimentally accessible.

IV. SUMMARY
As opposed to previous works that explore the influence of electronic states and energy band structures of uniform 2D systems on the carrier-mediated RKKY interaction, we have proposed a very different way to manipulate the RKKY interaction: tailoring the carrier propagation and interference via a well-known electron optics phenomenon in graphene P-N junctions -the anomalous focusing effect. This gives rise to rich spatial interference patterns and locally enhanced, symmetry-protected non-oscillatory RKKY interaction, which may pave the way towards long-range spin-spin interaction for scalable graphene-based spintronics devices.
The key physics leading to the focusing of the RKKY interaction is the focusing of the carrier spin fluctuation emanating from a local spin. In this context, we notice a very relevant work by Guimaraes et al. [64], which shows that a gate-defined curved boundary in graphene can focus the spin current emanating from a precessing magnetic moment onto a specific point. We expect that this spin current lens could be utilized as an alternative way to focus the RKKY interaction.
( ← − ∂ x ′ acting on the left) and the continuity conditions Then g 1D (x, x ′ , E) is obtain by first calculating the general solutions in the region x x ′ and then matching them using the boundary conditions.
To present the results in a physically intuitive way, we introduce the following concepts. Given the energy E and momentum q y , there is one right-going eigenstate |+, q N with q N ≡ (q N,x , q y ) and one left-going eigenstate |+,q N with q N ≡ (−q N,x , q y ) in the N region, as well as one right-going eigenstate |−, q P with q P ≡ (q P,x , q y ) and one left-going eigenstate |−,q P withq P ≡ (−q P,x , q y ) in the P region, where |s, q ∝ e iq·r |u s (q) is the eigenstate of uniform graphene in the conduction band (s = +) or valence band (s = −). For small |q y |, we choose q N,x > 0 and q P,x < 0, so that |+, q N and |−, q P (|+,q N and |−,q P ) propagate from the left (right) to the right (left) without decay. For large |q y |, we choose Im q N,x > 0 and Im q P,x > 0, so that |+, q N and |−, q P (|+,q N and |−,q P ) decays to zero at x → +∞ (x → −∞).
In terms of these left-going and right-going eigenstates with given energy E and q y , the 1D propagator from x ′ in the N region to x in the P region is where w(q y ) = 2 cos θ N /(e −iθ N + e iθ P ) is the transmission coefficient and v(q) ≡ v F q/|q| is the group velocity, with θ N (θ P ) is the incident (transmission) angle defined via v F (q N,x + iq y ) = (E +V 0 )e iθ N and v F (q P,x +iq y ) = −(V 0 −E)e iθ P . Substituting Eq. (A2) into Eq. (A1) gives the 2D propagator g(r, r ′ , E)| r ′ ∈N,r∈P in Eq. (6) of the main text. For x ′ < x < 0, the 2D propagator g(r, r ′ , E) is the sum of the direct, forward propagation from r ′ to r, dq y 2π |u + (q N ) u + (q N )| iv x (q N ) e iq N ·(r−r ′ ) , and the contribution from the reflected wave dq y 2π r(q y ) |u + (q N ) u + (q N )| iv x (q N ) e iq N ·r e −iq N ·r ′ via three steps: the propagation from r ′ to the PNJ, the reflection by the PNJ, and the propagation of the reflected wave from the PNJ to r. For x < x ′ < 0, g(r, r ′ , E) is obtained from g(r, r ′ , E)| x ′ <x<0 by replacing the direct, forward propagation by the direct, backward propagation: dq y 2π |u + (q N ) u + (q N )| iv x (q N ) e iq N ·(r−r ′ ) .
The above transformation properties of the propagator leads to the corresponding properties of the RKKY range function J, e.g., J is invariant upon t → −t. Similarly, upon (V 0 , E F ) → (−V 0 , −E F ), we have Using ∞ −∞ Im G 2 (R 2 , R 1 , E)dE = 0, we see that J is invariant upon (V 0 , E F ) → (−V 0 , −E F ). Also, when the two spins locate at mirror symmetric points R and R m , their RKKY interaction J is invariant upon E F → −E F .