# Graphene single electron transistor as a spin sensor for magnetic adsorbates

J. W. González(1),F. Delgado(1), and J. Fernández-Rossier(1,2) (1)
International Iberian Nanotechnology Laboratory (INL), Av. Mestre José Veiga,
4715-330 Braga, Portugal
(2) Departamento de Física Aplicada, Universidad de Alicante, 03690 San
Vicente del Raspeig, Spain

###### Abstract

We study single electron transport through a graphene quantum dot with
magnetic adsorbates. We focus on the relation between the spin order of the
adsorbates and the linear conductance of the device. The electronic structure
of the graphene dot with magnetic adsorbates is modeled through numerical
diagonalization of a tight-binding model with an exchange potential. We
consider several mechanisms by which the adsorbate magnetic state can
influence transport in a single electron transistor: by tuning the addition
energy, by changing the tunneling rate and, in the case of spin polarized
electrodes, through magnetoresistive effects. Whereas the first mechanism is
always present, the others require that the electrode has either an energy or
spin dependent density of states. We find that graphene dots are optimal
systems to detect the spin state of a few magnetic centers.

##  I Introduction

Graphene is a very promising candidate for high precision molecular sensing,
due to its extremely large surface to volume ratio, and its electrically
tunable large conductivity.Schedin07 ; Wernsdorfer_NL ; Pisana On the other
hand, being a zero-gap semiconductor with small mass and small density of
spinfull nuclei, makes graphene a material with potentially large spin
lifetime both, for carriers and host magnetic dopants. Pesin2012  Taken
together, these two ideas naturally lead to the use of graphene as a detector
of the spin state of extrinsic magnetic centers, in the form of magnetic
adatoms, vacancies and spinfull molecules. This connects with recently
reported Wernsdorfer_NL ; Candini-PRB-2011  experiments in which gated
graphene nanoconstrictions, operating in the single electron transport (SET)
regime, showed hysteresis in the linear conductance when a magnetic field is
ramped. This behavior was observed both when the molecular magnets were
intentionally deposited on graphene, as well as in carbon
nanotubes,Wernsdorfer_NatMat  but also in the case of bare graphene
nanojunctions,Candini-PRB-2011  where some type of graphene local moments
Yazyev07 ; Palacios2008 ; Soriano2010  is probably playing a role.

Figure 1: (Color online) (a) Scheme of a graphene constriction with randomly
distributed magnetic centers. (b) Diagram with the system energy levels and
graphene density of states.

The graphene spin sensor experiments of Ref. Wernsdorfer_NL, , are performed
in the Coulomb Blockade regime, showing a vanishing linear conductance except
in the neighborhood of specific values of the gate voltage $V_{g}$. This means
that the graphene nanoconstriction is weakly coupled to the electrodes and has
a charging energy larger than the thermal energy ($T\sim 100$ mK). The height
of the linear conductance peaks is significantly smaller than
$G_{0}=2e^{2}/h$, the quantum of conductance. These conditions imply that
transport takes place in the sequential regime.Beenakker  Thus, current flow
takes place due to sequential tunneling of electrons through the graphene
constriction, which we refer to as the central region in the rest of the
paper, and the entire device behaves like a single electron
transistor.Geim2008 ; review-dot

The aim of this work is to provide a theoretical background to understand how
the magnetic state of localized magnetic moments affects transport through the
graphene nanoconstriction in the SET regime. This is different from previous
works where the influence of the magnetic state of magnetic edges Fede2009
and adsorbed hydrogens Soriano2010  on the conductivity was studied in the
ballistic regime, with a central island strongly coupled to the electrodes,
and also in the diffusive regime,Castro-Neto  as well as SET through graphene
islands with magnetic zigzag edges.Ezawa

The paper is organized as follows. In Sec. II we discuss a tight-binding
Hamiltonian for the graphene island exchanged coupled to the spins of the
magnetic adsorbates. The results of this microscopic calculation justify the
use of a simple single-orbital spin-split model for the SET, discussed in
section III, together with the possible mechanisms that enable the magnetic
sensing. Finally, conclusions are presented in Sec. IV.

##  II Model for graphene island with magnetic adsorbates

###  II.1 Hamiltonian

Our starting point is a microscopic model for electrons confined in a graphene
nanoisland which are exchanged coupled to the magnetic centers. The graphene
central island is described with a tight-binding Hamiltonian for the honeycomb
lattice contained in a rectangular stripe of dimensions $L_{x}\times L_{y}$.
In order to avoid the spin-polarized states formed at the the zigzag
edges,Yazyev ; JP_JR  we impose periodic boundary conditions in one direction
so that the structure only has open edges of armchair type. The coupling to
the magnetic moments of the adsorbate molecules $\vec{m}(i)$ is then assumed
to be a local exchange $J$ or spin-dependent potential, affecting $N$ sites
randomly selected in the graphene central island. For simplicity, we consider
that the magnetic moments of the molecules are all oriented along the same
axis, which we choose as the spin quantization axis. These assumptions are
good approximations in the case of strongly uniaxial TbPC2
molecules.Wernsdorfer_NL  Hence, we can write the Hamiltonian of the graphene
and adsorbates as

$\mathcal{H}=\mathcal{H}_{0}+J\sum_{i}m_{z}(i)S_{z}(i)+eV_{g}\left(N_{TOT}-\hat{N}\right),$
(1)

where $\mathcal{H}_{0}$ is the tight-binding Hamiltonian for $\pi$-electrons
in graphene considering nearest neighbor interactions, $J$ is the strength of
the exchange coupling between the graphene electrons and the magnetic moment
of the molecules, which can take two values, $m_{z}(i)=\pm 1$. Finally, the
last term in the Hamiltonian describes the electrostatic coupling of the total
charge of the dot, which can be either $0$ or $e$, given by the difference in
the number of electrons $\hat{N}$ and the number of carbon atoms $N_{\rm TOT}$
in the central island. $S_{z}(i)$ is the local spin density of the $p_{z}$
electrons in graphene at site $i$,

$S_{z}(i)=\frac{1}{2}\left(c^{\dagger}_{i\uparrow}c_{i\uparrow}-c^{\dagger}_{i\downarrow}c_{i\downarrow}\right),$
(2)

where $c^{\dagger}_{i\uparrow}$ creates an $\uparrow$ electron at the $p_{z}$
orbital of site $i$ of graphene. In the following we assume that magnetic
fields controlling the spin orientation of the adsorbates, are applied along
the plane of graphene so that it is a good approximation to neglect the
diamagnetic coupling to the graphene electrons.

There are several independent microscopic mechanisms for spin dependent
interaction between magnetic adsorbates and the graphene $\pi$ electrons that
can be modeled with Eq. (1). In the case of magnetic molecule such as TbPC2,
used in in Ref. Wernsdorfer_NL, , the magnetic Tb atom is separated from the
graphene electrons by the non-magnetic atoms of the molecule, and the most
likely mechanism for spin coupling is kinetic exchange.Anderson61  This
coupling will generate a local Kondo-exchange between graphene electrons and
the molecules.Schrieffer-Wolff  More complicated scenarios, like coupling of
graphene electrons to unpaired electrons in the organic rings of the
molecules, would imply that every molecule affects several sites in graphene.
Direct dipolar coupling would also affect several sites per molecule, but the
average magnetic field created by a magnetic moment of 5 $\mu_{B}$ at 0.5 nm
on a disk with an area around 400 nm2, the graphene constriction area in Ref.
Wernsdorfer_NL, , is smaller than 1 $\mu$T, which would yield a negligible
maximal Zeeman coupling of neV per molecule.

###  II.2 Relevant energy scales

The reported dimensions of the central region, $L_{x}\simeq L_{y}\simeq 20$
nm, lead to an energy spacing $\Delta\epsilon$ of the single particle spectrum
much larger than the temperature and the charging energy.Geim2008 ; review-dot
We also assume that the exchange induced shifts are smaller than the single
particle splitting. As a result, the effect of exchange is to shift the bare
energy levels, without mixing them. Thereby, we can safely assume that
electrons tunnel through just one of the single particle levels, which might
be spin-split due to exchange with molecules. We assume that the charge of the
central island fluctuates between $q=0$ and $q=-|e|$, and that the transport
level is the lowest unoccupied level of the central island spectrum. The
energy of the transport level reads:

$\epsilon_{T}^{\sigma}=\epsilon_{0}+\sigma\frac{\delta}{2}-|e|V_{g},$ (3)

where $\epsilon_{0}$ is the single particle electron level, $\sigma=\pm 1$
denotes the spin direction and $\delta$ is the magnitude of the spin
splitting, which is a functional of the magnetic landscape $\\{m_{z}\\}$.

Within our model, a given magnetic landscape is defined by the location $i$
and the magnetic state $m_{z}(i)$ of the $N$ magnetic adsorbates. In Figs. 2
(a, b), we plot the value of $\delta$ for all the possible magnetic states of
a given arrangement of $N=10$ atoms, for two different values of $J$. This
choice corresponds to the estimated number of molecules in Ref.
Wernsdorfer_NL, . These figures show a correlation between the magnitude of
the splitting $\delta$ and the total magnetization $M_{T}$. The dispersion of
$\delta$ for a fixed total magnetization $M_{T}$ is the outcome of indirect
exchange couplingBrey-Das-Sarma

For comparison with the experiments, it is worth considering two extreme
magnetic landscapes. At large external field, all the magnetic moments are
aligned, i.e., $m_{z}(i)=+1$. We refer to this as the ferromagnetic (FM)
landscape. At magnetic fields smaller than the coercive field of the magnetic
molecules, their average magnetization should be zero and thus,
$\sum_{i}m_{z}(i)=0$. We refer to these cases as non-magnetic (NM). In order
to sample the positional disorder we perform an average over positional
configurations, both for NM and FM cases. For a fixed spin choice
$\\{m_{z}\\}$ with $M_{T}=0$, an average over positional configurations yields
$\langle\delta\rangle_{\rm pos}=0$. The reciprocal statement is also true: for
a fixed positional configuration, an average over all the magnetic landscapes
with $M_{T}=0$ also yields an average $\langle\delta\rangle_{\rm spins}=0$.

Figure 2: Spin splitting $\delta$ for a sample of dimensions $L_{y}=N_{y}a$
and $L_{x}=N_{x}\sqrt{3}a$, with $N_{x}=15$ and $N_{y}=17$, being $a=2.46\AA$
the lattice constant of graphene. $\delta$ versus total magnetization $M_{T}$
for a single spatial distribution of $N=10$ magnetic impurities with $J=5$ meV
(a) and $J=10$ (b). The red solid line indicates the average value
$\langle\delta\rangle_{\rm pos}$. Average $\langle\delta\rangle_{\rm pos}^{\rm
FM}$ versus number of magnetic molecules $N$ for $J=5$ meV (c) and exchange
energy $J$ for $N=4$ (d). The error bars correspond to the standard deviation.

In Fig. 2 we plot the average $\langle\delta\rangle_{\rm pos}^{\rm
FM}-\langle\delta\rangle_{\rm pos}^{\rm NM}=\langle\delta\rangle_{\rm
pos}^{\rm FM}$ over $500$ realizations as a function of the number of
molecules $N$ [Fig. 2(c)] and as a function of the molecule electron exchange
$J$ [Fig. 2(d)]. We have also calculated $\langle\delta\rangle_{\rm pos}^{\rm
FM}$ fixing the number of magnetic centers $N$, the strength of the coupling
$J$ and changing $N_{\rm TOT}=N_{x}N_{y}$, the total number of carbon atoms in
the island. We find that the results of all these simulations can be
summarized in the following equation:

$\left\langle\delta\right\rangle_{\rm pos}^{\rm FM}\approx\frac{N}{N_{\rm
TOT}}J.$ (4)

Whereas this result has been obtained from exact numerical diagonalization of
the Hamiltonian, this dependence can be rationalized using first order
perturbation theory, which yields the spin dependent shift of the transport
level:

$\Delta\epsilon_{T}^{\sigma}=\frac{\sigma}{2}J\sum_{i}^{N}|\phi_{T}(i)|^{2}m_{i},$
(5)

where $\phi_{T}(i)$ is the $J=0$ wave function of the transport orbital. We
now use $|\phi_{T}(i)|^{2}\simeq\frac{1}{N_{\rm TOT}}$ so that we can
approximate:

$\Delta\epsilon_{T}^{\sigma}\simeq\frac{\sigma}{2}\frac{J}{N_{\rm TOT}}M_{T}$
(6)

Using the fact that $M_{T}=N$ for the FM configurations and $0$ for the NM
ones, we arrive to Eq. (4).

##  III Spin-Split Single Orbital Model for SET

In this section we discuss SET across a central island with a single spin
split particle level. This is justified by the results of the previous
section. We obtain expressions for the current of the system and we discuss
the conditions under which the conductance depends on the magnetic state of
the single electron transport.

###  III.1 Single electron transistor with a spin-split single orbital model

We consider single electron transport though a spin split single transport
level,Recher2000  with energy $\epsilon_{T}^{\sigma}$. We assume that the
occupation of the transport level can be either 0 or 1, the doubly occupied
configuration being much higher in energy. Within these approximations, the
transport level has three relevant many body states: uncharged, and the two
charged with $\uparrow$ or $\downarrow$ spins. In the zero-applied bias limit,
each of these states will be occupied according to the thermal equilibrium
distribution, which we denoted as $P_{0}$, $P_{\uparrow}$ and $P_{\downarrow}$
respectively. In the SET regime we are interested and within the linear
response ($eV_{bias}\ll k_{B}T$), transport will be enabled only when the
addition energy lies within the thermally broadened transport window defined
by the applied bias.

Under these approximations, the current flowing from the left electrode to the
central island is given by:

$I=e\sum_{\sigma}\left\\{P_{0}W_{0\to\sigma}^{L}-P_{\sigma}W_{\sigma\to
0}^{L}\right\\},$ (7)

where $W_{0\to\sigma}^{L}$ and $W_{\sigma\to 0}^{L}$ are rates for electron
tunneling from the left electrode to the dot and vice-versa. Continuity
equation ensures that this current is identical to the current flowing towards
the right electrode and, thereby, equal to the net current flow. The tunneling
rates for electron tunneling out of and into the dot Haug_Jauho_book_1996  are
given respectively by

$W_{\sigma\to
0}^{L}=\frac{2\pi}{\hbar}|T_{L}|^{2}\rho_{L}(\epsilon_{T}^{\sigma})\left[1-f(\Delta_{\sigma}+V_{bias})\right],$
(8)

and

$W_{0\to\sigma}^{L}=\frac{2\pi}{\hbar}|T_{L}|^{2}\rho_{L}(\epsilon_{T}^{\sigma})f(\Delta_{\sigma}+V_{bias}),$
(9)

where $T_{L}$ is the strength of the dot - left electrode coupling, and

$\Delta_{\sigma}\equiv\epsilon_{T}^{\sigma}-E_{F},$ (10)

are the spin dependent addition energies. Importantly, both $\delta$ and
$V_{g}$ appear on equal footing, as additive quantities in this equation. The
density of states of the left electrode, evaluated at the spin-dependent
transport level energy, is denoted by $\rho_{L}(\epsilon_{T}^{\sigma})$, while
$f(\epsilon_{\sigma})=\left(e^{\beta\Delta_{\sigma}}+1\right)^{-1}$ denotes
the Fermi function. The electrode Fermi energy $E_{F}$ is taken to change
linearly with the bias voltage $V_{bias}$. In the zero bias limit, the linear
conductance reads:

$G=G_{0}\sum_{\sigma}\left(P_{0}+P_{\sigma}\right)\frac{\Gamma_{\sigma}}{k_{B}T}{\rm
Sech}^{2}\left(\frac{\beta\Delta_{\sigma}}{2}\right),$ (11)

where

$\Gamma_{\sigma}/\hbar=\frac{2\pi}{\hbar}|T_{L}|^{2}\rho_{L}(\epsilon_{T}^{\sigma}),$
(12)

is the single particle tunneling rate between the electrode and the transport
level.

###  III.2 Influence of magnetic state on conductance

From the above discussion it is apparent that, for a given gate potential
$V_{g}$ and temperature $T$, the linear conductance depends on the magnetic
landscape affecting the central island through two classes of independent
mechanisms, illustrated in Fig. 3:

  1. 1.

The change of the addition energies $\Delta_{\sigma}$ which, as we show below,
would result in a lateral shift of the $G(V_{g})$ resonance curve [Fig. 4(a)].

  2. 2.

The change of the electron lifetime $\Gamma=\sum_{\sigma}\Gamma_{\sigma}$,
that would result in a vertical resizing of the $G(V_{g})$ resonance curve
[Fig. 5(a)].

Figure 3: (Color online) (a) Scheme showing the transport level energy
splitting. Scheme of spin-dependence of transport due to: (b) detuning of the
transport level with respect to the electrode Fermi energy, (c)
magnetoresistance associated to spin polarized electrode(s), and (d) energy
dependent tunneling rates.

In the first mechanism the change in the magnetic state modifies the value of
$\delta$, which must have a similar effect than changing the gate potential.
It resembles the magneto-Coulomb effect,Ono ; Vanwees  by which the applied
magnetic field changes the Fermi energy of the electrode, shifting the
$G(V_{g})$ curves. However, this first mechanism necessarily implies a change
of sign of the variation of $G$ as the gate potential is scanned along the
resonance (see top panel of Fig. 4). Importantly, this is not observed in the
experiments with magnetic molecules,Wernsdorfer_NL  but it is observed in the
case of graphene nanoconstrictions.Candini-PRB-2011

Figure 4: (Color online) Conductance in units of $g_{0}=G_{0}\beta\Gamma$ as a
function of the gate voltage $eV_{g}$ and the level splitting $\delta$. Panel
(a) shows the normalized conductance as a function of the gate potential
$e\beta V_{g}$, where the labels corresponds to the different $\beta\delta$
values. Panel (b) shows the conductance as a function of $\beta\delta$ for
several values of the gate $e\beta V_{g}$. For the sake of clarity all curves
have been displaced by $0.1$. In (c, d) we present a contour plot of the dot
charge (defined as $Q=P_{\uparrow}+P_{\downarrow}$) and the dot magnetization
($m=P_{\uparrow}-P_{\downarrow}$) as a function of the gate voltage $eV_{g}$
and level splitting $\delta$.

Motivated by the behavior reported in Ref. Wernsdorfer_NL, , we pay attention
also to the second mechanism. For spin un-polarized transport, the change in
the transport energy level results in a change on tunneling rate
$\Gamma/\hbar$ only if the electrode density of states depends on energy,
which is exactly the case of graphene. For spin-polarized transport, the
relative orientation of the electrode and island magnetic moment give rise to
magneto-resistive effects that are accounted for by the changes in
$\Gamma_{\sigma}$

###  III.3 Transport for constant tunneling rates

We now discuss our transport simulations for the graphene single electron
transistor spin sensor. We focus on the first spin sensing mechanism in a
single electron transistor: changes in spin splitting of the transport level
produce changes in addition energies $\Delta_{\sigma}$ (figure (3)b). For that
matter, we neglect both the energy and spin dependence of the tunneling rates
$\Gamma_{\sigma}$. In Fig. 4(a) we show the linear conductance, in units of
$g_{0}=\beta\Gamma G_{0}$, as a function of gate voltage, for several values
of the transport level splitting $\delta$, in units of $k_{B}T$. It is
apparent that the Coulomb Blockade peaks undergoes a lateral shift, as
expected from the fact that $V_{g}$ and $\delta$ appear on equal footing on
the spin dependent addition energies. At $\Delta=0$ the two spin channels
contribute. Therefore, as we increase $\delta$, the height of the conductance
peaks decreases, because one of the two spin channels is removed from the
transport window of width $k_{B}T$

In figure (4)b we plot the variations in the linear conductances as a function
of the spin splitting $\delta$, for several values of $V_{g}$. We see two
types of curves. For values of $V_{g}$ such that the transport level is
occupied, as we increase $\delta$ the transport level is pushed downwards,
away from the elastic transport window, switching off the transistor
conductance. In contrast, for $V_{g}$ such that the transport level is empty
for $\delta=0$, lying above the elastic transport window, ramping $\delta$
makes one of the two spin states of the transport level enter the transport
window, giving rise to the double peak structure. The fact that $\delta$ and
$V_{g}$ play analogous roles is illustrated in Figs. 4(c, d), where we show
the average magnetization and occupation of the transport level in the phase
diagram defined by these two variables.

###  III.4 Transport with energy dependent tunneling rates

We now consider the second mechanism for spin sensing in a single electron
transistor: changes in spin splitting of the transport level produce changes
in the tunneling rates $\Gamma_{\sigma}$. This can happen for two reasons:

  1. 1.

One of the electrodes is spin polarized, so that
$\Gamma_{\downarrow}\neq\Gamma_{\uparrow}$. Spin polarized transport is
sensitive to the product of the magnetic moment of electrode and central
island. This type of effect has been thoroughly discussed in the case of SET
with ferromagnetic electrodes.Barnas98 ; Seneor2007

  2. 2.

The density of states of the electrode depends on energy. Thus, changes in the
transport level change $\Gamma_{\sigma}$, for both spins. This is a natural
scenario for graphene electrodes. dots_Ensslin ; Sols-Guinea

Let us consider first the case of spin polarized electrodes. We do the
assumption that the density of states are spin dependent but energy
independent: $\rho_{\sigma}=\rho_{0}(1+\sigma{\cal P})$, with ${\cal P}$ the
electrode polarization. In Fig. 5(a) we plot the linear conductance vs $V_{g}$
curves for several values of the splitting $\delta$, assuming a large value of
the electrode spin polarization ${\cal P}=0.9$. It is apparent that, on top of
the shift of the resonance curve whose origin was discussed in the previous
subsection, there is a change in the amplitude of the curve. Notice that
$G(\delta=k_{B}T)$ and $G(\delta=2k_{B}T)$ are smaller than $G(0)$ for all
values of $V_{g}$. In this specific sense, the gate-independent spin contrast
is similar to the experimental report with magnetic molecules.Wernsdorfer_NL
In Fig. 5(a) we show the linear conductance as a function of $\delta$ for
different values of $V_{g}$. It is apparent that, as opposed to the case of
non-magnetic electrodes shown in Fig. 5(b), the function $G(\delta$) is no
longer an even functions, reflecting the magneto-resistive behaviour.
Basically, transport is favored when the spin polarization of the electrode
and the central island are parallel.

Figure 5: (Color online) Normalized Conductance for ferromagnetic electrodes
as a function of the level splitting $\delta$ for several gate values. In (a,
c) the density of states depends of the polarization
$\rho=\rho_{0}(1+\sigma{\cal P})$, in this particular case we take a
polarization of ${\cal P}=0.9$. In (b, d) the electrode density of states is
linear with the energy,
$\rho(\Delta_{\sigma})=\rho(\epsilon_{0}+\sigma\delta/2-E_{D})$. For the sake
of clarity the conductance curves in (c, d) have been displaced by $0.2$ and
$4$ units respectively.

We now consider a non-spin polarized electrode with an energy dependent
density of states. This scenario occurs naturally in graphene. If we consider
idealized graphene electrodes, neglecting effects of interactions, disorder
and confinement, we have $\rho(\epsilon)=\rho_{0}|\epsilon-E_{D}|$, where
$E_{D}$ is the Dirac point. The $G(V_{g})$ curves, shown in Fig. 5(b) for
different values of $\delta$, shift and change amplitude. The shift is related
to the change of the addition energies, discussed in the previous subsection,
and the change in amplitude comes from the variation of the tunneling rate as
the transport level scans the energy dependent density of states of the
electrode.

In Fig. 5(d) we plot $G(\delta)$ for several values of $V_{g}$. The curves are
similar to the case with energy independent tunneling rates, except for the
dip at zero $\delta$ which occurs because we chose the bare transport level
right at Dirac point. This is the most favorable choice to maximize the effect
of energy dependence of $\Gamma$. From our results, and given the fact that
experimentally is not possible to put the Fermi energy arbitrarily close to
the Diract point,Geim-closest  we find it unlikely that this effect is playing
a role in the experiments.

###  III.5 Sensitivity of the single electron transistor spin sensor

We now discuss the sensitivity of the spin sensor based on the graphene single
electron transistor, as described by our model, neglecting changes in
$\Gamma$. From Fig. 4(a) we propose, as rule of thumb, that variations of
$\delta$ similar or larger than $k_{B}T$ can be resolved. Estimating $\delta$
from the case of fully spin polarized magnetic adsorbates, given in Eq. (4),
we find a relation between the minimal number of magnetic centers $N$ that can
be detected, and the temperature and exchange constant:

$\frac{N}{N_{\rm TOT}}>\frac{k_{B}T}{J}.$ (13)

It is apparent that decreasing the temperature, or increasing the spin-
graphene exchange coupling, increases the sensitivity of the device (makes it
possible to detect a smaller concentration of molecules). For instance, at
$100$ mK, and taking $N_{\rm TOT}\sim 15000$, which corresponds to an
approximate area of $400$ nm2, one could detect $10$ molecules for an exchange
coupling $J\gtrsim 15$ meV.

Recent reports have shown that it is possible to fabricate graphene nano
islands with lateral dimensions of 1nm.Barreiro2012  These dots have $N_{\rm
TOT}<100$. Thus, they would permit the detection of the spin of a single
magnetic adsorbate provided that $k_{B}T$ is kept hundred times smaller than
$J$. For $T=100$ mK this implies $J>1$ meV. Interestingly, in such a small dot
Coulomb Blockade persists even at room temperature, but increasing $k_{B}T$
keeping the sensitivity would require also to increase $J$.

##  IV  Discussion and conclusions

Because of its structural and electronic properties, graphene is optimal for a
spin sensor device. Being all surface, the influence of adsorbates on
transport should be larger than any other bulk material. Because of the linear
relation momentum and large Fermi velocity, energy level spacing in graphene
nano structures can easily be larger than the temperature, the tunneling
induced broadening, and the perturbations created by the adsorbates. One of
the consequences is that single electron transport takes place through a
single orbital level.

Our simulations show how the spin splitting $\delta$ of the transport level is
sensitive to the average magnetization of the magnetic adsorbates, which is
controlled by application of a magnetic field along the plane of graphene, to
avoid diamagnetic shifts. On the other hand, the linear conductance $G$ of the
single electron transistor depends on $\delta$, which accounts for the sensing
mechanism. More specifically, $G$ depends on $\delta$ due to either changes in
the spin dependent addition energies $\Delta_{\sigma}$ or changes in the
electrons lifetime $\Gamma_{\sigma}$. The first is independent of the nature
of the electrodes, whereas the second only happens if they are magnetic or
have an energy dependent density of states.

We have shown how, within an independent particle model and in the single
electron transport regime, the energy dependence of the graphene electrode
density of states can only be relevant if the transport energy level is fine
tuned to the Dirac point. However, this fine tuning is quite unlike to happen
in experimental conditionsGeim-closest . Still, the combined action of
disorder and Coulomb interaction could give rise to a so called Coulomb Gap in
the density of states of graphene, that might make the tunneling rates depend
on the energy.Coulomb-gap-graphene ; Coulomb-gap-graphene1 ; Coulomb-gap-
graphene2

Finally, we have assumed that both the edges of the graphene island and
graphene electrodes are non-magnetic. Our discussion of the effect of spin-
polarized electrodes on the transport properties of the device would be valid
for electrodes with ferromagnetic zigzag edges. A second possibility, out of
the scope of this work, is to consider a graphene single electron transistor
whose central island has ferromagnetic edges. This case has been already
studied. Ezawa

In conclusion, we have studied the mechanisms by which a graphene single
electron transistor could work as a sensor of the magnetic order of magnetic
atoms or molecules adsorbed on the graphene central region. Our work has been
motivated in part by recent experimental works, Wernsdorfer_NL ; Candini-
PRB-2011 . Whereas further work is still necessary to nail down the physical
mechanisms for the spin sensing principles underlying the experimental work,
our study provides a conceptual framework for graphene single electron
transistor spin sensors.

This work has been financially supported by MEC-Spain (Grant Nos.
FIS2010-21883-C02-01 and CONSOLIDER CSD2007-0010) as well as Generalitat
Valenciana, grant Prometeo 2012-11. We thank A. Candini for useful comments on
the manuscript.

## References

  * (1) F. Schedin, A. Geim, S. Morozov, E. Hill, P. Blake, M. Katsnelson, and K. Novoselov, Nature Materials 6, 652 (2007).
  * (2) A. Candini, S. Klyatskaya, M. Ruben, W. Wernsdorfer, and M. Affronte, Nano Letters 11, 2634 (2011).
  * (3) S. Pisana , P. M. Braganca , E. E. Marinero and B. A. Gurney, Nano Letters 10, 341 (2010).
  * (4) D. Pesin and A. MacDonald, Nature Materials 11, 409(2012).
  * (5) A. Candini, C. Alvino, W. Wernsdorfer, M. Affronte, Physical Review B 83, 121401 (2011).
  * (6) M. Urdampilleta, S. Klyatskaya, J. Cleuziou, M. Ruben, and W. Wernsdorfer, Nature Materials 10, 502 (2011).
  * (7) O. V. Yazyev and L. Helm, Physical Review B 75, 125408 (2007).
  * (8) J. J. Palacios, J. Fernández-Rossier, L. Brey, Physical Review B 77, 195428 (2008).
  * (9) D. Soriano, F. Muñoz-Rojas, J. Fernández-Rossier, J. J. Palacios, Physical Review B 81, 165409 (2010).
  * (10) C.W.J. Beenakker, Physical Review B 44, 1646 (1991).
  * (11) L. Ponomarenko, F. Schedin, M. Katsnelson, R. Yang, E. Hill, K. Novoselov and A. K. Geim, Science 320, 356 (2008).
  * (12) J. Güttinger, F. Molitor, C. Stampfer, S. Schnez, A. Jacobsen, S. Dröscher, T. Ihn and K. Ensslin, Rep. Prog. Phys. 75, 126502 (2012).
  * (13) F. Muñoz-Rojas, J. Fernández-Rossier, J. J. Palacios, Physical Review Letters 102, 136810 (2009).
  * (14) C. H. Lewenkopf, E. R. Mucciolo, and A. H. Castro Neto, Physical Review B 77, 081410(R) (2008).
  * (15) M. Ezawa, Physical Review B 77, 155411 (2008).
  * (16) J. Fernández-Rossier and J. J. Palacios, Phys. Rev. Lett. 99, 177204 (2007).
  * (17) O. V. Yazyev, Rep. Prog. Phys. 73, 056501 (2010).
  * (18) L. Brey, H. A. Fertig, and S. Das Sarma, Phys. Rev. Lett.99, 116802 (2007).
  * (19) P. W. Anderson Phys. Rev. 124, 41 (1961).
  * (20) J. R. Schrieffer and P. A. Wolff Phys. Rev. 149, 491 (1966).
  * (21) P. Recher, E. V. Sukhorukov, and Daniel Loss Phys. Rev. Lett. 85, 1962 (2000).
  * (22) H. Haug and A.-P. Jauho, Quantum kinetics in transport and optics of semi-conductors (Springer-Verlag, Berlin, 1996).
  * (23) K. Ono, H. Shimada, and Y. Ootuka, Journal of the Physical society of Japan 66, 1261 (1997).
  * (24) S. J. Van Der Molen, N. Tombros, and B. J. Van Wees, Physical Review B 73, 220406 (2006).
  * (25) J. Barnas, A. Fert Phys. Rev. Lett. 80 1058 (1998).
  * (26) P. Seneor, A. Bernand-Mantel and F. Petroff, Journal of Physics: Condensed Matter 19, 165222 (2007).
  * (27) A. Barreiro, H.S.J. van der Zant, L.M.K. Vandersypen, Nano Letters 12, 6096 (2012).
  * (28) A. Mayorov, D. C. Elias, I. S. Mukhin, S. V. Morozov, L. Ponomarenko, K. S. Novoselov, A. K. Geim, and R. V. Gorbachev, Nano Letters 12, 4629 (2012).
  * (29) S. Dröscher, H. Knowles, Y. Meir, K. Ensslin, and T. Ihn, Physical Review B 84, 073405 (2011).
  * (30) A. L. Efros, and B. I. Shklovskii, J. Phys. C: Solid State Phys. 8 49 (1971).
  * (31) B. Terrés, J. Dauber, C. Volk, S. Trellenkamp, U. Wichmann, and C. Stampfer, Applied Physics Letters 98, 032109 (2011).
  * (32) X. Liu, J. B. Oostinga, A. F. Morpurgo, and L. M. K. Vandersypen, Physical Review B 80, 121407(R) (2009).
  * (33) F. Sols, F. Guinea, and A. H. Castro Neto, Physical Review Letters 99, 166803 (2007).


††thanks: These authors contributed equally††thanks: These authors contributed
equally

# Structure, Stability and Superconductivity of N-doped Lutetium Hydrides at
kbar Pressures

Katerina P. Hilleke  Department of Chemistry, State University of New York at
Buffalo, Buffalo, NY 14260-3000, USA  Xiaoyu Wang  Department of Chemistry,
State University of New York at Buffalo, Buffalo, NY 14260-3000, USA  Dongbao
Luo  Department of Chemistry, State University of New York at Buffalo,
Buffalo, NY 14260-3000, USA  Nisha Geng  Department of Chemistry, State
University of New York at Buffalo, Buffalo, NY 14260-3000, USA  Busheng Wang
Department of Chemistry, State University of New York at Buffalo, Buffalo, NY
14260-3000, USA  Eva Zurek<EMAIL_ADDRESS>Department of Chemistry, State
University of New York at Buffalo, Buffalo, NY 14260-3000, USA

###### Abstract

The structure of the material responsible for the room temperature and near
ambient pressure superconductivity reported in an N-doped lutetium hydride
[Nature, 615, 244 (2023)] has not been conclusively determined. Herein,
density functional theory calculations are performed in an attempt to uncover
what it might be. Guided by a range of strategies including crystal structure
prediction and modifications of existing structure types, we present an array
of Lu-N-H phases that are dynamically stable at experimentally relevant
pressures. Although none of the structures found are thermodynamically stable,
and none are expected to remain superconducting above $\sim$17 K at 10 kbar, a
number of metallic compounds with _fcc_ Lu lattices – as suggested by the
experimental X-ray diffraction measurements of the majority phase – are
identified. The system whose calculated equation of states matches best with
that measured for the majority phase is fluorite-type LuH2, whose 10 kbar
superconducting critical temperature was estimated to be 0.09 K using the
Allen-Dynes modified McMillan equation.

##  I Introduction

Heike Kamerlingh Onnes’ 1911 discovery of mercury’s entrance into a
“new…superconductive state” at very low temperatures, where all electrical
resistance vanished [1], marked the beginning of a quest: could such a state
be observed at room temperature? Ever since, scientists have sought this “holy
grail”, steadily breaking through barriers such as the boiling point of liquid
nitrogen [2], 100 K [3], then near 200 K [4, 5, 6]. The latter breakthrough
can be directly traced back to Ashcroft’s proposal that hydrogen-rich alloys,
metallized at conditions of extreme pressure, albeit less extreme than those
required to metallize pure hydrogen, would be high-temperature phonon-mediated
superconductors [7]. It also marked a paradigm shift defined by a close
synergy between theory and experiment, with computations either predicting the
most promising superconducting phases or being instrumental in characterizing
the synthesized compounds [8, 9, 10].

For the materials with the highest superconducting critical temperatures,
$T_{c}$s, that were found two things were true: they featured high hydrogen
content and they required immense pressures – approaching those found in the
center of the Earth (350 GPa) – for stability. One prominent class of these
high-pressure high-temperature compounds are known as the “superhydrides”. All
of them are characterized by clathrate-like hydrogen-based lattices that
encapsulate an electropositive metal atom, typically an alkaline or rare
earth. Examples of compounds that have been both predicted and synthesized
include CaH6 ($T_{c}$ = 210-215 K, 160-172 GPa) [11, 12], LaH10 ($T_{c}$ = 260
K, 200 GPa) [13, 14], YH9 ($T_{c}$ = 262 K, 182 GPa) [15], YH6 ($T_{c}$ = 224
K, 166 GPa) [16], and mixed La/Y ternary hydrides [17, 18]

Clearly, the most prominent metal atoms in these phases are yttrium and
lanthanum, with supporting roles played by calcium, scandium, and other rare
earths. However, most of the heavier lanthanide hydrides are not expected to
be as promising because of the suppressive influence of the _f_ electrons on
superconductivity, with maximum $T_{c}$s decreasing rapidly once past La [19,
20]. As a result, the hydrides of lutetium received relatively little
attention despite the fact that the filled $4f$ shell of the metal is
chemically unreactive rendering its electronic properties similar to those of
Sc, Y and La… till now.

An early theoretical study generated a Lu-H convex hull using known
polyhydride structures, finding LuH4, LuH6, LuH9, and LuH10 as being
thermodynamically stable at various pressures up to 400 GPa [20]. Another
identified a unique $Immm$ structure for LuH8 with an estimated $T_{c}$ of
81-86 K at 300 GPa, based on a distorted version of the backbone of the
$Fm\bar{3}m$ LaH10 phase [21]. A theoretical comparison between the hydrides
of the rare earth elements with filled vs. unfilled $f$-states – Tm, Yb, and
Lu, found LuHn (_n_ =4-8, 10) phases either on or very near the Lu-H convex
hull at relatively low pressures (less than 200 GPa) [22]. Notably, LuH6, with
the same $Im\bar{3}m$ symmetry as CaH6, had an estimated $T_{c}$ of 273 K
(matching the melting point of ice) at 100 GPa. The filled $f$-shells of Lu
and Yb were suggested to confer a strong degree of phonon softening, thereby
resulting in a high electron-phonon coupling. Finally, a theoretical
investigation of trends in superconductivity proposed high-pressure $Cc$ LuH7
and $C222$ LuH12 phases, with the latter predicted to undergo a
superconducting transition below 6.7 K at 150 GPa [19].

On the experimental side, a recent work reported the synthesis of a Lu
hydride, suggested to be $Pm\bar{3}n$ Lu4H23, with a measured $T_{c}$ of 71 K
at 218 GPa [23]. This structure has previously been observed in experimental
studies in the La-H [24], Ba-H [25], and Eu-H [26] systems.

Thus, with reported $T_{c}$s of the superhydrides reaching temperatures not
uncommon for a typical winter-day in upstate New York, the focus of research
changed to predicting and synthesizing materials that could maintain high
$T_{c}$s, but at lower pressures, with the ultimate goal of realizing
superconductivity at ambient temperature and pressure. As the structures and
superconducting properties of the binary hydrides had been exhaustively
searched with no such candidate found, computations turned towards predicting
ternary hydrides that remained dynamically stable to pressures below 100 GPa
[27, 28, 29], or boron-carbon analogues of the superhydrides that were stable
at 1 atm [30].

It was therefore quite exciting when a recent experimental manuscript reported
superconductivity near room-temperature, $T_{c}$ = 294 K, at a very moderate
pressure of 10 kbar (1 GPa) in a nitrogen-doped lutetium hydride phase [31].
This pressure is low enough so that it becomes feasible to use pressure-
quenching [32] to stabilize the material to ambient conditions, or to use
careful strain engineering to achieve the desired superconductivity.
Unfortunately, though a variety of techniques including X-ray diffraction
(XRD), energy-dispersive X-ray measurements, elemental analysis and Raman
spectroscopy were used to characterize the superconducting material, its
composition and structure could not be fully resolved [31].

On the basis of the XRD and Raman analysis, the proposed room-temperature
superconducting material (referred to as compound A by the authors) was
indexed with space group $Fm\bar{3}m$, and both compound A and a minor
product, which was dubbed compound B, were suggested to consist of an _fcc_ Lu
network with additional N and H uptake [31]. At pressures above $\sim$30 kbar,
the superconducting compound A was found to undergo a pressure-induced
transition to a non-superconducting structure involving a symmetry reduction
of the Lu lattice to orthorhombic $Immm$ symmetry. The superconducting
compound was also observed to undergo a sequence of color changes
corresponding to structural transitions as pressure was applied, from blue to
pink (marking transition to the high-$T_{c}$ superconductor), to red.

Follow-up studies have, however, suggested that this color change is derived
in fact from pure LuH2 [33, 34]. Experiments reported no evidence for
superconductivity down to 1.5 K in LuH2 [33], or in LuH2±xNy from ambient
pressure to 6.3 GPa down to 10 K [34]. Moreover, DFT calculations [35]
concluded that LuH2 in the fluorite structure is the dominant phase of the
parent nitrogen-doped superconductor, based on its computed thermodynamic and
dynamic stability, optical properties and XRD pattern. A computational
exploration of the Lu-N-H phase diagram found no ternary phases on the convex
hull at pressures below 10 GPa, the binaries instead dominating, although a
few ternary phases (Lu20H2N17, Lu2H2N, LuH5N2, Lu3H6N, and Lu10HN8) were
within 100 meV/atom of the hull. A number of the identified phases were found
to be derived from either H vacancies or N-doping of LuH2 [36]. Another
computational study did not find any thermodynamically stable Lu-N-H phases at
1 GPa and the highest $T_{c}$ computed for N-doped $Fm\bar{3}m$-LuH3 did not
exceed 30 K [37].

Herein, we present a density functional theory (DFT) investigation of a series
of structures in the Lu-N-H system that were either constructed via
modification of known and theoretical prototype structures, via constrained
and unconstrained crystal structure prediction (CSP) searches, or by a
combination of these two methods. From the results of the unconstrained CSP
runs we obtain a baseline against which to measure the enthalpies of
constructed phases and to compare their properties. From constrained searches
and artificially-constructed structures we begin to understand the motifs that
contribute to dynamic stability at low pressures, and those which do not,
allowing us to narrow the range of possible structures for further
explorations into the Lu-N-H ternary system. The simulated X-ray diffraction
patterns of the optimized phases and calculated equations of states are
compared with available experimental data provided in Reference [31]. The
highest superconducting critical temperature we find – 17 K at 10 kbar – was
obtained for a CaF2-type LuNH phase that was far from thermodynamic stability.

##  II Computational Details

Precise geometry optimizations and electronic structure calculations were
performed using DFT in conjunction with the Perdew-Burke-Ernzerhof (PBE)
functional [38], as implemented in the Vienna ab initio simulation package
(VASP) [39, 40, 41]. The valence electrons of the hydrogen (H $1s^{1}$),
nitrogen (N $2s^{2}2p^{3}$), and lutetium (Lu $5p^{6}5d^{1}6s^{2}$) atoms were
simulated using plane wave basis sets with a cutoff energy of 600 eV. The core
electrons were treated with the projector augmented wave (PAW) method [42].
Detailed tests of the inclusion of the $4f$ electrons on the properties of
select structures, as well as the convergence of the plane wave basis were
performed and representative results are provided in the Supporting
Information. The reciprocal space was sampled using a $\Gamma$-centered
Monkhorst-Pack mesh [43], where the number of divisions along each reciprocal
lattice vector was chosen such that the product of this number with the real-
space lattice constant was 70 Å for density of states calculations and 50 Å
for static calculations. To interrogate the dynamic stability of promising
phases, phonon calculations were performed using the finite difference scheme,
as implemented in the Phonopy software package [44, 45].

The electron-phonon coupling (EPC) calculations were performed using the
Quantum Espresso (QE) package [46, 47] version 7.1 with the PBE functional. A
plane wave basis set with a cutoff energy of 80 Ry was used, along with a
charge density cutoff of 640 Ry for the valence electrons of hydrogen (H
$1s^{1}$), nitrogen (N $2s^{2}2p^{3}$), and lutetium (Lu
$5s^{2}5p^{6}6s^{2}5d^{1}$). The core electrons were treated with the PAW
pseudopotentials generated using the PSLibrary package [48]. The $k$-point and
$q$-point grids were selected to ensure the total electron-phonon-coupling
(EPC) constant, $\lambda$, was converged to within 0.05 at the desired
Gaussian broadening width for each structure, as summarized in the Supporting
Information.

The superconducting critical temperature ($T\textsubscript{c}$) was estimated
using the Allen-Dynes modified McMillan equation [49]:

$T\textsubscript{c}=\frac{\omega\textsubscript{ln}}{1.20}\exp\left[-\frac{1.04(1+\lambda)}{\lambda-\mu^{*}(1+0.62\lambda)}\right],$
(1)

in which the effective Coulomb potential, $\mu^{*}$, was set to 0.1, the
logarithmic average frequency $\omega\textsubscript{ln}$ was obtained by

$\omega\textsubscript{ln}=\exp\left(\frac{2}{\lambda}\int\frac{d\omega}{\omega}\alpha^{2}F(\omega)\ln\omega\right),$
(2)

and the electron phonon coupling constant, $\lambda$, was evaluated by

$\lambda=\int d\omega\alpha^{2}F(\omega)/\omega.$ (3)

The Eliashberg spectral function, $\alpha^{2}F(\omega)$, was obtained from the
QE calculations, and was also used to numerically solve the Eliashberg
equations [50].

The CSP searches were performed using the open-source evolutionary algorithm
(EA) XtalOpt [51, 52, 53] version 12 [54]. The initial generation consisted of
random symmetric structures created by the RandSpg algorithm [55]. Duplicate
structures were identified via the XtalComp algorithm [56] and discarded from
the breeding pool. Constrained XtalOpt searches were performed by determining
the symmetry of the Lu atoms using Pymatgen [57] and only keeping those
structures in the breeding pool that possessed an $Fm\bar{3}m$ symmetry Lu
sublattice. The parameters employed in the XtalOpt searches for the considered
stoichiometries (number of formula units, pressures at which the EA searches
were performed, and constraints employed) are provided in the Supporting
Information.

##  III Results

###  III.1 Known Ambient Pressure Phases

Before we begin our theoretical investigation of novel Lu-N-H combinations
that could be formed at mild pressures, let us review the structures and
properties of the known LuHx and LuN phases. Unlike the high-pressure
superhydrides, which bear little to no resemblance to the hydrides that are
known at ambient conditions, the 1 atm LuN and LuHx phases may provide the key
to the structure of Lu-N-H – or at least very good starting points – stemming
from the relatively low pressures required to stabilize this ternary phase.

At ambient pressure, LuN assumes the rock-salt, or $B1$, structure (Figure
1a), with the Lu atoms in the _fcc_ configuration. A transition to the $B2$ or
CsCl phase has been predicted near 250 GPa [58]. Our PBE calculations, which
likely underestimate the band gap, suggest semiconducting behavior at 1 atm
with a gap of 0.23 eV. In compounds, lutetium typically adopts the +3
oxidation state, and its hydrides can incorporate vacancies or extra hydrogen
atoms that go into the interstitial regions [59]. At 1 atm fluorite (CaF2)
type LuHx is adopted when $x=1.85-2.23$ (Figure 1(b)), usually resulting in a
metallic phase. Increasing the hydrogen content to $x=2.78-3$ yields a
hexagonal semiconducting phase [59]. This $P\bar{3}c1$ LuH3 transitions to a
cubic phase at $\sim$10 GPa (the AlFe3 or $D0_{3}$ structure type, Figure
1(c)), which can be stabilized at ambient via milling [60]. Recently,
superconductivity in $Fm\bar{3}m$-LuH3 was measured with a $T_{c}$ of 12.4 K
at 122 GPa [61].

Figure 1: Prototype Lu-N and Lu-H structures with _fcc_ Lu lattices: (a) NaCl-
type LuN, (b) CaF2-type LuH2, and (c) a high-pressure (hp) phase of LuH3.

To validate the computational settings used, we compared the lattice constants
of the known phases where the Lu atoms are found in the $fcc$ arrangement:
rock salt LuN (4.760 Å [62]) and fluorite type LuH2 (5.033 Å [63]) with those
of the optimized structures. The DFT lattice constants differed by only 0.17%
and 0.28% from experiment, further supporting the choice of our computational
parameters. These known ambient-pressure nitrides and hydrides of lutetium
provide a basis that could be used to build models of the high-$T_{c}$
superconducting phase reported in Ref. [31]. In fact, the similarity of the 1
atm lattice parameters of phase A (5.0289(4) Å) and the (presumably non-
superconducting) compound B (4.7529 Å) with the known dihydride and nitride of
lutetium, respectively, coupled with a comparison of the DFT-optimized unit
cell parameters of several hypothetical and selected partially-doped versions
of the known compounds were used to assign possible compositions [31]. Phase A
was tentatively assigned as LuH3-δNϵ, with partial N substitution onto H sites
in the cubic (high-pressure) LuH3, and phase B as LuN1-δHϵ, an H-doped variant
of rock-salt LuN [31]. On the other hand, a recent theoretical manuscript
proposed that CaF2-type LuH2 is the parent structure of the superconducting
phase, and compound B could be the rock-salt LuH structure, which is
dynamically stable at 0 GPa [35].

###  III.2 Newly Predicted Phases

The structures investigated herein were generated using a variety of
procedures including _ab initio_ CSP techniques, as well as modification of
known phases and compounds predicted using CSP. The advantage of CSP searches
is that they can, freed from structural preconceptions, locate the low-lying
configurations in a potential energy surface, whose complexity here is
heightened by the inclusion of three elements. Such searches can be
unconstrained, purely hunting down the lowest-enthalpy configurations given a
certain stoichiometry. Constraining a search to structures containing a
particular motif will narrow down the possible results, but could also miss
out on even lower-enthalpy alternatives that do not align with the
constraints.

To that end, a combination of both unconstrained and constrained CSP searches
were carried out for the Lu-N-H system using the XtalOpt EA. From the former
we can learn about the structural motifs that yield the most stability, and
comparison with the latter informs us of the enthalpic cost associated with a
specific structural feature. In addition, various structures were made “by
hand” via modification of known prototypes or CSP generated structures that
possess an _fcc_ Lu lattice. As we will soon see, a large structural variety
is present amongst the dynamically stable phases that we found, highlighting
the difficulties inherent in the computational prediction of metastable phases
that could potentially be synthesized.

Semiconductors:

Figure 2: Semiconducting Lu-N-H phases found using unconstrained evolutionary
crystal structure searches and prototype modification (Lu4NH9).

Unconstrained XtalOpt searches for the lowest-enthalpy structures were
performed for the Lu3NH11 and Lu4NH10,11 compositions at both 0 and 3 GPa, as
well as for Lu4NH6 and LuNH2 at 0 GPa. These EA runs located a number of
structurally diverse semiconducting phases with PBE band gaps that ranged from
1.1-2.1 eV; some are shown in Figure 2. A few of the predicted structures,
including $P2_{1}m$ LuNH2, and two Lu3NH11 phases – one with $P1$ symmetry at
0 GPa and one with $Cm$ symmetry at 3 GPa – possessed large empty regions.
$P2_{1}m$ LuNH2 (Figure 2(a)) is, in fact, a fully 2D compound. At 0 GPa $Pc$
Lu4NH11 (Figure 2(b)) was also identified; it consists of layers of trigonal
nets of Lu with H atoms in the resulting hexagonal channels, while the N atoms
are arranged in zigzag chains oriented along the _c_ -axis that weave into the
the Lu network (into the plane of the page).

Two of the structures from unconstrained searches – $P1$ Lu4NH10 (at 0 GPa)
and a second $Pc$ Lu4NH11 structure (at 3 GPa; Figure 2(c)) – possessed Lu
sublattices in slightly distorted _fcc_ arrangements. In $P1$ Lu4NH10, the N
atoms go into some of the sites octahedrally coordinated by Lu, while some H
atoms go into the tetrahedral interstices and the rest are scattered across
the unit cell, resulting in the very low symmetry. For $Pc$-I Lu4NH11 (Figure
2(c)), the N atoms go instead into the tetrahedral interstices and the
hydrogen atoms take the octahedral and most of the remaining tetrahedral
interstices. The _fcc_ Lu lattice is also preserved in a semiconducting $Amm2$
compound with Lu4NH9 stoichiometry (Figure 2(d)), which was produced not by
CSP but by modifying the geometry of the high-pressure AlFe3-type LuH3
compound. Here, H again partially occupies both tetrahedral and octahedral
interstices in _fcc_ Lu, leaving 1/4 of the tetrahedral interstices empty and
with 1/4 of the H atoms filling octahedral interstices being replaced by N.

From these results, it is clear that a variety of geometric motifs can be
found in the low-enthalpy Lu-N-H compounds, highlighting both the difficulty
of honing in on a single structure and the utility of guidance from
experimental data. The unit cell volumes of several of the systems identified
via unconstrained CSP searches were too large for them to stay as candidates
for the putative superconducting phase. Importantly, because all of the
aforementioned compounds were semiconducting it is impossible for any of them
to be superconductors. Our search continues, with inspiration taken from known
experimental phases or CSP searches guided via constraints towards desired
structural features – or both.

Structures from Prototype Modification:

The relatively low pressures needed to stabilize the putative room-temperature
superconducting phase highlight the importance of – and inspiration to be
gleaned from – examining the ambient- and low-pressure compounds formed
between Lu and either N or H. Notably, within most of these, the Lu atoms
adopt the _fcc_ arrangement that has been suggested for the superconducting
phase.

In addition to the ambient pressure $B1$ mononitride, LuN (Figure 1(a)), we
considered a hypothetical rock-salt monohydride, LuH (Figure 3(a)), and
hypothetical zinc-blende (or $B3$) LuN and LuH phases (Figure 3(b,c)). To
explore the potential of a solid solution between the two rock-salt phases
calculations were carried out on the unit cells shown in Figure 3(d). From
these, only LuN0.25H0.75 and LuN0.5H0.5 were dynamically stable at 0 GPa.
Similarly, solid solutions of zinc-blende LuN and LuH were optimized (Figure
1(e)) and from these LuH, LuN0.5H0.5, LuN0.75H0.25 and LuN were 0 GPa
dynamically stable.

N/H substitution into the fluorite-type, or $C1$, LuH2 phase (Figure 1(b)),
yielded another set of potential candidates (Figure 3(f)), with LuN0.5H1.5 and
LuNH being dynamically stable at 0 GPa. LuNH is a half-Heusler-like compound
with equal amounts of N and H occupying the tetrahedral interstices of the Lu
lattice. From the dynamically stable phases identified in this section, $C1$
LuN0.5H1.5 is weakly metallic under PBE-DFT, and thus likely in actuality to
be a non-metal. The rest are metallic. Below, we will compare the pressure-
volume relation calculated for the phases discussed in this section with the
experimental results obtained for compounds A and B, and discuss the
thermodynamic stability, electronic structure and potential for
superconductivity in these prototype-based Lu-N-H phases.

Figure 3: Illustrations of hypothetical (a) rock-salt ($B1$) LuH, and zinc-
blende ($B3$) (b) LuH and (c) LuN phases. (d) Rock-salt and (e) zinc-blende
LuNxH(1-x), and (f) fluorite ($C1$) Lu(NxH(1-x))2 solid solution models that
were considered.

Structures Inspired by Evolutionary Searches:

Figure 4 illustrates a number of 0 GPa dynamically stable, metallic phases
with _fcc_ Lu sublattices that were found in a variety of ways. The
$Fd\bar{3}m$ Lu4NH7 phase (Figure 4(a)) was found in an unconstrained
evolutionary search performed at 1 GPa. It can be constructed from a modified
2$\times$2$\times$2 supercell of CaF2-type LuH2, in which 1/8 of the
tetrahedral interstices of the Lu lattice are occupied by N rather than H. The
distribution of the N atoms throughout the unit cell is in a diamond-like
lattice. In this structure, the octahedral interstices of the Lu lattice are
left empty. This structure belongs to the same family of phases illustrated in
Figure 3(f), representing another N-substituted CaF2-type LuH2 derivative.
However, rather than being derived from prototype modification, it was located
in an XtalOpt search and then served as a template to construct additional
metastable phases. One of these, $Fd\bar{3}m$ Lu2NH5 (Figure 4(b)), was
generated by placing H2 units into some of the empty octahedral interstices of
the Lu lattice, and replacing an additional H atom from Lu4NH7 by N, so that
the N atoms now trace out a _bcc_ network within the structure, leaving H2
molecules lying along only half of the N-N contacts.

Figure 4: Crystal structures of various dynamically stable Lu-N-H phases
obtained from a combination of CSP searches – some constrained – and
subsequent modification. Lu2NH5 and Lu2NH3 had PBE band-gaps of 1.09 and 0.06
eV at 10 kbar.

Another (incomplete, or prematurely terminated) XtalOpt search at 1 GPa
identified $P\bar{4}3m$ Lu4NH6 (Figure 4(c)), which was chosen for further
analysis and modification due to its dynamic stability, and the good match
between its simulated XRD pattern with experiment. Like $Fd\bar{3}m$ Lu4NH7,
$P4\bar{3}m$ Lu4NH6 is similarly a variant of the CaF2-type LuH2 structure, in
which 1/8 of the tetrahedral interstices of the Lu lattice are occupied by N
rather than H and an additional 1/8 of the tetrahedral interstices are left
empty. Rather than the diamond-like distribution of N atoms found in Lu4NH7,
the substituting N atoms and vacancies are arranged in a CsCl-type framework.
Filling the vacancies in Lu4NH6 with N atoms yields the $Pn\bar{3}m$ Lu2NH3
structure (Figure 4(d)).

In the above phases, N atoms were positioned in the tetrahedral interstices of
an _fcc_ Lu framework, whereas in Lu2NH5 the octahedral interstices were
partially occupied by H2 molecular units. In $R3m$ Lu4NH4, which was
identified using an XtalOpt search carried out at 6 GPa where the Lu
sublattice was constrained to maintain the $Fm\bar{3}m$ space group, the N
atoms are not found within the tetrahedral holes but instead lie on 1/4 of the
octahedral holes of the Lu lattice (Figure 4(e)). The N atoms in Lu4NH4 trace
out a simple cubic arrangement, with their positions shifted slightly off of
the center of the surrounding Lu6 octahedra, while the tetrahedral interstices
are half occupied by H and half are left empty. The remaining H atoms can be
grouped into H@H6 vertex-sharing octahedra.

The unconstrained 0 GPa XtalOpt searches that mainly uncovered the
semiconducting compounds shown in Figure 2 also produced the metallic $R3m$
Lu4NH6 phase (Figure 4(f)). In this phase, the Lu-N and Lu-H interactions
become separated, with layers of N@Lu6 octahedra – in essence, slabs of $B1$
LuN – interrupting a CaF2-type packing of Lu and H. Because this phase was
found using an EA search that generated sufficient structures to explore the
potential energy landscape, it was 136.1 meV/atom lower in enthalpy than the
previously discussed $P\bar{4}3m$ Lu4NH6, and at 5 GPa this difference
increased to 175 meV/atom. Perhaps this structure, with N-rich layers
intercalated into a LuH2 matrix, could hint at a strategy for inducing
epitaxial strain on simple LuHn structures, thereby altering their electronic
and mechanical properties from those of their parent.

##  IV Properties: Stability, Equation of States, Electronic Structure,
Superconductivity

Thermodynamics:

The thermodynamic stability of the new structures was investigated by
calculating their formation enthalpies relative to the solid elemental phases
as a function of pressure. The reference phases employed were Lu: $\alpha$-Sm
(0-8 GPa [64]) and the hexagonal phase (9-10 GPa [65]); H2: $P6_{3}/m$ phase
(0-10 GPa [66]); and N2: $\alpha$-N2 phase (0-7 GPa [67]) and $\epsilon$-N2
phase (8-10 GPa [68]). Known experimental phases including fluorite type LuH2,
$B1$ LuN, $P\bar{3}c1$ LuH3, and $P2_{1}3$ NH3 [69] were also considered.

The 0 GPa convex hull shown in Figure 5 illustrates that only the known
structures are thermodynamically stable, and all of the previously discussed
Lu-N-H compounds are thermodynamically unstable within the static lattice
approximation. Up to 10 GPa only the known phases lie on the hull, while all
others lie above it. A structure’s thermodynamic stability can be
characterized by its distance to the convex hull, which is listed in the
Supporting Information, where it is also plotted as a function of pressure
(for all compounds, regardless of their dynamic stability). Herein, we employ
70 meV/atom, which corresponds to the 90th percentile of the DFT-calculated
metastability of all of the known inorganic crystalline materials [70], as a
gauge to identify those structures that could potentially be synthesized. At 0
GPa only five structures – all found using our unconstrained crystal structure
search – fall in this range. From these only $R3m$ Lu4NH6 was metallic.

Figure 5: Convex hull at 0 GPa. Only dynamically stable structures within 300
meV/atom above the hull are shown. If multiple structures exist for the same
stoichiometry, only the most stable structure is listed. Black dots represent
thermodynamically stable phases on the hull, and the colored points are
colored by their distance from the hull in meV/atom. Triangles: structures
generated via evolutionary search. Boxes: structures from prototype
modification. Circles: structures generated by inserting atoms into structures
derived from EA searches.

Let us now turn to the metallic phases with _fcc_ Lu lattices. For the rock-
salt solid-solution family, hydrogen concentrations ranging from 25-100% were
roughly within 150-250 meV/atom from the hull, with $B1$ LuH as the lower
boundary. For the zinc-blende solid-solution system this range expanded to
100-500 meV/atom, with $B3$ LuH corresponding to the lower boundary as well.
Doping fluorite LuH2 causes its energy to explode quickly: 25% nitrogen
content results in an increase of energy by $\sim$200 meV/atom above the
convex hull, which rises to $\sim$550 meV/atom for a 50% composition, and 1.1
eV/atom for 75% nitrogen content. The 0 GPa ternary convex hull plot shows
that most of the low-enthalpy metastable structures are found at the bottom
left hand corner. The reason for this is that these are the only regions where
full unconstrained CSP searches were performed, and the survivor bias makes us
think that this region is where stable structures might appear. It should be
noted, however, that these stoichiometries were chosen because exploratory
calculations suggested their volumes were likely to provide the best match
with the experimental equation of states of compound A. This will be explored
shortly below.

Assuming linear behavior of the enthalpy-pressure relation, we were able to
estimate the pressure where the considered phases may become thermodynamically
stable if the slope of the distance from the hull versus pressure is negative.
This estimate does not take into account the dynamic stability, nor does it
include temperature or effects arising from the zero point motion of the
nuclei. The results suggest that $B1$ LuN-LuH mixtures become favored at high-
pressures: LuH by $\sim$30 GPa, $\sim$50 GPa for LuN0.25H0.75, about 70 GPa
for LuN0.5H0.5, and 80 GPa for LuN0.75H0.25. The higher the hydrogen
concentration, the lower the predicted stabilization pressure, with a lower
boundary of 30 GPa for $B1$ LuH. The slope of the other two solid solutions
considered, $B3$ and $C1$ type, is positive suggesting they will never be
stabilized. Two further phases that could potentially be stabilized within the
megabar range are $R3m$ Lu4NH6 (16 GPa) and $P_{1}$ Lu4NH10 (34 GPa) because
they are very close to the hull. The rest of the structures either possess a
positive slope, or cannot be stabilized until at least 140 GPa.

Equation of States and X-ray Diffraction Patterns:

One of the key experimental observables guiding our choice of stoichiometries
was the pressure-volume relation, or equation of states (EoS), of the majority
phase presented in Reference [31], which was assigned tentatively as an
$Fm\bar{3}m$ structure with a LuH3-δNϵ stoichiometry (or compound A). Above
$\sim$30 kbar a first-order structural phase transition with a $\sim$0.3%
volume discontinuity was observed suggesting that the metal lattice of the
resulting non-superconducting phase distorted to the $Immm$ spacegroup. In
Figure 6 we plot the EoS fits from Reference [31] for the majority phase,
which were obtained for two pressure ranges. Choosing stoichiometries whose
volumes matched well with experiment was initially non-intuitive because the
effective radius of the metal atom changes substantially with its oxidation
state being largest for Lu and smallest for Lu3+.

Figure 6: The DFT calculated pressure-volume relationship or equation of
states (EoS) of the Lu-N-H phases considered in this study. The colored
squares correspond to the specified structures and the open diamond, triangle
and circles to various structures comprising the $B1$, $B3$ and $C1$ solid
solution series (see Figure 3), except for $C1$ LuH2 and $B3$ LuH. The black
lines represent the EoS fitted using the Birch-Murnaghan method for compound A
using data from the pressure ranges $0<P<40$ kbar (solid) and $P>42.7$ kbar
(dashed) [31].

From all of the phases we considered, both fluorite LuH2 and zinc-blende LuH
presented the best match with the experimental data below 40 kbar. At higher
pressures, however, the volume of $B3$ LuH was computed to become
progressively smaller than the measured value for compound A. The good
agreement with $C1$ LuH2, on the other hand, remained up to at least 80 GPa.
At 0 GPa $B3$ LuH was slightly larger than $C1$ LuH2, in-line with the general
notion that the effective radius of Lu+ is larger than that of Lu2+. However,
the volume of $B3$ LuH shrinks much faster (with a slope that is similar to
that of $B1$ LuH) with increasing pressure as compared to that of $C1$ LuH2,
while the volume of cubic LuH3 shrinks at an even slower rate. Thus, the
compressibility in these compounds appears to be dependent upon the repulsion
exhibited between the ionic cores, with a larger number of H- anions resulting
in a higher resistance to compression. Due to its larger ionic radius, N3- is
less compressible than H-. Since the computed EoS of cubic LuH3 has a smaller
slope than the EoS derived from experiment, and introduction of nitrogen will
decrease the slope further, it could be expected that a compound with the
LuH3-δNϵ stoichiometry that was proposed for compound A would not have a slope
that coincides with the experimentally derived EoS.

Because the calculated EoS of fluorite LuH2 across the whole pressure range
yielded the best fit with the experimentally reported EoS, we employed the
quasiharmonic approximation to obtain a temperature-dependent EoS. Fitting the
resulting EoS using the Birch-Murnaghan method at 300 K yielded $V_{0}$
(reference volume at $P=0$) of 31.85 Å3, $K_{0}$ (bulk modulus at $P=0$) of
922.8 kbar, $K_{0}^{\prime}$ ($dK_{0}/dP$ at $P=0$, dimensionless) of 3.7
(data for 0 and 100 K can be found in the SI). This compares well with the
values presented in Ref [31] obtained using fits to data collected below
(above) 40 kbar of 31.74 (31.6) Å3, 886.4(900) kbar, and 4, respectively.

To determine if the structures discussed here could yield XRD patterns similar
to those observed in experiment, their simulated 0 GPa XRD patterns were
generated, as was an XRD pattern for a model $Fm\bar{3}m$ Lu cell whose
lattice constant ($a=5.029$ Å) was in-line with the refined unit cell
suggested for superconducting compound A at 0 GPa (plots are provided in the
Supporting Information). The PyXtal XRD Similarity tool [71] was used to
assess the similarity between the simulated powder XRD patterns of the
proposed structures and that of the model $Fm\bar{3}m$ Lu cell. The strongest
matches came from the experimental phases CaF2-type LuH2 (0.9848), AlFe3-type
LuH3 (0.9316), and from ZnS-type LuH (0.9962) – in-line with the volume of
$B3$ LuH adhering closely to the experimental EoS near 0 GPa. Of the N/H-doped
NaCl, ZnS, and CaF2-type structures, the best XRD matches could be attributed
to the ZnS-based phases, with the NaCl-based phases comparing most poorly. Of
the phases directly obtained from XtalOpt searches or based on modifications
of XtalOpt results, $Fd\bar{3}m$ Lu4NH7 and $R3m$ Lu4NH4 provided the best
matches, although their enthalpies place them well above the convex hull in
the pressure range of interest.

Electronic Structure and Superconductivity:

Superconductivity has been measured in elemental Lu at pressures above
$\sim$100 kbar, with $T_{c}$ rising to $\sim$0.6 K near 160 kbar [72]. Adding
hydrogen and mild pressure does not improve the superconducting properties
much or at all: superconductivity in LuH2 was not observed down to 1.5 K at
pressures as high as 7.7 GPa [33]. These recent experimental results are in
agreement with our computed values at 10 kbar, obtained via the Allen-Dynes
modified McMillan equation, which is thought to be appropriate for phonon-
mediated superconductors whose $\lambda<\sim$1-1.5. As shown in Table 1, we
found that the $T_{c}$ of fluorite-type LuH2 was $\sim$0.1 K, owing to a small
$\omega_{\text{ln}}$ combined with a modest $\lambda=0.29$.

To study the potential for superconductivity in ternary Lu-N-H compounds we
performed EPC calculations for the previously discussed metallic phases that
were dynamically stable at 10 kbar – the pressure at which the maximum $T_{c}$
was observed in Reference [31]. Table 1 shows that though the $T_{c}$s of most
of these phases (with the exception of LuN0.5H1.5) were predicted to surpass
that of $C1$ LuH2, they do not even reach the boiling point of liquid
nitrogen, in agreement with recent theoretical calculations that did not find
any Lu-N-H phases with room temperature superconductivity [37].

Structure | $\lambda$ |  $\omega\textsubscript{ln}$ (K) |  $T\textsubscript{c}$ (K)
---|---|---|---
CaF2-type LuH2 | 0.29 | 302 | 0.09
CaF2-type LuNH | 0.78 | 377 | 16.9 (18.3)
CaF2-type LuN0.5H1.5 | 0.11 | 680 | 0.0
$R3m$ Lu4NH4 | 0.64 | 151 | 4.2
$P\bar{4}3m$ Lu4NH6 | 0.48 | 291 | 2.9
$Fd\bar{3}m$ Lu4NH7 | 0.47 | 435 | 4.2
$R3m$ Lu4NH6 | 0.29 | 306 | 0.12
Table 1: The electron phonon coupling, $\lambda$, logarithmic average
frequency, $\omega\textsubscript{ln}$, and superconducting critical
temperature, $T_{c}$, estimated using the Allen-Dynes modified McMillan
equation with $\mu^{*}=0.1$ at 10 kbar for select Lu-N-H compounds. For LuNH
the numerical solution of the Eliashberg equations was employed to obtain the
value in brackets.

The highest $T_{c}$ compound we found, fluorite type LuNH, can be derived from
LuH2 by replacing 50% of the hydrogen atoms by nitrogen (Figure 3(f)). This
chemical substitution dramatically increased the EPC, placing it in the realm
of the ambient pressure conventional superconductor with the highest confirmed
$T_{c}$, MgB2. However, the larger $\lambda$ of 0.78 was attained at a cost of
the thermodynamic stability: while $C1$ LuH2 fell on the 10 kbar hull, LuNH
was 564 meV/atom above the hull, suggesting it could never be made. The
$T_{c}$ of LuNH ($\sim$17 K) was estimated to be a factor of two smaller than
that of MgB2 with its strong covalent B-B bonds, whose motions, with
frequencies around 600 cm-1, yield an $\omega\textsubscript{ln}$ of 504 cm-1
(or 725 K) [73]. As we shall soon see, in LuNH the EPC are relatively evenly
distributed from the high frequency motions of the hydrogen vibrations to the
very low frequency acoustic modes. Their $\alpha^{2}F$-weighted logarithmic
average yields an $\omega\textsubscript{ln}$ of 257.2 cm-1 (or $\sim$370 K).
Numerical solution of the Eliashberg equations raised the $T_{c}$ of LuNH only
slightly to $\sim$18 K.

Let us examine the electronic structure of CaF2-type LuNH and its
contributions to the EPC to better understand how these factors influence the
$T_{c}$. Replacing H by N in LuH2 increases the density of states (DOS) at the
Fermi level ($E_{F}$) by around 50% from 0.019 states/eV/Å3 to 0.027
states/eV/Å3, concomitantly increasing the $T_{c}$. As shown in Figure 7, the
major contributions to the DOS at $E_{F}$ are the H $1s$ and N $2p$ states,
with a negligible amount from the metal, indicative of a +3 oxidation state.
The primitive cell of LuNH contains one formula unit, and as a result its
conduction band is half-filled. The reaction
$\text{LuNH}+\frac{1}{2}\text{H}_{2}\rightarrow\text{LuN}+\text{H}_{2}$ is
exothermic by 400 meV/atom; we would therefore expect the products of this
reaction to be found in a CSP search for unit cells whose sizes approach
infinity.

Figure 7: PBE band structure and projected densities of states of fluorite-
type LuNH at 10 kbar.

Pivoting to the phonon band structure in Figure 8, we observe that the large
differences in the mass between the three elements splits their bands nicely
into three regions. The vibrational modes of lutetium are mainly below 140
cm-1 (acoustic region), nitrogen are between 380-470 cm-1 and hydrogen above
660 cm-1. It should be noted that due to the extremely heavy mass of lutetium
versus nitrogen and hydrogen (175 vs. 14 and 1 a.u.), lutetium moves roughly
ten times slower than the hydrogen, and four times slower than the nitrogen.
As a result, the atomic displacements of the nitrogen and hydrogen atoms along
the low frequency acoustic modes are still significant.

Because of the separation of these vibrational modes, their contribution to
the total EPC can be obtained: motions from the acoustic modes contribute 41%,
23% from the nitrogen active region, and 35% from the hydrogen active region.
The largest Lu-based contribution originates from the lower two acoustic
phonon branches around the middle of the $\Gamma$-$K$ path, and also around
the $L$ point. Visualization of these motions show they result in the
formation of N-Lu-H molecular fragments and a hexagonal-like Lu lattice. In
the nitrogen-active region the largest EPC is found at the $\Gamma$ point,
resulting from the nitrogen atoms approaching lutetium to form N-H motifs. In
the hydrogen-active region, the entire bands exhibit moderate EPC, especially
at several points where the modes are softened; visualization shows that these
correspond to the motion of hydrogen atoms closer to lutetium to form H-Lu
units.

Figure 8: Phonon dispersion curve and projected EPC constant
($\lambda_{\mathbf{q}\nu}$). Blue color indicates $\lambda_{\mathbf{q}\nu}$
approaches 0, and red indicates $\lambda_{\mathbf{q}\nu}$ approaches the
maximum value of 0.36. The atom projected phonon density of states is
illustrated, along with the total $\lambda$ and the integral of
$\lambda(\omega))$ separated into regions comprising the Lu, N and H-based
modes.

The $\omega\textsubscript{ln}$ of our Lu-N-H compounds ranged from $\sim$150 K
($R3m$ Lu4NH4) to 680 K (CaF2-type LuN0.5H1.5). The absence of high frequency
vibrations in these compounds, resulting from the low pressure and the absence
of covalent bonds, suggests that higher $\omega\textsubscript{ln}$ are
unlikely to be found in other Lu-N-H compounds at 10 kbar with _fcc_ Lu
lattices. Generally speaking, the $\omega\textsubscript{ln}$ calculated for
hydrogen and the high-$T_{c}$ hydrides at extreme pressures is significantly
higher, with values of 1200-1800 K not being uncommon. For those hydrides
where comparable $\omega\textsubscript{ln}$ values have been calculated, room
temperature superconductivity has only been predicted in phases with a very
large EPC (e.g. $Fmmm$ ThH18 at 400 GPa, $\omega\textsubscript{ln}=$568 K,
$\lambda=$3.39, $T_{c}=$296 K) [74]. Therefore, we speculate that similar EPC
constants are required for a Lu-N-H compound to be superconducting near room
temperature, provided the mechanism is phonon-mediated.

##  V Conclusions

Density functional theory calculations were performed to explore Lu-N-H
containing compounds that could be (meta)stable in a pressure range of about
0-100 kbar (10 GPa). The computations were biased towards systems where the Lu
atoms adopt an _fcc_ arrangement, because it was recently suggested that a
compound with this structural feature could be responsible for the near-
ambient superconducting critical temperature, $T_{c}$, of 294 K reported at 10
kbar [31]. Based on the results of our calculations we conclude that:

  * •

The Lu-N-H potential energy landscape, within the static lattice approximation
and neglecting quantum nuclear and anharmonic effects, contains many local
minima with _fcc_ Lu lattices. Other geometries, not explicitly considered
here, could be generated via altering the N/H ratio of the solid-solution
prototypes we discuss. Which of these structures are synthesizable, and which
are kinetically and/or thermally stable and relatively chemically inert is
currently unknown.

  * •

None of the ternary compounds studied here are thermodynamically stable (e.g.
they do not lie on the convex hull) up to 10 GPa at 0 K. Only the known
binaries, LuH2, LuH3, LuN and NH3, comprise the convex hull. Thermal effects
and the role of configurational entropy on the thermodynamic stability are not
known.

  * •

From all of the phases considered here, the one whose equation of states (EoS)
had the closest match with the fits to experimental data obtained for compound
A was fluorite-type LuH2, with errors smaller than 0.3% up to 80 kbar. EoS
calculations on model compounds suggest that the previously proposed formula
for compound A, LuH3-δNϵ, would not have the same slope as what was observed
experimentally.

  * •

XRD similarity indices for the compounds studied here compared to a pure _fcc_
Lu lattice with the experimental lattice constant indicated a fair match at 0
GPa for binaries LuH2, LuH3, and ZnS-type LuH, N-substituted ZnS-type LuH,
$Fd\bar{3}m$ Lu4NH7 and $R3m$ Lu4NH4.

  * •

Many, though not all, of the investigated phases exhibited metallic behavior,
and their density of states at the Fermi level (DOS at $E_{\text{F}}$) varied
greatly. For example, the main contributions to the DOS at $E_{\text{F}}$ for
$B1$ and $B3$ LuH were the Lu $d$ states; for LuH2 the DOS at $E_{\text{F}}$
was very small and mainly lutetium $p$-like, and for LuNH the main components
arose from hydrogen $s$ and nitrogen $p$ states.

  * •

The logarithmic average frequency, $\omega_{\text{ln}}$, of the Lu-N-H
compounds whose superconducting properties we studied ranged from
$\sim$150-680 K, and the electron phonon coupling constants, $\lambda$, varied
between 0.1-0.8. Assuming conventional superconductivity, the $T_{c}$s of such
compounds can be estimated using the modified McMillan Allen Dynes equation.
Under this approximation, and with $\mu^{*}=0.1$ at 10 kbar we obtain a
$T_{c}$ of 0.09 K for fluorite LuH2. The highest $T_{c}$ compound we found was
fluorite-type LuNH with a $T_{c}$ of 17 K.

Though we have not uncovered an Lu-N-H containing phase with a superconducting
critical temperature near what was recently reported in Reference [31], we
believe our computations shed light on the structures that contain these
elemental combinations at mild pressures. Future work will ascertain if our
choice of standard DFT parameters (gradient corrected exchange functional,
neglect of spin polarization and strong electron correlations, and inclusion
of $f$ electrons in the core) impact our conclusions. Our work also highlights
the complexity inherent in the computational search for phases that may be
metastable with desired structural and property characteristics in multi-
element _ab initio_ (or even machine-learning-assisted) crystal structure
prediction.

##  VI Acknowledgments

We are grateful to R. Dias for sharing experimental data, as well as G.W.
Collins and R.J. Hemley for useful discussions. K.H. acknowledges the
Chicago/DOE Alliance Center under Cooperative Agreement Grant No. DE-
NA0003975, and N.G. the U.S. National Science Foundation (DMR-2132491). This
material is based upon work supported by the U.S. Department of Energy, Office
of Science, Fusion Energy Sciences funding the award entitled High Energy
Density Quantum Matter under Award Number DE-SC0020340. Computations were
carried out at the Center for Computational Research at the University at
Buffalo (http://hdl.handle.net/10477/79221).

## References

  * Onnes [1911] H. K. Onnes, Proc. K. Ned. Akad. Wet. B 14, 113 (1911).
  * Wu  _et al._ [1987] M. K. Wu, J. R. Ashburn, C. J. Torng, P. H. Hor, R. L. Meng, L. Gao, Z. J. Huang, Y. Q. Wang, and C. W. Chu, Phys. Rev. Lett. 58, 908 (1987).
  * Sheng and Hermann [1988] Z. Z. Sheng and A. M. Hermann, Nature 332, 138 (1988).
  * Drozdov  _et al._ [2015] A. P. Drozdov, M. I. Eremets, I. A. Troyan, V. Ksenofontov, and S. I. Shylin, Nature 525, 73 (2015).
  * Osmond  _et al._ [2022] I. Osmond, O. Moulding, S. Cross, T. Muramatsu, A. Brooks, O. Lord, T. Fedotenko, J. Buhot, and S. Friedemann, Phys. Rev. B 105, L220502 (2022).
  * Minkov  _et al._ [2020] V. S. Minkov, V. B. Prakapenka, E. Greenberg, and M. I. Eremets, Angew. Chem. Int. Ed. 59, 18970 (2020).
  * Ashcroft [2004] N. W. Ashcroft, Phys. Rev. Lett. 92, 187002 (2004).
  * Hilleke and Zurek [2022] K. P. Hilleke and E. Zurek, J. Appl. Phys. 131, 070901 (2022).
  * Hilleke  _et al._ [2022] K. P. Hilleke, T. Bi, and E. Zurek, J. Appl. Phys. 128, 441 (2022).
  * Zurek and Bi [2019] E. Zurek and T. Bi, J. Chem. Phys. 150, 050901 (1 (2019).
  * Ma  _et al._ [2022] L. Ma, K. Wang, Y. Xie, X. Yang, Y. Wang, M. Zhou, H. Liu, X. Yu, Y. Zhao, H. Wang, G. Liu, and Y. Ma, Phys. Rev. Lett. 128, 167001 (2022).
  * Li  _et al._ [2022] Z. Li, C. Zhang, X. Wang, S. Zhang, Y. Jia, S. Feng, K. Lu, J. Zhao, J. Zhang, B. Min, Y. Long, R. Yu, L. Wang, M. Ye, Z. Zhang, V. Prakapenka, S. Chariton, P. A. Ginsberg, J. Bass, S. Yuan, H. Liu, and C. Jin, Nature Commun. 13, 2863 (2022).
  * Somayazulu  _et al._ [2019] M. Somayazulu, M. Ahart, A. K. Mishra, Z. M. Geballe, M. Baldini, Y. Meng, V. Struzhkin, and R. J. Hemley, Phys. Rev. Lett. 122, 027001 (2019).
  * Drozdov  _et al._ [2019] A. Drozdov, P. Kong, V. Minkov, S. Besedin, M. Kuzovnikov, S. Mozaffari, L. Balicas, F. Balakirev, D. Graf, V. Prakapenka, E. Greenberg, D. Knyazev, M. Tkacz, and M. Eremets, Nature 569, 528 (2019).
  * Snider  _et al._ [2021] E. Snider, N. Dasenbrock-Gammon, R. McBride, X. Wang, N. Meyers, K. V. Lawler, E. Zurek, A. Salamat, and R. P. Dias, Phys. Rev. Lett. 126, 117003 (2021).
  * Troyan  _et al._ [2021] I. A. Troyan, D. V. Semenok, A. G. Kvashnin, A. V. Sadakov, O. A. Sobolevskiy, V. M. Pudalov, A. G. Ivanova, V. B. Prakapenka, E. Greenberg, A. G. Gavriliuk, I. S. Lyubutin, V. V. Struzhkin, A. Bergara, I. Errea, R. Bianco, M. Calandra, F. Mauri, L. Monacelli, R. Akashi, and A. R. Oganov, Adv. Mater. 33, 2006832 (1 (2021).
  * Semenok  _et al._ [2021a] D. V. Semenok, I. A. Troyan, A. G. Ivanova, A. G. Kvashnin, I. A. Kruglov, M. Hanfland, A. V. Sadakov, O. A. Sobolevskiy, K. S. Pervakov, I. S. Lyubutin, K. V. Glazyrin, N. Giordano, D. N. Karimov, A. L. Vasiliev, R. Akashi, V. M. Pudalov, and A. R. Oganov, Mater. Today 48, 18 (2021a).
  * Song  _et al._ [2021a] P. Song, Z. Hou, P. Baptista de Castro, K. Nakano, K. Hongo, Y. Takano, and R. Maezono, Chem. Mater. 33, 9501 (2021a).
  * Semenok  _et al._ [2020] D. V. Semenok, I. A. Kruglov, I. A. Savkin, A. G. Kvashnin, and A. R. Oganov, Current Opinion in Solid State and Materials Science 24, 100808 (2020).
  * Peng  _et al._ [2017] F. Peng, Y. Sun, C. J. Pickard, R. J. Needs, Q. Wu, and Y. Ma, Phys. Rev. Lett. 119, 107001 (2017).
  * Sun  _et al._ [2020] W. Sun, X. Kuang, H. D. J. Keen, C. Lu, and A. Hermann, Phys. Rev. B 102, 144524 (2020).
  * Song  _et al._ [2021b] H. Song, Z. Zhang, T. Cui, C. J. Pickard, V. Z. Kresin, and D. Duan, Chin. Phys. Lett. 38, 107401 (2021b).
  * Li  _et al._ [2023] Z. Li, X. He, C. Zhang, K. Lu, B. Bin, J. Zhang, S. Zhang, J. Zhao, L. Shi, S. Feng, X. Wang, Y. Peng, R. Yu, L. Wang, Y. Li, J. Bass, V. Prakapenka, S. Chariton, H. Liu, and C. Jin, arXiv:05117 (2023).
  * Laniel  _et al._ [2022] D. Laniel, F. Trybel, B. Winkler, F. Knoop, T. Fedotenko, S. Khandarkhaeva, A. Aslandukova, T. Meier, S. Chariton, K. Glazyrin, V. Milman, V. Prakapenka, I. A. Abrikosov, L. Dubrovinsky, and N. Dubrovinskaia, Nature Commun. 13, 6987 (2022).
  * Peña Alvarez  _et al._ [2021] M. Peña Alvarez, J. Binns, M. Martinez-Canales, B. Monserrat, G. J. Ackland, P. Dalladay-Simpson, R. T. Howie, C. J. Pickard, and E. Gregoryanz, J. Phys. Chem. Lett. 12, 4910 (2021).
  * Semenok  _et al._ [2021b] D. V. Semenok, D. Zhou, A. G. Kvashnin, X. Huang, M. Galasso, I. A. Kruglov, A. G. Ivanova, A. G. Gavriliuk, W. Chen, N. V. Tkachenko, A. I. Boldyrev, I. Troyan, A. R. Oganov, and T. Cui, J. Phys. Chem. Lett. 12, 32 (2021b).
  * Di Cataldo  _et al._ [2021] S. Di Cataldo, C. Heil, W. von der Linden, and L. Boeri, Phys. Rev. B 104, L020511 (2021).
  * Durajski and Szczȩśniak [2021] A. P. Durajski and R. Szczȩśniak, Phys. Chem. Chem. Phys. 23, 25070 (2021).
  * Zhao  _et al._ [2022] Y. Zhao, X. Zhang, X. Li, S. Ding, Y. Liu, and G. Yang, J. Mater. Chem. C 39, 14626 (2022).
  * Geng  _et al._ [2023] N. Geng, K. P. Hilleke, L. Zhu, X. Wang, T. A. Strobel, and E. Zurek, J. Am. Chem. Soc. 145, 1696 (2023).
  * Dasenbrock-Gammon  _et al._ [2023] N. Dasenbrock-Gammon, E. Snider, R. McBride, H. Pasan, D. Durkee, N. Khalvashi-Sutter, S. Munasinghe, S. Dissanayake, K. V. Lawler, A. Salamat, and R. P. Dias, Nature 615, 244 (2023).
  * Deng  _et al._ [2021] L. Deng, T. Bontke, R. Dahal, Y. Xie, B. Gao, X. Li, K. Yin, M. Gooch, D. Rolston, T. Chen, Z. Wu, Y. Ma, P. Dai, and C. Chu, PNAS 118, e21089381138 (2021).
  * Shan  _et al._ [2023] P. Shan, N. Wang, X. Zheng, Q. Qiu, Y. Peng, and J. Cheng, arXiv:076718 (2023).
  * Ming  _et al._ [2023] X. Ming, Y. J. Zhang, X. Zhu, Q. Li, C. He, Y. Liu, B. Zheng, H. Yang, and H. H. Wen, arXiv:08759 (2023).
  * Liu  _et al._ [2023] M. Liu, X. Liu, J. Li, J. Liu, Y. Sun, X.-Q. Chen, and P. Liu, arXiv:06554 (2023).
  * Xie  _et al._ [2023] F. Xie, T. Lu, Z. Yu, Y. Wang, Z. Wang, S. Meng, and M. Liu, arXiv:11683 (2023).
  * Huo  _et al._ [2023] Z. Huo, D. Duan, T. Ma, Q. Jiang, Z. Zhang, D. An, F. Tian, and T. Cui, arXiv:12575 (2023).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Kresse and Furthmüller [1996a] G. Kresse and J. Furthmüller, Comput. Mater. Sci. 6, 15 (1996a).
  * Kresse and Furthmüller [1996b] G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996b).
  * Kresse and Joubert [1999] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
  * Blöchl [1994] P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Monkhorst and Pack [1976] H. J. Monkhorst and J. D. Pack, Phys. Rev. B 13, 5188 (1976).
  * Togo and Tanaka [2015] A. Togo and I. Tanaka, Scr. Mater. 108, 1 (2015).
  * Togo [2023] A. Togo, J. Phys. Soc. Jpn. 92, 012001 (2023).
  * Giannozzi  _et al._ [2009] P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21, 395502 (2009).
  * Giannozzi  _et al._ [2017] P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. B. Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. D. Corso, S. de Gironcoli, P. Delugas, R. A. DiStasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H.-Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H.-V. Nguyen, A. O. de-la Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics: Condensed Matter 29, 465901 (2017).
  * Dal Corso [2014] A. Dal Corso, Comput. Mater. Sci. 95, 337 (2014).
  * Allen and Dynes [1975] P. B. Allen and R. C. Dynes, Phys. Rev. B 12, 905 (1975).
  * Eliashberg [1960] G. Eliashberg, Sov. Phys. JETP 11, 696 (1960).
  * Lonie and Zurek [2011] D. C. Lonie and E. Zurek, Comput. Phys. Commun. 182, 372 (2011).
  * [52] Http://xtalopt.github.io/.
  * Falls  _et al._ [2021] Z. Falls, P. Avery, X. Wang, K. P. Hilleke, and E. Zurek, J. Phys. Chem. C 125, 1601 (2021).
  * Avery  _et al._ [2019] P. Avery, C. Toher, S. Curtarolo, and E. Zurek, Comput. Phys. Commun. 237, 274 (2019).
  * Avery and Zurek [2017] P. Avery and E. Zurek, Comput. Phys. Commun. 213, 208 (2017).
  * Lonie and Zurek [2012] D. C. Lonie and E. Zurek, Comput. Phys. Commun. 183, 690 (2012).
  * Ong  _et al._ [2013] S. P. Ong, W. D. Richards, A. Jain, G. Hautier, M. Kocher, S. Cholia, D. Gunte r, V. L. Chevrier, K. A. Persson, and G. Ceder, Computational Materials Science 68, 314 (2013).
  * Singh and Verma [2015] S. K. Singh and U. Verma, in  _J. Phys.: Conf. Ser._ , Vol. 640 (IOP Publishing, 2015) p. 012029.
  * Libowitz [1972] G. G. Libowitz, Berichte der Bunsengesellschaft für physikalische Chemie 76, 837 (1972).
  * Kataoka  _et al._ [2021] R. Kataoka, T. Kojima, K. Tada, M. Kitta, N. Takeichi, K. Sakaki, M. Nozaki, T. Kimura, and A. Kamegawa, Materialia 15, 100956 (2021).
  * Shao  _et al._ [2021] M. Shao, S. Chen, W. Chen, K. Zhang, X. Huang, and T. Cui, Inorg. Chem. 60, 15330 (2021).
  * Pebler and Wallace [1962] A. Pebler and W. Wallace, J. Phys. Chem. 66, 148 (1962).
  * Olcese [1979] G. L. Olcese, J. Phys. F: Met. Phys. 9, 569 (1979).
  * Liu [1975] L.-g. Liu, J. Phys. Chem. Solids 36, 31 (1975).
  * Spedding  _et al._ [1961] F. Spedding, J. Hanak, and A. Daane, J. Less Common Metals 3, 110 (1961).
  * Pickard and Needs [2007] C. J. Pickard and R. J. Needs, Nat. Phys. 3, 473 (2007).
  * Donohue [1961] J. Donohue, Acta Crystallographica 14, 1000 (1961).
  * Mills  _et al._ [1986] R. Mills, B. Olinger, and D. Cromer, J. Chem. Phys. 84, 2837 (1986).
  * Pickard and Needs [2008] C. J. Pickard and R. Needs, Nat. Mater. 7, 775 (2008).
  * Sun  _et al._ [2016] W. Sun, S. T. Dacek, S. P. Ong, G. Hautier, A. Jain, W. D. Richards, A. C. Gamst, K. A. Persson, and G. Ceder, Sci. Adv. 2, e1600225 (2016).
  * Fredericks  _et al._ [2021] S. Fredericks, K. Parrish, D. Sayre, and Q. Zhu, Computer Physics Communications 261, 107810 (2021).
  * Wittig  _et al._ [1974] J. Wittig, C. Probst, and W. Wiedemann, Low Temperature Physics-LT 13: Volume 3: Superconductivity , 490 (1974).
  * Kong  _et al._ [2001] Y. Kong, O. Dolgov, O. Jepsen, and O. Andersen, Physical Review B 64, 020501 (2001).
  * Zhong  _et al._ [2022] X. Zhong, Y. Sun, T. Iitaka, M. Xu, H. Liu, R. J. Hemley, C. Chen, and Y. Ma, J. Am. Chem. Soc. 144, 13394 (2022).