# Single-Molecule Magnet Mn12 on GaAs-supported Graphene:
Gate Field Effects from First Principles

Shuanglong Liu,1,2,3 Maher Yazback,1,2,3 James N. Fry,1
Xiao-Guang Zhang,1,2,3 Hai-Ping Cheng1,2,3<EMAIL_ADDRESS>1Department of
Physics, University of Florida, Gainesville, Florida 32611, USA
2Quantum Theory Project, University of Florida, Gainesville, Florida 32611,
USA
3Center for Molecular Magnetic Quantum Materials, University of Florida,
Gainesville, Florida 32611, USA

###### Abstract

We study gate field effects on the heterostructure Mn12O12(COOH)16(H2O)4 $|$
graphene $|$ GaAs via first-principles calculations. We find that under
moderate doping levels electrons can be added to but not taken from the
single-molecule magnet Mn12O12(COOH)16(H2O)4 (Mn12). The magnetic anisotropy
energy (MAE) of Mn12 decreases as the electron doping level increases, due to
electron transfer from graphene to Mn12 and change in the band alignment
between Mn12 and graphene. At an electron doping level of $-5.00\times
10^{13}\,\textrm{cm}{}^{-2}$, the MAE decreases by about 18% compared with
zero doping. The band alignment between graphene and GaAs is more sensitive to
electron doping than to hole doping, since the valence band of GaAs is close
to the Fermi level. The GaAs substrate induces a small band gap in the
supported graphene under zero gate field and a nearly strain-free
configuration. Finally, we propose a vertical tunnel junction for probing the
gate dependence of MAE via electron transport measurements.

##  I Introduction

Mn12O12(COOR)16(H2O)4, where R represents -CH3 or other ligands, is a
prototypical single-molecule magnet (SMM) [1] whose magnetic and electronic
properties have been studied since the late 90s [2, 3, 4, 5, 6]. This molecule
is also interesting as a spin system because its total magnetic moment puts it
near the boundary between classical and quantum regimes. Tunneling magnetism
measurements show its quantum nature [7], but its big magnetic moment ($S=10$)
makes it almost classical. The magnetic anisotropy energy (MAE) of
Mn12O12(COOR)16(H2O)4, which forms the barrier for magnetic tunneling, depends
on the type of ligand [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32] as well as the charge state of the
molecule. In experiments, one or two electrons reduced Mn12O12(COOR)16(H2O)4
have been synthesized by adding PPh4^+ or other cations to the molecular
crystal [33, 34, 35, 36, 13, 37]. The resulting negatively charged
Mn12O12(COOR)16(H2O)4 with integer number of electrons yields a decreased MAE
for -CHCl2 [35] or -C6F5 [13] ligands. Based on computations, a decreased MAE
has also been reported for negatively charged Mn12O12(COOR)16(H2O)4 with -H
[6] or -C6H5 [38] ligands. When a magnetic molecule is adsorbed on a surface,
a fractional number of electrons may be transferred to/from the molecule,
modifying its MAE [39, 40, 41, 42].

Most earlier experiments are performed either in solution or in molecular
crystals. Recently, Hebard’s [43] group experimentally investigated the gate-
voltage dependence of transport properties of Mn12O12(COOR)16(H2O)4 on
graphene surfaces (supported by Si) with different ligands. The applied gate
voltage is believed to affect transport properties by modulating the charge
transfer between graphene and Mn12O12(COOR)16(H2O)4. However, such experiments
were unable to determine whether the MAE of Mn12O12(COOR)16(H2O)4 is also
tuned by the gate voltage and the charge transfer induced by it. It is also
not clear how the coupling to a semiconductor substrate changes the magnetic
and electronic properties of the molecule. In this study, we aim to answer
these questions through a first-principles calculation using
Mn12O12(COOH)16(H2O)4 (Mn12) as an example.

We choose GaAs-supported graphene ($\textrm{Gr}|\textrm{GaAs}$) as substrate.
Compared to a Si substrate, GaAs provides free carriers at lower temperatures
[44, 45]. This allows transport measurements at lower temperatures, where the
magnetic properties of SMMs are better observed. Common stable GaAs surfaces
include (100) [46, 47] and (111) [48, 49] with various surface
reconstructions. The GaAs (111) surface has closest lattice match with the
graphene lattice. Munshi et al. reported the growth of GaAs nanorods on few-
layer graphene, where a GaAs(111) surface is in contact with the top graphene
layer [50]. Several epitaxial atomic structures of the interface between
GaAs(111) and graphene have been proposed in literature but all with
relatively large strain in graphene [50, 51]. Here, we propose a different
atomic structure where the strain in graphene is close to zero. Using this
structure, we simulate gate field effects on the $\textrm{Gr}|\textrm{GaAs}$
interface with and without the adsorption of Mn12. For brevity, we denote the
heterostructure with Mn12 by $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$.

The rest of the paper is organized as follows. We describe the computational
details in Section II. We present the atomic structure, the electronic
structure, and the magnetic anisotropy of the
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ heterostructure in Sections III.1–III.3.
In Section III.4, we propose a vertical tunnel junction based on the
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ heterostructure. Finally, we conclude in
Section IV.

##  II Method

All calculations are based on density functional theory (DFT) [52, 53] as
implemented in the Vienna Ab initio Simulation Package (VASP) [54, 55] and the
SIESTA package [56]. We use VASP to relax atomic structures with no applied
electric field and SIESTA to calculate electronic structures in the presence
of a gate electric field. The atomic structure is kept at the relaxed
configuration with no electric field.

In VASP calculations, we apply an energy cutoff of $500\,\textrm{eV}$ for
plane waves and projector augmented wave (PAW) pseudopotentials [57]. We adopt
the exchange correlation energy functional proposed by Klimes et al., optB86b
[58], to include the van der Waals interaction between graphene and Mn12 (or
GaAs). Given the large supercell size ($a=b\approx 3.2\,\textrm{nm}$,
$c=5.2\,\textrm{nm}$), only the $\Gamma$ point is sampled in reciprocal space.
The energy tolerance for electronic self-consistency and the force tolerance
for ionic relaxation are set to $1\times 10^{-6}\,\textrm{eV}$ and
$0.02\,\textrm{eV}/\textrm{\AA}$ respectively. We use Gaussian smearing with a
smearing parameter of $0.2\,\textrm{eV}$ to facilitate the electronic self-
consistent procedure. In order to eliminate interaction between periodic
images in the perpendicular direction (with respect to the GaAs slab), a
vacuum layer of at least $18\,\textrm{\AA}$ is added and electric dipole
corrections (for both energy and force) are enabled. On-site Coulomb
interaction ($U$) within the DFT+$U$ method and spin-orbit interactions are
not considered for ionic relaxations.

In SIESTA calculations, we apply double-$\zeta$ polarized (DZP) basis
functions [56] for Mn and O atoms and single-$\zeta$ polarized (SZP) basis
functions for C, H, Ga, and As atoms. Such a mixed basis set allows us to
describe the magnetic properties (due to Mn atoms) accurately with less
computational load. A SZP basis set is often not sufficiently accurate for
structural relaxation, but it usually produces reasonably good electronic
structure for a fixed atomic structure. The basis functions are optimized for
Mn12, graphene, and the GaAs slab separately. Detailed specifications of the
basis functions are presented in Appendix A. Since the optB86b functional is
not available in SIESTA and the atomic structure is fixed, we use the Perdew-
Burke-Ernzerhof (PBE) exchange correlation energy functional [59] instead. We
apply norm-conserving pseudopotentials as generated by the Troullier-Martins
scheme [60] and a mesh cutoff of $200\,\textrm{Ry}$ for real space sampling.
To accurately determine the Fermi energy, we adopt a $6\times 6$ $k$-grid for
sampling reciprocal space [61] and the $4^{\textrm{th}}$ order Methfessel-
Paxton smearing method [62] with a smearing temperature of $200\,\textrm{K}$.
Results of convergence tests with respect to the smearing temperature and with
respect to the $k$-grid are given in Tables B1 and B2 respectively (see
Appendix B). The effects of a single back gate are modeled via the effective
screening medium (ESM) method [63]. In the single gate configuration, the
boundary condition for the back-gate (vacuum) side is constant electrostatic
potential (vanishing first derivative of the potential). Such a non-periodic
boundary condition for the electrostatic potential is imposed by the
corresponding Green’s function in the ESM method. In our simulations, the
back-gate (vacuum) boundary is $15\,\textrm{\AA}$ below (above) the bottom
(the top) of the system under study. In order to improve the numerical results
for band alignment within the heterostructure
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$, we treat the semicore Ga $d$ electrons
as valence electrons and adopt the DFT+1/2 approach for GaAs [64, 65]. A
DFT+1/2 cutoff radius of $3.8\,\textrm{Bohr}$ is applied to limit the range of
self-energy potential for As $4p$ orbitals [66]. The DFT+1/2 band gap of bulk
GaAs is calculated to be $1.526\,\textrm{eV}$, in good agreement with the
experimental value of $1.52\,\textrm{eV}$ at low temperatures [44], as opposed
to $0.236\,\textrm{eV}$ without this approach. For the same purpose of
improving the band alignment, we apply the DFT+$U$ method proposed by Dudarev
et al. [67] and set the on-site Coulomb interaction $U$ to $4\,\textrm{eV}$
for Mn atoms. This value for $U$ yields good agreement of the density of
states of Mn12 compared with X-ray photoemission spectra (XPS) and X-ray
emission spectra (XES) measurements [68].

We also use the SIESTA package to calculate the MAE of Mn12 on graphene
without the GaAs substrate (to be justified later). Spin-orbit interactions
are included via the pseudopotentials [69, 70] and evaluated in an on-site
approximation [71]. At high electron doping levels, the self-consistent
calculation fails to converge with the Methfessel-Paxton smearing method. This
convergence problem is solved by using the Fermi-Dirac smearing method. To
improve the numerical accuracy, we apply a $24\times 24$ $k$-grid together
with a smearing parameter of $0.041\,\textrm{eV}$ and an energy tolerance of
$1\times 10^{-6}\,\textrm{eV}$ for electronic self-consistency. We set the
spin-orbit coupling strength parameter to be 1.34, such that the calculated
MAE ($5.2\,\textrm{meV}$, or $61\,\textrm{K}$) of Mn12O12(COOCH3)16(H2O)4 is
close to the experimental value [72]. The DFT+$U$ method is not applied for
calculating the MAE.

##  III Results

###  III.1 Atomic structure

Fig. 1 shows the atomic structure of the heterostructure
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$. The GaAs(111) surface is modeled by a
slab consisting of six Ga atomic layers and six As atomic layers. The top two
atomic layers are stabilized by a $2\times 2$ reconstruction with Ga vacancies
[49]. Each As atom at the bottom is terminated by a pseudo hydrogen atom with
0.75 electrons to avoid fictitious surface bands. The lattice constants of
bulk GaAs and graphene are 5.653 and $2.461{\,\textrm{\AA}}$, respectively. We
found a good lattice match for the supercell shown in Fig. 1, which contains
$4\times 4$ GaAs unit cells (with surface reconstruction) and $13\times 13$
graphene unit cells. This allows graphene to match the GaAs lattice with only
a 0.04% compression.

Figure 1:  Side view (left) and top view (right) of the
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ heterostructure. In the side view,
$d_{1}\approx 2.46\,\textrm{\AA}$ and $d_{2}\approx 3.42\,\textrm{\AA}$. In
the top view, only graphene and the atoms adjacent to graphene are shown. The
hydrogen atoms of Mn12 are magnified for visibility. Purple: manganese, red:
oxygen, gray: carbon, white: hydrogen, green: gallium, and blue: arsenic.

During atomic relaxation, all atoms are relaxed except the bottom three Ga
(As) atomic layers in order to mimic the bulk environment of GaAs. Without the
Mn12 molecule, graphene already buckles slightly due to the nonuniform
interaction with GaAs. The maximal out-of-plane displacement of graphene
(carbon atoms) is about $\pm 0.08\,\textrm{\AA}$, which is one order of
magnitude larger than in-plane displacements. After the adsorption of Mn12,
graphene is further distorted, with a maximal out-of-plane displacement of
about $\pm 0.03\,\textrm{\AA}{}$. Due to relatively weak chemical bonds
between the core of the Mn12 molecule and the surrounding ligands, the
molecule is prone to distortion and losing parts of ligands. For example, Mn12
loses its structural integrity when deposited on a Au(111) surface [73]. In
contrast, previous DFT calculations suggest that Mn12 remains intact on
graphene [4]. In the current study, the structure of Mn12 also remains intact
with slight structural distortion when it is adsorbed on
$\textrm{Gr}|\textrm{GaAs}$. Figs. 2a and 2b show histograms of the number of
chemical bond lengths and bond angles versus change in the bond length or bond
angle. All bond lengths or bond angles of Mn12 change by no more than
$0.03\,\textrm{\AA}$ or $4^{\circ}$. A fraction 69.5% of all chemical bond
lengths and 60.4% of bond angles of Mn12 change within the range of
$[-0.005:0.005]\,\textrm{\AA}{}$ and $[-0.25:0.25]\,\deg$. However, there are
8 Mn-O bonds which change by more than $0.015\,\textrm{\AA}$ in length. Each
of these Mn-O bonds stems from a Mn^3+ ion of the outer Mn8O8 ring, and all of
them point towards graphene (from Mn to O). Also, there are two Mn-O-H bond
angles that change by more than $3.5\,\deg$. For both of these Mn-O-H
entities, the O-H bond belongs to a H_2O unit and the Mn-O bond length changes
by more than $0.015\,\textrm{\AA}$.

Figure 2:  Statistics of (a) bond length changes and (b) bond angle changes of
Mn12. The comparison is between the Mn12 molecule adsorbed on
$\textrm{Gr}|\textrm{GaAs}$ and an isolated free Mn12 molecule.

###  III.2 Electronic structure

Now, we turn to the electronic structure of the
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ heterostructure under both zero and
finite gate electric fields. Without spin-orbit coupling, isolated graphene is
a semi-metal, without an energy gap at the Fermi level. When graphene is
supported on the GaAs(111) surface, calculation using SIESTA shows an energy
gap of $2.2\,\textrm{meV}$ at the Dirac point. There are two possible factors
that can induce such an energy gap in graphene, 1) a structural distortion in
graphene itself, and 2) the non-uniform potential due to the GaAs substrate.
To identify which factor is responsible for the gap, we remove the GaAs
substrate and compute the band structure for isolated graphene but with the
same structural distortion. The resulting band structure is gapless to within
our numerical precision ($\sim 0.1\,\textrm{meV}$). From this we conclude that
the $2.2\,\textrm{meV}$ energy gap is due to the non-uniform potential of the
GaAs substrate. The band gap of graphene does not change after adsorption of
Mn12.

Figure 3:  Band structure of $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$. Squares,
circles, and triangles represent, respectively, Mn12, graphene, and GaAs
states. The size of a symbol (square, circle, triangle) is proportional to the
projected density of states of the Mn12 molecule, graphene layer, or 1st GaAs
layer). The Fermi level is set to zero.

In the $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ heterostructure, both graphene
and GaAs are nonmagnetic, and according to our spin-polarized DFT
calculations, both the highest occupied molecular orbital (HOMO) and the
lowest unoccupied molecular orbital (LUMO) of Mn12 have the same spin. We
define this spin to be spin up. Fig. 3 shows the spin-up energy bands of this
heterostructure. The Fermi level lies within the $2.2\,\textrm{meV}$ energy
gap of graphene. Relative to the Fermi energy, the HOMO orbital of Mn12, the
LUMO orbital of Mn12, the valence band maximum (VBM) of GaAs and the
conduction band minimum (CBM) of GaAs lie at
$E_{\textrm{HOMO}}=-0.681\,\textrm{eV}$,
$E_{\textrm{LUMO}}=0.192\,\textrm{eV}$,
$E_{\textrm{VBM}}=-0.035\,\textrm{eV}$, and
$E_{\textrm{CBM}}=1.496\,\textrm{eV}$. Compared with these four typical
energies, the band gap of graphene is much smaller. For this reason, and for
convenience of discussion, we will denote the Dirac point as the middle of the
apexes of the upper and lower Dirac cones, even though the two cones are not
quite connected. The LUMO energy $E_{\textrm{LUMO}}$, the energy of the Dirac
point $E_{\textrm{Dirac}}$, and the valence band maximum $E_{\textrm{VBM}}$
together dictate the band alignment between Mn12, graphene, and GaAs.

Figure 4:  The typical energies $E_{\textrm{LUMO}}$, $E_{\textrm{Dirac}}$, and
$E_{\textrm{VBM}}$ of $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ versus charge
doping level.

Fig. 4 shows how $E_{\textrm{LUMO}}$, $E_{\textrm{VBM}}$, and
$E_{\textrm{Dirac}}$ are affected by carrier density, i.e., the dependence of
the band alignment on electrostatic doping. A negative (positive) carrier
density means electron (hole) doping. All energies are measured relative to
the Fermi energy, which is set at zero. Overall, all these three typical
energies decrease with electron doping and increase with hole doping level.
$E_{\textrm{Dirac}}$ in particular is more sensitive to electron doping than
hole doping. For example, $E_{\textrm{Dirac}}$ is $-0.195\,\textrm{eV}$ at a
charge density of $\rho=-0.50\times 10^{13}\,\textrm{cm}^{-2}$, but only
$+0.008\,\textrm{eV}$ at $+0.50\times 10^{13}\,\textrm{cm}^{-2}$. The
asymmetric response of $E_{\textrm{Dirac}}$ to charge doping can be understood
from two aspects. First, electrons are mainly added to (taken from) graphene
before the HOMO of Mn12 (the VBM of GaAs) is brought to the Fermi level,
whence Mn12 (GaAs) becomes charged (doped). However, the VBM of GaAs is much
closer to the Fermi level than the HOMO of Mn12. Second, GaAs has a much
higher density of states than graphene. Therefore, the same amount of charge
causes a smaller shift in the energy bands when GaAs is doped than when
graphene is doped. At $\rho=-1.00\times 10^{13}\,\textrm{cm}^{-2}$, the LUMO
orbital of Mn12 is only a few meV above the Fermi level. According to our
calculations, the LUMO of Mn12 becomes partially occupied, and thus Mn12 is
negatively charged. It is noteworthy that the Dirac point of graphene becomes
lower in energy than the VBM of GaAs at a point between $\rho=0.50\times
10^{13}\,\textrm{cm}^{-2}$ and $\rho=1.00\times 10^{13}\,\textrm{cm}^{-2}$. We
also simulate gate field effects on the band alignment of
$\textrm{Gr}|\textrm{GaAs}$ without Mn12. By comparing with
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$, we find that the influence of Mn12 on
both $E_{\textrm{Dirac}}$ and $E_{\textrm{VBM}}$ is smaller than
$1\,\textrm{meV}$ within the charge density range of $[-0.50:1.00]\times
10^{13}\,\textrm{cm}^{-2}$. At $\rho=-1.00\times 10^{13}\,\textrm{cm}^{-2}$,
$E_{\textrm{Dirac}}$ ($E_{\textrm{VBM}}$) increases by 3.3 (2.0) meV upon the
adsorption of Mn12.

Table 1:  Mulliken charge analysis for $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ at different charge doping levels. Each entry represents the amount of excess charge per atom. Positive (negative) values indicate gain (loss) of electrons relative to the case of zero doping. The number of atoms for Mn12, graphene, the 1st Ga layer, and the 1st As layer are 100, 338, 48, and 64 respectively.  Doping level ($\times 10^{13}\;\textrm{cm}^{-2}$) |  Mn12 $(\times 10^{-3})$ |  Graphene $(\times 10^{-3})$ |  1st Ga layer $(\times 10^{-3})$ |  1st As layer $(\times 10^{-3})$ |  2nd Ga layer $(\times 10^{-3})$ |  2nd As layer $(\times 10^{-3})$ |  3rd Ga layer $(\times 10^{-3})$ |  3rd As layer $(\times 10^{-3})$
---|---|---|---|---|---|---|---|---
$-1.00$ | 0.61 | 1.67 | $-6.36$ | $-3.38$ | 0.26 | 0.03 | 0.17 | $-0.14$
$-0.50$ | 0.00 | 0.91 | $-3.15$ | $-1.66$ | 0.13 | 0.02 | 0.08 | $-0.06$
$-0.10$ | 0.00 | 0.17 | $-0.59$ | $-0.31$ | 0.03 | 0.01 | 0.01 | $-0.01$
$-0.05$ | 0.00 | 0.09 | $-0.30$ | $-0.15$ | 0.01 | 0.00 | 0.01 | $-0.01$
0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00
0.05 | 0.00 | $-0.07$ | 0.25 | 0.10 | 0.00 | $-0.01$ | $-0.01$ | 0.00
0.10 | 0.00 | $-0.18$ | 0.62 | 0.30 | $-0.02$ | $-0.01$ | $-0.01$ | 0.01
0.50 | 0.00 | $-0.19$ | 0.56 | 0.13 | $-0.13$ | $-0.20$ | $-0.06$ | -0.07
1.00 | 0.00 | $-0.17$ | 0.41 | $-0.03$ | $-0.21$ | -0.36 | $-0.12$ | -0.20

Table 1 shows the distribution of excess charges at the interface of
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ based on Mulliken charge analysis. The
electron distribution for $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ under zero
doping level is subtracted as a reference. As seen from the table, the number
of electrons on Mn12 barely changes except at $\rho=-1.00\times
10^{13}\,\textrm{cm}^{-2}$, where 0.06 ($\approx 0.61\times 10^{-3}\times
100$) electrons are added to Mn12. Graphene gains 0.56 ($\approx 1.67\times
10^{-3}\times 338$) at a doping level of $\rho=-1.00\times
10^{13}\,\textrm{cm}^{-2}$, and it loses 0.06 ($\approx\lvert-0.17\rvert\times
10^{-3}\times 338$) electrons at a doping level of $+1.00\times
10^{13}\,\textrm{cm}^{-2}$. This is consistent with our previous observation
that electron doping results in a larger shift of the Dirac point than hole
doping. The 1st GaAs layer tends to lose (gain) electrons when graphene gains
(loses) electrons. This is likely due to the Coulomb repulsion between these
two adjacent atomic layers. The number of electrons gained by the 1st GaAs
layer does not change monotonically with the hole doping level, but reaches a
maximum at $\rho=+0.10\times 10^{13}\,\textrm{cm}^{-2}$. This observation
supports that GaAs is doped with holes at a doping level between
$\rho=+0.10\times 10^{13}\,\textrm{cm}^{-2}$ and $\rho=+0.50\times
10^{13}\,\textrm{cm}^{-2}$.

Figure 5:  Mulliken charge analysis for Mn12 within
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ at a doping level of $\rho=-1.00\times
10^{13}\,\textrm{cm}^{-2}$. The charge distribution at zero doping is
subtracted as a reference. $\textrm{O}_{\textrm{Mn}^{\textrm{3+}}}$
($\textrm{O}_{\textrm{Mn}^{\textrm{4+}}}$) denotes oxygen atoms in the Mn8O8
ring (Mn4O4 center).

Compared with zero doping, the Mn8O8 ring, the Mn4O4 core, the COOH groups,
and the H2O units of Mn12 gain $8.18\times 10^{-3}$, $1.84\times 10^{-3}$,
$47.70\times 10^{-3}$, and $3.59\times 10^{-3}$ electrons respectively at the
charge doping level of $\rho=-1.00\times 10^{13}\,\textrm{cm}^{-2}$. Although
the COOH groups gain much more charge than the Mn8O8 ring in total, one Mn^3+
ion (in the Mn8O8 ring) gains more electrons than one atom in the COOH groups
on average (see Fig. 5). The average amount of electrons gained by one Mn^3+
ion is about four times larger than that gained by one Mn^4+ ion (in the Mn4O4
core). For both Mn^3+ and Mn^4+ ions, the spin-up channel gains and the spin-
down channel loses electrons. Since the highest occupied spin down molecular
orbital is well below the Fermi level and thus remains fully occupied, the
loss of electrons in the spin down channel is due to the deformation of
occupied states. According to Boukhvalov et al. [74], redistribution of
electron density affects the magnetic exchange interactions between Mn atoms.

Figure 6:  Hartree potential difference $\delta
V_{\textrm{H}}=V_{\textrm{H}}(\rho,z)-V_{\textrm{H}}(0,z)$ of
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$. The horizontal arrows indicate the
three regions indexed by 1, 2, and 3. H_Mn12 denotes a -H ligand of Mn12. The
inset zooms in on the region shaded in gray.

Let $V_{\textrm{H}}(\rho,z)$ be the plane-averaged (over the $x$-$y$ plane)
Hartree potential of $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ at position $z$ for
charge doping density $\rho$. Fig. 6 shows $\delta
V_{\textrm{H}}=V_{\textrm{H}}(\rho,z)-V_{\textrm{H}}(0,z)$, which can be
viewed as the gate potential. We focus on three regions: 1) the region between
GaAs and graphene, 2) the region between graphene and Mn12, and 3) the region
beyond Mn12. Let regions $i$ ($i=1,2,$ and $3$) begin at $z_{1}^{i}$ and end
at $z_{2}^{i}$. We define the gate potential buildup across region $i$ to be
$\Delta_{i}=\delta V_{\textrm{H}}(z_{2}^{i})-\delta
V_{\textrm{H}}(z_{1}^{i})$, where the argument $\rho$ is dropped for brevity.
Within the charge density range $[-0.50:1.00]\times
10^{13}\,\textrm{cm}^{-2}$, $\Delta_{1}$ is effectively tuned, whereas
$\Delta_{2}$ is not. This is consistent with Fig. 4, in the sense that the
band alignment between graphene and GaAs (Mn12) is (is not) effectively tuned.
Within the same charge density range, $\Delta_{3}$ remains almost zero, which
signifies that the gate electric field does not extend across Mn12. At
$\rho=-1.00\times 10^{13}\,\textrm{cm}^{-2}$, both $\Delta_{2}$ and
$\Delta_{3}$ are significantly tuned as a result of Mn12 being doped with
electrons. The gate electric field $\mathcal{E}_{i}$ in region $i$ is
determined by $e\mathcal{E}_{i}=\Delta_{i}/(z_{2}^{i}-z_{1}^{i})$ where $e$ is
the unit charge (a positive value). $\mathcal{E}_{1}$, $\mathcal{E}_{2}$, and
$\mathcal{E}_{3}$ are 0.69, 0.08, and $0.05\,\textrm{V}/\textrm{nm}$
respectively at $\rho=-1.00\times 10^{13}\,\textrm{cm}^{-2}$. Recall that GaAs
becomes doped with holes at a doping level between $\rho=0.10\times
10^{13}\,\textrm{cm}^{-2}$ and $\rho=0.50\times 10^{13}\,\textrm{cm}^{-2}$. As
consequence, $\mathcal{E}_{1}=-0.06\,\textrm{V}/\textrm{nm}{}$ at
$\rho=0.10\times 10^{13}\,\textrm{cm}^{-2}$ is larger in magnitude than
$\mathcal{E}_{1}=-0.05\,\textrm{V}/\textrm{nm}{}$ at $\rho=0.50\times
10^{13}\,\textrm{cm}^{-2}$, although the doping level for the latter is
higher. Note that the adsorption of Mn12 does not affect the Hartree potential
between GaAs and graphene by much (less than $1\,\textrm{meV}$) within the
charge density range $[-1.00:1.00]\times 10^{13}\,\textrm{cm}^{-2}$. This
implies that the electron tunneling rate between GaAs and graphene in
transport measurements is similar with and without the presence of Mn12.

###  III.3 Magnetic anisotropy

Next, we examine how electrostatic doping affects the magnetic anisotropy
energy of Mn12. We have shown in Section III.2 that electrons can be added to
but cannot be taken from Mn12 at moderate doping levels. Therefore, we
consider the effects of electron doping on MAE in particular. When the charge
doping density is higher in magnitude than $-1.00\times
10^{13}\,\textrm{cm}{}^{-2}$, there are only graphene and Mn12 states near the
Fermi level. According to second order perturbation theory [75], the MAE is
dominated by pairs of occupied and unoccupied states around the Fermi level.
Therefore, it should be a good approximation to calculate MAE without GaAs.
Based on the heterostructure $\ce{Mn12}|\textrm{Gr}$ without GaAs, we
calculate the MAE as $\mathrm{MAE}=E_{\perp}-E_{\parallel}$, where $E_{\perp}$
($E_{\parallel}$) is the DFT total energy for the spin of Mn12 perpendicular
(parallel) to the magnetic easy axis. The calculated MAE of
$\ce{Mn12}|\textrm{Gr}$ versus the charge doping density is shown in Fig. 7.
As seen from the figure, the MAE decreases as the electron doping level
increases. The MAE at $\rho=-5.0\times 10^{13}\,\textrm{cm}{}^{-2}$ decreases
by about 18% compared with the value at zero doping. It is noteworthy that the
MAE of $\ce{Mn12}|\textrm{Gr}$ at zero doping is 1% smaller than the MAE of an
isolated Mn12 [4].

Figure 7:  Magnetic anisotropy energy (MAE) versus charge doping density for
$\ce{Mn12}|\textrm{Gr}$.

We can relate the decrease in MAE in terms of electron transfer from graphene
to Mn12. Fig. 8 shows the number of electrons transferred from graphene to
Mn12 versus the charge doping density. The number of electrons added to Mn12
increases with the electron doping density, and about 0.18 electrons are added
to Mn12 at $\rho=-5.0\times 10^{13}\,\textrm{cm}{}^{-2}$. Along with the
electron transfer, the magnetic moment of Mn12, in units of Bohr magneton,
increases by nearly the same amount as the number of electrons transferred.
The magnetic moment of Mn12 increases rather than decreases since the spin
polarization of the LUMO orbital, which receives the added electrons, is
parallel to the total spin of the whole Mn12 molecule. Previously, Park and
Pederson considered potassium addition to introduce extra electrons on Mn12
[6]. They found that the MAE of Mn12 decreases by 15%, 37%, and 56% with 1, 2,
and 4 extra electrons respectively. As a numerical experiment, we add extra
electrons to a single isolated Mn12 molecule which has the same atomic
positions as the Mn12 molecule of the $\ce{Mn12}|\textrm{Gr}$ heterostructure.
The MAE of the isolated Mn12 molecule decreases as the number of added
electrons increases from zero to one. The MAE of isolated Mn12 decreases by
15% with one additional electron, coincident with Park and Pederson’s
findings. This confirms that electron transfer from graphene to Mn12 is
responsible for the decrease in the MAE of $\ce{Mn12}|\textrm{Gr}$ under
electron doping.

However, with 0.18 additional electrons, the MAE of isolated Mn12 decreases by
2%, which is much smaller than the decrease of 18% in the MAE of
$\ce{Mn12}|\textrm{Gr}$. This can be understood within second order
perturbation theory, where pairs of occupied and unoccupied states can be
classified into three types, Mn12-Mn12, graphene-graphene, and Mn12-graphene.
The 2% decrease in the MAE of the isolated Mn12 can be understood as a
decrease in the contribution from the Mn12-Mn12 pairs. Since contributions to
the MAE from graphene-graphene pairs are negligibly small (as seen from
isolated graphene), we owe the extra reduction in the MAE of
$\ce{Mn12}|\textrm{Gr}$ to a decrease in the contribution from the
Mn12-graphene pairs. Such an argument is supported by the observation that
electron doping dramatically changes the band alignment between Mn12 and
graphene and thus the pairs of occupied graphene (Mn12) states and unoccupied
Mn12 (graphene) states around the Fermi level. Note that the change in the
band alignment between Mn12 and graphene is mainly caused by electrostatic
doping rather than by electron transfer between Mn12 and graphene. Thus, the
change in the band alignment between Mn12 and graphene is another factor that
is responsible for the decrease in the MAE of $\ce{Mn12}|\textrm{Gr}$ under
electron doping.

In order to investigate the possible role of a local electric field between
Mn12 and graphene, we take the isolated Mn12 molecule again and move it close
to ($2.46\,\textrm{\AA}{}$ away from) the plane where the boundary condition
of constant electrostatic potential is enforced. In this way, we retain a
strong local electric field below Mn12 when the Mn12 molecule is charged.
Otherwise, the local electric field is not as strong in our earlier numerical
experiment using the isolated Mn12 molecule. It turns out that the decrease in
MAE is still 2% with 0.18 additional electrons even though a strong local
electric field is present below Mn12. Therefore, we exclude the local electric
field between Mn12 and graphene and the associated electric field across Mn12
as a separate reason for the decrease in the MAE of $\ce{Mn12}|\textrm{Gr}$
under electron doping. From perturbation theory [76], magnetic anisotropy of
molecules is closely related to a molecular quadrupole moment, which should
change with gate electric field. Since we determine molecular orbitals and
thus electron density self-consistently, the change in molecular quadrupole
moment and its influence on MAE have been captured implicitly by our
calculations. However, the effect of electric field on MAE via spin-orbit
coupling amplitude is still unknown.

Figure 8:  Number of electrons transferred from graphene to Mn12 (left-hand
scale) and magnetic moment of Mn12 (right-hand scale) versus electrostatic
doping level. A negative charge density means electron doping.

###  III.4 A vertical tunnel junction

In this section, we propose a vertical tunnel junction for detecting the gate-
dependence of MAE discussed in the previous section. The junction is
illustrated in Fig. 9a where a Mn12 molecule is sandwiched between two
graphene layers. The bottom graphene layer should be supported by $n$-type
GaAs, which has conducting charge carriers at low temperature (several
Kelvins) [44, 45]. A gate voltage $V_{\textrm{g}}$is applied between the
bottom graphene layer and the GaAs and a bias voltage $V_{\textrm{b}}$ between
the bottom and top graphene layers. In the proposed vertical tunnel junction,
GaAs is electron-doped to bring the conduction band near the Fermi level. Such
a band alignment may change the behavior of the system substantially compared
to the results presented in previous sections. Assuming a similar response of
MAE to a gate electric field, we argue that the gate-dependence of MAE can be
probed by measuring the electron tunneling current $I$ through Mn12 as a
function of $V_{\textrm{g}}$ and $V_{\textrm{b}}$.

Figure 9:  (a) Illustration of a
$\textrm{Gr}|\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$ tunnel junction. $V_{g}$ is
gate voltage, and $V_{b}$ is bias voltage. The two graphene layers are
separated by about $1.42\,\textrm{nm}$. (b) Calculated plane-averaged
electrostatic potential of the junction. The inset shows schematic band
alignment between the Mn12 molecule and the two graphene layers. (c) Spin
excitations with zero gate and with a finite gate. $\lvert
m_{\textrm{s}}\rangle$ are spin states with $m_{\textrm{s}}=\pm 10,\,\pm
9,\,\dots,\,0$. (d) Band alignment when a spin excitation enters the bias
window.

First, we construct an atomic structure of the tunnel junction by adding a
flat graphene layer above the $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$
heterostructure, and calculate from first principles the electronic structure
of the junction at zero gate and zero bias. Fig. 9b shows the plane-averaged
electrostatic potential (Hartree potential and ionic potential) across the
tunnel junction. The electrostatic potential around the bottom graphene layer
is similar to that around the top graphene layer. The inset of Fig. 9b depicts
a schematic band alignment of the junction. According to our calculations, the
Dirac points of either graphene layer are at the Fermi level and the LUMO of
Mn12 is $0.230\,\textrm{eV}$ above the Fermi level. To estimate the structural
change in Mn12 caused by a finite bias voltage, we apply static electric
fields on a single isolated Mn12 molecule. At $0.2\,\textrm{V/\AA}$ ($\approx
2.84\,\textrm{V}$ across the junction), the maximal change in Mn-O bonds is
$0.032\,\textrm{\AA}$ and the change in all other chemical bonds is below
$0.003\,\textrm{\AA}$. At $0.1\,\textrm{V/\AA}$, the change in all bond
lengths is within $0.018\,\textrm{\AA}$. Because the $S_{z}$ degeneracy is
broken by the spin-orbit interaction, below the LUMO energy there should be
spin states of Mn12 on an energy scale of $1\,\textrm{meV}$ relative to the
ground state energy of Mn12. Such spin states are detectable via transport
measurements in single-molecule break junctions based on previous studies [77,
78, 79].

The mechanism for detecting spin excitations in a tunneling current is similar
to the approach for magnetic tunnel junctions [80, 81]. The spin excitations
can also be probed in the vertical tunnel junction shown in Fig. 9a. As
established in Section III.3, the MAE of Mn12 can be reduced by electron
doping. Consequently, the energy spacings between spin excitations are
reduced, as illustrated in Fig. 9c. This is likely the case even though Mn12
is fractionally charged and the spin projection $m_{\textrm{s}}$ is no longer
a good quantum index. At finite bias voltages, the Fermi level of the bottom
graphene layer $E_{\textrm{F}}$ is different from that of the top graphene
$E_{\textrm{F}}^{\prime}$. We call the energy range between $E_{\textrm{F}}$
and $E_{\textrm{F}}^{\prime}$ the bias window. Whenever a spin excitation
enters the bias window (see an illustration in Fig. 9d), an additional
transmission path through the excitation of the spin state is opened, leading
to a sudden increase in the tunneling current $I$. The record of multiple spin
excitations allows determination of the energy spacings between these spin
excitations, and thus the MAE of Mn12. As gate voltage varies, the gate-
dependence of MAE can then be probed.

In general, one can expect sensitive change of the tunneling current due to
the spin state of the SMM in the junction. Therefore tunneling measurements
can be used as a more general probe of the spin state of the SMM than just the
MAE. We comment that the usual DFT plus non-equilibrium Green’s function
(DFT$+$NEGF) approach [82] cannot capture the transport signals due to spin
excitations. See Appendix C for DFT$+$NEGF results at zero bias voltage.

##  IV Conclusion and discussion

In conclusion, we have simulated gate field effects on the heterostructures
$\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$, $\textrm{Gr}|\textrm{GaAs}$, and
$\ce{Mn12}|\textrm{Gr}$. In the $\ce{Mn12}|\textrm{Gr}|\textrm{GaAs}$
heterostructure, Mn12 can gain but not lose electrons at moderate doping
levels. The MAE of Mn12 adsorbed on graphene decreases by about 18% at an
electron doping level of $-5.00\times 10^{13}\,\textrm{cm}{}^{-2}$. Such a
decrease in the MAE is due to electron transfer from graphene to Mn12 as well
as the change in the band alignment between Mn12 and graphene. The band
alignment between graphene and GaAs is more sensitive to electron doping than
to hole doping. At an electron doping level of $-1.00\times
10^{13}\,\textrm{cm}{}^{-2}$, the Dirac point of graphene is lifted by about
$0.24\,\textrm{eV}$ relative to the valence band maximum of GaAs. Compared
with electrostatic doping, the adsorption of Mn12 does not have much effect on
the band alignment between graphene and GaAs. The GaAs substrate induces a
band gap of about $2.2\,\textrm{meV}$ in graphene due to the interaction
between GaAs and graphene. It remains to study gate field effects on similar
heterostructures with other ligands around the Mn12O12(COOR)16(H2O)4 molecule.

###### Acknowledgements.

This work was supported by the Center for Molecular Magnetic Quantum
Materials, an Energy Frontier Research Center funded by the U.S. Department of
Energy, Office of Science, Basic Energy Sciences under Award No. DE-SC0019330.
Computations were done using the utilities of the National Energy Research
Scientific Computing Center and the University of Florida Research Computing.

##  Appendix A Basis set in SIESTA calculations

Our input specifications for basis set in SIESTA are given below. For accurate
description of magnetic properties, we applied DZP basis functions for Mn and
O atoms. The remaining atoms (C, H, Ga, and As) employ SZP basis functions.
The basis set for C atoms of Mn12 is different from those for graphene due to
the different chemical environment, and the two basis sets are optimized for
Mn12 and graphene separately.

⬇

%block PAO.Basis

Mn 2

n=4 0 2 P 1

6.4382730 5.1834909

1.000 1.000

n=3 2 2

4.8440717 2.8316085

1.000 1.000

O 2

n=2 0 2

4.3035057 2.9518194

1.000 1.000

n=2 1 2 P 1

4.9096245 3.3500338

1.000 1.000

H 1

n=1 0 1 P 1

5.2026590

1.000

C_Mn12 2

n=2 0 1

4.8900415

1.000

n=2 1 1 P 1

5.8351229

1.000

C_Gr 2

n=2 0 1

4.7137437

1.000

n=2 1 1 P 1

4.7446991

1.000

Ga 3

n=3 2 1

4.9967221

1.000

n=4 0 1

5.9850501

1.000

n=4 1 1 P 1

8.4542189

1.000

As 2

n=4 0 1

5.0625027

1.000

n=4 1 1 P 1

5.2670527

1.000

H.75 1

n=1 0 1 P 1

5.5154366

1.000

%endblock PAO.Basis

##  Appendix B Convergence tests

Table B1 shows test results on $\ce{Mn12}|\textrm{Gr}$ for different smearing
temperatures. When the smearing temperature is reduced from $400\,\textrm{K}$
to $50\,\textrm{K}$, the number of electrons on Mn12 (Mulliken charge) varies
by less than $1\times{10}^{-5}$, the total magnetic moment of the hetero-
structure varies by less than $0.01\ \mu_{B}$, and the position (in energy) of
the Dirac point varies by less than $0.1\ \,\textrm{meV}$. This indicates that
a smearing temperature of $200\,\textrm{K}$ yields reliable electronic and
magnetic properties. Table B2 shows test results on $\ce{Mn12}|\textrm{Gr}$
for different $k$-point grids. All $k$-point grids contain the K point where
the Dirac point lies, and they yield the same energy of the Dirac point. The
variation in both the number of electrons and the total magnetic moment are
insignificant when the $k$-point grid is denser than $3\times 3$. These test
results support that a $6\times 6$ $k$-point grid yields reliable electronic
and magnetic properties.

Table B1:  Convergence test on $\ce{Mn12}|\textrm{Gr}$ with respect to the smearing temperature $T$. The $k$-point grid is fixed at $6\times 6$. $N_{e}$ is the number of electrons on Mn12 (Mulliken charge). $M$ is total magnetic moment. $E_{\textrm{Dirac}}$ is the position (in energy) of the Dirac point relative to the Fermi energy. $T$ (K) | $N_{e}$ |  $M$ ($\mu_{B}$) |  $E_{\textrm{Dirac}}$ (meV)
---|---|---|---
400 | 460.036436 | 19.9992 | 0.1
300 | 460.036438 | 19.9989 | 0.1
200 | 460.036437 | 19.9983 | 0.1
100 | 460.036434 | 19.9964 | 0.1
50 | 460.036437 | 19.9918 | 0.1
Table B2:  Convergence test on $\ce{Mn12}|\textrm{Gr}$ with respect to the $k$-point grid. The smearing temperature $T$ is fixed at 200 K. The quantities are the same with those in Table B1. $k$-grid | $N_{e}$ |  $M$ ($\mu_{B}$) |  $E_{\textrm{Dirac}}$ (meV)
---|---|---|---
$3\times 3$ | 460.036385 | 19.9917 | 0.1
$6\times 6$ | 460.036437 | 19.9983 | 0.1
$9\times 9$ | 460.036448 | 19.9994 | 0.1
$12\times 12$ | 460.036450 | 19.9997 | 0.1

##  Appendix C Electron Transmission

We calculated electron transmission for a $\textrm{Gr}|\ce{Mn12}|\textrm{Gr}$
junction under zero bias by Caroli’s formula [83]. The Green’s functions in
Caroli’s formula are based on the Kohn-Sham Hamiltonian with spin-orbit
coupling. Due to the huge system size (1148 atoms), we included only the
$\Gamma$ point for both SCF and transport calculations. Fig. C1 shows the
electron transmission as a function of energy above the Fermi level (set at
zero). There are four peaks, corresponding to the LUMO (at
$0.2\,\textrm{eV}$), LUMO+1, LUMO+2, and LUMO+3 orbitals respectively. There
are no peaks due to spin excitations, which are a few meV above the Fermi
level in Fig. C1. It requires theoretical development beyond the DFT+NEGF
method to obtain electron transport signals due to spin excitations. It is
noteworthy that both absolute and relative peak positions are tunable by
electron doping as well as by spin alignment.

Figure C1:  Electron transmission of a $\textrm{Gr}|\ce{Mn12}|\textrm{Gr}$
junction under (a) zero doping and (b) an electron doping level of
$\rho=-5\times 10^{13}\,\textrm{cm}^{-2}$. The black (red) line is for the
spin aligned in the easy (a hard) direction. The Fermi level is set to zero.

## References

  * Sessoli  _et al._ [1993a] R. Sessoli, H. L. Tsai, A. R. Schake, S. Y. Wang, J. B. Vincent, K. Folting, D. Gatteschi, G. Christou, and D. N. Hendrickson, Journal of the American Chemical Society 115, 1804 (1993a).
  * Barra  _et al._ [1997] A. L. Barra, D. Gatteschi, and R. Sessoli, Physical Review B 56, 8192 (1997).
  * Pederson and Khanna [1999] M. R. Pederson and S. N. Khanna, Physical Review B 59, R693 (1999).
  * Li  _et al._ [2014] X. G. Li, J. N. Fry, and H. P. Cheng, Physical Review B 90, 125447 (2014).
  * Michalak  _et al._ [2010] L. Michalak, C. M. Canali, M. R. Pederson, M. Paulsson, and V. G. Benza, Physical Review Letters 104, 017202 (2010).
  * Park and Pederson [2004] K. Park and M. R. Pederson, Physical Review B 70, 054414 (2004).
  * Hill [2013] S. Hill, Polyhedron 64, 128 (2013).
  * Artus  _et al._ [2001] P. Artus, C. Boskovic, J. Yoo, W. E. Streib, L. C. Brunel, D. N. Hendrickson, and G. Christou, Inorganic Chemistry 40, 4199 (2001).
  * Barra  _et al._ [2007] A. L. Barra, A. Caneschi, A. Cornia, D. Gatteschi, L. Gorini, L. P. Heiniger, R. Sessoli, and L. Sorace, Journal of the American Chemical Society 129, 10754 (2007).
  * Bian  _et al._ [2005] G. Q. Bian, T. Kuroda-Sowa, N. Gunjima, M. Maekawa, and M. Munakata, Inorganic Chemistry Communications 8, 208 (2005).
  * Bian  _et al._ [2004] G. Q. Bian, T. Kuroda-Sowa, H. Konaka, M. Hatano, M. Maekawa, M. Munakata, H. Miyasaka, and M. Yamashita, Inorganic Chemistry 43, 4790 (2004).
  * Chakov  _et al._ [2006] N. E. Chakov, S. C. Lee, A. G. Harter, P. L. Kuhns, A. P. Reyes, S. O. Hill, N. S. Dalal, W. Wernsdorfer, K. A. Abboud, and G. Christou, Journal of the American Chemical Society 128, 6975 (2006).
  * Chakov  _et al._ [2005a] N. E. Chakov, M. Soler, W. Wernsdorfer, K. A. Abboud, and G. Christou, Inorganic Chemistry 44, 5304 (2005a).
  * Chakov  _et al._ [2003] N. E. Chakov, W. Wernsdorfer, K. A. Abboud, D. N. Hendrickson, and G. Christou, Dalton Transactions 2003, 2243 (2003).
  * Chakov  _et al._ [2005b] N. E. Chakov, L. N. Zakharov, A. L. Rheingold, K. A. Abboud, and G. Christou, Inorganic Chemistry 44, 4555 (2005b).
  * Fonin  _et al._ [2009] M. Fonin, S. Voss, S. Herr, G. de Loubens, A. D. Kent, M. Burgert, U. Groth, and U. Rudiger, Polyhedron 28, 1977 (2009).
  * George  _et al._ [2010] S. M. George, S. I. Mamun, and J. Kim, Inorganic Chemistry Communications 13, 429 (2010).
  * Hoshino  _et al._ [2013] N. Hoshino, T. Shiga, H. Oshio, and T. Akutagawa, Dalton Transactions 42, 4377 (2013).
  * Kuroda-Sowa  _et al._ [2003] T. Kuroda-Sowa, T. Nogami, H. Konaka, M. Maekawa, M. Munakata, H. Miyasaka, and M. Yamashita, Polyhedron 22, 1795 (2003).
  * Kushch  _et al._ [2012] L. A. Kushch, V. D. Sasnovskaya, A. I. Dmitriev, E. B. Yagubskii, O. V. Koplak, L. V. Zorina, and D. W. Boukhvalov, Dalton Transactions 41, 13747 (2012).
  * Lampropoulos  _et al._ [2013] C. Lampropoulos, M. Murugesu, A. G. Harter, W. Wernsdofer, S. Hill, N. S. Dalal, A. P. Reyes, P. L. Kuhns, K. A. Abboud, and G. Christou, Inorganic Chemistry 52, 258 (2013).
  * Lampropoulos  _et al._ [2010] C. Lampropoulos, G. Redler, S. Data, K. A. Abboud, S. Hill, and G. Christou, Inorganic Chemistry 49, 1325 (2010).
  * Lim  _et al._ [2006] J. M. Lim, Y. Do, and J. Kim, European Journal of Inorganic Chemistry 2006, 711 (2006).
  * Ruiz-Molina  _et al._ [2002] D. Ruiz-Molina, P. Gerbier, E. Rumberger, D. B. Amabilino, I. A. Guzei, K. Folting, J. C. Huffman, A. Rheingold, G. Christou, J. Veciana, and D. N. Hendrickson, Journal of Materials Chemistry 12, 1152 (2002).
  * Soler  _et al._ [2001] M. Soler, P. Artus, K. Folting, J. C. Huffman, D. N. Hendrickson, and G. Christou, Inorganic Chemistry 40, 4902 (2001).
  * Tsai  _et al._ [2001] H. L. Tsai, D. M. Chen, C. I. Yang, T. Y. Jwo, C. S. Wur, G. H. Lee, and Y. Wang, Inorganic Chemistry Communications 4, 511 (2001).
  * Tsai  _et al._ [2005] H. L. Tsai, H. A. Shiao, T. Y. Jwo, C. I. Yang, C. S. Wur, and G. H. Lee, Polyhedron 24, 2205 (2005).
  * Verma  _et al._ [2016] S. Verma, A. Verma, A. K. Srivastava, A. Gupta, S. P. Singh, and P. Singh, Materials Chemistry and Physics 177, 140 (2016).
  * Yoon  _et al._ [2003] S. W. Yoon, M. Heu, W. S. Jeon, D. Y. Jung, B. J. Suh, and S. Yoon, Physical Review B 67, 052402 (2003).
  * Zagaynova  _et al._ [2011] V. S. Zagaynova, T. L. Makarova, N. G. Spitsina, and D. W. Boukhvalov, Journal of Superconductivity and Novel Magnetism 24, 855 (2011).
  * Zhao  _et al._ [2004] H. H. Zhao, C. P. Berlinguette, J. Bacsa, A. V. Prosvirin, J. K. Bera, S. E. Tichy, E. J. Schelter, and K. R. Dunbar, Inorganic Chemistry 43, 1359 (2004).
  * Zobbi  _et al._ [2005] L. Zobbi, M. Mannini, M. Pacchioni, G. Chastanet, D. Bonacchi, C. Zanardi, R. Biagi, U. Del Pennino, D. Gatteschi, A. Cornia, and R. Sessoli, Chemical Communications 2005, 1640 (2005).
  * Aubin  _et al._ [1999] S. M. J. Aubin, Z. M. Sun, L. Pardi, J. Krzystek, K. Folting, L. C. Brunel, A. L. Rheingold, G. Christou, and D. N. Hendrickson, Inorganic Chemistry 38, 5329 (1999).
  * Eppley  _et al._ [1995] H. J. Eppley, H. L. Tsai, N. Devries, K. Folting, G. Christou, and D. N. Hendrickson, Journal of the American Chemical Society 117, 301 (1995).
  * Soler  _et al._ [2000] M. Soler, S. K. Chandra, D. Ruiz, E. R. Davidson, D. N. Hendrickson, and G. Christou, Chemical Communications 2000, 2417 (2000).
  * Tsai  _et al._ [2000] H.-L. Tsai, T.-Y. Jwo, G.-H. Lee, and Y. Wang, Chemistry Letters 29, 346 (2000).
  * Coronado  _et al._ [2004] E. Coronado, A. Forment-Aliaga, A. Gaita-Arino, C. Gimenez-Saiz, F. M. Romero, and W. Wernsdorfer, Angewandte Chemie-International Edition 43, 6152 (2004).
  * Rostamzadeh Renani and Kirczenow [2012] F. Rostamzadeh Renani and G. Kirczenow, Physical Review B 85, 245415 (2012).
  * Hu and Wu [2013] J. Hu and R. Q. Wu, Physical Review Letters 110, 097202 (2013).
  * Tsukahara  _et al._ [2009] N. Tsukahara, K. Noto, M. Ohara, S. Shiraki, N. Takagi, Y. Takata, J. Miyawaki, M. Taguchi, A. Chainani, S. Shin, and M. Kawai, Phys Rev Lett 102, 167203 (2009).
  * Zhang [2017] Y. C. Zhang, Journal of Chemical Physics 146, 194705 (2017).
  * Barraza-Lopez  _et al._ [2007] S. Barraza-Lopez, M. C. Avery, and K. Park, Physical Review B 76, 224413 (2007).
  * Zhu  _et al._ [2020] X. C. Zhu, A. Hale, G. Christou, and A. F. Hebard, Journal of Applied Physics 127, 064303 (2020).
  * Blakemore [1982] J. S. Blakemore, Journal of Applied Physics 53, R123 (1982).
  * Tongay  _et al._ [2009] S. Tongay, A. F. Hebard, Y. Hikita, and H. Y. Hwang, Phys. Rev. B 80, 205324 (2009).
  * Bachrach  _et al._ [1981] R. Z. Bachrach, R. S. Bauer, P. Chiaradia, and G. V. Hansson, Journal of Vacuum Science & Technology 18, 797 (1981).
  * Biegelsen  _et al._ [1990a] D. K. Biegelsen, R. D. Bringans, J. E. Northrup, and L. E. Swartz, Physical Review B 41, 5701 (1990a).
  * Biegelsen  _et al._ [1990b] D. K. Biegelsen, R. D. Bringans, J. E. Northrup, and L. Swartz, Phys Rev Lett 65, 452 (1990b).
  * Ohtake  _et al._ [2001] A. Ohtake, J. Nakamura, T. Komura, T. Hanada, T. Yao, H. Kuramochi, and M. Ozeki, Phys. Rev. B 64, 045318 (2001).
  * Munshi  _et al._ [2012] A. M. Munshi, D. L. Dheeraj, V. T. Fauske, D. C. Kim, A. T. J. van Helvoort, B. O. Fimland, and H. Weman, Nano Letters 12, 4570 (2012).
  * Liu  _et al._ [2019] Q. L. Liu, Z. Y. Zhao, and J. H. Yi, Journal of Alloys and Compounds 795, 351 (2019).
  * Hohenberg and Kohn [1964] P. Hohenberg and W. Kohn, Physical Review 136, B864 (1964).
  * Kohn and Sham [1965] W. Kohn and L. J. Sham, Physical Review 140, A1133 (1965).
  * Kresse [1995] G. Kresse, Journal of Non-Crystalline Solids 193, 222 (1995).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Physical Review B 59, 1758 (1999).
  * Soler  _et al._ [2002] J. M. Soler, E. Artacho, J. D. Gale, A. Garcia, J. Junquera, P. Ordejon, and D. Sanchez-Portal, Journal of Physics-Condensed Matter 14, 2745 (2002).
  * Blochl [1994] P. E. Blochl, Physical Review B 50, 17953 (1994).
  * Klimes  _et al._ [2011] J. Klimes, D. R. Bowler, and A. Michaelides, Physical Review B 83, 195131 (2011).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Physical Review Letters 77, 3865 (1996).
  * Troullier and Martins [1991] N. Troullier and J. L. Martins, Physical Review B 43, 1993 (1991).
  * Pack and Monkhorst [1977] J. D. Pack and H. J. Monkhorst, Physical Review B 16, 1748 (1977).
  * Methfessel and Paxton [1989] M. Methfessel and A. T. Paxton, Physical Review B 40, 3616 (1989).
  * Otani and Sugino [2006] M. Otani and O. Sugino, Physical Review B 73, 115407 (2006).
  * Ferreira  _et al._ [2008] L. G. Ferreira, M. Marques, and L. K. Teles, Physical Review B 78, 125116 (2008).
  * Yuan  _et al._ [2018] J. H. Yuan, Q. Chen, L. R. C. Fonseca, M. Xu, K. H. Xue, and X. S. Miao, Journal of Physics Communications 2, 105005 (2018).
  * Ribeiro  _et al._ [2011] M. Ribeiro, L. R. C. Fonseca, and L. G. Ferreira, EPL 94, 27001 (2011).
  * Dudarev  _et al._ [1998] S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Physical Review B 57, 1505 (1998).
  * Boukhvalov  _et al._ [2007] D. W. Boukhvalov, M. Al-Saqer, E. Z. Kurmaev, A. Moewes, V. R. Galakhov, L. D. Finkelstein, S. Chiuzbaian, M. Neumann, V. V. Dobrovitski, M. I. Katsnelson, A. I. Lichtenstein, B. N. Harmon, K. Endo, J. M. North, and N. S. Dalal, Physical Review B 75, 014419 (2007).
  * Bachelet and Schlüter [1982] G. B. Bachelet and M. Schlüter, Physical Review B 25, 2103 (1982).
  * Kleinman [1980] L. Kleinman, Physical Review B 21, 2630 (1980).
  * Fernández-Seivane  _et al._ [2006] L. Fernández-Seivane, M. A. Oliveira, S. Sanvito, and J. Ferrer, Journal of Physics: Condensed Matter 18, 7999 (2006).
  * Sessoli  _et al._ [1993b] R. Sessoli, D. Gatteschi, A. Caneschi, and M. A. Novak, Nature 365, 141 (1993b).
  * Voss  _et al._ [2008] S. Voss, M. Burgert, M. Fonin, U. Groth, and U. Rudiger, Dalton Trans , 499 (2008).
  * Boukhvalov  _et al._ [2010] D. W. Boukhvalov, V. V. Dobrovitski, P. Kogerler, M. Al-Saqer, M. I. Katsnelson, A. I. Lichtenstein, and B. N. Harmon, Inorganic Chemistry 49, 10902 (2010).
  * Wang  _et al._ [1993] D. S. Wang, R. Q. Wu, and A. J. Freeman, Physical Review B 47, 14932 (1993).
  * Brereton  _et al._ [1981] M. P. Brereton, M. K. Cooper, G. R. Dennis, and G. L. D. Ritchie, Australian Journal of Chemistry 34, 2253 (1981).
  * Heersche  _et al._ [2006] H. B. Heersche, Z. de Groot, J. A. Folk, H. S. van der Zant, C. Romeike, M. R. Wegewijs, L. Zobbi, D. Barreca, E. Tondello, and A. Cornia, Phys Rev Lett 96, 206801 (2006).
  * Jo  _et al._ [2006] M. H. Jo, J. E. Grose, K. Baheti, M. M. Deshmukh, J. J. Sokol, E. M. Rumberger, D. N. Hendrickson, J. R. Long, H. Park, and D. C. Ralph, Nano Letters 6, 2014 (2006).
  * Zyazin  _et al._ [2010] A. S. Zyazin, J. W. G. van den Berg, E. A. Osorio, H. S. J. van der Zant, N. P. Konstantinidis, M. Leijnse, M. R. Wegewijs, F. May, W. Hofstetter, C. Danieli, and A. Cornia, Nano Letters 10, 3307 (2010).
  * Zhang  _et al._ [1997] S. Zhang, P. M. Levy, A. C. Marley, and S. S. P. Parkin, Phys. Rev. Lett. 79, 3744 (1997).
  * Du  _et al._ [2010] G. X. Du, S. G. Wang, Q. L. Ma, Y. Wang, R. C. C. Ward, X.-G. Zhang, C. Wang, A. Kohn, and X. F. Han, Phys. Rev. B 81, 064438 (2010).
  * Taylor  _et al._ [2001] J. Taylor, H. Guo, and J. Wang, Physical Review B 63, 121104(R) (2001).
  * Caroli  _et al._ [1971] C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James, Journal of Physics C: Solid State Physics 4, 916 (1971).