Nonequilibrium RKKY interaction in irradiated graphene

We demonstrate that the Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction in graphene can be strongly modiﬁed by a time-periodic driving ﬁeld even in the weak drive regime. This effect is due to the opening of a dynamical band gap at the Dirac points when graphene is exposed to circularly polarized light. Using Keldysh-Floquet Green’s functions, we develop a theoretical framework to calculate the time-averaged RKKY coupling under weak periodic drives, and we show that its magnitude in undoped graphene can be decreased controllably by increasing the driving strength, while mostly maintaining its ferromagnetic or antiferromagnetic character. In doped graphene, we ﬁnd RKKY oscillations with a period that is tunable by the driving ﬁeld. When a sufﬁciently strong drive is turned on that brings the Fermi level completely within the dynamically opened gap, the behavior of the RKKY coupling changes qualitatively from that of doped to undoped irradiated graphene.


I. INTRODUCTION
The Ruderman-Kittel-Kasuga-Yosida (RKKY) coupling is the indirect exchange interaction in a magnetically doped system, in which the coupling between the localized impurity spins is mediated via the conduction electrons of the host material [1][2][3].As a function of the impurity separation, it oscillates between positive and negative values indicating ferromagnetic and antiferromagnetic spin couplings with a period determined by the host's Fermi level [4].The envelop of these oscillations decays in a power law that depends on the dimensionality and the specific band structure of the host material [5][6][7][8][9][10][11][12].For systems with parabolic bands, the envelope of the RKKY oscillations decays as R −3 in three dimension (3D) and as R −2 in two dimension (2D) [1][2][3].
The 2D nature and Dirac energy dispersion of graphene have spawned the discovery of many unique and intriguing electronic properties [13,14].In particular, the RKKY interaction between magnetic impurities deposited on graphene exhibits peculiar properties that are quite different from conventional systems.For undoped graphene, the power law of the RKKY oscillations is R −3 in contrast to R −2 expected in 2D.In addition, the RKKY oscillations do not alternate between ferromagnetic and antiferromagnetic values but can instead display short-ranged fluctuations depending on the locations of the impurity atoms on the honeycomb lattice [5][6][7][8]15].This dramatic departure from the normal 2D RKKY behavior is not unexpected since the RKKY oscillation is fundamentally a Fermi surface phenomenon, and undoped graphene is characterized by a pair of Fermi points at the Dirac nodes.In doped graphene, some of the conventional 2D RKKY behavior is recovered due to the presence of a finite Fermi surface.Hence, on top of possible short-ranged fluctuations due to the lattice, its RKKY interaction oscillates between ferromagnetic and antiferromagnetic values, with the envelope decreasing as R −2 .This makes the RKKY coupling of doped * These two authors contributed equally.
graphene qualitatively similar to that of a 2D systems with a parabolic dispersion [8,10].
The long-range nature of the indirect exchange interaction provides an important mechanism in layered magnetic structures allowing for the controlled storage and transfer of spin-encoded information [16,17].The capability to further manipulate the indirect exchange interaction would add another dimension towards the possibility of engineering new material behaviors and functionalities.Over the past two decades, there has been much interest and progress in the studies of ultrafast optical control of spin dynamics and magnetic order [18].An emerging strategy to achieve this end in magnetically ordered materials is through optical control of direct exchange coupling [19,20].In magnetically doped or layered systems, RKKY interaction is an important type of exchange coupling to consider.Optical manipulation of RKKY interaction was first studied in semiconductor quantum dot systems and II-VI diluted magnetic semiconductors [21][22][23][24], where the exchange interaction can be mediated by virtual electron-hole pairs or excitons generated by a sub-band gap optical excitation.Manipulation of RKKY interaction is also possible through materials engineering of magnetic multilayers, where the interlayer interaction is predominantly due to RKKY coupling [25,26].An applied static electric field [27][28][29][30] could also provide a means to control the RKKY interaction in materials and structures with considerable spin-orbit coupling [11,12,31,32].
Recently, there has been a burgeoning interest in the study of periodically driven Floquet systems and the tantalizing possibility to engineer and control such non-equilibrium quantum states [33].Periodic driving through monochromatic irradiation of a solid renormalizes the band structure through non-perturbative lightmatter coupling into Floquet-Bloch bands.By tuning the intensity and frequency of the radiation one can dynamically tune these photon-dressed bands and hence control the properties of the irradiated system.Besides the control of direct exchange coupling already mentioned, Floquet dynamics also strongly influence transport properties [34][35][36], optical properties [37][38][39], topological phases [40][41][42] and interlayer tunneling [43,44].

arXiv:2004.11337v1 [cond-mat.mes-hall] 23 Apr 2020
The dynamic tunability and exquisite level of control offered by optical driving can provide a significant advantage over other proposed mechanisms for controlling the RKKY interaction.In this work, we study the nonequilibrium RKKY interaction in magnetically doped graphene under illumination of monochromatic circularly polarized light.Due to their linear dispersion and valley degrees of freedom, chiral Dirac fermions in graphene couple distinctively to electromagnetic field compared to Schrödinger fermions with a parabolic dispersion.Indirect exchange interaction mediated by irradiated Dirac fermions are therefore expected to display unconventional properties not found in equilibrium.Our theory is developed using the Keldysh-Floquet formalism in order to capture the non-perturbative Floquet dynamics of the RKKY interaction.We derive a formula for the timeaveraged RKKY coupling for the weak drive, off-resonant regime, and perform analytical and numerical studies as a function of impurity locations on the lattice sites and Fermi levels, for various driving strengths.Our results reveal that for undoped graphene a weak periodic driving primarily decreases the RKKY exchange interaction in proportion to the driving strength, whereas for doped graphene increasing the driving strength increases the period of RKKY oscillation allowing for a optical tuning of indirect exchange coupling.We are also able to identify a parameter regime where the exchange coupling can be tuned between ferromagnetic and antiferromagnetic.
Our paper is organized as follows.In Sec.II we introduce the theoretical model of graphene with magnetic adatoms under monochromatic illumination.We derive the Floquet Hamiltonian for irradiated graphene and then use the Keldysh-Floquet formalism to derive a general formula for the time-averaged non-equilibrium RKKY interaction.We then specialize to the weak drive limit and obtain an approximate Floquet Hamiltonian under the photon-number representation and the corresponding expression of the irradiated RKKY coupling.We next proceed to obtain the explicit analytical formulas of real-space non-equilibrium Green's functions in Sec.III.In Sec.IV we present and analyze our results for the time-averaged RKKY coupling for the cases of undoped graphene (Sec.IV A) and doped graphene (Sec.IV B) under different impurity configurations.These numerical results are elucidated with analytical studies where we obtain approximate analytic expressions for the irradiated RKKY coupling and exhibit their dominant dependence on the impurity separation.Before concluding in Sec.VI, we present a few relevant remarks on our results and an outlook for future direction in Sec.V.

II. FORMULATION
The low-energy electrons in single layer graphene are described by two inequivalent high-symmetry points (K and K ) in the Brillouin zone known as the Dirac points.
Electrons in the vicinity of these points disperse linearly and can be described by a Dirac Hamiltonian [13,14] where v F = 10 6 ms −1 is the band velocity and τ = 1 (τ = −1) corresponds to the K (K ) Dirac points.v F is related to the graphene tight-binding model parameters by v F = 3γa/2 = 6.6 eV Å with γ = 3 eV being the nearest-neighbor hopping amplitude and a = 1.4 Å the carbon-carbon distance.In the presence of magnetic impurity adatoms, electrons in graphene couple with the impurity spins through a short-range interaction.In the spirit of standard RKKY theory, such an interaction is modeled as a local potential at the impurity site and the single-particle Hamiltonian takes the form where S is the electron spin operator, S i (i = 1, 2) are the spin angular momenta of the magnetic impurities assumed to be located at the origin and at R, and λ is the coupling between the itinerant and localized spins.
The honeycomb graphene lattice can be considered as two superposing sublattices A and B. The location of the two impurity spins can be either at the sites of the same sublattice A or different sublattices A and B.
We consider normal incidence on the graphene plane with a circularly polarized (CP) light with frequency Ω and field amplitude E 0 .The electric field couples to the Hamiltonian in Eq. ( 1) via the minimal coupling scheme and the resulting timedependent Hamiltonian of the irradiated graphene system becomes where A = v F eE 0 /( √ 2Ω).Due to the time-periodicity of the Hamiltonian Eq. ( 3), the Floquet-Bloch theorem is satisfied granting a solution of the following form to the time-dependent Dirac equation where l ∈ Z labels the Floquet modes and l,k is the quasienergy modulo Ω.The time periodic nature of the Floquet states, |u l,k (t + T ) = |u l,k (t) , where T = 2π/Ω is the driving period, further allows a Fourier series representation of these states, Substituting |ψ l,k (t) in Eq. ( 3) maps the time-dependent Dirac equation into the Floquet Hamiltonian, H F = H 0 (t) − i ∂ t in the Fourier domain [45], such that where H F,mn = H mn − n Ωδ m,n is the Floquet Hamiltonian and is known as the Floquet matrix [46,47].For our system, the Floquet matrix takes a block-tridiagonal form and can be explicitly written as where τ = + (τ = −) corresponds to the K (K ) point.
Here we should note that the adoption of the continuum description of graphene is essential to obtaining the block-tridiagonal form of the Floquet Hamiltonian of irradiated graphene.The adoption of this description, however, sets an upper bound to the light-matter coupling strength A = A/( Ω), where A is given in Eq. ( 3).This upper bound can be seen from the following consideration, starting first with the tight-binding model of graphene.The Floquet Hamiltonian can be derived from this tight-binding description, which consists of a Floquet matrix with every elements being generally non-zero.If we require this to recover the continuum description of the Floquet Hamiltonian in Eq. ( 8), all elements of the block-pentadiagonal and higer-ordered bands of this matrix must be set to zero (here the terminology "bands" refer to the bands of the matrix).The matrix elements in those bands are proportional to , where J n is the Bessel function of order n.The condition J n 2 [ √ 2aA/( v F )] ≈ 0 gives aA/ v F 2 and thus A 3γ/ Ω.Here a and γ are defined under Eq.(1).For our numerical calculations in Sec.IV, we choose Ω = 6.6γ, which puts the driving field in off resonance, and A 0.06 so that both the above condition and the weak drive condition are satisfied with A 3γ/ Ω < 1.The above consideration was first discussed in the context of irradiated two-dimensional electron gases in Ref. [48] and can be seen to apply quite generally for other systems whenever one tries to obtain a Floquet Hamiltonian within a continuum low-energy description.
Turning on the driving field at t = 0 starts pumping energy into the system.Crucially, the electron subsystem in a solid is always an open system, with energy relaxation pathways through various inelastic scattering processes (electron-electron and electron-phonon scattering) as well as coupling to extrinsic degrees of freedom in the ambient environment.In the presence of both energy injection and relaxation, naturally heat does not build up infinitely over time.We adopt a specific model to account for energy relaxation by coupling graphene to an external heat reservoir, which will be specified in more details in a later section, Sec.III.After initial transients has subsided, the system dynamics settles into a non-equilibrium steady state (NESS).In NESS, time periodicity is restored and the system's Green's functions become periodic in time, G R,A,< k (r, t + T ; r , t + T ) = G R,A,< k (r, t; r , t ), where G R,A,< denotes the retarded, advanced, and lesser Green's function.Periodicity allows these Green's functions to be expressed in the Floquet representation [47], [G(r, r , ω)] mn = (9) dt rel e i(ω+mΩ)t−i(ω+nΩ)t G(r, t; r , t ), where t av = (t + t )/2 and t rel = t − t are the average time and relative time.The frequency variable ω in the above is defined in the reduced zone ω ∈ (−Ω/2, Ω/2].We will use the notation ω for reduced-zone frequency and ω ∈ (−∞, ∞) for extended-zone frequency, which will be relevant for the Green's functions we will discuss further below.
The RKKY interaction coupling is mediated by the spin polarization of the graphene electrons induced by the impurity spins.We use perturbation theory to obtain the interaction energy up to leading order in S 1 • S 2 .The exchange energy between impurity spins S 1 and S 2 is given by the interaction energy between one impurity spin and the induced electron spin density due to the other impurity spin, (10) where S(r, t) = −iTr[SδG < (r, t; r, t)] is the expectation value of the induced spin density, with δG < (r, t; r, t) denoting the change in G < due to the perturbation of S 1 .Dyson's equation up to first order in the local potential V (r) = λδ(r)S • S 1 gives δS(r, t) = −iTr{SG R (r, t; r , t )V (r )G < (r , t ; r, t) +SG < (r, t; r , t )V (r )G A (r , t ; r, t)}, (11) where G R,A,< are the non-interacting Green's functions of the irradiated graphene sheet, the trace is over r , t and the spin degrees of freedom.Let the impurity spin S 1 be sitting on the sublattice site β and S 2 be sitting on α, where α, β = {A, B}.Then, substituting Eq. (11) in Eq. ( 10), expressing the time-dependent Green's functions in the Floquet representation and averaging over time, we obtain the Floquet representation of the timeaveraged exchange energy for impurity spins located at sublattice sites α and β: + G < αβ (R, ω)G A βα (−R, ω)}, where µ, ν are the x, y, z-projections of the impurity spin, a, b label the electron spin degrees of freedom, ω is the reduced-zone frequency, and the trace here is over the Floquet index of the Floquet Green's functions.In the above, it should be emphasized that all summations are written out explicitly and α and β are not summed over.Since the graphene Hamiltonian Eq. ( 1) is spinindependent, the trace over spins is rendered only on the electron spin operators S in the above equation.We have assumed that the concentration of the magnetic adatoms (e.g. from transition metal elements) is dilute enough not to affect the band structure of graphene significantly.In the case of sufficient adatom concentration, graphene electrons could acquire spin-orbit coupling from the adatoms [49][50][51][52] and the calculation here would need to be generalized to account for spin-orbit coupling.
The explicit summation over spins can be carried out on the electron spin operator S µ = ( /2)σ µ , by using the trace relation Tr{σ µ σ ν } = 2δ µν for Pauli matrices σ µ , yielding so that the exchange energy becomes , where the exchange coupling strength can be written as (see Appendix D and Eq.(A9) for details),

III. QUASIENERGY SPECTRUM AND REAL-SPACE FLOQUET GREEN'S FUNCTIONS
To compute the RKKY coupling from Eq. ( 14), in this section we derive the real-space Floquet Green's functions of irradiated graphene.The Floquet Green's functions can be obtained directly from the Floquet Hamiltonian.For an arbitrary driving strength, this generally needs to be done numerically and the calculation of Eq. ( 14) would require extensive computational effort, and increasingly so towards stronger driving.In this work, we focus only on the weak drive regime A 1, in which the theory becomes analytically tractable.To this end, we first project the Floquet Hamiltonian into a new set of basis, which will be called photon-number representation, that makes it more physically transparent.Its main idea can be motivated from the observation that the graphene Floquet Hamiltonian Eq. ( 8) under CP light breaks into a set of decoupled two-level systems at k = 0.For instance, at the K-point, those two levels comprise the A sublattice component of the n − 1 Floquet mode and the B sublattice component of the n Floquet mode, they are coupled via A = v F eE 0 /( √ 2Ω) but are decoupled from the all other Floquet modes.Hence, we can project the entire Floquet Hamiltonian for all k values into the set of basis that diagonalize its block diagonal representation at k = 0.This method was first introduced in Ref. [53] and we can express it in the following form consisting of three consecutive unitary transformations: where τ = ±1 and the operators with τ = + (τ = −) act on the Hamiltonian at the K (K ) valley.The definitions of the transformations U 1,τ , U 2,τ , U 3 are relegated to Appendix B. The resulting Hamiltonian HF is in the basis , where the ↑, ↓ symbols label a new "isospin" degrees of freedom: for the K point and for the K point, where One notices that the explicit matrix forms of HF for the two valleys happen to be the same when expressed in the bases Eqs. ( 16)- (17).It is also useful to note that the basis of the K point can be transformed to the basis of the K point via where Σ y = iσ y ⊗ I ∞ , I ∞ is the identity matrix in the Floquet space, and Φ K ↑,↓ and Φ K ↑,↓ correspond to Φ ↑,↓ in Eqs. ( 16) and (17), respectively.
Up to this point, there is no approximation and the transformed Hamiltonian HF captures the same physics as the original Floquet Hamiltonian in Eq. ( 6), but it provides a more convenient picture of the phenomena at play.The expression of HF contains the following parameters Each of these parameters captures a successive photon process mimicking a tight-binding Hamiltonian in the Floquet-ladder space; i.e.F 0 captures the zero-photon process that opens a gap at k = 0 resulting from timereversal symmetry breaking ("on-site").F 1 captures the one-photon resonance resulting from the hybridization of the first nearest-neighbor (|n − m| = 1) Floquet-ladder bands at k = Ω/(2v F ), while F 2 captures the two-photon resonance resulting from the hybridization of the second nearest-neighbor (|n − m| = 2) Floquet-ladder bands at k = Ω/v F .Similarly, there are other coefficients F N >2 describing higher-order photon processes, but wiill not be needed for our present theory.The values of F N =0,1,2 in the Floquet-ladder tightbinding Hamiltonian HF depends on the driving strength A = A/( Ω).For consistency with the regime of validity of the continuum Floquet Hamiltonian that sets an upper limit on the value of A (see Eq. ( 8) and the relevant discussion in Sec.II), in the following we will focus on the weak drive regime A 1 which renders F 1 and F 2 negligible compared to F 0 because they are smaller by a factor of A and A 2 .This approximation is known as the F 0 approximation [53].Hence, HF becomes dependent only on F 0 and takes a block-diagonal form HF = ∞ n=−∞ hn , where stands for the matrix direct sum over the Floquet space, with the following 2 × 2 block Hamiltonian hn : where ∆ = √ 4A 2 + 2 Ω 2 − Ω and θ k is the angle between k and the K − K direction.Within the F 0approximation the Hamiltonian in Eq. ( 20) captures the photon-induced band gap ∆ at the K and K points and the photon-induced renormalization of the graphene Fermi velocity by F 0 .The quasienergies for K or K are simply gapped Dirac dispersions shifted by integer multiples of Ω, as shown for one of the valleys in Fig. 1.
Here we can make a connection to the leading-order approximate Floquet Hamiltonian H FM obtained under the Floquet-Magnus expansion [54] that is commonly used as starting point for calculation of Floquet dynamics in the high-frequency off-resonant regime.For irradiated Dirac system such as graphene, it gives a static gapped Dirac Hamiltonian with a gap ∆ FM = A 2 /( Ω) describing only the n = 0 Floquet mode [35].On the other hand, the F 0 Hamiltonian captures not only the induced gap but also the Fermi velocity renormalization for all Floquet modes n ∈ Z.For the case of graphene, one can obtain the leading-order H FM from the F 0 approximation by expanding Eq. ( 20) in ∆ and F 0 up to second order in A and then setting n = 0.
Having obtained the Floquet Hamiltonian, we can now proceed to obtain the real-space Green's functions in the photon-number representation under the F 0 approximation.The simple form of the F 0 Hamiltonian allows us to obtain analytically tractable Green's functions.We start with the momentum-space Green's functions and obtain the corresponding real-space counterparts through a Fourier transformation.The retarded Green's function in the momentum space is GR = ∞ n=−∞ gR n , where the 2 × 2 block Green's function corresponding to the n th Floquet mode is given by gR n (k, ω) = [(ω + iη)I σ − hn ] −1 .In the following, we write the Green's function in terms of the physical, extended-zone frequency ω = ω − n Ω with gR (k, ω) ≡ gR n=0 (k, ω): where The color intensity in Fig. 1 shows the spectral function of the Floquet states calculated from Eq. ( 22), where the trace is over is the isospin degree of freedom (↑ , ↓).The real-space Green's functions for each valley are obtained by their corresponding Fourier transformation of the momentum-space Green's functions Substituting Eq. ( 22) into the above, we obtain the retarded Green's function for each Floquet mode in the photon-number representation (see Appendix C for details of the derivation), [ where + (−) on the right-hand side of Eqs. ( 25)-( 26) corresponds to ↑↑ (↓↓) and ↑↓ (↓↑) configurations of the isospins in the left-hand side of these equations [see Eqs. ( 16) and ( 17) for isospin definition], is the angle between R and the x-axis, and where H (1) 0 and H (1) 1 are the zeroth-order and first-order Hankel functions of the first kind, respectively.
We assume thermalization of the irradiated system is achieved in NESS through a coupling to an external fermion bath (e.g., a metallic lead).The Floquet-Keldysh formalism enables one to determine the Floquet lesser Green's function in the momentum space as where Σ < = −Σ R F + F Σ A is the lesser self-energy, and In the wide-band approximation for the fermion bath, the retarded self-energy takes the momentum-independent form Σ R = −iηI σ ⊗ I ∞ , and the lesser self-energy becomes Σ < = 2iη Transforming to the photon-number representation, we find the lesser Green's function The real-space lesser Green's function is then obtained from the Fourier transformation of Eq. (31) where With all the required Green's functions obtained, we return to the expression of the RKKY coupling in Eq. ( 14).The main contribution to the RKKY coupling in graphene comes from the K and K points.Therefore the total real-space Green's function consists of contributions from both valleys, where G + (G − ) denotes the Green's function associated with the K (K ) Dirac points.The above Green's function G(R, ω), which enters Eq. ( 14), is written in the original Floquet representation.A crucial observation is that the unitary transformation U T,± in Eq. ( 15) for the Hamiltonian is k-independent; this allows the same transformation to carry over to the transformation for real-space Green's function also.Through this transformation we can formally express all the Green's functions in Eq. ( 14) from the original Floquet representation into the photon-number representation.We choose the basis for the K point, Eq. ( 17), as the common basis.Then, applying the transformation U T,− [Eq.(15)] and defining where we have used Eq. ( 18) and have written G(k, ω) ≡ G+ (k, ω) = G− (k, ω), since G± (k, ω) have the same explicit matrix form in the photon-number representation.With Eq. ( 35) and the F 0 approximation, Eq. ( 14) takes on the simplified form (see Appendix D for details), Im Tr GR αβ (R, ω) G< βα (−R, ω) .
Under the F 0 approximation, the Green's function in Eq. ( 36) are block diagonal in the Floquet index n consisting of the 2 × 2 Green's functions gn (R, ω).This allows us to further express Eq.( 36) in terms of the extended zone frequency ω ∈ (−∞, ∞), (38) with g(R, ω) given by Eq. ( 24).In the above, we have used a common numerical representation {+1, −1} for the pseudospins {A, B} and the K isospins {↑, ↓}.In the absence of irradiation, Eq. ( 37) can be reduced to its equilibrium limit [5,7] as shown in Appendix D.

IV. NUMERICAL AND ANALYTICAL RESULTS
Eq. (37) allows for efficient numerical and analytical evaluation of the non-equilibrium RKKY coupling.In this section we present numerical results for J(R) and perform additional approximations to obtain analytic results.We consider both undoped and doped graphene for different values of A, for impurities located at the same sublattice sites or at different sublattice sites, with the separation R along the zigzag or armchair direction.In all cases we also include results for the undriven case with A = 0, which are consistent with the equilibrium results [5][6][7][8] in the literature.
A. Undoped Graphene, EF = 0 The panels of Fig. 3 show our numerically evaluated RKKY coupling J for different cases of impurity locations and separation directions for Fermi level located at the Dirac point.The choice of the range of driving strength A is informed by its upper limit A 3γ/ Ω < 1 (see Sec. II) for the particular frequency value Ω = 6.6γ we have chosen.For this large frequency value the driving field is off-resonant and the system is in the weak drive regime.
When the magnetic impurities are located on the same sublattice (AA or BB) and connected through the zigzag direction [Fig.3(a)], J displays short-scale oscillations as a function of the impurity separation R for all the driving strengths A studied.This is a general feature in systems with a multi-valley band structure [4,55] and arise from intervalley scattering.Graphene's band structure contains two valleys K and K that originate from the honeycomb lattice structure.As such, these short-scale oscillations persist with and without irradiation and have an intrinsic period that is independent of Fermi level.They have a period 2π/|K − K | = 3 √ 3a/2, where |K − K | is the distance between two adjacent valleys in the Brillouin zone and a is given in Eq. (1).It is convenient to take the k x -axis along an adjacent pair of K and K with the origin chose at the M -point between them, such that K − K = 2K.One may notice that these short-range oscillations can only manifest for certain choices of impurity separation in graphene, namely, when the phase difference (K − K ) • ∆R = 2π, where ∆R is the increment that gives the next point of impurity separation on the lattice, and is given by the distance between two adjacent impurity sites along the zigzag or armchair direction.Along one of these directions, it is the vector joining adjacent A sublattice sites for impurity spins on the same sublattice, and the vector joining adjacent A and B sublattice sites for impurity spins on different sublattices.As a result, these sharp oscillations appear in the case of the AA (BB) impurities separated along the zigzag direction, which has (K − K ) • ∆R = 2π/3, in addition to the overall decrease with increasing R.For weak driving strengths A, we observe that irradiation suppresses the RKKY coupling for all values of R while maintaining its sign, hence remaining ferromagnetic in character J(R) < 0 like in equilibrium.Interestingly however, we find that for large enough driving J(R) can flip sign and display antiferromagnetic character within a rather large range of R, as seen for A = 0.06 in Fig. 3(a).
To gain further insights into this behavior, we proceed to derive approximate analytic results.Our results in the following are valid up to O(∆ 2 ) in the weak drive A 1 regime when the gap ∆/ Ω A 2 1. Taking η → 0 in Eq. ( 37) with the expressions of the Green's functions [g R (R, ω)] ↑↑,↓↓ and [g < (R, ω)] ↑↑,↓↓ from Eqs. ( 25) and (32), and defining z = −α |(2ω/∆) 2 − 1|, the RKKY coupling for α = β = A can be written as In Eqs. ( 39)-( 40), J 0 , Y 0 are the zeroth-order Bessel and Neumann functions, I 0 , K 0 are the zeroth-order modified Bessel functions of the first and second kind.We have also made use of K − K = 2K and defined In Eq. ( 39), since z < α 1, we can use the small argument expansions of the K 0 and I 0 and expand the integrands up to second order in z.After integration we find both integrals are ∼ (∆/2) 3 and are negligible, hence J(R) ≈ J0 (R) .In Eq. ( 40) for J0 (R), since ∆/ Ω is small, α 1 will be satisfied for a large range of R. We can therefore extend the lower integration limit from α to 0, which incur an error ∼ (∆/2) 3 that is negligible.Now J(R) can be written as which is analytically integrable leading to When ∆ = 0 and F 0 = 1, Eq. ( 43) recovers the the equilibrium RKKY coupling, with J(R) decaying as 1/R 3 modulated by the short-range oscillation factor [1 + cos (2K • R)] that arises from intervalley scattering.Notice that these oscillations do not cause any sign change in J as a function of R and the equilibrium RKKY coupling retains its ferromagnetic character for all R [5][6][7].Irradiation modifies this behavior in two interesting ways.In Eq. ( 43), the first term, which is inherited from the equilibrium contribution, is renormalized by F 2 0 which decreases with increasing driving strength A. The second term is a new contribution due to irradiation.It decays more slowly as 1/R and increases as ∆ 2 with the driving strength A. It displays the same oscillation period since this aspect is uniquely determined by the graphene lattice and the impurity orientation.The interplay between these two terms in Eq. ( 43) explains the numerically observed RKKY coupling in Fig. 3(a) as follows.We notice that the second term has an opposite sign and partly cancels the first, ferromagnetic-like term.Increasing A thus suppresses the RKKY coupling, as we have found numerically.If A is large enough, the second term becomes more dominant and the RKKY coupling switches to an antiferromagnetic character.The persistence of this behavior for large R values as observed in Fig. 3(a) is due to the much slower 1/R dependence.
We now consider the case for two impurities located at opposite sublattice sites (AB or BA) connected through the zigzag direction.Fig. 4(b) shows the numerically calculated J(R) for the AB zigzag case, which displays short-range oscillations similar to the AA zigzag case since (K − K ) • ∆R = 2π also for this configuration.As opposed to the AA zigzag case, the exchange coupling for the AB case is antiferromagnetic in character in equilibrium, and remains so under weak driving A 1. Following similar steps as in the AA zigzag case, we obtain the following analytical result for the AB zigzag RKKY coupling The above recovers the equilibrium case [5][6][7] in the absence of irradiation (∆ = 0 and F 0 = 1), showing an antiferromagnetic coupling.Similar to the AA case, Eq. ( 44) consists of two terms, with the second term ∼ 1/R arising purely from irradiation.It is opposite in sign and increases as ∆ 2 with the driving strength, causing a suppression of J(R) with increasing A. For weak driving in the range of A considered, this second term does not become large enough to dominate the first term for all values of R. As a result, J(R) does not change sign as observed in Fig. 4(b), and the irradiated case resembles qualitatively the equilibrium case.The effect of the relative orientation between the two impurities on the graphene lattice is brought out in the case when the separation is along the armchair direction.Figs.4(c)-(d) show the cases when the impurities are located at the same sublattice sites AA and different sublattice sites AB.The separation ∆R between the two impurities satisfy (K − K ) • ∆R = 2π, and as a result J(R) does not display any short-scale oscillations.Its smooth decay profile resembles the RKKY coupling in conventional insulators [56][57][58][59].The RKKY coupling remains ferromagnetic for the AA case and antiferromagnetic for the AB case in the weak drive regime, with mag-nitudes that decrease with increasing A. We also obtain the following analytic results for the AA case and the AB case Except for the absence of the short-range oscillation factors, the above expressions for the armchair case resemble those in the zigzag cases, with similar dependence on R, F 0 and ∆ in each of the two constituent terms.To summarize the above analyses of different impurity configurations on undoped graphene, we have found that the main effect of irradiation under weak driving A 1 is to decrease the magnitude of the RKKY coupling from its equilibrium value.Interestingly, for the particular AA zigzag case, we find that a moderately strong drive (still within A 1) can result in a sign change of the RKKY interaction.In the next section, we will present our results for the doped case with a finite Fermi energy.

B. Doped Graphene, EF = 0
In our calculations we take the Fermi level to be fixed by coupling to the fermion bath (e.g. through attachment to a metallic lead).In addition to the short-scale oscillations resulting from intervalley scattering for the zigzag case, the RKKY interaction displays an additional, longer-range oscillation arising from the intravalley scattering processes near the Fermi surface of each valley.This long-range RKKY oscillation is familiar to the usual case as it fluctuates between both signs as a function of R. Figs.4(a)-(b) show the AA and AB zigzag cases respectively with a Fermi level above the photon-induced band gap, E F > ∆/2.The longerrange oscillation is due to the q = 2k F intravalley scattering and has a period given by π/k F , where k F is the Fermi wave vector.While the short-range oscillation period remains unchanged under irradiation, we find that the longer-range oscillation period increases with A. This is due to the increased photon-induced gap ∆ with A, resulting in a smaller Fermi surface and a smaller 4(c)-(d) show the RKKY coupling respectively for impurities at AA and AB along the armchair direction with a Fermi energy E F > ∆/2.The oscillation is marked by the absence of the shorter-range oscillations as in the undoped case [Figs.3(c)-(d)] and a prevailing period π/k F that increases with A.
In the following we derive the approximate analytical results for the RKKY coupling up to O(∆ 2 ) for different impurity configurations when 2E F > ∆.The total exchange coupling can be written as a sum of the intrinsic contribution already obtained for undoped graphene and a second contribution due to a finite Fermi energy, We start with the AA zigzag configuration.The intrinsic contribution J int (R) is given by Eq. ( 43).To evaluate the extrinsic contribution, for convenience we define y = α |(2ω/∆) 2 − 1| and write J ext (R) = Ĵ1 (R) + J1 (R), where Ĵ1 (R) and J1 (R) take the following form and Similar arguments for Eqs. ( 39)-( 40) applies here.Expansions of the integrands in Eq. ( 47) up to second order in y < α 1 show that both integrals are ∼ (∆/2) 3 and are negligible.Eq. ( 48) for J(R) can be evaluated by extending the lower integration limit to 0 with an negligible error ∼ (∆/2) 3 .Then J ext (R) ≈ J1 (R) can be written as The above integrals can be performed analytically, yielding J(R) = J int (R) + J ext (R) as where and G m,n p,q is the Meijer G-function [60].A more familiar form of the RKKY interaction can be obtained by considering the asymptotic regime k F R 1 when the impurity separation is long compared to the Fermi wavelength.From the asymptotic forms of the Meijer G-functions we find from which we obtain At equilibrium, the dominant contribution in J(R) for finite Fermi energy goes as 1/R 2 similar to 2D metals [5][6][7], in contrast to the 1/R 3 dependence in undoped graphene.Under irradiation, the doped case has a notable distinction from the undoped case, the dominant R-dependence in the second purely non-equilibrium term of Eq. ( 54) decays more rapidly going as 1/R 2 , as compared to 1/R in Eq. ( 43).As a result, the effect of this term in suppressing the exchange coupling becomes less important in the doped case.From Eq. ( 54) we can clearly distinguish the lattice, Fermi surface, and irradiation effects at play in the system.The lattice effect is reflected by the period of short-range oscillations π/|K| and is not influenced by irradiation.On the other hand, the Fermi surface effect is manifested via the longer period π/k F , which is dependent on ∆ and increases with the driving strength A.
If we keep increasing the driving strength for a fixed R we can leave the asymptotic regime and enter into a regime with k F R becomes smaller than 1, since k F decreases with A. This can be attained while keeping A 1. Thus we can consider the limit k F R 1 and use the small argument expansions of the Meijer G-functions to obtain, where γ is the Euler-Mascheroni constant.Using the above, we find that Eq. ( 50) reduces to the same form as Eq. ( 43) for irradiated undoped graphene.Physically, when the electrochemical potential is maintained constant in the system, increasing the driving strength A allows to "drain out" the Fermi sea via the increase in gap size ∆.To shed light further on this point, we draw a comparison between the two cases with the Fermi level far above the band gap [Fig.4(a)] and just above the band gap [Fig.5(a)].Fig. 4(a) shows, for a fixed E F , the nonequilibrium and the equilibrium cases.Their difference is relatively small since the Fermi momentum k F under irradiation is still close to the equilibrium value when E F is large compared to ∆/2.The difference between the irradiated and equilibrium cases is most dramatically illustrated when E F is fixed at a value just above the photon-induced gap under irradiation.Fig. 5(a) shows J(R) for such a value of E F in and out of equilibrium.The significant difference is due to the qualitatively distinct scenarios realized by increasing A. When undriven, the RKKY behavior is that of doped graphene in equilibrium, which displays long-range RKKY oscillations with a period π/k F .Under a driving field that induces a value of ∆ 2E F , the Fermi surface is drastically reduced and the RKKY coupling approaches that of irradiated undoped graphene [Fig.3(a)].This drastic change in the behavior of the RKKY exchange interaction under irradiation reveals a strong qualitative difference between the equilibrium and the irradiated cases; furthermore, it implies that by tuning the laser amplitude and frequency one can switch between the two qualitatively distinct regimes of RKKY interactions for doped and undoped graphene under irradiation.
One can follow similar calculation steps as outlined above for the AA zigzag case to derive the analytical results for the AB zigzag, and the AA and AB armchair configurations.In the asymptotic regime find, for the AB zigzag case, whereas for the AA armchair case, and for the AB armchair case, For all these cases, we also find from Eq. ( 50) that the k F R 1 behavior of J(R) reduces to their corresponding results in the undoped regime under irradiation [ABzigzag, Eq. ( 44); AA-armchair, Eq. ( 45); AB-armchair, Eq. ( 46)].The underlying reason is already explained above in the case for AA zigzag case.
The effect of the Fermi energy on the non-equilibrium RKKY interaction can be further elucidated by plotting J(R) as a function of E F at a fixed R, as shown in Fig. 6(a).It becomes clear that the RKKY coupling under irradiation approaches the equilibrium case when the Fermi energy is far from the conduction band edge E F ∆/2. Departure from the equilibrium case occurs most considerably in the region of low Fermi energies, which is enlarged as shown in Fig. 6(b).In particular, one notices that J(R) in the irradiated case remains constant throughout E F ∈ [0, ∆/2], because the Fermi level stays within the photon-induced gap.Thus at a fixed E F > 0, as one increases the driving strength A from zero, the gap ∆ increases, and the RKKY coupling changes from having a predominantly intravalley scattering character to having a purely intervalley scattering character.For large enough A when the dynamical gap encloses the Fermi level, the RKKY coupling remains constant behaving as irradiated undoped graphene.
V. DISCUSSION AND OUTLOOK Our theory applies for dilute concentration of magnetic impurities located on the graphene sublattice sites.One interesting direction for future consideration is the inclusion of spin-orbit coupling effects.For sufficient concentration of magnetic adatoms (e.g.transition metal elements), the band structure of graphene becomes strongly modified due to the induced spin-orbit coupling and Zeeman splitting [49][50][51][52].Our formalism can be further generalized to account for the induced spin-orbit coupling in the graphene electrons.
In equilibrium, it has been known that the RKKY interaction in undoped graphene is antiferromagnetic for impurity spins located at the same sublattice sites and ferromagnetic at different lattice sites.This statement is general for bipartite lattices and is satisfied due to the particle-hole symmetry in undoped graphene [5].In the presence of a driving field, one surprising finding from our study is the sign change of J(R) observed for the AA zigzag case when A = 0.06 [Fig.3(a)], resulting in a ferromagnetic coupling over a wide range of R.This is also depicted in Fig. 5(a), in which the A = 0.06 case exhibits essentially the same behavior as undoped graphene due to the small Fermi surface resulting from a large induced dynamical gap.Within our analytic approach, we have shown that this is due to competition between the equilibrium-like antiferromagnetic contribution and a purely light-induced ferromagnetic contribution [Eq.(43)].When A is large enough, the latter term wins out.An open question for future investigation would be to fully examine this sign change under larger values of driving strengths beyond the weak drive regime.
In doped system, we have demonstrated that the oscillation period of the RKKY coupling can be tuned by the driving field.This tunability is due to the dynamical band gap opened at the Dirac point under circularly polarized light and is more enhanced in gapless systems like graphene.One can compare this to the RKKY interaction in intrinsically gapped system such as doped III-V semiconductors.Because of its large intrinsic gap, the band gap increase due to renormalization by optical Stark effect under a subgap frequency irradiation will be relatively small, and the tunability of its RKKY interaction will be less dramatic.The qualitatively distinct behaviors of the RKKY coupling for doped and undoped irradiated graphene, marked by the presence and absence of long-range oscillations respectively, pose a tantalizing prospect for optical switching of ferromagnetic/antiferromagnetic coupling.Starting from slightly doped graphene at equilibrium, turning on the laser can switch the coupling from antiferromagnetic to ferromagnetic for the AA zigzag case [Fig.5(a .This observation suggests further investigation into Floquet control of the magnetic ordering in magnetically doped graphene would be promising.

VI. CONCLUSION
This work has presented a formalism for calculating the RKKY interaction in graphene irradiated by a circularly polarized light.We have obtained a generic expression of the time-averaged RKKY coupling in terms of Keldysh-Floquet Green's functions.By transforming the original Floquet Hamiltonian into a new basis, we arrived at a new Hamiltonian that carries a transparent physical meaning as a tight-binding-like Hamiltonian in the Floquet space, allowing for a systematic order-byorder inclusion of different photon processes.In this work we focus on the off-resonant, weak drive regime where only virtual photon absorptions and emissions happen, which allow us to truncate terms corresponding to first and higher order photon processes.In this framework we have obtained the numerical and approximate analytical results for the time-averaged RKKY coupling under a combinations of different impurity locations and separation directions.When the Fermi level is at the Dirac point (undoped graphene), we find that the magnitude of the RKKY coupling is decreased by increasing irradiation strength while maintaining its ferromagnetic or antiferromagnetic character.When the driving field becomes strong enough, our results suggest that the ferromagnetic coupling in the case of impurities at the same sublattice sites and separated along the zigzag direction can be switched to antiferromagnetic coupling for a wide range of impurity separation.When the Fermi level is above the Dirac point (doped graphene), the RKKY coupling is characterized by long-range oscillations with a period that can be tuned by the driving field.For large values of driving strength, the RKKY behavior of doped graphene turns into that of undoped graphene through the optically induced metal-insulator transition.and finally we take the time average, and then sum over m to arrive at the exchange in the Floquet representation in Eq. ( 12).
From the definitions of the real-space Green's functions [61], one can show that . These relations can also be seen explicitly generically from the momentum-space Green's functions, e.g., With these relations we can write the integrand in Eq. (12) as which leads to Eq. ( 14).

Appendix B: Unitary transformations for the photon-number representation
In this appendix we describe the set of unitary transformations used to obtain the Hamiltonian in the photon number representation HF in Eq. ( 15) [53].The first transformation rotates the original Floquet Hamiltonian H F at k = 0 so that it becomes block diagonal.At the K point this unitary transformation takes the form while for the K point this transformation takes the form U 1,τ =−1 = I ∞ , where I ∞ is identity operator in the Floquet space.
The second transformation is constructed as the transformation that diagnonalizes the transformed Floquet Hamiltonian U † 1,τ H F U 1,τ at k = 0.If we denote the spinor eigenvectors corresponding to the eigenvalues n,± = ±∆/2 + n Ω as φ τ n,± , then U 2,τ = [. . ., φ τ n−1,+ , φ τ n−1,− , φ τ n,+ , φ τ n,− , φ τ n+1,+ , φ τ n+1,− , . ..],where ∆ is given in Eq. ( 20) and τ = + (τ = −) corresponds to the K (K ) point.This transformation takes the explicit form where we have defined the projection operator P A,B for projecting onto the A or B sublattice Then we apply this transformation on the Hamiltonian U † 1,τ H F U 1,τ for all k to arrive at U † 2,τ U † 1,τ H F U 1,τ U 2,τ .Finally, in order to establish a scalar product structure between the isospin and momentum in the transformed Hamiltonian that mimics the equilibrium Hamiltonian of graphene, we apply the valley-independent unitary transformation U 3 , Eq. (B3).It satisfies the following properties, P α P ᾱ = 0, P α P α = P α , (D1) where ᾱ stands for the complement of α, i.e., ᾱ = B if α = A, and vice versa.Then the integrand in Eq. ( 14) can be written as where in the last line we made use of the cyclic property of the trace.We can express the Green's functions in the Floquet representation in terms of the K photon-number basis: Then Eq. (D2) can be expressed as Using the projection operator property P α P α = P α in Eq. (D1) and the cyclic properties of the trace we can rewrite Eq. (D6) as By replacing the integrand in Eq. ( 14) by Eq. (D7) we obtain at Eq. (36) in the main text.The Green's functions in Eq. ( 36) are block-diagonal in the Floquet index n and composed of the 2 × 2 Green's functions gn (R, ω), hence the trace in Eq. ( 36) results in a simple sum over n.By writing gn (R, ω) ≡ g(R, ω − n Ω) for each valley and using Eq. ( 38), Eq. ( 36) can be written as Making the change of variable ω = ω − n Ω, the Green's functions in Eq. (D8) become independent of n and the limits of integration become [−(1/2+n) Ω, (1/2−n) Ω].Then, by performing the sum over the Floquet index n we arrive at Eq. (37) where ω ∈ (−∞, ∞) in this equation is the extended zone frequency.Next, we discuss how the exchange energy in Eq. ( 37) reduces to its equilibrium counterpart in the absence of periodic driving.In the equilibrium limit A → 0, ∆ = 0, we have F 0 = 1 and n = 0, so that the Hamiltonian in Eq. ( 20) takes the form of H0 = v F (σ x k x − σ y k y ) for both the K and K points.Meanwhile, Eq. ( 18) gives h − = g + = 1 and h + = g − = 0, the isospin bases Eqs. ( 16) and ( 17) become (Φ ↑ , Φ ↓ ) T = σ x (φ A , φ B ) T for K and (Φ ↑ , Φ ↓ ) T = σ z (φ A , φ B ) T for K .Therefore, the equilibrium Hamiltonian H0 in the photon-number representation is related to that in the original pseudospin basis [Eq.(1)] by the following unitary transformations H 0 = σ x H0 σ x for K and H 0 = σ z H0 σ z for K .(D9) It follows that the equilibrium Green's functions in the photon-number representation are also related to those in the original pseudospin basis by g K = σ x gσ x for K and g K = σ z gσ z for K .(D10) In order to express Eq. ( 37) in terms of the Green's functions g K,K in the original pseudospin basis, we start by recalling that the Green's function in Eq. ( 38) is written in K isospin basis.Hence in equilibrium this Green's function transforms to the pseudospin basis as σ z gF0 (R, ω)σ z = g(R, ω), such that g(R, ω) = e iK•R σ z σ y g(R, ω)σ y σ z +e iK •R σ z g(R, ω)σ z . (D11) Since σ y σ z = −σ y σ z = iσ x , the first term on the right in the above equation becomes σ x g F0 σ x .Using Eq. (D10), Eq. (D11) immediately reduces to its expected equilibrium form g(R, ω) = e iK•R g K (R, ω) + e iK •R g K (R, ω) .(D12) Hence Eq. ( 37) in equilibrium becomes

FIG. 1 .
FIG. 1. Spectral function ρ(k, ω) [Eq.(23)] of Floquet states as a function of k and ω, obtained within the F0 approximation.The side color bar indicates ρ(k, ω) in logarithmic scale.Here the system is driven by a field with strength A = A/ Ω = 0.2 and frequency Ω = 2.2γ.

FIG. 2 .
FIG. 2. Honeycomb lattice structure of graphene.The filled circles represent the atoms of sublattice A while the unfilled ones represent the atoms of sublattice B. The zigzag and armchair directions in the honeycomb structure are indicated with arrows.

FIG. 3 .
FIG.3.Time-averaged RKKY coupling J(R)/J0 for undoped graphene under irradiation.The driving field is taken to have a frequency Ω = 6.6γ and strength A ∈ [0, 0.06], with A = 0 being the equilibrium case.The panels correspond to cases with the impurities located at different sublattice sites (A or B) and separated along different directions (zigzag or armchair): (a) AA zigzag; (b) AB zigzag; (c) AA armchair; (d) AB armchair.The inset in (a) shows the negative values of J(R)/J0 over the full range of R down to the smallest separation.J0 is given in Eq.(41).

FIG. 4 .
FIG.4.Time-averaged RKKY coupling J(R)/J0 for doped graphene under irradiation with a Fermi energy EF = 0.065γ.The driving field is taken to have a frequency Ω = 6.6γ and strength A = 0 and 0.06.For the latter value of A the photon-induced gap ∆ = 0.71EF .The panels correspond to cases with the impurities located at different sublattice sites (A or B) and separated along different directions (zigzag or armchair): (a) AA zigzag; (b) AB zigzag; (c) AA armchair; (d) AB armchair.The insets in (c) and (d) show the negative and positive values of J(R)/J0 respectively over the full range of R down to the smallest separation.

FIG. 5 .
FIG. 5. Time-averaged RKKY coupling J(R)/J0 for doped graphene under irradiation with a Fermi energy EF = 0.0236γ.The driving field has the same frequency and strengths as in Fig. 4. For A = 0.06 the photon-induced gap ∆ = 1.99EF .The panels correspond to cases with the impurities located at different sublattice sites (A or B) and separated along different directions (zigzag or armchair): (a) AA zigzag; (b) AB zigzag; (c) AA armchair; (d) AB armchair.The insets in (c) and (d) show the negative and positive values of J(R)/J0 respectively over the full range of R down to the smallest separation.

FIG. 6 .
FIG. 6.(a) Dependence of J(R)/J0 on the Fermi energy.In this case we the impurity spins to be separated by fixed distance R = 70 √ 3a on the same sublattice sites AA along the zigzag direction.Panel (b) shows an enlarged view of the small Fermi energy region marked by the rectangle in (a).The position of the conduction band edge at EF = ∆/2 is indicated by vertical line.The driving frequency is the same as in Figs.