PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017)

Received: 19 September 2017, Accepted: 12 October 2017
Edited by: K. Hallberg
Licence: Creative Commons Attribution 3.0
DOI: http://dx.doi.org/10.4279/PIP.090008

www.papersinphysics.org

ISSN 1852-4249

Dilute antiferromagnetism in magnetically doped phosphorene

A. Allerdt,1 A. E. Feiguin1∗

We study the competition between Kondo physics and indirect exchange on monolayer black phos-
phorous using a realistic description of the band structure in combination with the density matrix
renormalization group (DMRG) method. The Hamiltonian is reduced to a one-dimensional problem
via an exact canonical transformation that makes it amenable to DMRG calculations, yielding exact
results that fully incorporate the many-body physics. We find that a perturbative description of the
problem is not appropriate and cannot account for the slow decay of the correlations and the complete
lack of ferromagnetism. In addition, at some particular distances, the impurities decouple forming
their own independent Kondo states. This can be predicted from the nodes of the Lindhard function.
Our results indicate a possible route toward realizing dilute anti-ferromagnetism in phosphorene.

I. Introduction

Phosphorene, a single layer of black phosphorus, is one
of the many allotropes of the element. Others include
red, white, and velvet phosphorus. Of these, black is
the most thermodynamically stable and least reactive
[1]. Its crystalline structure resembles that of graphite,
whose 2D form is famously known as graphene, and can
similarly be fabricated into 2D layers. The main qual-
itative distinction lies in the fact that phosphorene has
a “puckered” hexagonal structure which is responsible
for opening a gap in the band structure. In addition, the
bonding structure is vastly different, since phosphorene
does not have sp2 bonds, but is composed of 3p orbitals
[2, 3].

Successful production of monolayer black phospho-
rus has been achieved only in the last few years, with the
first publications concerning its single layer properties
appearing in 2014 [4–7]. While being a semi-conductor
with a small direct band gap in bulk form (∼ 0.3 eV)
[8], the gap increases as the number of layers is de-

∗a.feiguin@northeastern.edu

1 Department of Physics, Northeastern University, Boston, Mas-
sachusetts MA 02115, USA.

creased.

Since its appearance, there have been many proposed
promising applications and exotic properties. Besides
being a semi-conductor with a band gap in the optical
range, it has ample flexibility and high carrier mobil-
ity [4, 9]. Further, it also possesses interesting opti-
cal properties such as absorbing light only polarized in
the armchair direction, indicating a possible future as a
linear polarizer [10]. Its stable excitons present possi-
ble applications in optically driven quantum computing
[11]. Interest on its applications continues to grow. For
a comprehensive review, we refer to Ref. [11].

Our understanding of magnetic doping in phospho-
rene is still in its infancy. A thorough study of metal
adatoms adsorbed on phosphorene was performed in
Ref. [12] where a variety of structural, electronic and
magnetic properties emerge. The authors show the
binding energies are twice the amount in graphene. A
DFT study has proposed possible chemical doping by
means of adsorption of different atoms, ranging from
n-type to p-type, as well as transition metals with finite
magnetic moment [13]. Another conceivable method of
moving the Fermi level is to induce strain on the lat-
tice, causing the level to move into the conduction band
[14]. The case of magnetic impurities is considerably

090008-1



PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017) / A. Allerdt et al.

Figure 1: Phosphorene lattice showing the directions
chosen to place the impurities for calculations shown in
Figs. 3, 4, and 5. The numbering of lattice sites refers
to the impurity separations in the same figures. Also
shown are the five hopping parameters borrowed from
Ref. [5]. The labels A, B, C and D refer to the four
sublattices.

non-trivial, since one has to account for the Kondo ef-
fect [15] with the impurity being screened by the con-
duction spins in the metallic substrate, and the RKKY
interaction [16–18], an effective indirect exchange be-
tween impurities mediated by the conduction electrons.
These two phenomena are expected to be present and
compete in phosphorene when the Fermi level is not
sitting in the gap. The RKKY interaction in phospho-
rene with the inclusion of mechanical strain is investi-
gated in Ref. [14]. However, all their calculations are
based on second order perturbation theory and ignore
the effects that Kondo physics can induce. In this work,
we present numerical results for two Kondo impurities
in phosphorene that capture the full many-body physics
and we discuss the competition between Kondo and in-
direct exchange.

II. Model and numerical method

The Hamiltonian studied in this work is the two-
impurity Kondo model, generically written as:

H = Hband + JK

(
~S1 ·~sr1 + ~S2 ·~sr2

)
, (1)

where Hband is the non-interacting part describing the
band structure of the material while the impurities are
locally coupled to the substrate at positions r1 and r2
via an interaction JK. Here, ~Si represent the impurities,
and ~sri is the spin of the conduction electrons at site
ri. This problem has been theoretically studied in other

Figure 2: Local density of states of phosphorene. Verti-
cal lines show positions of the different Fermi energies
used in the calculations.

materials using a variety of approaches. This work,
however, will take a non-perturbative approach reveal-
ing important subtle points. Traditionally, the RKKY
interaction is characterized by the Lindhard function in
second order perturbation theory, which takes the form:

χij = 2<


 ∑
Eα>Ef>Eβ

〈i|α〉〈α|j〉〈j|β〉〈β|i〉
Eα −Eβ


 , (2)

where |i,j〉 represent lattice sites and |α,β〉 are eigen-
states of Hband. This function oscillates and changes
sign with a period that depends on the density of the
electrons. When the full geometry of the problem is
taken into account, the picture becomes more complex,
since a system could have a complex Fermi surface even
with more than one band crossing the Fermi level in
some cases [19]. However, it has been shown by the
authors in previous studies [19–21] that this does not
uncover the full picture.

We adapt the four-band model introduced in Ref. [5]
to obtain a tight-binding description of the phospho-
rene band structure. Five hopping parameters are em-
ployed to closely reproduce the bands near the Fermi
level with a gap of approximately 1.6 eV. These are
shown schematically in Fig. 1, and their values are re-
stated in Table 1. A plot of the density of states is shown
in Fig. 2. All energies in the calculations are in units of
eV. To see the band dispersions, we refer to Ref. [5].

In order to solve the interacting problem we first
perform a unitary transformation to map the quadratic
part of the Hamiltonian Hband onto an equivalent one-
dimensional one, as described in detail in Refs. [20, 22].

090008-2



PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017) / A. Allerdt et al.

t1 −1.220
t2 3.665
t3 −0.205
t4 −0.105
t5 −0.055

Table 1: Tight binding hopping parameters in eV used
in the calculations throughout this work.

The full many-body calculation can in turn be car-
ried out using the density matrix renormalization group
(DMRG) algorithm [23–25] with high accuracy and
without approximations. For all DMRG calculations,
the total system size is L = 124 (including impuri-
ties), and we fix the truncation error to be smaller than
10−7. It is known that phosphorene nano-ribbons can
host quasi-flat edge states whose emergence is topologi-
cally similar to graphene [26]. A study of their even/odd
properties under application of an electric field has re-
cently been performed [27], revealing a gap opening in
these edge states. Edge states can introduce notable
finite-size effects [21]. However, due to the geometry
of the lattices produced by the Lanczos transformation
employed in this work, these edge states are not present.
In addition, we shall focus our study in a regime away
from half-filling.

III. Results

We start by calculating the Lindhard function as a ref-
erence using Eq. (2). Figure 3 shows results for two
different Fermi energies, one in the conduction band,
and one in the valence band. Notice the different pe-
riod of oscillation, which is to be expected from conven-
tional RKKY theory [16–18]. We also show, in differ-
ent panels, the DMRG result for the impurity spin-spin
correlations 〈Sz1Sz2〉. Since the system is SU(2) sym-
metric, only the z-component is displayed. For small
values of JK we find that he RKKY correlations satu-
rate at -1/4. This is an indication that the two impuri-
ties are decoupled forming “free moments”. When this
happens, the spins can be pointing in either direction
yielding a four-fold degenerate ground state. Since we
force Sz = 0 in our calculations, the impurities have no
choice but to align anti-ferromagnetically. This is a fi-
nite size effect that is expected when the interaction JK
is of the order of the level spacing in the non-interacting
bands. Apart from short distances, the spin correla-
tions roughly follow the oscillation patterns predicted

Figure 3: (a) Lindhard function and (b) spin-spin corre-
lations along the zig-zag direction with the Fermi level
at position A in Fig. 2. (c) and (d) are the same quanti-
ties with the Fermi level at position C.

by the Lindhard function with one striking difference:
ferromagnetism is non-existent. This lack of ferromag-
netism is also seen in previous works [19–21]. How-
ever, in doped graphene, for example [21], the correla-
tions depart drastically from the Lindhard function, and
Kondo physics dominates after just a few lattice spac-
ings. Here, we find that RKKY physics plays the dom-
inant role, with correlations persisting even for quite
large coupling (JK=2.0 eV), which is greater than the
band gap.

Figure 4 again the Lindhard function and spin cor-
relations, but for a different Fermi energy and in an
extended range. This is to highlight the fact that the
RKKY oscillations survive at large distances with lit-
tle signs of decay. Ferromagnetism is again absent in
contrast with the predictions of perturbation theory, ex-
emplifying the need for exact numerical techniques that
can capture the many-body physics. The impurities,
however, do appear to form their independent Kondo
singlets at some particular distances. This results in the
impurities being completely uncorrelated 〈Sz1Sz2〉 = 0
and typically occurs at positions where the Lindhard
function has a node, as observed in Ref. [20].

To show that these effects are not due to the par-
ticular directions chosen, calculations were done along
other paths with qualitatively similar results. As an ex-
ample, we show in Fig. 5 the correlations along the
perpendicular direction with the Fermi level at posi-

090008-3



PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017) / A. Allerdt et al.

Figure 4: Lindhard function (top) and spin-spin cor-
relations (bottom) along the zig-zag direction with the
Fermi level at position B in Fig. 2. Results are shown
with an extended range, highlighting the persistent os-
cillations even for large JK.

tion B. As can be seen from the lattice structure, sin-
gle layer phosphorene contains four sublattices, as la-
beled in Fig. 1. Assuming the first impurity is placed
on an ‘A’ site, along the zig-zag direction the impurities
follow an “A − A,A − B,A − A,A − B...” pattern.
Alternatively, the second impurity could be placed on
the other “layer” along the zig-zag direction resulting
in a “A − C,A − D,A − C,A − D...” pattern. It is
clear that A−D is equivalent to B −C. The armchair
direction, as it has been defined, consists of impurities
always on the same sublattice. In all cases (not shown),
ferromagnetism is absent and we find dominant anti-
ferromagnetic RKKY interactions.

IV. Conclusions

By means of a canonical transformation and exact nu-
merical calculations using the DMRG method, we have
studied the competition between Kondo and RKKY
physics on phosphorene. The method is numerically ex-
act, and even though it is limited to finite systems, these
can be very large, of the order of a hundred lattice sites
and more. Our results highlight the non-perturbative
nature of the RKKY interaction and the non-trivial ab-

Figure 5: Lindhard function (top) and spin-spin corre-
lations (bottom) along the armchair direction with the
Fermi level at position B in Fig. 2.

sence of ferromagnetism. This remains an outstanding
question and should stimulate more research in this di-
rection. It is possible that by adding a repulsive inter-
action between conduction electrons, the system may
acquire a net magnetic moment, a behavior of this type
has been observed in graphene doped with hydrogen de-
fects [28–30]. According to Lieb’s theorem [31], this is
expected to occur for bi-partite lattices when two impu-
rities are on the same sublattice. Even though phospho-
rene has four sublattices (and hence, the theorem does
not rigorously apply), the system may still realize simi-
lar physics. Unfortunately, our approach can only de-
scribe non-interacting/quadratic Hamiltonians and we
cannot prove this conjecture. On the other hand, the
dominant anti-ferromagnetism at all dopings and dis-
tances (more robust than in graphene [21]) indicates
a route toward realizing dilute 2D anti-ferromagnetism
with phosphorene.

Acknowledgements - The authors are grateful to the
U.S. Department of Energy, Office of Basic Energy Sci-
ences, for support under grant DE-SC0014407.

[1] R B Jacobs, Phosphorus at high temperatures and
pressures, J. Chem. Phys. 59, 945 (1937).

090008-4



PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017) / A. Allerdt et al.

[2] L Pauling, M Simonetta, Bond orbitals and bond
energy in elementary phosphorus, J. Chem. Phys.
20, 29 (1952).

[3] R R Hart, M B Robin, N A Kuebler, 3p orbitals,
bent bonds, and the electronic spectrum of the p4
molecule, J. Chem. Phys. 42, 3631 (1965).

[4] H Liu et al., Phosphorene: An unexplored 2d
semiconductor with a high hole, ACS Nano 8,
4033 (2014).

[5] A N Rudenko, M I Katsnelson, Quasiparticle
band structure and tight-binding model for single-
and bilayer black phosphorus, Phys. Rev. B 89,
201408 (2014).

[6] A S Rodin, A Carvalho, A H Castro Neto, Strain-
induced gap modification in black phosphorus,
Phys. Rev. Lett. 112, 176801 (2014).

[7] L Li et al., Black phosphorus field-effect transis-
tors, Nat. Nanotechnol. 9, 372 (2014).

[8] R W Keyes, The electrical properties of black
phosphorus, Phys. Rev. 92, 580 (1953).

[9] X Ling, H Wang, S Huang, F Xia, M S Dres-
selhaus, The renaissance of black phosphorus, P.
Natl. Acad. Sci. USA 112, 4523 (2015).

[10] V Tran, R Soklaski, Y Liang, L Yang, Layer-
controlled band gap and anisotropic excitons in
few-layer black phosphorus, Phys. Rev. B 89,
235319 (2014).

[11] A Carvalho et al., Phosphorene: From theory to
applications, Nat. Rev. Mater. 1, 16061 (2016).

[12] V V Kulish, O I Malyi, C Persson, P Wu, Adsorp-
tion of metal adatoms on single-layer phospho-
rene, Phys. Chem. Chem. Phys. 17, 992 (2015).

[13] P Rastogi, S Kumar, S Bhowmick, A Agarwal, Y
S Chauhan, Effective doping of monolayer phos-
phorene by surface adsorption of atoms for elec-
tronic and spintronic applications, IETE J. Res.
63, 205 (2017).

[14] H-J Duan et al., Anisotropic RKKY interaction
and modulation with mechanical strain in phos-
phorene, New J. Phys. 19, 103010 (2017).

[15] A C Hewson, The Kondo Problem to Heavy
Fermions, Cambridge University Press, New York
(1983).

[16] K Yosida, Magnetic properties of Cu-Mn alloys,
Phys. Rev. 106, 893 (1957).

[17] M A Ruderman, C Kittel, Indirect exchange cou-
pling of nuclear magnetic moments by conduction
electrons, Phys. Rev. 96, 99 (1954).

[18] T Kasuya, A theory of metallic ferro- and anti-
ferromagnetism on Zener’s model, Prog. Theor.
Phys. 16, 45 (1956).

[19] A Allerdt, R Žitko, A E Feiguin, Nonperturba-
tive effects and indirect exchange interaction be-
tween quantum impurities on metallic (111) sur-
faces, Phys. Rev. B 95, 235416 (2017).

[20] A Allerdt, C A Büsser, G B Martins, A E Feiguin,
Kondo versus indirect exchange: Role of lattice
and actual range of RKKY interactions in real ma-
terials, Phys. Rev. B 91, 085101 (2015).

[21] A Allerdt, A E Feiguin, S D Sarma, Competi-
tion between Kondo effect and RKKY physics in
graphene magnetism, Phys. Rev. B 95, 104402
(2017).

[22] C A Büsser, G B Martins, A E Feiguin, Lanczos
transformation for quantum impurity problems in
d-dimensional lattices: Application to graphene
nanoribbons, Phys. Rev. B 88, 245113 (2013).

[23] S R White, Density matrix formulation for quan-
tum renormalization groups, Phys. Rev. Lett. 69,
2863 (1992).

[24] S R White, Density-matrix algorithms for quan-
tum renormalization groups, Phys. Rev. B 48,
10345 (1993).

[25] A E Feiguin, The density matrix renormalization
group method and its time-dependent variants,
AIP Conf. Proc. 1419, 5 (2011).

[26] M Ezawa, Topological origin of quasi-flat edge
band in phosphorene, New J. Phys. 16, 115004
(2014).

[27] B Zhou, B Zhou, X Zhou, G Zhou, Even–odd
effect on the edge states for zigzag phosphorene
nanoribbons under a perpendicular electric field,
J. Phys. D: Appl. Phys. 50, 045106 (2017).

090008-5



PAPERS IN PHYSICS, VOL. 9, ART. 090008 (2017) / A. Allerdt et al.

[28] O V Zazyev, L Helm, Defect-induced magnetism
in graphene, Phys. Rev. B 75, 125408 (2007).

[29] S Casolo, O M Løvvik, R Martinazzo, G F Tan-
tardini, Understanding adsorption of hydrogen
atoms on graphene, J. Chem. Phys. 130, 054704
(2009).

[30] H González-Herrero et al., Atomic-scale control
of graphene magnetism by using hydrogen atoms,
Science 352, 437 (2016).

[31] E H Lieb, Two theorems on the Hubbard model,
Phys. Rev. Lett. 62, 1201 (1989).

090008-6