RKKY Interactions in Graphene Landau Levels

TL;DR

This paper investigates how Landau levels in graphene influence RKKY interactions between magnetic impurities, revealing enhanced interactions at the lowest Landau level and exploring resulting magnetic phases and transitions.

Contribution

It introduces a non-perturbative analysis of RKKY interactions in graphene Landau levels, highlighting the role of the lowest Landau level and developing mean-field theories for strained and real magnetic field scenarios.

Findings

RKKY interactions are enhanced by the lowest Landau level.

Impurity interactions decay as 1/R^3 beyond the lowest Landau level.

Strained graphene exhibits a ferrimagnetic state with small net moment.

Abstract

We study RKKY interactions for magnetic impurities on graphene in situations where the electronic spectrum is in the form of Landau levels. Two such situations are considered: non-uniformly strained graphene, and graphene in a real magnetic field. RKKY interactions are enhanced by the lowest Landau level, which is shown to form electron states binding with the spin impurities and add a strong non-perturbative contribution to pairwise impurity spin interactions when their separation $R$ no more than the magnetic length. Beyond this interactions are found to fall off as $1/R^3$ due to perturbative effects of the negative energy Landau levels. Based on these results, we develop simple mean-field theories for both systems, taking into account the fact that typically the density of states in the lowest Landau level is much smaller than the density of spin impurities. For the strain field case, we find that the system is formally ferrimagnetic, but with very small net moment due to the relatively low density of impurities binding electrons. The transition temperature is nevertheless enhanced by them. For real fields, the system forms a canted antiferromagnet if the field is not so strong as to pin the impurity spins along the field. The possibility that the system in this latter case supports a Kosterlitz-Thouless transition is discussed.

Tables

Table 1. Table 1: Numerical values of spatially averaged RKKY coupling strength for a strain magnetic field of strength B = 10 𝐵 10 B=10 T, assuming parameters for Co as discussed in the Introduction. Top line denotes different impurity concentrations, and AA,AB,BB represent sublattice site locations of the impurities.

${\overline{J}}_{\text{RKKY}}^{\mu \nu }$ in eV ($\times {10}^{-5}$)	1%	2%	3%	4%	5%
AA	-0.8545	-2.2997	-4.0877	-6.1355	-8.3965
AB	1.5404	4.0948	7.2181	10.7650	14.6553
BB	-0.8273	-2.2618	-4.0441	-6.0888	-8.3484

Equations

$$H_0^{\tau }=v_F(\tau q_x{\sigma }_x+q_y{\sigma }_y),$$

$$\psi^{(\tau=+1)}_{n,k}=\frac{1}{\sqrt{2}}e^{+i\vec{K}\cdot\vec{r}}\left(\begin{array}[]{c}\phi_{n-1,k}\\ \text{sgn}(n)\phi_{n,k}\\ \end{array}\right),$$

$$\psi^{(\tau=+1)}_{n=0,k}=e^{+i\vec{K}\cdot\vec{r}}\left(\begin{array}[]{c}0\\ \phi_{n=0,k}\\ \end{array}\right)$$

$$\psi^{(\tau=-1)}_{n,k}=\frac{1}{\sqrt{2}}e^{-i\vec{K}\cdot\vec{r}}\left(\begin{array}[]{c}-\text{sgn}(n)\phi_{n,k}\\ \phi_{n-1,k}\\ \end{array}\right),$$

$$\psi^{(\tau=-1)}_{n=0,k}=e^{-i\vec{K}\cdot\vec{r}}\left(\begin{array}[]{c}\phi_{n=0,k}\\ 0\\ \end{array}\right)$$

$${\varphi }_{n,k}=\frac{1}{\sqrt{L_x}}{e}^{ikx}{e}^{-\frac{(y-y_0{)}^2}{2{\ell }^2}}{N}_{|n|}{H}_{|n|}(\frac{y-y_0}{\ell }),$$

$$e\vec{A}(\vec{x})=\frac{\beta}{a}\left(\begin{array}[]{c}u_{xx}-u_{yy}\\ -2u_{xy}\\ \end{array}\right)$$

$$\psi^{(\tau=\pm 1)}_{n,k}=\frac{1}{\sqrt{2}}\left(\begin{array}[]{c}\phi_{n-1,\pm k}\\ \pm\text{sgn}(n)\phi_{n,\pm k}\\ \end{array}\right),$$

$$\psi^{(\tau=\pm 1)}_{n=0,k}=\left(\begin{array}[]{c}0\\ \phi_{n=0,\pm k}\\ \end{array}\right).$$

$$V^{({\mu }_1,{\mu }_2)}\equiv V^{({\mu }_1)}+V^{({\mu }_2)}=J\left[{\vec{S}}_1\cdot \vec{s}\delta (\vec{r}-{\vec{R}}_1)P_{{\mu }_1}+{\vec{S}}_2\cdot \vec{s}\delta (\vec{r}-{\vec{R}}_2)P_{{\mu }_2}\right].$$

$$H_Z=H_Z^{(e)}+H_Z^{(imp)}=g_0{\mu }_{B}B{s}_z+g_{imp}^{(0)}{\mu }_B\sum _i\hat{z}\cdot {\vec{S}}_i.$$

$${\varphi }_{n=0,m}(z)\propto z^m{e}^{-\frac{|z{|}^2}{4}},$$

$${\tilde{\varphi }}_m(z)\propto (z^2-{\eta }^2)z^m{e}^{-\frac{|z{|}^2}{4}},m=0,1,....,M-3$$

$${\xi }_1(z)=\frac{1}{\sqrt{2\pi \sinh (\frac{{\eta }^2}{2})}}{e}^{-\frac{|z{|}^2}{4}}\sinh (\frac{\eta z}{2}),$$

$${\xi }_2(z)=\frac{1}{\sqrt{2\pi \cosh (\frac{{\eta }^2}{2})}}{e}^{-\frac{|z{|}^2}{4}}\cosh (\frac{\eta z}{2}).$$

$$\bar{H}$$

$$\displaystyle\;\left\{\left(\begin{array}[]{cccc}n^{(1)}_{z}\sqrt{\tanh(\frac{\eta^{2}}{2})}&n^{(1)}_{z}&\left(n^{(1)}_{x}-in^{(1)}_{y}\right)\sqrt{\tanh(\frac{\eta^{2}}{2})}&n^{(1)}_{x}-in^{(1)}_{y}\\ n^{(1)}_{z}&\frac{n^{(1)}_{z}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}&n^{(1)}_{x}-in^{(1)}_{y}&\frac{n^{(1)}_{x}-in^{(1)}_{y}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}\\ \left(n^{(1)}_{x}+in^{(1)}_{y}\right)\sqrt{\tanh(\frac{\eta^{2}}{2})}&n^{(1)}_{x}+in^{(1)}_{y}&-n^{(1)}_{z}\sqrt{\tanh(\frac{\eta^{2}}{2})}&-n^{(1)}_{z}\\ n^{(1)}_{x}+in^{(1)}_{y}&\frac{n^{(1)}_{x}+in^{(1)}_{y}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}&-n^{(1)}_{z}&-\frac{n^{(1)}_{z}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}\\ \end{array}\right)\right.$$

$$\displaystyle\left.+\left(\begin{array}[]{cccc}n^{(2)}_{z}\sqrt{\tanh(\frac{\eta^{2}}{2})}&-n^{(2)}_{z}&\left(n^{(2)}_{x}-in^{(2)}_{y}\right)\sqrt{\tanh(\frac{\eta^{2}}{2})}&in^{(2)}_{y}-n^{(2)}_{x}\\ -n^{(2)}_{z}&\frac{n^{(2)}_{z}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}&in^{(2)}_{y}-n^{(2)}_{x}&\frac{n^{(2)}_{x}-in^{(2)}_{y}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}\\ \left(n^{(2)}_{x}+in^{(2)}_{y}\right)\sqrt{\tanh(\frac{\eta^{2}}{2})}&-n^{(2)}_{x}-in^{(2)}_{y}&-n^{(2)}_{z}\sqrt{\tanh(\frac{\eta^{2}}{2})}&n^{(2)}_{z}\\ -n^{(2)}_{x}-in^{(2)}_{y}&\frac{n^{(2)}_{x}+in^{(2)}_{y}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}&n^{(2)}_{z}&-\frac{n^{(2)}_{z}}{\sqrt{\tanh(\frac{\eta^{2}}{2})}}\\ \end{array}\right)\right\}.$$

$$\frac{2\pi E_1}{JS}$$

$$\frac{2\pi E_2}{JS}$$

$$\frac{2\pi E_3}{JS}$$

$$\frac{2\pi E_4}{JS}$$

$$E_1=$$

$$+\frac{\pi JS}{(h^2+1){\ell }^2}(\frac{2{\left(h^2-1\right)}^{2}JS{g}_0{\mu }_{B}B(X+1)\left((n_z^{(1)})+(n_z^{(2)})\right)}{\pi {\ell }^2}$$

$$+\frac{(h^2-1{)}^2{J}^2{S}^2(X+1)\left(h^4(X+1)-2h^2(X-3)+X+1\right)}{{\pi }^2(h^2+1{)}^2{\ell }^4}$$

$${{(g_0{\mu }_{B}B{)}^2\left({\left(h^2+1\right)}^2\left((n_z^{(1)}{)}^2+(n_z^{(2)}{)}^2\right)+2\left(h^4-6h^2+1\right)n_z^{(2)}{n}_z^{(1)}\right))}^{1/2}]}^{1/2}$$

$$E_2=$$

$$-\frac{\pi JS}{(h^2+1){\ell }^2}(\frac{2{\left(h^2-1\right)}^{2}JS{g}_0{\mu }_{B}B(X+1)\left(n_z^{(1)}+n_z^{(2)}\right)}{\pi {\ell }^2}$$

$$+\frac{(h^2-1{)}^2{J}^2{S}^2(X+1)\left(h^4(X+1)-2h^2(X-3)+X+1\right)}{{\pi }^2(h^2+1{)}^2{\ell }^4}$$

$${{+(g_0{\mu }_{B}B{)}^2\left({\left(h^2+1\right)}^2\left((n_z^{(1)}{)}^2+(n_z^{(2)}{)}^2\right)+2\left(h^4-6h^2+1\right)n_z^{(2)}{n}_z^{(1)}\right))}^{1/2}]}^{1/2}.$$

$$E_1+E_2\approx$$

$$-\frac{JS\left(n_z^{(1)}+n_z^{(2)}\right)}{2\pi {\ell }^2}-g_0{\mu }_{B}B.$$

$$(E_1+E_2){|}_{\eta \to \infty }=\frac{J^2{S}^2\left((n_z^{(1)}{)}^2+(n_z^{(2)}{)}^2-2\right)}{4{\pi }^2{\ell }^4{g}_0{\mu }_{B}B}-\frac{JS\left(n_z^{(1)}+n_z^{(2)}\right)}{2\pi {\ell }^2}-g_0{\mu }_{B}B.$$

$$(E_1+E_2)-(E_1+E_2){|}_{\eta \to \infty }=-\frac{J^2{S}^2}{2{\pi }^2{\ell }^4{g}_0{\mu }_{B}B}{e}^{-2{\eta }^2}\left(n_x^{(1)}{n}_x^{(2)}+n_y^{(1)}{n}_y^{(2)}\right).$$

$$\Delta E^{(2)}=\sum _{\begin{array}{c}n<0,\\ n^{\prime }\ne n\end{array}}\sum _{k,k^{\prime }}\sum _{s,s^{\prime }}\sum _{\tau ,{\tau }^{\prime }}\frac{|V_{n^{\prime }{k}^{\prime }{s}^{\prime }{\tau }^{\prime },nks\tau }{|}^2}{E_n^{(0)}-E_{n^{\prime }}^{(0)}},$$

$$I_{n_1,n_2}({\vec{R}}_1,{\vec{R}}_2)=\frac{L_x}{2\pi }{\int }_{-\infty }^{+\infty }dk{\varphi }_{n_1,k}({\vec{R}}_1){\varphi }_{n_2,k}^{\ast }({\vec{R}}_2),$$

$$I_{n_1,n_2}({\vec{R}}_1,{\vec{R}}_2)=\frac{1}{2\pi {\ell }^2}\frac{1}{\sqrt{2^{n_1-n_2}(n_2!)/(n_1!)}}{e}^{-{\eta }^2+i{y}_m\Delta x/{\ell }^2}(-1{)}^{n_1-n_2}{\left(\frac{\Delta y+i\Delta x}{2\ell }\right)}^{n_2-n_1}{L}_{n_1}^{n_2-n_1}(2{\eta }^2).$$

Full text

I Introduction

Graphene is one of the most interesting platforms for the two-dimensional electron gas (2DEG) to have become available in the laboratory in recent years, both for its fundamental physics and for the potential applications it offers Castro Neto et al. (2009); Kotov et al. (2012); Aoki and Dresselhaus (2014). Among its many unique characteristics, the possibility that it can sustain magnetic order has been an ongoing subject of investigation. There seems to be little experimental evidence that pristine graphene has such order Sepioni et al. (2010), but theoretical studies strongly suggest that antiferromagnetic order can be sustained at ribbon edges Son et al. (2006) or among moments forming on vacancy defects in the structure Castro Neto et al. (2009); Palacios et al. (2010). To date there is no convincing observation of such order, and whether it is realized in the real material is unclear. One strategy that has been pursued to enhance magnetism in graphene is to combine it with magnetic impurities Dugaev et al. (2006); Brey et al. (2007); Saremi (2007); Black-Schaffer (2010); Fabritius et al. (2010); Sherafati and Satpathy (2011); Kogan (2011); Lee et al. (2012); Roslyak et al. (2013); Gorman et al. (2013); Crook et al. (2015); Min et al. (2017). In these systems, impurity magnetic moments locally couple to the electron spin density of the 2DEG, effectively coupling the impurity moments magnetically via RKKY interactions Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). When the graphene is doped, this leads to Heisenberg coupling $J_{RKKY}^{\mu,\nu}\vec{S}_{i}\cdot\vec{S}_{j}$ between impurity spins $i,j$, with the effective exchange constant behaving as $J_{RKKY}^{\mu,\nu}\sim\sin(k_{F}R_{ij})/R_{ij}^{2}$, where $\mu,\nu$ are the sublattice upon which the impurities at $i$, $j$ reside, $R_{ij}$ is the impurity separation, and $k_{F}$ the Fermi wavevector. In behavior analogous to that of vacancies Palacios et al. (2010), the sign of the coupling changes depending on relative sublattice, such that $J_{RKKY}^{A,A}=J_{RKKY}^{B,B}=-J_{RKKY}^{A,B}$ Brey et al. (2007).

The spatial oscillatory behavior reflects the presence of a Fermi surface, and allows impurities on the same (opposite) sublattice to couple (anti-)ferromagnetically up to a separation of order $R_{ij}\sim\pi/k_{F}$. By contrast, for undoped graphene, the point-like form of the Fermi surface Castro Neto et al. (2009) leads to a non-oscillatory form, $J_{RKKY}\sim 1/R_{ij}^{3}$ Brey et al. (2007), which is still opposite for opposite sublattices. The slow fall-off of this interaction has interesting consequences for spin stiffness in the system, for example introducing non-analyticity into the effective energy functional for spin gradients Reja et al. (2019). In principle this behavior can be modified by shifting the chemical potential of the system via an electric gate Reja et al. (2017, 2019). The possibility of controlling the magnetic properties in a single graphene system is one of the reasons it is of such intrinsic interest.

In this paper, we consider alternative strategies to modifying and controlling magnetism in graphene: application of non-uniform strain Levy et al. (2010); Jiang et al. (2017); Liu et al. (2018), or a magnetic field applied perpendicular to the system. These seemingly different modifications of the graphene system have in common the restructuring of the electronic spectrum into Landau levels. While uniform strain has quantitative but not qualitative effects Gorman et al. (2013) on RKKY interactions in graphene, non-uniform strain if applied appropriately can introduce an effective “pseudo”-magnetic field into the low-energy Hamiltonian Guinea et al. (2009); Low and Guinea (2010), with the effective field directed oppositely for the two valleys of the graphene band structure. Such fields have indeed been created in graphene bubbles Lu et al. (2012); Klimov et al. (2012), and are possible in artificially structured graphene analogs Gomes et al. (2012); Tian et al. (2015). By its nature, a strain-induced pseudofield couples to the electronic orbital degrees of freedom, but to neither the spin of the electrons nor of the impurities. A real field, by contrast, couples to spin as well orbital degrees of freedom, but for low fields this does not fully polarize the impurity spin density, so that non-trivial order can set in, as we explain below.

The RKKY analysis for these systems differs in important ways from that of graphene in the absence of a field Brey et al. (2007). The key reason for this is the existence of a zero-energy Landau level Gusynin and Sharapov (2005); Goerbig (2011) in which the Fermi level resides if the system is not strongly doped. The standard RKKY analysis, which depends on second order perturbation theory Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957), becomes invalid because of the large degeneracy associated with the Landau level. We demonstrate, however, that the analytic properties of a Landau level allow one to compute the re-organization of the energy states in the lowest Landau level due to the presence of two spin impurities essentially exactly, introducing four bound states that separate off from the degenerate Landau level, two of which are filled at the electron densities we consider. The bound states introduce a spin coupling between the two impurities which scales linearly with the $sd$ coupling constant $J$ between an impurity spin and the 2DEG, and so is formally considerably stronger than standard RKKY interactions, which scale as $J^{2}$. However, the range of this coupling is limited, falling off as a Gaussian with length scale $\ell=\sqrt{\hbar c/eB}$, with $B$ the effective field. Importantly, the remaining Landau levels in the spectrum may be taken into account perturbatively, yielding an interaction of the same Heisenberg form $\vec{S}_{i}\cdot\vec{S}_{j}$ as in zero field, which is of magnitude $J^{2}$, but falls off much more slowly, as $1/R^{3}$.

An interesting difference between the Landau levels in the strain case and those of the applied field case is that, in the former, the non-analytic contribution only applies to the spins on the same sublattice, so that the RKKY couplings on one sublattice are stronger than on the other. This raises the possibility that the system could sustain a net magnetic moment, so that the order (at least at the mean-field level) is ferrimagnetic. In principle this net magnetic moment would make detection of magnetism via magnetization measurements in these graphene systems considerably easier than the antiferromagnetism expected of a perfect graphene lattice Brey et al. (2007).

To investigate this last effect, we perform a mean-field analysis for strained graphene that supports Landau levels. For reasonably size fields, we find that the density of impurities is actually quite large compared to the density of states in the lowest Landau levels, so that a model with effective pairwise interactions between individual impurities becomes inappropriate. To deal with this situation, our mean-field theory is developed in terms of two sets of spin impurities, ones that bind electrons in the lowest Landau level, and ones that do not. Within this model we compute the net magnetization as a function of temperature on each of the sublattices, and find that for reasonable impurity densities and couplings to the electron system, their difference is quite small, so that direct detection of magnetic order remains a challenge for these systems.

We also analyze the mean-field phases for the real magnetic field case. The presence of Zeeman coupling in both the impurity and electron spins always induces a magnetization component along the direction of the field. However, because of the antiferromagnetic coupling between impurities on opposite sublattices, at weak applied field the mean-field state is a canted antiferromagnet, with broken U(1) symmetry, due to spontaneous ordering of the magnetization component of impurity spins transverse to the applied field. Formally such long-range order is unstable to thermal fluctuations at any finite temperature Mermin and Wagner (1966), but the system can nevertheless support a true thermodynamic phase transition due to vortex excitations of the U(1) degree of freedom, which support a Kosterlitz-Thouless transition Brink et al. (2018). This is a distinguishing feature of the graphene - impurity spin system when it is embedded in a real magnetic field, and we present estimates for the Kosterlitz-Thouless transition temperature below for reasonable impurity densities and couplings.

This article is organized as follows. In Section II we introduce the basic model used for our analysis of RKKY couplings. Section III is focused on an analysis of the coupling strength for a single pair of spin impurities coupled to graphene electrons with a Landau level spectrum. In Section IV we present our mean-field theory for the impurity magnetization for the Landau levels produced by strain, and in Section V we present the corresponding analysis for Landau levels produced by a real field. Finally in Section VI we present a summary and a discussion of implications of and speculations about our results.

II Continuum Hamiltonian

We begin with a description of the models we adopt to anaylze RKKY interactions of spin impurities coupled to graphene Landau levels. We adopt a continuum model for the Hamiltonian in the vicinity of a Dirac point at the center of a valley,

[TABLE]

where $\vec{q}$ is the momentum relative to the $K$ ($\tau=1$) or $K^{\prime}$ ($\tau=-1$) points, $v_{F}$ is the speed of the electrons in their vicinity, and $\sigma_{x,y}$ are Pauli matrices acting on spinors whose entries encode the wavefunction amplitude on the $A$ and $B$ sublattices. In these equations and what follows, we have set $\hbar=1$. Magnetic fields, be they effective fields due to strain or a real magnetic field, are introduced into the orbital Hamiltonian by the Peierl’s substitution, $\vec{q}\rightarrow\vec{q}\pm e\vec{A}$, where $\vec{A}$ is the vector potential, which for our purposes corresponds to one associated with a uniform magnetic field. Because $\vec{A}$ is position dependent, the momentum must now be regarded as an operator, $\vec{q}\rightarrow-i\vec{\nabla}$. In the case of a real field, $\vec{q}\rightarrow\vec{q}+e\vec{A}$, and choosing $\vec{A}=-By\hat{x}$, eigenstates of $H_{0}^{\tau}$ have the form Goerbig (2011)

[TABLE]

for $n\neq 0$, and

[TABLE]

for $n=0$, in the $K$ valley. For the $K^{\prime}$ valley, the corresponding expressions are

[TABLE]

for $n\neq 0$, and

[TABLE]

for $n=0$. In these expressions, the wavevector $\vec{K}\equiv\frac{4\pi}{3\sqrt{3}a}\hat{x}$ denotes the position of the $K$ point relative to the $\Gamma$ point in the Brillouin zone, with $a=0.142$ nm the nearest neighbor carbon bond length, while the $K^{\prime}$ point is located at $-\vec{K}$. The functions $\phi_{n,k}$ are Landau level states localized near a guiding center coordinate $y_{0}=k\ell^{2}$,

[TABLE]

with $L_{x}$ the system size in the $\hat{x}$ direction, and $N_{|n|}\equiv\sqrt{1/{\pi^{1/2}\ell 2^{n}|n|!}}$ a normalization constant. $H_{|n|}$ are Hermite polynomials. The eigenvalues of $H_{0}^{\tau}$ associated with these wavefunctions are $E_{n,k}^{\tau}\equiv E_{n}^{(0)}=\text{sgn}(n)v_{F}\sqrt{2\hbar eB_{s}n}.$

For the case of strain-induced magnetic fields, we follow the approach developed in Ref. Guinea et al., 2009. Briefly, this involves encoding lattice distortions that vary very slowly on the scale of the graphene lattice constant in a vector potential given by a strain tensor $u_{ij}$, with

[TABLE]

and $\beta=-\frac{\partial\ln t}{\partial\ln a}\approx 2$ specifying the change in the nearest neighbor tunneling parameter $t$ when the lattice constant $a$ changes. The effective gauge field couples with ${\it opposite}$ signs for the two valleys: $\vec{q}\rightarrow\vec{q}+e\vec{A}$ for the $K$ valley ($\tau=1$), while $\vec{q}\rightarrow\vec{q}-e\vec{A}$ for the $K^{\prime}$ valley ($\tau=-1$). The eigenstates of $H_{0}^{\tau}$ in this case are

[TABLE]

for $n\neq 0$, while for $n=0$,

[TABLE]

Note the important distinction from the real magnetic field case: the support for the two valleys is essentially the same for each sublattice, whereas in the real field case the roles of the sublattices is switched. The difference reflects the fact that the strain breaks the inversion symmetry of the graphene lattice, so that the two sublattices come in asymmetrically in the Hamiltonian. This ultimately opens the possibility that the magnetization can have different magnitudes on each of the sublattices, with a net magnetization resulting.

To analyze RKKY coupling between impurities carrying spin by the electrons, we will consider Hamiltonians with two impurities at locations $\vec{R}_{1,2}$ on specified sublattices. For $\mu_{1,2}=A,B$, the coupling takes the form

[TABLE]

In this expression, $J$ is an assumed $sd$ coupling constant between the impurity spins and the electron gas, $\vec{S}_{1,2}$ are the impurity spins (assumed classical, as is standard in RKKY analyses Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957)), $\vec{s}$ is the electron gas spin operator, and $P_{\mu}$ is a projection operator onto the $\mu$ sublattice. Note that the exchange constant $J$ can vary widely depending on the type of impurity adsorbed on the surface. For quantitative estimates given below we will adopt a value appropriate for Co when bound to individual carbon atoms (and so to a particular sublattice), $Ja_{C}\approx 1$eV, where $a_{C}$ is the area per carbon atom in the graphene lattice, with effective spin $S=3/2$ Fritz and Vojta (2013).

For strained graphene, the terms above are sufficient, and the effective Hamiltonian for a two impurity system is $H_{strain}=\sum_{\tau}H_{0}^{\tau}+V^{(\mu_{1},\mu_{2})}$, with the vector potential properly substituted into $H_{0}^{\tau}$. In a real magnetic field, one must also account for the Zeeman coupling between the field and the electrons, as well as the impurities. This introduces terms of the form

[TABLE]

Note that $g_{0}\approx 2$ for electrons in graphene, and $g_{imp}^{(0)}\approx 2$ for Co adatoms Fritz and Vojta (2013). For the system in a real magnetic field, the Hamiltonian for a pair of impurities adsorbed on graphene becomes $H_{field}=\sum_{\tau}H_{0}^{\tau}+V^{(\mu_{1},\mu_{2})}+H_{Z}$.

We now turn to the computation of the effective RKKY coupling between two impurity spin degrees of freedom adsorbed on graphene in these situations.

III RKKY Interaction

As described in the introduction, the computation of RKKY interactions in this problem necessarily involves a non-perturbative contribution, due to the high degeneracy of a Landau level. In particular, if we assume the system to be only moderately doped, so that the Fermi energy lies in the $n=0$ Landau level states – which we will from hereon refer to as the lowest Landau level (LLL) – then one must understand how these levels become energetically organized in the presence of the impurities. We begin by showing how this can be done.

III.1 Lowest Landau Level Energies: Exact Solution

We begin with the case of a vector potential induced by strain. From the form of Eq. 8, it is apparent that these states can only couple impurities together when they are both on the $B$ sublattice, for the specific form of strain we consider, and we focus for the moment on this case. With a change to circular gauge, states in the LLL can be written in the form Yoshioka (2002) $(0,\phi_{n=0,m}(z))^{{\dagger}}$, where $z=(x-iy)/\ell$, and

[TABLE]

where $m$ is an angular momentum index. These states have the interesting property that they are peaked at a distance $r_{m}=\sqrt{2m}\ell$, and have a width of $\ell$. A generic state in the LLL takes the form $(0,f(z)e^{-\frac{|z|^{2}}{4}})^{{\dagger}}$, with $f$ an analytic function in $z$. For a finite system, a natural requirement is that a power law expansion of $f$ contains terms of order smaller than some (large) integer $M$. The dimension of the LLL with this condition is $M$.

Now suppose we place the two impurities at positions $\vec{R}_{1,2}=\pm\eta\vec{e}_{x}$. Then states of the form

[TABLE]

are completely decoupled from the impurities, and will have the same energy as in their absence. It is easy to see that the dimension of this subspace is $M-2$; this means that the set of states affected by the impurities within the LLL can be reduced to a single pair, which must be orthogonal to the states in Eq. 12. Remarkably, in the limit $M\rightarrow\infty$ these states may be written explicitly, and take the form

[TABLE]

Note the states above are orthonormal.

These states as written do not include spin. When taken into account, we have reduced the problem of finding the energy spectrum for electrons in the presence of the two impurities, when projected into the LLL, into a $4\times 4$ matrix diagonalization problem. The energies of these 4 states within the LLL will be sensitive to the relative orientation of the two spins, and this determines the LLL contribution to the RKKY interaction between the spins. Since the LLL states have zero energy in the absence of the impurities, the Hamiltonian for this case can be taken to simply be $V^{(B,B)}$ (see Eq. 9). Writing the single particle states in the order $(|\xi_{1},\uparrow\rangle,|\xi_{1},\downarrow\rangle,|\xi_{2},\uparrow\rangle,|\xi_{2},\downarrow\rangle)$, the projected Hamiltonian becomes

[TABLE]

In this expression, we have written the two classical spin vectors in the form $\vec{S}_{i}=S\vec{n}^{(i)}$ with $i=1,2$. Notice the result is linearly proportional to $J$ and $S$; the contribution to the RKKY interaction from these states will be linear in these quantities, whereas standard RKKY interactions usually contain only contributions quadratic in these.

If we choose our spin $\hat{z}$ axis to be along the direction $\hat{n}^{(1)}$, and write $\hat{n}^{(2)}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$, the energy eigenvalues of $\bar{H}$ can be written explicitly in the form

[TABLE]

Fig. 1 illustrates the resulting LLL energy structure resulting from this analysis. Most states remain at zero energy in spite of the two impurities, but two break away to negative energy, and two to positive energy. If the Fermi energy $E_{F}$ is in the main band of zero energy states, then the two states at $E_{1}$ and $E_{2}$ lower the total electronic energy of the system. As can be seen from the explicit forms in Eqs. 15, the amount by which this energy is lowered depends on the relative orientation of the two spins.

Fig. 2 illustrates the behavior of the four energy levels as a function of separation $\eta$ and relative orientation angle $\theta$. Notice that the total lowering of the energy is always maximized when the two spins are aligned, so that this contribution to the spin-spin coupling is ferromagnetic. This is the same sign of coupling for spins on the same sublattice in the absence of a magnetic field Brey et al. (2007).

The situation for electrons in a real magnetic field is similar, but one must include the Zeeman term $H_{Z}^{(e)}$ (Eq. 10) in the Hamiltonian. In this case the induced interaction by the LLL is the same for two impurities on the same sublattice, with one of the two valleys inducing the interaction in each case (see Eqs. 3 and 5.) (Spin interactions from the LLL on opposite sublattices should be negligibly small due the rapidly oscillating relative phase factors in the real space wavefunctions, $\sim e^{\pm 2i\vec{K}\cdot\vec{r}}$.) The resulting energies for the two filled negative energy states for impurities on the same sublattice are

[TABLE]

and

[TABLE]

In these expressions, $X\equiv\cos\theta$ and $h=\sqrt{\tanh(\frac{\eta^{2}}{2\ell^{2}})}$. These rather complicated expressions can be considerably simplified for most physically relevant situations, for which the ratio $J/\pi\ell^{2}g_{0}\mu_{B}B$ is small; for example, for our Co estimate we obtain $\sim 0.1$ for this ratio. (Note this ratio is independent of the field strength $B$.) In this case an expansion in this ratio yields for the energy of the two filled states, to second order,

[TABLE]

For large separations ($\eta\rightarrow\infty$), the two spins decouple, but the energy remains dependent on the individual spin orientations:

[TABLE]

Neglecting the constant term, we see that the ${\cal O}(J)$ term effectively renormalizes the impurity gyromagnetic ratio, $g_{imp}^{(0)}\rightarrow g_{imp}=g_{imp}^{(0)}+g_{1}$, with $g_{1}=-J/\pi\ell^{2}g_{0}\mu_{B}B$. The ${\cal O}(J^{2})$ term creates a spin anisotropy favoring an in-plane spin orientation. The form for the interaction when the asymptotic energy ($\eta\rightarrow\infty$) is removed is particularly simple:

[TABLE]

Our results for the LLL contribution to the RKKY interaction are summarized in Fig. 3. For effective magnetic fields generated by strain, a ferromagnetic interaction between spins on one of the two sublattices is generated. This interaction scales linearly in $J$ and so is relatively strong for small $J$ at length scales below $\ell$, but falls off rapidly above this length scale. At such larger distances, the RKKY interaction becomes dominated by by the contributions from other Landau levels, as we will discuss in the next subsection. In the case of real magnetic field, the effect of the electron Zeeman coupling simplifies the behavior from the LLL, introducing an effective renormalization of impurity spin $g$-factor at linear order in $J$, and inducing anisotropy in the spin-spin interaction at quadratic order. Again, this contribution falls off rapidly for impurity separations large than $\ell$.

In both cases, for large impurity separations the RKKY interaction is dominated by contributions from $n\neq 0$ Landau levels. These can be handled in perturbation theory, as we discuss in the next subsection.

III.2 $n\neq 0$ Levels: Perturbation Theory

The underlying physics of RKKY interactions between impurity spins is that the energies of electrons in the system are modified in a way that depends upon their relative spin orientation. In many systems this can be handled at second order in perturbation theory Ruderman and Kittel (1954); Kasuya (1956); Yosida (1957). Because of the degeneracy of the partially filled $n=0$ Landau level, these states had to be handled carefully, as discussed above. The remaining levels can be handled in the more standard fashion. For $n<0$, it is possible to show that contributions to the energy at linear order in $J$ from $V^{(\mu_{1},\mu_{2})}$ (Eq. 9) vanish from these levels since they are completely filled, and are thus singlets in the electron spin. The first non-vanishing contributions come at second order, with the total change in energy of these filled levels given by

[TABLE]

with $V_{n^{\prime}k^{\prime}s^{\prime}\tau^{\prime},nks\tau}=\langle n^{\prime}k^{\prime}s^{\prime}\tau^{\prime}|V^{(\mu_{1},\mu_{2})}|nks\tau\rangle$ the matrix element of the perturbation to the electron gas of the two impurity spins, and $n$,$k$,$s$, and $\tau$ are respectively the Landau level, wavevector (proportional to guiding center coordinate $y_{0}$), spin, and valley indices of a state. The energies $E_{n}^{(0)}$ are as given in Sec. II. Note that we do not include the electron Zeeman contribution in the energy denominator, as this is in practice quite small in comparison with $|E_{n}^{(0)}-E_{n^{\prime}}^{(0)}|$, and so may also be handled perturbatively for the situation of real magnetic field. We begin first with the case of Landau levels produced by non-uniform strain.

III.2.1 Strain-Induced Landau Levels

We begin with the case where the impurities are both on the $B$ sublattice. To compute $\Delta E^{(2)}$ for this case, which we will call $E^{(2)}_{BB}$, we need to compute sums over $k$, which in the thermodynamic limit ($L_{x}\rightarrow\infty$) become integrals. The relevant integral has the form

[TABLE]

and can be evaluated to yield, for $n_{1}\leq n_{2}$,

[TABLE]

In this expression, $y_{m}=(\vec{R}_{1}+\vec{R}_{2})\cdot\hat{y}$, $\Delta x=(\vec{R}_{1}-\vec{R}_{2})\cdot\hat{x}$, $\Delta y=(\vec{R}_{1}-\vec{R}_{2})\cdot\hat{y}$, and $L_{n}^{m}$ is an associated Laguerre polynomial. For $n_{2}<n_{1}$, the result is the same , with $n_{1}\leftrightarrow n_{2}$, and $y_{1}\leftrightarrow y_{2}$. In the case $n_{1}=n_{2}$, $I_{n_{1},n_{2}}(\vec{R}_{1},\vec{R}_{2})$ depends only on $|\vec{R}_{1}-\vec{R}_{2}|=2\eta$, and to simplify the notation we write $I_{n,n}(\vec{R}_{1},\vec{R}_{2})=I_{n}(\eta)$. Summing over Landau level index, spin, and valley, after considerable algebra we arrive at the expression

[TABLE]

In light of Eq. 23, this can be cast in the form

[TABLE]

where $L_{n}(2\eta^{2})\equiv L_{n}^{0}(2\eta^{2})$ is a Laguerre polynomial. The Landau index sums appearing in Eq. 24 can be computed numerically, and must be cut off at a maximum value that is determined by the density of electrons in $p_{z}$ orbitals in graphene, which in turn is the density of carbon atoms. A similar calculation for impurities both on $A$ sites yields

[TABLE]

while for one impurity on an $A$ site and the other on a $B$ site we obtain

[TABLE]

Writing $E^{(2)}_{\mu\nu}\equiv J_{RKKY}^{\mu\nu}\vec{S}_{1}\cdot\vec{S}_{2}$, we show representative results for $J_{RKKY}$ in Fig. 4 as a function of $\eta=|\vec{R}_{1}-\vec{R}_{2}|/2\ell$. Several points are worth noting. As in the case of zero magnetic field, the coupling between spins on the same sublattice is ferromagnetic, while on opposite ones they are antiferromagnetically coupled. The presence of rapid oscillations can be traced back to the need for a cutoff in the sum over Landau levels, with the largest (negative) Landau index retained determined to give the correct overall electron density. For the most part the oscillations do not change the sign of the coupling; in the few cases where it does this leads neither to a change in the average sign of the coupling or the sign at very large distances. Finally, the overall scale of the interaction falls off as $1/\eta^{3}$, as illustrated in 5.

III.2.2 Magnetic Field-Induced Landau Levels

The case of Landau level states from a real magnetic field differs from the strain-induced ones in that there is inversion symmetry, so that the wavefunctions are symmetric under an interchange of $K$ and $K^{\prime}$ and the $A$ and $B$ sublattices. Moreover, the energies of the electron states in this case are spin-dependent due to the Zeeman coupling. While the latter effect can in principle be accounted for in $E^{(2)}$ without further approximations beyond the perturbation theory we are using in $V^{(\mu_{1},\mu_{2})}$, in practice the Zeeman splitting is very small compared to the Landau level energy differences without it at any achievable laboratory magnetic field. Thus it is sufficient and simplifying to treat the electron Zeeman energy perturbatively as well.

Our perturbation now takes the form

[TABLE]

and, working to linear order in the Zeeman coupling, we use

[TABLE]

in the numerator of Eq. 21. (Note that retaining the second order term in $V^{(z)}$ adds a contribution to the energy that is independent of the relative orientations of $\vec{S}_{1}$ and $\vec{S}_{2}$, and so is irrelevant for our current purpose.) The computation for the first term of Eq. 27 runs very similarly to that of the last subsection, and yields the results

[TABLE]

and

[TABLE]

where $\Delta\theta$ is the angle between the relative position vector $\vec{R}_{1}-\vec{R}_{2}$ and the $\hat{x}$-direction. Again writing $E^{(2)}_{\mu\nu}\equiv J_{RKKY}^{\mu\nu}\vec{S}_{1}\cdot\vec{S}_{2}$, Fig. 6 illustrates representative results.

Finally, the electron Zeeman terms yields a contribution to the energy of the form

[TABLE]

with $\varepsilon_{0}\equiv\sqrt{2e\hbar B}v_{F}$. Summing through the discrete indices produces a result that is independent of which sublattice the impurity spin resides upon, and has the simple form

[TABLE]

which is effectively a further renormalization of the impurity $g$-factor. Note the sum has an upper cutoff determined by the electron density.

IV Magnetic Order for Strain-Induced Landau Levels: Mean-Field Theory

To assess the effect of the electronic Landau level structure on magnetic order in the impurity spins, we consider simple mean-field theories of the magnetization. We begin in this section with the case of Landau levels induced by strain. The simplest approach to this would be treat the system as a system of classical spins, using the RKKY interactions computed in the previous section as a model for pairwise spin interactions. However, in doing this one assumes that the number of states in the LLL is sufficiently large to provide one bound electron state for every spin impurity on one of the two sublattices, which we take to be the $B$ sublattice. In practical situations this turns out not to be the case. Equating the density of spin impurities on one sublattice, $n_{imp}$, to the degeneracy of a Landau level per unit area (including spin), $1/\pi\ell^{2}$, produces a minimum magnetic field scale (in Tesla) $B_{c}\approx 7.9\times 10^{4}\tilde{n}_{imp}$(T), where $\tilde{n}_{imp}$ is the ratio of impurity atom density on one of the sublattices to graphene carbon atom density. With $\tilde{n}_{imp}$ typically being of order a few percent, we see that the effective field would need to exceed $\sim 1000$T to reach this limit. To date, strain-induced fields Levy et al. (2010) are at most of order several hundred Tesla, so that we need to consider situations with fewer LLL states than impurities – typically, much fewer.

To handle this, we consider a mean field theory in which only a fraction of impurities, $f_{b}=1/\pi\ell^{2}n_{imp}$, residing on the $B$ sublattice actually form bound states from the lowest Landau level, while all the spin impurities on both sublattices interact with one another through the perturbative contributions of the $n\neq 0$ Landau levels. Denoting the impurity spin degree of freedom at such a bound site on the $B$ sublattice as $\vec{S}_{b}$ and a spin degree of freedom on one of the remaining $B$ sites as $\vec{S}_{g}$, the total average spin direction for moments on the $B$ sites becomes

[TABLE]

where $\langle\cdot\cdot\cdot\rangle$ here denotes an average over sites and thermal fluctuations. The corresponding average magnetization (normalized to unity) for spins on the $A$ sublattice is denoted by $\vec{M}_{A}$. Referring to Eqs. 25 and 26, we see the pairwise interactions of a spin on an $A$ site at position $\vec{R}_{i}$ and either an $A$ or $B$ site at position $\vec{R}_{j}$ may be written in the forms $E^{(2)}_{AA}=J_{RKKY}^{AA}(|\vec{R}_{i}-\vec{R}_{j}|)\vec{S}_{i}\cdot\vec{S}_{j}$ and $E^{(2)}_{AB}=J_{RKKY}^{AB}(|\vec{R}_{i}-\vec{R}_{j}|)\vec{S}_{i}\cdot\vec{S}_{j}$, respectively. To form a single spin average, we adopt a simple model pair distribution function

[TABLE]

of finding an impurity on one of the sublattices at a displacement $\vec{R}$ within an area $d^{2}R$, given that there is an impurity at the origin. Assuming the impurity at the origin is on the $A$ sublattice, we can then write an average energy functional for its spin of the form

[TABLE]

where $\hat{S}_{A}$ denotes the orientation of the spin at the origin, and

[TABLE]

Similarly, for a $B$ site lacking bound electrons, the effective energy functional is

[TABLE]

The $\bar{J}_{\text{RKKY}}^{\mu\nu}$ coefficients can be computed numerically using our results from the previous section and our model $P_{\text{imp}}$. Results of these calculation are presented in Table 1.

To write an effective energy functional for a spin at a site binding an electron, we need to revisit the quantum problem yielding the bound state energy. Recalling in deriving the RKKY interaction for a single pair of spin impurities, we computed the bound state energies for two electrons interacting (through the $sd$ Hamiltonian) with classical spins localized at two sites. For our mean-field estimate, we will consider the microscopic potential due to a single spin in this collection, taken to be at the origin, and the average potential from all other sites on the $B$ sublattice, which carries average magnetic moment per unit area $n_{\text{imp}}S\vec{M}_{B}$. The potential due to the impurity at the origin has the form

[TABLE]

which, when projected into the LLL, couples only to the $m=0$ angular momentum state (see Eq. 11). Projecting into this one spatial state, the mean-field Hamiltonian for the electron spin becomes

[TABLE]

where we have defined the quantity

[TABLE]

Eigenstates of $H_{e}$ are easily seen to have energy

[TABLE]

With a bound electron occupying the lower energy state, and adding in the perturbative contribution to the RKKY interaction, we arrive at a mean-field energy functional for impurity spins on the $B$ sublattice with bound electrons of the form

[TABLE]

where $g_{v}$=2 is the number of valleys.

Our (classical) mean-field theory now proceeds by computing the average normalized magnetization, which without loss of generality we can assume to lie in the $\hat{z}$ direction, $M_{A}=\langle\hat{z}\cdot\hat{S}_{A}\rangle$, $M_{g}=\langle\hat{z}\cdot\hat{S}_{g}\rangle$, and $M_{b}=\langle\hat{z}\cdot\hat{S}_{b}\rangle$, using Boltzmann probability distributions at temperature $T$ proportional to $e^{-E_{A}/k_{B}T}$, $e^{-E_{g}/k_{B}T}$, and $e^{-E_{b}/k_{B}T}$, respectively, with $M_{B}$ related to $M_{b}$ and $M_{g}$ by Eq. 31. This set of equations can be straightforwardly solved numerically.

Figs. 7 and 8 illustrate typical results for effective fields of $10T$ and $2.5T$, respectively. Two features of the curves are of particular note. First, $M_{A}$ and $M_{B}$ nearly cancel, so that the correlations between sublattices are largely antiferromagnetic, as is the case for unstrained graphene Brey et al. (2007). However, $M_{b}$, the normalized magnetization on the bound electron sites, remains non-zero to larger temperatures than would be the case without the LLL contribution. The effect survives to increasingly high temperature as the field increases, and significantly increases the mean-field transition temperature. Unfortunately, because the number of bound electron sites is small, direct observation of this effect is challenging. While for any $T>0$, $M_{A}-M_{B}\neq 0$, indicating the system is formally a ferrimagnet, the difference is very small, yielding only a small net magnetic moment, as illustrated in Fig. 9. This occurs because $B$ is rather small compared to $B_{c}$, so that there is a relatively small number of electrons available in the LLL to bind to the impurity spins. An interesting possibility to circumvent this might be to consider non-uniform strain with small local regions of very high effective field, where locally an appreciable net magnetization could set in. Nevertheless, it is interesting that even in the small (relative to $B_{c}$) uniform fields we consider, the enhancement of $T_{c}$ is considerable. This may indicate that the effect of $M_{b}$ could be seen in thermodynamic measurements which are sensitive to spin fluctuations.

V Magnetic Order for Field-Induced Landau Levels: Mean-Field Theory

We next turn to the formulation of mean-field theory for the case of a real magnetic field. As discussed in Section III, the main complication in this situation is the introduction of Zeeman coupling to the spins, both that of the electrons and those of the impurity spins. Moreover, in this more symmetric situation, bound electrons may appear on either sublattice. The modifications for the basic energy functionals from the last section are nevertheless straightforward. For $A$ and $B$ sites without bound electrons we have

[TABLE]

and

[TABLE]

respectively. For an electrons bound to an impurity on sublattice $\mu$ the Hamiltonian (cf. Eq. 35) is modified to

[TABLE]

with eigenenergies

[TABLE]

in which

[TABLE]

and

[TABLE]

The resulting energy functional for sites binding an electron takes the form

[TABLE]

where $\bar{\mu}=B(A)$ if $\mu=A(B)$. In analogy with the previous section, $\hat{S}_{\mu,b}$ and $\hat{S}_{\mu,g}$ are thermally averaged with these energy functionals, and we search for self-consistent solutions satisying

[TABLE]

Because of the antiferromagnetic coupling between impurity spins on the $A$ and $B$ sublattices, in the presence of the Zeeman coupling to the spins we do not expect them to be collinear; the mean-field state should be a canted antiferromagnet, in which $\vec{M}_{A}$ and $\vec{M}_{B}$ have parallel components along $\vec{B}$, and antiparallel components perpendicular to it. It is important to recognize that this is a broken symmetry state, with the planar angle of the latter components in the groundstate determined arbitrarily; i.e., the state has broken U(1) symmetry. Without loss of generality, for the purposes of the mean-field solutions we can take the magnetizations to lie in the $\hat{x}-\hat{z}$ plane, assuming $\vec{B}=B\hat{z}$. Note that because of the explicit Zeeman coupling, there will always be non-zero components of $\vec{M}_{A}$ and $\vec{M}_{B}$ along the $\hat{z}$ direction at any temperature. The spontaneous ordering is captured by the non-zero $\hat{x}$ components of these.

Figs. 10 and 11 illustrate typical results at $B=0.401$T. At low temperature $\hat{z}\cdot\vec{M}_{A}=\hat{z}\cdot\vec{M}_{b}$ becomes pinned to a value less than one at these fields, indicating that the spins have become canted. Above a crossover temperature this value begins to fall, showing that polar fluctuations of the spin direction begin to become important. By contrast, $\hat{x}\cdot\vec{M}_{A}=-\hat{x}\cdot\vec{M}_{b}$ falls continuously with temperature, reflecting the behavior of in-plane spin fluctuations. This component truly drops to zero at a mean-field transition temperature, and the broken U(1) symmetry of the spin ordering is restored. Fig. 12 illustrates representative mean-field phase diagrams at different magnetic fields.

There are well-known limitations to the use of mean-field theory results for systems with broken continuous symmetries in two dimensions, such as in our analyses in this and the previous section. We now turn to a discussion of these, as well as other implications and speculations related to our study.

VI Discussion and Summary

In the analyses above, we considered pairwise RKKY interactions between spins in two situations where the electronic states are organized as Landau levels, in one case generated by non-uniform strain in the system, in the other by a magnetic field. We found that when the Fermi energy is in the lowest Landau level (LLL), the large degeneracy of states leads to a situation in which two pairs of states break off from the LLL, two above and two below the Fermi energy, whose precise energies depend on the relative orientations of the impurities. This means that the impurity pair binds a pair of electrons, which leads to a relatively strong interaction between them, up to a distance of order the magnetic length. Beyond this, the RKKY interactions are dominated by perturbative changes to the negative energy states, which on a coarse-grained scale leads to ferromagnetic interactions on the same sublattice, and antiferromagnetic interactions across sublattices. Rapid oscillations around these behaviors can be traced back to the effect of the negative energy cutoff, which is determined by the total density of electrons in the $p_{z}$ orbitals of the carbon atoms.

In practical situations, it is not sufficient to treat dilute spin impurities on the graphene as simply coupled by the RKKY interaction formally computed above, because the density of electrons in the LLL at field scales that are practically realizable is quite small compared to the density of impurities. In this situation we expect most impurities will not capture bound electrons. We considered the implications of this in a simple mean-field theory, in which a fraction of impurity spins $f_{b}$ binds electrons, of sufficient number to precisely deplete the LLL. We found that the effect of this can lead to a significant increase in the mean-field $T_{c}$, but that the magnetization above $T_{c}$ for the situation with $f_{b}=0$ is very small, because $f_{b}$ is typically quite small. Interestingly, in the case of magnetic fields due to non-uniform strain, the imbalance between sublattices leads to a stronger magnetization on one than the other at any finite temperature, so the ordered state is formally ferrimagnetic. Again, because of the relative smallness of $f_{b}$, the net moment is small.

The nature of the broken symmetries in these two-dimensional magnets is such that we do not truly expect spontaneous long-range magnetic order to set in at any finite temperature Mermin and Wagner (1966). The mean-field phase diagram instead indicates a cross-over temperature at which measurements of the spin-spin correlation length begins to exceed the typical distance between neighboring magnetic moments. In principle this effect could be observed by slowly lowering the temperature of the system to some value below $T_{c}$, and then rapidly quenching the system to very low temperature, so that thermal magnetic fluctuations are frozen out. In principle, a local measurement of the magnetization at this low temperature would reveal domains of size the correlation length at the temperature from which the system was quenched.

In the case of Landau levels generated by real fields, the presence of Zeeman coupling has interesting implications. The symmetry of the effective Hamiltonian for the spins is lowered from SU(2) to U(1), and the mean-field groundstate when ordered is a canted antiferromagnet. Even in the absence of long-range magnetic order, the system should still exhibit a true thermodynamic, Kosterlitz-Thouless (KT) phase transition, in which vortices of the in-plane component of the magnetization unbind at some temperature $T_{KT}$. To estimate this, we consider a very simple model, in which the U(1) spin degrees of freedom are on a square lattice, with nearest neighbor coupling constant $J_{eff}$. By computing the mean-field $T_{c}$ of this model, and matching this to the simplest estimate of the (KT) transition, $k_{B}T_{KT}=\pi J_{eff}/2$ Nelson (2002), one can obtain an estimate for $T_{KT},$ which turns out to be simply proportional to the mean-field $T_{c}$, as illustrated in Fig. 12. An important caveat with respect to this estimate is that it does not include any renormalization of $J_{eff}$ due to vortex-antivortex pairs, which always decreases the transition temperature. In particular this renormalization should become quite large as $B\rightarrow 0$, and SU(2) symmetry is restored, in which case there will be no KT transition, so that $T_{KT}\rightarrow 0$. The low barrier to spins tilting into the $\hat{z}$ direction in this limit suggests that the core energy of vortices becomes small, so that our simple estimate becomes increasingly unreliable. Fig. 12 schematically shows the form we expect $T_{KT}$ to take in a more sophisticated treatment.

We conclude with some speculations about further interesting behaviors this system might host. In the case of non-uniform strain, as noted above, the relatively large energy scale associated with binding of LLL electrons to sites leads to two temperature scales, a lower one at which the majority of spins lose most of their collective order, and higher one at which those binding electrons do so. While direct magnetization measurements are unlikely to detect this, an interesting possibility is that it might be visible in transport due to scattering of electrons from these spin degrees of freedom, or via na anomalous Hall effect. Beyond this, the behavior of the system with greater levels of doping may show interesting effects, for example yielding a mean-field $T_{c}$ that drops as the LLL is depleted. Further doping may yield similar physics to that of the LLL when the Fermi energy reaches negative Landau levels, which may yield oscillations in $T_{c}$. Finally, in our mean-field analysis, collective behavior among the degrees of freedom was not considered. For example, correlations among nearby impurity spins could lead to electrons binding to multiple impurities rather than individual ones. In principle our mean-field theory could accommodate this phenomenologically by adopting larger values of $f_{b}$; since $\ell$ is typically much larger than the average distance between impurities, this renormalization could be considerable. If such effects are important, we expect ferrimagnetism in a strain-induced field could be notably larger than our estimates above. Beyond this, interactions among the electrons themselves may induce a spin stiffness, which could lessen the correlation between electron and impurity spins in the system, and would tend to work against the induced ferrimagnetism. Thus the net effect of correlated behaviors in this system is unclear. We leave their consideration for future research.

Acknowledgements – This work was supported by the NSF through Grant Nos. DMR-1506263 and DMR-1506460, by the US-Israel Binational Science Foundation, and by Indiana University through an FRSP grant. HAF thanks the Aspen Center for Physics, where part of this work was done.