# Broad Band Mott Localization is all you need for Hot Superconductivity:
Atom Mott Insulator Theory for Cu-Pb Apatite

G. Baskaran  The Institute of Mathematical Sciences, C.I.T. Campus, Chennai
600 113
Department of Physics and QuCenDiEM, Indian Institute of Technology, Chennai
600 036, India &
Perimeter Institute for Theoretical Physics, Waterloo, ON, Canada

###### Abstract

A hypothetical non-dimerized Cu chain in equilibrium is a spin-$\frac{1}{2}$
atom Mott insulator (AMI), eventhough its band width is high $\sim$ 10 eV.
This RVB reservoir has a large exchange coupling J $\sim$ 2 eV. This idea of,
broad band Mott localization was used by us in our earlier works, including
prediction of high Tc superconductivity in doped graphene, silicene and a
theory for hot superconductivity reported in Ag-Au nanostructures (TP 2008).
In the present work we identify possible random AMI subsystems in Cu-$\rm
Pb_{10}(PO_{4})_{6}O$ and develop a model for reported hot superconductivity
(LKK 2023). In apatite structure, network of interstitial columnar spaces run
parallel to c-axis and ab-plane. They accomodate excess copper, as $\rm
Cu^{0}$ clusters, chains and planar segments. They are our emergent AMI’s.
Electron transfer from AMI’s to insulating host, generates strong local
superconducting correlation, via phyics of doped Mott insulator. Josephson
coupling between doped AMI’s, establishes hot superconductivity. A major
Challenge to superconducting order in real material is competing insulating
phases - valence bond solid (spin-Peirels)-lattice distortions etc. AMI theory
points to ways of making the elusive superconductivity palpable. We recommend
exploration of hot superconductivity in the rich world of minerals and
insulators, via metal atom inclusion.

## Introduction

A milestone discovery of high Tc superconductivity in cuprates, by Bednorz and
Muller [1] in 1986 heralded serious attempts to discover new superconductors,
possibly with higher Tc’s. Resonating valence bond theory of high Tc
superconductivity, a response to understand cuprates, pioneered by Anderson
[2, 3, 4, 5, 6, 7] created a fertile ground for new physics, variety of
quantum spin liquids [8, 9] etc. Doped spin-$\frac{1}{2}$ Mott insulators
became a playground for high Tc superconductivity and offered new hope,
inspite of competing orders [10].

Announcement of hot superconductivity in certain hydrides [11, 12], Ag-Au
nanostrucutes [13] and more recently, Cu-$\rm Pb_{10}(PO_{4})_{6}O$ [14] has
excited the community. Historically there have been claims of hot
superconductivity, named unstable and elusive superconductors[15].

In the past, inspired by possibility of broad band Mott localization we have
developed theories of superconductivity in doped graphene-graphite [16, 17,
18], silicene [19], B doped diamond [20], K3 p-Terphenyl [21], Hydride [22]
and Ag-Au nanostructures [24, 23]. Here atom Mott insulators (AMI)and near
Mott localization in low dimensional monoatom assemblies play special role.

In what follows, we propose a model for possibile hot superconductivity in
Cu-$\rm Pb_{10}(PO_{4})_{6}O$ system, based on our AMI theory [23, 24] applied
to reported superconductivity in Ag-Au nanostructures. Our arguments are
rather qualitative. However, they enumerate rich possibilities, likely to have
implications for future studies.

In $\rm Pb_{10}(PO_{4})_{6}O$ , interstitial sites of Pb octahedral columns
are half empty. In addition, a network of c-axis interstitial columns and ab-
plane interstital columns are present in the ionic-covalent insulating matrix.
Our hypothesis is that excess Cu gets accommodated as a random network of
weakly connected clusters, chains and 2D patches of $\rm Cu^{0}$ in the
network of interstitial columns. These are the weakly coupled random network
of emergent AMI’s

Difference in electronegativity induces electron transfer from Mott insulating
$\rm Cu^{0}$ AMI subsystem to the host insulator. It is a good assumption that
transferred electrons get localized in the random envoronment of the disturbed
host. A random network of hole doped Mott insulator emerges. In doped Mott
insulator, resonating singlets generate strong Cooper pair correlation. In the
pair tunneling mechanism of superconductivity [25], because of spin-charge
decoupling present in the Hubbard chains, single electron tunneling between
chains is blocked at low energy scales. However, pair tunneling is not
blocked. This Josephson type pair tunneling selfconsistently builds a 3D
superconducting state[25, 26].

This weakly coupled doped reference Mott insulator, a coupled t-J model of
chains and clusters explains observed hot superconductivity, for our estimates
of model parameter for Cu-Pb apatite. A major Challenge to our mechanism of
superconductivity is competing (insulating) orders such as valence bond solid
(spin-Peirels) and structural distortions, within AMI subsystem.

Interestingly anomalous diamagnetism and ferromagnetism at room temperatures
have been reported in monoatomic Cu, Ag and Au nonosystems [27, 28, 29, 30,
31, 32, 33] This underscores importance of electron electron Coulomb
interaction effects in broad band monovalent atom systems in low dimensions.

In the present paper we focuss on possible doped Mott insulator physics that
has proven validity in the context of single band system such as cuprates and
other strongly correlated superconductors, consistent with our earlier
proposal of ‘Five fold way to new superconductors’ [10] and emergent Mott
insulators in effectively atomic Mott systems [16, 18, 19, 20, 21].

Our article is organized as follows. First we present arguments for how
isolated bare narrow bands (for example in the $\rm Cu^{2+}$ substitution
scenerio) will not give us scale of Tc reported and a need for broad band Mott
localization. After a summary of phenomenology of superconducting Cu-Pb
Apatite, we present a structural-electronic model for Cu subsystem in Cu-Pb
Apatite. Emergence of weakly coupled effectively Mott insulating cluster,
chains and patches of Cu0 AMI’s is presented (4s orbital are weakly hybridized
with orbitals of surrounding ions). Then we present inter cluster, chain, 2D
patches pair tunneling mechanism, coupled t-J chain model relevant for weakly
coupled doped Mott insulating chains and estimates of superconducting Tc’s.
Obstacles to reach hot superconductivity, via competing orders prevalent in
low dimensional Mott insulating systems and electron lattice coupling is
discussed next.

## Why Broad Band Mott Localization ?

Many mechanisms for high Tc superconductivity exists. One of them is the doped
spin-$\frac{1}{2}$ Mott insulator, RVB mechanism. In this route an effective
single band close to half filling plays a key role. Known doped Mott
insulators, for example high Tc cuprates have a bare band width of about 2 eV
and singlet exchange coupling is about 0.15 eV. Maximum Tc possible in cuprate
system (from theory point of view) is about 200 K - it is never realized in
practice, because of unavoidable competing orders. In organic superconductors
(which could be viewed as self doped Mott insulators) bands are about 0.5 eV
wide; experimentlly seen maximum Tc is about 12 K.

In cuprates the physics is not one of emergence of heavy quasiparticles at the
Fermi leve, but that of celebrated spin-charge decoupling, failure of Fermi
liquid state and absence of Fermi liquid quasi particles at the Fermi level.
It is interesting to recall emergence of heavy effective mass of Fermi liquid
quaiparticles and an emergent flat band at the Fermi level in heavy Fermions.
Here bare band width of conduction electrons is very wide. However, Kondo
exchange with localized spins generates, at low energy scale an effective
narrow quasi particle band (built out of Kondo resonances) at the Fermi level.
This strongly correlated state has superconductivity. However, scale of Tc is
limited by Kondo temperature in an exponential fashion.

If one assumes, as the authors of LK-99 have suggested, replacement of one
chain of Pb2+ by $\rm Cu^{2+}$ , one should get a narrow and isolated bare
half filled band at the Fermi level. This is because of large Cu-Cu
separation, and an unfriendly envoronment (as explained below)for $\rm
Cu^{2+}$ . Conequently, as recent ab-initio calculations [34] show, one gets
isolated bare narrow bands. Correlation effects will further narrow the band.
When we estimate Tc, using doped Mott insulator physics, we get a Tc not
exceeding 10K. Even this Tc will be reduced further by disorder effects and
singlet localization via lattice distortion.

So we are forced to think about an alternative. Here, broad band Mott
localization comes in handy and could push Tc’s to the scales observed.

## Superconducting Phenomenology of LK99

$\rm Pb_{10}(PO_{4})_{6}O$ , an ionic insulator is a member of a large family
of apatite minerals, with a nominal charge state
Pb${}^{2+}_{10}$(PO${}_{4})^{3-}_{6}$O2-. The authors [14] state that
replacement of about a tenth of Pb2+ ions by Cu2+ ions results in an insulator
to metal transition. They report Meissner signal, vanishing resistivity,
levitation. and Tc exceeding 400 K. It is important to point out that their
structural analysis using XRays, shows a minor phase of Cu2S. Defect
concentrations, nature of structural irregularies, size of crystallites are
not reported.

They also sketch a mechanism of superconductivity based on substitution of
$\rm Pb^{2+}$ by $\rm Cu^{2+}$ in the c-axis chains of Pb. Very recent ab-
initio calculation[34] shows presence of well isolated narrow 3d band of $\rm
Cu^{2+}$ . In our estimate the entire isolated band is too narrow and could
support only very low Tc superconductivity, via Mott localization physics, as
briefly explained in the next section.

Recent works [35, 36] report succesfull synthesis of LK99; their search for
superconductivity continues.

## A structural model and reference electronic state for LK99

What are the charge state of added Cu in LK-99 ? In general, local solid state
and quantum chemistry and energetics dictate stable and possible metastable
charge states of added atom. Cu is known to become Cu2+ ion, when there is
covalent mixing with oxygen atoms, as in copper complexes, solid CuSO4, CuO
and cuprates.

Cu2+ in principle can replace Pb2+ in the Pb chain. However, the trigonal
prismatic local envoronment provided by corners of 6 surround phosphate
tetrahedra, is not conducive for Cu2+ ion replacement. Ionization energy of 4s
electron of Cu is high, about 7.7 eV. It is also likely that lone pair
electrons of Pb2+ has an important role in not selecting Cu2+ state. When one
replaces $\rm Pb^{2+}$ ions by $\rm Cu^{2+}$ , according to very recent ab-
initio studies [34] one gets narrow and well isolated Cu 3d band at the Fermi
velvel (even without correlation effects). As mentioned earlier
superconducting Tc’s given by electronic mechanism in this bare and isolated
narrow band is very small.

In what follows we present possibility of excess Cu atoms getting accommodated
as clusters, chains and 2D patches of $\rm Cu^{0}$ and becoming a source of
atom Mott insulator (AMI) and broad band Mott localization. To appreciate
this, we discuss the structure first. Apatite minerals have complex structures
and large unit cells. Nature has used millions of years to arrive at a pretty
periodic structure in $\rm Pb_{10}(PO_{4})_{6}O$ .

Pairs of $\rm(PO_{4})^{3-}$ tetrahedra stack to form columns, which run
parallel to c-axis and form a Kagome lattice. $\rm Pb^{2+}$ ions create c-axis
chains and form a hexagonal lattice. Pb${}^{2+}_{6}$ octahedra share faces and
create c-axis columns of triangular lattice. Half of the interstitial sies of
the octahedra are occupied by $\rm O^{2-}$ ions. (Triangular and Kagome
lattices are formed respectively by centers of elementary hexagon and centers
of nearest neighbor bond of honeycomb lattice). Asymmetric pairing of
Phosphate tetrahedra adds possibility of chirality orders in apatite.

The half empty interstitials in Pb octahedral columns in real systems are
likely to have disordered O2- occupantion, leaving pieces of empty
interstitial columns. In addition network of c-axis interstitial columns and
ab-plane interstital columns are present in the ionic-covalent netowrk.
Further minerals, including apatites with tetrahedral moieties are known to be
good ionic (for example O2-) conductors [37]. Since ionic radii of O2- ion
$\rm Cu^{0}$ are very close, we expect our proposed channel network to
accomodate $\rm Cu^{0}$ atoms.

Our hypothesis is that excess Cu gets accommodated as a random network of
clusters, chains and 2D patches of $\rm Cu^{0}$ . In what follows we discuss
how the $\rm Cu^{0}$ network becomes reservoirs of resonating and strong
singlets, via broad band Mott localization.

## Emergence of singlet reservoirs in $\rm Cu^{0}$ subsystem

Monovalent metals are special - they are half filled, very broad tight binding
bands made of s like Wannier orbitals. Copper in bulk is a broad band 4s
metal, well described by a half filled band of repulsive Hubbard model.
However, 3 dimensionality makes Hubbard repulsive U irrelevant at low energy
scales. A robust Fermi liquid state (which resits superconductivity, even at
the lowest temperatures) wins and we get a jellium metal. Any electron
correlation effects go off the energy shell [24, 23]

What is fascinating is that in low, namely zero and one dimensions low energy
physics changes qualitatively and gives way to broad band Mott localization.

Nano Clusters of $\rm Cu^{0}$ as Singlet Reservoirs

Let us consider zero dimensions or small clusters of monovalent atoms for
emergent strong pairing correlations. Monovalent metal atom clusters have been
extensively studied from the point of view of electron correlation effects and
magnetism, using Hubbard model for example [38, 39, 40, 41, 42, 43, 44].
Nearest neighbor hoping t $\approx$ 3 eV and onsite repulsion U $\approx$ 7\.
For example Cu2 exists as a stable molecule in free space as a tightly bound
spin singlet with a binding energy of $\approx$ 2 eV.

In small Cu clusters, large antiferromagnetic exchange J $\sim$ 2 eV gives
rise to strong singlet bonds and resonances. In comparison to cuprates, where
J $\approx$ 0.15 eV, we have an order of magnitude larger singlet pairing
energy scale. Consequently local pairing temperature scale we get, on doping
are several hundreds of Kelvin.

Effective Hamiltonian for a monovalent atom cluster is the Heisenberg
antiferromagnetic Hamiltonian:

$H_{c}=\sum_{ij}J_{ij}({\bf S}_{i}\cdot{\bf
S}_{j}-1/4)\equiv-\sum_{ij}J_{ij}b^{\dagger}_{ij}b_{ij}$ (1)

Spin operator ${\bf S}_{i}\equiv
d^{\dagger}_{i\alpha}{\hat{\sigma}}_{\alpha,\beta}d_{j\beta}$. Spin-spin
interaction term has been rewritten using bond singlet operator
$b^{\dagger}_{ij}\equiv(1/\sqrt{2}(d^{\dagger}_{i\uparrow}d^{\dagger}_{j\downarrow}-d^{\dagger}_{i\downarrow}d^{\dagger}_{j\uparrow})$,
(d’s are electron operators) to emphasize the fact that strongly coupled bond
singlets are natural basic units (rather than spin moments) in clusters. It is
the resonance of these bond singlets (metal cluster aromaticity) which is at
the heart of our mechanism for superconductivity. More importantly resonance
energies of bond singlets are comparable to cohesion energy of monovalent
clusters [40].

Consider a simple tetrahedral cluster, an elementary closely packed unit in
fcc lattice. Two singlets resonate and stabilize a singlet ground state.
Larger clusters have been studied from the point of view of equilibrium shapes
of clusters. For us what is important is the singlet content and their
resonance. The tetrahedral RVB system has a doubly degenerate singlet ground
state. The symmetric tetrahedra distorts via valence bond localization - this
is an analogue of spin-Peierls distortion in the finite cluster.

We have suggested presence of a distribution of clusters in Cu Pb Apatite. As
we go to larger clusters there are more complex resonances, depending on size
and shape of a cluster. Singlets and their resonances continue to be present
till a critical size. Beyond the critical size singlets and their resonances
go off the energy shell [24, 23] and the bulk of the cluster behaves as a
jellium metal (filled shells) at low energies. Surfaces of larger clusters are
likely to have interesting singlet correlations and resonances.

It is clear that small clusters of various kinds are important as singlet
reservoirs, as they provide a large singlet pairing energy scale J $\approx$ 2
eV.

$\rm Cu^{0}$ Chains as Singlet Reservoirs

In an ideal Cu chain in equilibrium, without dimerization, half filled band of
4s electrons get Mott localized [45, 46, 47], as a consequence of low
dimensionality, even though the band is broad ($\sim 8$ eV wide) band. This
spin-half Mott insulating chain becomes a reservoir of resonating singlets and
large exchange energy scale J $\sim$ 2 eV. If we ignore lattice dimerazation
instabilities, 1D Hubbard model predicts a Mott insulator with a Mott Hubbard
gap of $\sim$ 1 eV (Lieb and Wu [45]) and superexchange J $\sim$ 2 eV for Cu.
Effective Hamiltonian of neutral Cu chain is given by an antiferromagnetic
Heisenberg Hamiltonian similar to equation 1.

When we move away from half filling we get, Luttinger liquid with spin-charge
decoupling, enhanced local pairing and a very large fluctuationg pairing gap
in ideal doped Cu chains. An organization of Cu chains into weakly interacting
chain and doping easily brings in ambient Tc superconductivity as discussed
below.

For illustration, consider a 2D array of weakly coupled (for a range of
interchain hopping t⟂) Cu chains. To avoid competing antiferromagnetic 2
sublattice order, we stagger the chains, so that when t⟂ = t we get isotropic
triangular lattice. Using known studies we find doped Mott insulator physics
and ambient Tc superconductivity, provided competing instabilities such as
dimerization and valence bond order are taken care of.

In the above system pair tunneling plays an important role in stabilizing 2D
or 3D superconducting order. We have used Wheatley-Hsu-Anderson’s mechanism
[25] of interlayer pair tunneling and similar interchain pair tunneling theory
of Zachar, Kivelson and Emery [26] to estimate superconducting Tc’s.

Scale of Tc is given by kBTc $\approx$ t${}^{2}_{\perp}$/$\Delta_{s}$. Here
$\Delta_{s}$ is a local spin gap scale. It can be approximated by
superexchange J, in the present context. It is easily seen that for reasonable
choice of t⟂, in our model system of Ag and Au we get 2D superconducting Tc at
ambient temperature scales. Thus we conclude that nanoscale Cu segments are
also important singlet reservoirs, with a high pairing temperature scale.

2D Planar $\rm Cu^{0}$ net as Singlet Reservoirs

Two dimensions is equally interesting. We considered hypothetical 2D
triangular lattices of Cu. Firstly, U is not strong enough to result in Mott
localization. However, a form of (self) doped Mott insulator physics
(dynamical Mott localization) survives at and close to half filling. Using
results from existing body of literature [48, 49, 50, 51, 52, 53, 54] we have
estimated the Tc. After doping, results for Tc are not very different for
triangular lattice and square lattice, even though symmetry of order
parameters are different (d-wave for square lattice[5], and d + id for
triangular lattice and honeycomb lattices [55, 17]).

Interestingly, in these model studies kBTc is a finite fraction of hoping
parameter t, even for intermediate repulsions. In view of this non-BCS
dependence we obtain ambient temperature (and higher) pairing temperature
scales for Cu subsystem at and close to half filling.

What is important for LK99 is that the nanoscale 2D patches will support
strong singlet pairing, which could eventually contribute to observed 3D hot
Tc granular superconductivity.

## Randomly Coupled Charged AMI’s and Emergence of Granular Superconductivity

As we mentioned earlier, weakly coupled reservoir of neutral singlets
(clusters, chand and 2D patches) could transfer a fraction of 4s valence
electrons to the insulating host. Transferred electrons will in general get
trapped and become immobile, via disorder and correlation effects. In the
charged AMI’s resonating neutral singlets become charged Cooper pairs and get
delocalized within the cluster, via physics of doped Mott insulator. Our
clusters and chains, become Cooper pair boxes and are ready for Josephson
coupling with neighbors.

We have a model of randomly Josephson coupled Cooper pair boxes. In view of
the strong local pairing scale, our estimate of superconducting Tc easily
reach ambient temperature superconducting scales. We hope to discuss details
in a later article.

## Competing Orders

Competing orders are ubiquitous in strongly correlated electron systems such
as cuprates. These orders are typically charge and spin stripe orders
accompanied by strong lattice distortions. Competing orders are encouraged
when one has certain structural chemistry. For example, in LSCO, a monolayer
cuprate, octahedral rotation strongly encourages stripe order and
superconducting Tc comes down close to 30 K. Whereas in YBCO, a bilayer
cuprate, the bilayer nature prevents easy octahedral rotation and Tc is more
than 90 K. Trilayer cuprates, the middle layers are protected from competing
orders and Tc’s exceed 120 K.

In clusters and chains there is a strong tendency for dimerization type of
instability. In chain it is the spin-Peierls (valence bond) order accommanied
by lattice dimerization

It was our prediction [24] that monovalent single layer atomic net is a
potential hot superconductor. We also emphasize that in real 2D monoatomic
monovalent metal system [56], when synthesized, competing orders such as
valence bond order and lattice dimerization will compete and supress
superconductivity. Free hanging 2D layers are good bet, as they are not
influenced by 2D layer grown on substrates.

A deeper understanding of structural chemistry will enable structural design
with reduced competing order.

## Discussion and Conclusion

Excited by the report of observation of hot superconductivity in the Cu-Pb
Apatite system, we have presented a couple of scenarios, keeping solid state
and quantum chemical constraints in mind and suggested possibility of hot
superconductivity via electronic mechanism of doped Spin$\frac{1}{2}$ Mott
insulator, utilizing broad band Mott localization.

Our theory makes use of doped spin-$\frac{1}{2}$ doped Mott insulator route to
superconductivity, with the emphasis on the fact that monovalent metal atoms
under suitable conditions (such as reduced dimensions) readily undergo Mott
localization. Over years we have used this AMI idea to develop tangible models
to understand superconductivity in complex solid state systems such as
hydrides, K3p-Terphenyl, B-doped diamond and Ag-Au nanostructures, doped
graphene and silicene.

In the present article we have proposed a compelling scenario, without full
justification from quantum chemical calculations etc. In the same spirit, we
briefly discuss possibility of a related, but somewhat different scenario.
Here Cu atom get accommodated as a neutral impurity and creates an s-like
donor impurity state in the band gap of $\rm Pb_{10}(PO_{4})_{6}O$ . As
impurity concentration is increased impurity state wave functin will overlap
and we could have a Anderson-Mott insulator to superconductor transition, as
outlined in our theory for B doped diamond [20]. Here also, depending on the
binding energy of the donor state we could in principle get Tc reaching the
scales claimed.

As a key idea and notion, superconductivity in doped Mott insulators began
while attempting to understand high Tc cuprates by Anderson in 1987. It is
becoming increasingly clear that this idea has a wide applicability. We once
again emphasize, excited by our theory, experimental exploration of the rich
field of low dimensional monovalent metals, monovalent metal doped minerals,
insulators and in metallurgy, where variety of nanostructures are norm of the
day.

We also wish to reiterate that while $\rm Cu^{2+}$ doping is an attractive
idea, correponding isolated bare band widths are so low that Cooper pairs will
get localized strongly by lattice distortion and disorder. For the scale of
reported Tc, we need broad band Mott localization emphasized in the present
article, rather than narrow band Mott localization. It is important to recall
that even in high Tc cuprates an effective single band (which is not isolated)
width is of the order of 2 eV, which is capable of giving superconductivity.
at temperatures close to 100 K.

Our broad band Mott localization and AMI is likely to play a key role in other
systems perhaps. We conclude our article by reproducing a paragraph from our
2019 arXiv article on broad band Mott localization and emergent AMI (reference
[23], page 6, 2nd column, 3rd paragraph in discussion session):

‘Does ambient temperature granular superconductivity exist already, in some
form or other, in the rich world of minerology and metallurgy - natural ones
and man made ? An exploration into minerology and metallurgy is likely to
offer surprises’.

## Acknowledgement

I thank Veer Awana for sharing his very recent results on synthesis and
measurements of LK-99 samples. Discussions with V N Muthukumar, after several
years, on issues of hot superconductivity and renewal of our collaboration is
acknowledge. Brief but encouraging discussions with Venky Venkatesan and
Rajaram Nityananda are acknowledged. Valuable help from Bala Yogendra for
arXiv uploading is acknowledged. I thank the Institute of Mathematical
Sciences, Chennai, IITMadras and Perimeter Institute of Theoretical Physics,
Waterloo, Canada for continuing hospitality. Perimeter Institute for
Theoretical Physics, Waterloo, Canada is supported by the Government of Canada
through Industry Canada and by the Province of Ontariothrough the Ministry of
Research and Innovation.

P.W. Anderson showed us a fertile direction to achieve hot supercondtivity. A
phenomenologist par excellence, who listened to nature always, Phil was
skeptical about existence of hot superconductors, before cuprates. In his
article [57] Superconductivity in the Past and Future (1969) he states, ‘Room-
temperature magnets or superconducting skis ? Probably not. For one thing, so
many diverse types of materials exist in nature or have been made in the
laboratory that large, obvious effects would probably already have been seen.
But room-temperature interferometers ? Who knows ?’

My obsession with possible hot superconductivity, via broad band Mott
localization, ever since I got entangled with RVB theory from January 1987,
did bring constant appreciation from Phil. I dedicate this article to the
memory of P.W. Anderson, my mentor, whose 100th birth anniversary is on the
horizon.

## References

  * [1] G. Bednorz and K.A. Muller, Z. Phys. B 64 189 (1986)
  * [2] P.W. Anderson, Science 235 1196 (1987)
  * [3] G. Baskaran, Z. Zou and P.W. Anderson, Solid St. Com. 63 973 (1987); P.W. Anderson, G. Baskaran, Z. Zou and T. Hsu, Phys.Rev.Lett., 58 2790 (1987); G. Baskaran and P.W. Andersdon, Phys. Rev., B 37 580 (1988)
  * [4] S.A. Kivelson, D.S. Rokhsar, and J. P. Sethna Phys. Rev., B 35 8865(R) (1987)
  * [5] G. Kotliar, Phys. Rev. B 37 3664 (1988)
  * [6] P.W. Anderson et al., J. Physics: Condensed Matter, 16 R755 (2004)
  * [7] G. Baskaran, Iranian J. Phys. Res. 6 163 (2006); arXiv:cond-mat:0611548
  * [8] X.G. Wen, Quantum Field Theory of Many-Particle Systems (Oxford University Press 2010)
  * [9] Y. Zhou, K. Kanoda T.K. Ng, Rev. Mod. Phys., 89 025003 (2017)
  * [10] G. Baskaran, Pramana, 73 61 (2009);
https://www.ias.ac.in/article/fulltext/pram/073/01/0061-0112

  * [11] N. Dasenbrock-Gammon et al., Nature, 615 244 (2023)
  * [12] A.D. Grockowiak et al., Frontiers in Electronic Materials, 2 2.837651 (2022)
  * [13] D.K.Thapa arXiv:1807.08572; S.K. Saha et al 2022 Supercond. Sci. Technol. 35 084001 (2022)
  * [14] Sukbae Lee, Ji-Hoon Kim and Young-Wan Kwon, arXiv:2307.12008; Sukbae Lee et al., arXiv:2307.12037
  * [15] Y. Kopelevich, R.R. da Silva and B.C. Camargo, Physica C 514 237 (2015); Y. Kopelevich and P.D. Esquinazi, J. Low Temp. Phys. 146 629 (2007); P.D. Esquinazi et al., Quantum Studies: Mathematics and Foundations, 5 41 (2018)
  * [16] G. Baskaran, Phys. Rev. B65 212505 (2002)
  * [17] Annica M. Black-Schaffer and Sebastian Doniach, Phys. Rev.B75 134512 (2007)
  * [18] S. Pathak, V.B. Shenoy and G. Baskaran, Phys. Rev. B 81 085431 (2010)
  * [19] G. Baskaran, arXiv:1309.2242 (2013); Book Chapter, Many-body Approaches in Different Scales (Eds. G.G.N. Angilella and C. Amovilli; Springer Nature 2018)
  * [20] G. Baskaran, arXiv:cond-mat:0404286 (2004); Sci. Technol. Adv. Mater. 7 S49 (2006); Sci. Technol. Adv. Mater. 9 044104 (2008)
  * [21] G Baskaran, arXiv:1704.08153 (2017)
  * [22] G. Baskaran, arXiv:1507.03921 (2015)
  * [23] G. Baskaran, arXiv:1906.02143
  * [24] G. Baskaran, Int. J. Mod. Phys. B36 2250184 (2022); arXiv:1808.02005 (2018)
  * [25] J. Wheatley, T. Hsu and P.W. Anderson, Nature, 333 121 (1988); V.N. Muthukumar and G. Baskaran, Mod. Phys. Lett. B8 699 (1994); G. Baskaran, Phil. Mag., B76 119 (1997)
  * [26] Z. Zachar, S.A. Kivelson and V. J. Emery, Journal of Superconductivity, 10 373 (1997)
  * [27] S. Ramachandra Rao, Nature, 136 436 (1935)
  * [28] J. Imry, Phys. Rev., B 91 104503 (2015)
  * [29] A Hernando et al., New Journal of Physics 16 073043 (2014)
  * [30] P.G. van Rhee et al., Phys. Rev. Lett., 111 127202 (2013)
  * [31] G. Frommeye and K.F. Schneider, Phys. Stat. Sol:(a)41 653 (1977)
  * [32] R. Deblock et al., Phys. Rev. Lett., 89 206803 (2002)
  * [33] W. Luo et al., Nano Letters, 7 3134 (2007); V Tuboltsev et al., ACS Nano, 7-8 6691 (2013)
  * [34] Junwen Lai et al., arXiv:2307.16040; S.M. Griffin, arXiv:2307.16892; Liang Si and Karsten Held, arXiv:2308.00676; Rafal Kurleto et al., arXiv:2308.00698
  * [35] Kapil Kumar, N.K. Karn and V.P.S. Awana, arXiv:2307.16402
  * [36] Li Lu et al., arXiv:2307.16802
  * [37] Xiaoyan Yang, Alberto J. Fernández-Carrión and Xiaojun Kuang, Chem. Rev., (24th July 2023), online: https://doi.org/10.1021/acs.chemrev.2c00913
  * [38] M Torrini and E Zanazzi, J. Phys. C: Solid State Phys.9 63 (1976)
  * [39] G. M. Pastor, R. Hirsch and B. Muhlschlegel, Phys. Rev. B 53 10382 (1996)
  * [40] G.M. Pastor in Atomic clusters and nanoparticles (Ed.C. Guet et al.) Les Houces Session LXXIII (Springer 2000)
  * [41] Eva M. Fernández et al., Phys. Rev. B 70 165403 (2004)
  * [42] W.A. de Heer, Rev. Mod. Phys., 65 611 (1993)
  * [43] B. K. Rao and P. Jena, Phys. Rev., B32 2058 (1985)
  * [44] M.N. Huda and A.K. Ray, Eur. Phys. J., D 22 217 (2003)
  * [45] E. Lieb and F.Y. Wu, Phys. Rev. Lett., 20 1445 (1968)
  * [46] M. H. McAdon and W.A. Goddard III, The J. of Chem. Phys., 88 277 (1988)
  * [47] Metal Chain/Chains of Metal, Ed. M. Springborg and Yi Dong (Elsevier 2007)
  * [48] M. Laubach et al., Phys. Rev., B 91 245125 (2015)
  * [49] J. Wang et al., Int. J. Mod. Phys. C28 1750127 (2017)
  * [50] C-D. Hebert, Semon and A.-M.S. Tremblay, Phys. Rev.,B 92 195112 (2015); S. Acheche et al., Phys. Rev., B 94 245133 (2016)
  * [51] K. Misumi, T. Kaneko, and Y. Ohta, Phys. Rev., B 95 075124 (2017)
  * [52] W.W.De Silva et al., Phys. Rev., B 93 205111 (2016)
  * [53] Tanmoy Das, R.S. Markiewicz and A. Bansil (2014) Advances in Physics, 63 151 (2014) DOI: 10.1080/00018732.2014.940227
  * [54] K.-S. Chen et al., Phys. Rev. B 88 245110 (2013)
  * [55] G. Baskaran, Phys. Rev. Lett. 91 097003 (2003)
  * [56] Janne Nevalaita and Pekka Koskinen, Phys. Rev. B97 035411 (2018)
  * [57] P.W. Anderson, in Superconductivity, pp. 1343-1358 Part 2, edited by R. Parks (Marcel Dekker, New York, 1969); this article is reproduced in P.W. Anderson’s book ‘A Career in Theoretical Physics’ (World Scientific, 2004)


# The structural stability of tungsten nanoparticles

L. Pizzagalli  S. Brochard  J. Godet  J. Durinck  Laurent.Pizzagalli@univ-
poitiers.fr Institut P’ , CNRS UPR 3346, Université de Poitiers, SP2MI,
Boulevard Marie et Pierre Curie, TSA 41123, 86073 Poitiers Cedex 9, France

###### Abstract

Motivated by contradicting reports in the literature, we have investigated the
structural stability of tungsten nanoparticles using density functional theory
calculations. The comparison of BCC, FCC, A15, disordered, and icosahedral
configurations unequivocally shows that BCC is the energetically most stable
structure when the number of atoms is greater than 40. A disordered structure
is more stable for smaller sizes. This result conflicts with an earlier
theoretical study on transition metal nanoparticles, based on a semi-empirical
modeling of nanoparticles energetics [D. Tománek et al., Phys. Rev. B 28, 665
(1983)]. Examining this latter work in the light of our results suggests that
an erroneous description of clusters geometry is the source of the
discrepancy. Finally, we improve the accuracy of the semi-empirical model
proposed in this work, which will be useful to calculate nanoparticle energies
for larger sizes.

##  I Introduction

Due to its high melting and evaporation temperatures and excellent mechanical
properties, tungsten is the selected wall material to be used in fusion
reactor divertors like ITER. Despite a high thermal stability and a good
resistance to sputtering, the interaction with the plasma during operation
leads to the formation of dust composed of tungsten nanoparticles of various
sizes [1]. This aspect motivated extensive investigations on the properties of
W nanoparticles, and in particular their environmental and biological impacts
in a fusion-like environment [2]. Tungsten-based nanomaterials are also
increasingly used in biomedicine applications [3].

The most stable phase of bulk tungsten is named $\alpha$-W and has a body-
centered structure (BCC). Two metastable allotropes are known, which can form
in specific conditions. The first one is named $\beta$-W and has a A15 cubic
structure. $\beta$-W thin films have been extensively studied owing to the
report of a giant spin Hall effect [4]. The second one is $\gamma$-W, with a
face-centered cubic (FCC) structure [5]. Naturally, one would assume that
tungsten nanoparticles are made of the most stable bulk phase, i.e.
$\alpha$-W. However, in an early theoretical paper, Tománek and co-workers
predicted that a FCC structure would be favored in small nanoparticles for
several transition metals, including tungsten [6]. Their conclusions were
based on the argument that surface energies associated with the FCC structure
are typically lower than the ones for the BCC structure, and that this effect
will dominate at small sizes when the surface-to-bulk ratio increases. In the
specific case of tungsten, they calculated a nanoparticle size threshold of
5-6 nm, below which FCC should be the lowest energy structure. Later X-ray
diffraction experiments on nanometer sized W clusters confirmed the
prediction, and determined the BCC-FCC structural transition at 7 nm [7].
Similar conclusions in support of this argument were reached for Mo and Cr
clusters [8, 9]. But other experiments led to different observations. For
instance Iwama and Hayakawa found that 3–20 nm Mo and W nanoparticles show
crystalline structures being either BCC or A15 for Mo, and A15 for W, but not
FCC [10]. More recently, Schöttle and co-workers reported transmission
electron microscopy observations of 1-2 nm W nanoparticles with a BCC
structure [11].

On the numerical side, classical molecular dynamics simulations were carried
out for investigating this possible BCC-FCC transition as a function of the
size and shape in W nanoparticles. Hence Chen and co-workers reported that the
BCC structure is significantly more stable than the FCC, except for a singular
high energy shape [12]. This work confirms an earlier study by Marville and
Andreoni, who found that tungsten nanoparticles were more stable in a BCC
structure than in a FCC or icosahedral arrangement [13]. These simulations
then support the recent microscopy measurements against the early predictions
and experiments. However one has to be cautious about definite conclusions,
since interatomic potentials are not always accurate in the description of
under-coordinated atoms at the surfaces and edges of nanoparticles. Another
aspect is the potential stabilization of the A15 phase in the W nanoparticles,
which was overlooked in published numerical works. This motivates us to
perform a thorough investigation of the structural stability of small tungsten
nanoparticles using first-principles calculations. In particular we consider
numerous systems with BCC, FCC, A15, disordered, and icosahedral structures.

##  II Methods

###  II.1 Electronic structure calculations

Table 1: Computed data for tungsten in BCC, FCC, and A15 structures, compared to experimental results when available, and DFT computed data from the literature. The lattice parameter $a_{0}$ is expressed in angstrom, the bulk modulus $B$ in GPa, and the surface energies $\gamma$ in J m-2. The bulk energy per atom $\varepsilon_{0}$ are given relatively to the BCC phase (i.e. $\varepsilon_{0}(\mathrm{BCC})=0$), and are expressed in eV at-1. Note that the energies for {110} and {111} surfaces of FCC W (marked with a *) correspond to the unrelaxed slab configuration (see text for details). |  |  | DFT
---|---|---|---
|  | Exp. | This work | Others
BCC | $a_{0}$ | 3.165111Reference 14 | 3.1753 | 3.1741222Reference 15, 3.172333Reference 16
$B$ | 310444Reference 17 | 311 | 309555Reference 18
$\gamma_{100}$ |  | 4.071 | 3.954666Reference 19
$\gamma_{110}$ | 3.265777Reference 20, 3.675888Reference 21 | 3.302 | 3.23066footnotemark: 6
$\gamma_{111}$ |  | 3.569 | 3.46666footnotemark: 6
FCC | $a_{0}$ | 4.1311footnotemark: 1 | 4.0256 | 4.02533footnotemark: 3, 3.960999Reference 22
$\varepsilon_{0}$ |  | 0.494 | 0.4933footnotemark: 3, 0.5099footnotemark: 9
$B$ |  | 282 | 286101010Reference 23
$\gamma_{100}$ |  | 3.284 |
$\gamma_{110}$ |  | 3.736∗ |
$\gamma_{111}$ |  | 2.637∗ |
A15 | $a_{0}$ | 5.0511footnotemark: 1 | 5.0665 | 5.05933footnotemark: 3
$\varepsilon_{0}$ | 0.082111111Reference 24 | 0.093 | 0.0933footnotemark: 3
$B$ |  | 298 | 29866footnotemark: 6

We perform first principles calculations in the framework of density
functional theory (DFT) for determining the relaxed structure and the
associated energy of various nanoparticle models. The code pw.x in the Quantum
Espresso software [25] is used for this purpose. A well converged electronic
structure is achieved by using a plane wave energy cutoff of 40 Ry and a
charge density cutoff of 320 Ry. Exchange-correlation contributions are
described with the Perdew-Burke-Ernzerhof functional [26]. We use the
Projector Augmented-Wave method [27] for ion-electron interactions, with the
valence electron configuration 5s25p65d46s2. Finally, the cold smearing method
[28] of Marzari and co-workers is applied to improve the convergence of the
electronic structure.

We first compute the lattice parameter and bulk modulus for the three tungsten
allotropes, and their energie differences, using supercells and a very dense
grid of k-points for Brillouin zone sampling. Our results in Tab. 1 are in
excellent agreement with other recent calculations and available experiments.
The table also includes the energies of surfaces with low Miller indexes for
BCC and FCC W. Those are calculated using a slab configuration, a large number
of layers, and a k-point grid of $20\times 20\times 1$, in order to obtain
converged results. The relaxation is achieved when all components of all ionic
forces are lower than $2.6\times 10^{-4}$ eV Å-1. Calculated surface energies
for BCC W are in good agreement with the literature. We also consider the
canonical surfaces of FCC W, for which no information seems to be available in
the literature. In the case of the $(110)$ and $(111)$ surfaces, we observe
distortions occurring in the center of the slab during force relaxation. We
have not pursued the analysis of this issue any further, since it is not the
focus of the study. For these two surfaces, we report in Tab. 1 the energies
of the unrelaxed initial configuration.

For the relaxation of nanoparticles, we employ supercells large enough to
allow for at least 10 Å between periodic replicas in all dimensions. The
k-point sampling is made at the $\Gamma$-point as is usual for 0D systems.
Finally, forces are relaxed until all components for all ions are lower than
$2.6\times 10^{-4}$ eV Å-1.

###  II.2 Nanoparticles selection

The shape of 1-2 nm nanoparticles, i.e. including about thirty to a few
hundred atoms, is often poorly documented, and tungsten is no exception.
Therefore we follow a well proven, standard methodology [29], in order to
generate low energy initial configurations. For BCC and FCC for which the
energies of surfaces are known, we use a Wulff construction as a first option.
Alternatively, we also carve nanoparticles out of bulk, with spherical or
smoothed cubic shapes. A coordination analysis is next performed to select
configurations with as few low coordination atoms as possible, since the
latter are associated with a high energy in metals. These two techniques are
both used to obtain initial configurations of BCC and FCC nanoparticles. For
A15 for which surfaces are not known, only spherical nanoparticles are
generated.

The carving strategy is more difficult to implement in the case of a
disordered structure, since bulk tungsten does not exist in an amorphous phase
to our knowledge. Recently, Jana and Caro performed an extensive search for
the most stable structures of iron nanoparticles including less than 200 atoms
[30]. Except for a few crystalline configurations at magic numbers, they found
that most of these structures are amorphous. We select several of these
configurations from their freely available structures database, which can be
used as initial configurations after a small rescaling. These systems are
labeled D (for disordered) in the remainder of the paper.

Finally, we also consider icosahedral nanoparticles, identified with label
’I’. The two possible candidates in the investigated size range contain 55 and
147 atoms [31]. Overall we select 27 different nanoparticles, with atoms
numbers ranging from 55 to 169, to be used as starting configurations in
first-principles calculations. The set includes 8 BCC, 7 FCC, 5 A15, 5
disordered (D), and 2 icosahedral (I) nanoparticles.

##  III Stability

Figure 1:  Examples of DFT relaxed tungsten nanoparticles, with various sizes
(given as a number of atom in the label) and different atomic structures (BCC,
FCC, A15, D for disordered, I for icosahedral). Figure 2:  Nanoparticle energy
per atom (eV) as a function of the number of atoms $N$, the reference being
the BCC bulk energy, after DFT relaxation (BCC: blue circles; FCC, orange
squares; A15: magenta diamonds; D: green triangles; I: golden stars). Dashed
lines are obtained by fitting all data points for a given set with Eq. 3 (fit
parameters reported in Tab. 2).

Selected examples of DFT relaxed nanoparticles are represented in Fig. 1. In
most cases we observe limited structural changes compared to the initial
configurations. In particular only a slight surface relaxation is observed for
BCC and I nanoparticles, with an energy gain of 0.15 eV/at in average. For A15
and FCC systems, the surface relaxation is greater, with an average relaxation
energy of 0.26 eV/at and 0.32 eV/at, respectively. The larger energy gain is
obtained for D nanoparticles, for which significant surface atom displacements
are identified. In average, the relaxation of D nanoparticles yields an energy
gain of 0.83 eV/at.

The energy $E$ (with respect to the perfect BCC bulk) of a relaxed
nanoparticle made of $N$ atoms can be written

$E=N\varepsilon_{0}+E_{s}(N),$ (1)

with $\varepsilon_{0}$ the bulk energy per atom relatively to the ground state
BCC structure ($\varepsilon_{0}(\mathrm{BCC})=0$). $E_{s}$ is an excess
energy, akin to a surface energy, although it can in principle includes
further contributions associated with Laplace pressure, surface stress
relaxation, or quantum confinement effects. Assuming that only surface atoms
contributes to $E_{s}$, and that the number of surface atoms is roughly
$N^{2/3}$, Eq. 1 becomes

$E=N\varepsilon_{0}+\beta N^{2/3},$ (2)

with $\beta$ the average energy contribution to $E_{s}$ from surface atoms.
The energy per atom of a nanoparticle can then be written

$E/N=\varepsilon_{0}+\beta N^{-1/3}.$ (3)

Figure 2 shows the nanoparticle energies calculated by DFT for all
configurations. As predicted by Eq. 3, $E/N$ values increase for $N\rightarrow
0$ due to the growing surface contribution, and converge to $\varepsilon_{0}$
for $N\rightarrow\infty$. The most striking result is that for a given $N$,
BCC nanoparticles appear to be always more stable than the others. Fitting the
Eq. 3 separately for each kind of structure better highlights this finding
(dashed lines in Fig. 2). We observe overall
$E(\mathrm{BCC})<E(\mathrm{A15})<E(\mathrm{D})<E(\mathrm{I})<E(\mathrm{FCC})$
in the investigated $N$ range. Therefore our results disagree with predictions
made in Ref. 6. In fact our calculations clearly show that FCC nanoparticles
are not energetically more stable than BCC nanoparticles. They also reveal
that BCC is the lowest energy structure and that A15, disordered and
icosahedral nanoparticles are more stable than FCC nanoparticles.

Table 2:  $\beta$ parameters obtained by fitting Eq. 3 with DFT calculations (shown in Fig. 2). For BCC, FCC, and A15 nanoparticles, $\varepsilon_{0}$ values calculated for bulk systems are used in the fit. For disordered (D) and icosahedral (I) nanoparticles, $\varepsilon_{0}$ is computed together with $\beta$ by fit. The table also includes the corresponding surface energies calculated using the spherical model described in the appendix, and the number of atoms $N_{c}$ below which each phase becomes more stable than BCC. |  $\beta$ (eV) |  $\varepsilon_{0}$ (eV) |  $\gamma$ (J m-2) | $N_{c}$
---|---|---|---|---
BCC | 7.149 | 0.000 | 3.729 |
FCC | 5.861 | 0.494 | 3.019 | 17
A15 | 7.000 | 0.093 | 3.614 | 4
D | 6.099 | 0.307 | 3.142 | 40
I | 5.836 | 0.433 | 3.006 | 27

The $\beta$ values determined in the fitting process are reported in Tab. 2.
We adjust $\beta$ while using $\varepsilon_{0}$ from bulk calculations (Tab.
1), except for disordered and icosahedral nanoparticles for which
$\varepsilon_{0}$ is also adjusted. Note that only minor changes are observed
when both $\beta$ and $\varepsilon_{0}$ are adjusted on DFT data for BCC, FCC,
and A15 nanoparticles.

$\beta$ in Eq. 3 corresponds to an energy per surface atom. Conversion into
the usual surface energies $\gamma$, i.e. energies per surface area, is
straightforward if one assumes that nanoparticles are spherical (see the
Appendix). Computed values are included in Tab. 2. We find that the lowest
$\gamma$ values are obtained for I and FCC, followed by D, A15, and finally
BCC. This is in agreement with arguments based on a lower surface to volume
ratio for icosahedral and FCC [6, 31, 8]. We also observe that $\gamma$ values
for BCC and FCC are in the range of surface energies calculated for well
defined orientations as reported in Tab. 1. However, it is difficult to push
further the comparison. Swaminarayan and co-workers proposed that the surface
energy of spherical nanoparticles made of FCC metals is close to
$\gamma_{110}$ [32]. This is clearly not valid in the present work. Our value
of 3.019 J m-2 is also larger than predictions made using an analytical model
[33].

The largest surface energy for BCC systems necessarily implies that other
structures will be energetically favored below a given size. Using Eq. 3, the
critical transition between BCC and phase $X$ is predicted to occur at

$N_{c}=\left[\frac{\beta(\mathrm{BCC})-\beta(\mathrm{X})}{\varepsilon_{0}(\mathrm{X})}\right]^{3},$
(4)

since $\varepsilon_{0}(\mathrm{BCC})=0$. Computed values are reported in Tab.
2. We find that FCC nanoparticles become more stable than BCC nanoparticles
for $N\leq 17$. This is dramatically smaller than previously reported sizes of
$N_{c}=5660$ (Ref. 6) and $N_{c}=10470$ (Ref. 7). In addition, such a
transition is not relevant here because our data also suggest that a
disordered state is favored for $N\leq 40$. It is well known that a non-
crystalline molecular configuration should prevail for all metals at the
smallest scales. Our computed threshold is close to a measured value of
$N=30$, but this good agreement may be fortuitous given the large experimental
uncertainty [7] and the limited set of disordered configurations considered in
our simulations.

To conclude this section, our DFT calculations unambiguously show that a BCC
structure is favored for small tungsten nanoparticles of 1-1.6 nm.
Nanoparticles with FCC or A15 structures are always higher in energies. When
the number of atoms is lower than about 40, a disordered state becomes
favorable. These results are at odds with a theoretical analysis [6] and X-ray
diffraction experiments [7]. They are however in agreement with microscopy
measurements [11] and classical molecular dynamics calculations [13, 12]. The
main factor explaining the BCC stability over FCC despite the higher surface
energy contribution, is the large bulk energy difference in favor of BCC. It
is then obvious that the BCC structure will be favored for nanoparticles
larger than those investigated here. Finally, it is important to keep in mind
that we compute 0 K energy differences, and an evaluation of nanoparticle free
energies would be needed for more definite conclusions. As a first hint, the
bulk entropy difference between FCC and BCC is estimated to be approximately 1
$k_{B}$ per atom [34], thus at least one order of magnitude lower than the
energy differences in Tab. 1 at 300 K.

##  IV Semi-empirical modeling

Figure 3:  Surface energy $E_{s}$ computed with the Eq. 5 plotted against $E_{s}$ calculated by DFT, for all nanoparticles (BCC: blue circles; FCC, orange squares; A15: magenta diamonds; D: green triangles; I: golden stars). Open symbols show values obtained using the initial parameters as given in Ref. 6 and reported in Tab. 3 (Opt0, see text for details), while filled symbols correspond to those obtained using optimized parameters (Opt2). The dashed black line marks the perfect agreement. The inset graph shows the root mean square deviations (RMSD) for each set (lines with same colors than the points in the main graph) and different optimization levels (Tab. 3). Table 3: Root mean square deviations (RMSD, in eV) of $E_{s}$(DFT)-$E_{s}$(Eq. 5) for each set and different optimization levels. Opt0: $E_{s}$(Eq. 5) is calculated using the original parameters ($\lambda$,$\eta$) given in Ref. 6. Opt1: $\eta$ is optimized. Opt2: both $\lambda$ and $\eta$ are optimized. Opt |  | BCC | FCC | A15 | D | I
---|---|---|---|---|---|---
0 | $\eta$ | 0.4 | 0.08 | 0.4 | 0.08 | 0.08
$\lambda$ | 0.5 | 0.5 | 0.5 | 0.5 | 0.5
RMSD | 23.17 | 13.06 | 29.46 | 22.01 | 9.84
1 | $\eta$ | 1.00 | 0.20 | 0.42 | 0.00 | 0.33
$\lambda$ | 0.5 | 0.5 | 0.5 | 0.5 | 0.5
RMSD | 12.91 | 11.19 | 27.76 | 17.94 | 4.15
2 | $\eta$ | 1.00 | 0.20 | 0.77 | 0.04 | 0.12
$\lambda$ | 0.546 | 0.461 | 0.579 | 0.619 | 0.540
RMSD | 4.35 | 7.17 | 2.20 | 2.76 | 0.00

In this section, we focus on the semi-empirical model proposed by Tománek and
co-workers [6]. We first aim at comparing and improving predictions from this
model with our DFT calculated nanoparticle energies. In a second step, we try
to understand why using this model leads to overestimated sizes for the BCC-
FCC transition.

We summarize the basics of this model in the first place. Using the second
moment approximation, Tománek and co-workers proposed that $E_{s}$ in Eq. 1
can be approximated by

$E_{s}=E_{c}\sum_{i=1}^{N_{s}}\left[(Z_{i}/Z_{b})^{\lambda}-1\right],$ (5)

where $N_{s}$ is the number of surface atoms and $\lambda=1/2$. $E_{c}$ is the
cohesive energy of the bulk state. $Z_{i}$ is an effective coordination number
including both first and second neighbors: $Z_{i}=Z_{i}^{1}+\eta Z_{i}^{2}$
where $Z_{i}^{1}$ and $Z_{i}^{2}$ refer to nearest neighbors and next-nearest
neighbors, respectively. In their paper Tománek and co-workers reported that
appropriate values for parameter $\eta$ are 0.08 for FCC and 0.4 for BCC [6].
They also defined surface atoms as those with $Z_{i}$ lower than 10. Finally,
$Z_{b}$ is the coordination for bulk atoms, which depends on the atomic
structure.

We first compute $Z_{i}$ for all BCC and FCC nanoparticles studied in the
present work, from the DFT relaxed configurations and following the rules
mentioned above. Using $E_{c}(\mathrm{BCC})=-8.90$ eV/at [35] and
$E_{c}(\mathrm{FCC})=-8.90+\varepsilon_{0}=-8.406$ eV/at, $E_{s}$ is
calculated according to Eq. 5. The results are plotted against DFT data in
Fig. 3, and labeled as ’Opt0’. It appears that compared to DFT $E_{s}$ for BCC
systems is systematically underestimated, whereas it is overestimated for FCC
systems. The root mean square deviations (RMSD) are reported in Tab. 3. One
can see that the error is significant, in particular for BCC nanoparticles.

Nanoparticles with A15, disordered and icosahedral structures are not
considered in Ref. 6, and we have to determine appropriate guess for model
parameters. For disordered D and icosahedral I nanoparticles, we use the same
parameters than for FCC ones. We know that $E_{c}(\mathrm{D})=-8.593$ eV/at
and $E_{c}(\mathrm{I})=-8.467$ eV/at, using $\varepsilon_{0}$ from Tab. 2. In
the A15 structure, there are 8 atoms in the elementary cell, with two
different atomic environments. 2 of them have 12 neighbors at 2.83 Å, and the
remaining 6 have 2 neighbors at 2.53 Å, 4 at 2.83 Å, and 8 at 3.102 Å. We
assume that the first neighbors distance is lower than 2.95 Å, i.e.
$Z_{i}^{1}=7.5$ on average. Neighbors atoms at 3.102 Å are considered second
neighbors, i.e $Z_{i}^{2}=6$ on average. We use $\eta=0.4$ as initial guess by
analogy with BCC, and $E_{c}(\mathrm{A15})=-8.807$ eV/at. In all cases, atoms
are identified as belonging to the surface if their effective coordination
number is lower than 10, as in Ref. 6. $E_{s}$ values computed with the model
for A15, D and I nanoparticles are included in Fig. 3. In most cases they are
underestimated compared to the DFT results. The RMSD for D and A15 are similar
to the BCC value (Tab.3). Note that the low RMSD value for I is due to the
limited size of the set (two nanoparticles).

In order to improve the accuracy of the semi-empirical model, an optimization
method is applied with $\eta$ as a variable parameter (Opt1). The new $\eta$
and RMSD values are shown in Tab. 3. The best improvements are obtained for
BCC and I, but only a moderate RMSD reduction is observed for FCC, A15 and D
nanoparticles. Finally, we allow both $\eta$ and $\lambda$ to vary during the
optimization (Opt2). Final results are plotted in Fig. 3. A remarkable
refinement is achieved for A15, D and BCC systems, and to a lesser extent for
FCC. This is clearly demonstrated by the RMSD values which become lower than 3
eV for D and A15. Overall, we find that the model is highly sensitive to
$\lambda$, and much less to $\eta$. Except for FCC, a $\lambda$ value slightly
greater than 0.5 greatly increases the model accuracy.

Figure 4:  $\delta(N)$ as a function of $N$, the number of atoms: data from
Fig. 3 in Ref. 6 (blue line), and computed from DFT data using Eq. 9 (orange
line). The dotted blue line shows the best fit of $\delta(N)$ data in Ref. 6
with the expression $\mu N^{-1/3}$. The dashed green line is positioned at
$(E_{c}(\mathrm{BCC})-E_{c}(\mathrm{FCC}))/E_{c}(\mathrm{BCC})=0.0555$ (using
$E_{c}(\mathrm{BCC})=-8.90$ eV and $E_{c}(\mathrm{FCC})=-8.406$ eV).

Surface energies calculated using the original model in Ref. 6 are therefore
not very accurate, but it is not clear whether this is the main reason behind
the difference in prediction for the BCC-FCC transition. In the same paper,
the following expression is defined:

$\delta(N)=\frac{1}{N}\left[\sum_{\begin{subarray}{c}i=1\\\
\mathrm{FCC}\end{subarray}}^{N_{s}}\left[(Z_{i}/Z_{b})^{1/2}-1\right]-\sum_{\begin{subarray}{c}i=1\\\
\mathrm{BCC}\end{subarray}}^{N_{s}}\left[(Z_{i}/Z_{b})^{1/2}-1\right]\right].$
(6)

$\delta(N)$ depends only on the cluster geometry and is represented in the
Fig. 3 of Ref. 6 and reproduced in our Fig. 4. Tománek and co-workers proposed
that the FCC-BCC transition occurs when the following condition is met:

$\delta(N)\stackrel{{\scriptstyle
N=N_{c}}}{{=}}\frac{E_{c}(\mathrm{BCC})-E_{c}(\mathrm{FCC})}{E_{c}(\mathrm{BCC})},$
(7)

with the right-hand side a constant, plotted as a dashed line in Fig. 4. They
found that $N_{c}=5660$ using their $\delta(N)$ values.

We remark that $\delta(N)$ can also be computed from the $\beta$ values in
Tab. 2 determined from DFT results. In fact, using Eq. 5, we can write

$\delta(N)=\frac{1}{N}\left[\frac{E_{s}(FCC)}{E_{c}(\mathrm{FCC})}-\frac{E_{s}(BCC)}{E_{c}(\mathrm{BCC})}\right].$
(8)

Since $E_{s}=\beta N^{2/3}$, we finally obtain

$\delta(N)=N^{-1/3}\left[\frac{\beta(\mathrm{FCC})}{E_{c}(\mathrm{FCC})}-\frac{\beta(\mathrm{BCC})}{E_{c}(\mathrm{BCC})}\right].$
(9)

Figure 4 shows $\delta(N)$ computed with Eq. 9, $\beta$ values in Tab. 2,
$E_{c}(\mathrm{BCC})=-8.90$ eV and $E_{c}(\mathrm{FCC})=-8.406$ eV. There is
clearly a large difference with $\delta(N)$ as given in Ref. 6, which likely
explains the disagreement concerning the critical size of the BCC-FCC
transition. With our data and using Eq. 7, we find a transition at $N_{c}=7$.
The small difference with $N_{c}=17$ in Tab. 2 can be explained by the fact
that an additional approximation is made in Ref. 6 to derive the Eq. 7.

It is also noteworthy that the $\delta(N)$ curve provided in Ref. 6 does not
seem physically correct. In fact, it should mainly obey a $N^{-1/3}$
variation. The best fit of $\delta(N)$ with the expression $\mu N^{-1/3}$ is
represented as a dotted blue line in Fig. 4, with $\mu=0.6646$. The agreement
is obviously not satisfactory. This cast some doubts about the validity of
these $\delta(N)$ data. In addition, according to Eq. 9 one can write

$\mu=\frac{\beta(\mathrm{FCC})}{E_{c}(\mathrm{FCC})}-\frac{\beta(\mathrm{BCC})}{E_{c}(\mathrm{BCC})}.$
(10)

This expression can be employed to compute $\beta(FCC)$ assuming that
$\beta(BCC)=7.149$ eV, $E_{c}(\mathrm{BCC})=-8.90$ eV,
$E_{c}(\mathrm{FCC})=-8.406$ eV, and $\mu=0.6646$. We find $\beta(FCC)=1.166$
eV, which corresponds to $\gamma=0.601$ J m-2 (see the Appendix). Such a
surface energy value is too low for W to be physically meaningful. This
confirms that the main source of discrepancy between our calculations and
predictions made in Ref. 6 is the $\delta(N)$ curve. Since the latter is not a
material dependent quantity, this unfortunately questions the validity of
their predictions for other metals.

##  V Conclusions

In this paper we report investigations on the structure and stability of
tungsten nanoparticles, based on first-principles DFT calculations. In
particular, various nanoparticles with BCC, FCC, A15, disordered and
icosahedral structures, are considered. These models include 55 to 169 atoms,
equivalent to sizes of about 1–1.6 nm. Our first conclusion is that BCC
nanoparticles are the most stable energetically, followed by A15, disordered,
icosahedral, and FCC in this order. Variations of nanoparticle energy as a
function of size reveal a BCC-disordered transition at small sizes (at 40
atoms), and no BCC-FFC transition. It contradicts an earlier benchmark
theoretical study on transition metal nanoparticles [6]. Our investigations
suggest that the discrepancy could be explained by an erroneous description of
the cluster geometry in this study. Finally, we also analyze the proposed
semi-empirical model based on the second moment approximation, and improve its
accuracy by adjusting the model parameters with respect to our DFT data. This
will be useful for accurately calculating the energy of tungsten nanoparticles
for all investigated structures, at much larger sizes and without the need of
explicit atomistic calculations.

As it stands, our calculations do not explain why A15 and FCC tungsten
nanoparticles were observed [10, 7], and not only BCC ones [11]. The most
likely rationale is that these structures can be formed in specific growth
conditions, with a possible influence of substrate, impurities, or
surfactants. Although they are metastable compared to BCC, once formed they
are probably stable at ambient conditions or during moderate annealings.

In perspective to this study, in the light of the results presented here, it
seems worthwhile to perform DFT calculations of the structure and stability of
nanoparticles made of other transition metals, in particular those for which a
BCC-FCC transition was predicted like Mo, Ta, Nb, Cr and V.

###### Acknowledgements.

Computer time for this study was partially provided by the MCIA (Mésocentre de
Calcul Intensif Aquitain). This work pertains to the French Government program
“Investissements d’Avenir” (EUR INTREE, reference ANR-18-EURE-0010, and LABEX
INTERACTIFS, reference ANR-11-LABX-0017-01).

## Author Declarations

### Conflict of interest

The authors have no conflicts to disclose.

### Author contributions

Laurent Pizzagalli: Conceptualization (lead); Investigation (lead); Writing –
original draft (lead); Writing – review & editing (equal). Sandrine Brochard:
Conceptualization (supporting); Writing – review & editing (equal). Julien
Durinck: Conceptualization (supporting); Writing – review & editing (equal).
Julien Godet: Conceptualization (supporting); Writing – review & editing
(equal).

## Data Availability Statement

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

*

##  Appendix A Spherical nanoparticle model

A general expression of the energy $E$ of a nanoparticle of $N$ atoms is

$E=N\varepsilon_{0}+\gamma S(N),$ (11)

with $S(N)$ the surface of the nanoparticle, $\gamma$ the surface energy, and
$\varepsilon_{0}$ the energy per atom of the corresponding bulk phase. This is
equivalent to Eq. 1 if $E_{s}(N)=\gamma S(N)$. Assuming a spherical shape for
the nanoparticle and a radius $r$, the nanoparticle surface is

$S=4\pi r^{2},$ (12)

and the nanoparticle volume is

$V=\frac{4\pi}{3}r^{3}=Nv_{0},$ (13)

with $v_{0}$ the volume of one atom. Combining Eqs. 12 and 13, we obtain

$r=\left(\frac{3v_{0}N}{4\pi}\right)^{1/3},$ (14)

and

$S=4\pi\left(\frac{3v_{0}}{4\pi}\right)^{2/3}N^{2/3}.$ (15)

Finally, Eq. 11 can be written

$E=N\varepsilon_{0}+4\gamma\pi\left(\frac{3v_{0}}{4\pi}\right)^{2/3}N^{2/3},$
(16)

and

$E/N=\varepsilon_{0}+4\gamma\pi\left(\frac{3v_{0}}{4\pi}\right)^{2/3}N^{-1/3}.$
(17)

In comparison with Eq. 3, we finally get

$\gamma=\frac{\beta}{4\pi}\left(\frac{4\pi}{3v_{0}}\right)^{2/3}.$ (18)

$v_{0}$ is taken to be equal to the bulk atomic volume, assuming that the
nanoparticle relaxation is small, and that $v_{0}$ is the same for all
nanoparticle atoms. $v_{0}$ is then easily calculated for each crystalline
structure using DFT determined $a_{0}$ (given in Tab. 1). For the disordered D
and icosahedral I nanoparticles, we use the same $v_{0}$ than for the FCC
structure.

## References

## References

  * Rubel  _et al._ [2018] M. Rubel, A. Widdowson, J. Grzonka, E. Fortuna-Zalesna, S. Moon, P. Petersson, N. Ashikawa, N. Asakura, D. Hamaguchi, Y. Hatano, K. Isobe, S. Masuzaki, H. Kurotaki, Y. Oya, M. Oyaidzu, and M. Tokitani, Dust generation in tokamaks: Overview of beryllium and tungsten dust characterisation in jet with the iter-like wall, Fusion Engineering and Design 136, 579 (2018), special Issue: Proceedings of the 13th International Symposium on Fusion Nuclear Technology (ISFNT-13).
  * Acsente  _et al._ [2020] T. Acsente, L. Gabriela Carpen, E. Matei, B. Bita, R. Negrea, E. Bernard, C. Grisolia, and G. Dinescu, Tungsten nanoparticles produced by magnetron sputtering gas aggregation: Process characterization and particle properties, in  _Progress in Fine Particle Plasmas_ (IntechOpen, 2020).
  * Gu  _et al._ [2022] C. Gu, Z. Wang, Y. Pan, S. Zhu, and Z. Gu, Tungsten-based nanomaterials in the biomedical field: A bibliometric analysis of research progress and prospects, Advanced Materials 35, 10.1002/adma.202204397 (2022).
  * Pai  _et al._ [2012] C.-F. Pai, L. Liu, Y. Li, H. W. Tseng, D. C. Ralph, and R. A. Buhrman, Spin transfer torque devices utilizing the giant spin Hall effect of tungsten, Applied Physics Letters 101, 122404 (2012).
  * Chopra  _et al._ [2004] K. L. Chopra, M. R. Randlett, and R. H. Duff, FACE-CENTERED-CUBIC TUNGSTEN FILMS OBTAINED BY, Applied Physics Letters 9, 402 (2004).
  * Tománek  _et al._ [1983] D. Tománek, S. Mukherjee, and K. H. Bennemann, Simple theory for the electronic and atomic structure of small clusters, Phys. Rev. B 28, 665 (1983).
  * Oh  _et al._ [1999] S. J. Oh, S. H. Huh, H. K. Kim, J. W. Park, and G. H. Lee, Structural evolution of w nano clusters with increasing cluster size, J. Chem. Phys. 111, 7402 (1999).
  * Huh  _et al._ [2000] S. H. Huh, H. K. Kim, J. W. Park, and G. H. Lee, Critical cluster size of metallic cr and mo nanoclusters, Phys. Rev. B 62, 2937 (2000).
  * Vystavel  _et al._ [2005] T. Vystavel, S. A. Koch, G. Palasantzas, and J. T. M. De Hosson, Structural dynamics of gas-phase molybdenum nanoclusters: A transmission electron microscopy study, Appl. Phys. Lett. 86, 10.1063/1.1886261 (2005).
  * Iwama and Hayakawa [1985] S. Iwama and K. Hayakawa, Preparation of ultrafine mo and w particles by the gas evaporation technique with electron beam heating, Surf. Sci. 156, 85 (1985).
  * Schöttle  _et al._ [2014] C. Schöttle, P. Bockstaller, D. Gerthsen, and C. Feldmann, Tungsten nanoparticles from liquid-ammonia-based synthesis, ChemComm 50, 4547 (2014).
  * Chen  _et al._ [2017] L. Chen, J. L. Fan, and H. R. Gong, Phase transition and mechanical properties of tungsten nanomaterials from molecular dynamic simulation, J. Nanopart. Res. 19, 10.1007/s11051-017-3812-z (2017).
  * Marville and Andreoni [1987] L. Marville and W. Andreoni, Size dependence of the structural properties of transition-metal aggregates from an empirical interatomic potential scheme, J. Phys. Chem. 91, 2645 (1987).
  * Lassner and Schubert [1999] E. Lassner and W.-D. Schubert,  _Tungsten: Properties, Chemistry, Technology of the Element, Alloys, and Chemical Compounds_ (Springer US, 1999).
  * Becquart and Domain [2006] C. Becquart and C. Domain, Migration energy of he in w revisited by ab initio calculations, Phys. Rev. Lett. 97, 196402 (2006).
  * Heinola and Ahlgren [2010] K. Heinola and T. Ahlgren, Diffusion of hydrogen in bcc tungsten studied with first principle calculations, J. Appl. Phys. 107, 10.1063/1.3386515 (2010).
  * Lowrie and Gonas [1967] R. Lowrie and A. M. Gonas, Single-crystal elastic properties of tungsten from 24° to 1800°c, Journal of Applied Physics 38, 4505 (1967).
  * Jiang  _et al._ [2020] D. Jiang, S. Zhong, W. Xiao, D. Liu, M. Wu, and S. Liu, Structural, mechanical, electronic, and thermodynamic properties of pure tungsten metal under different pressures: A first-principles study, International Journal of Quantum Chemistry 120, 10.1002/qua.26231 (2020).
  * Tran  _et al._ [2016] R. Tran, Z. Xu, B. Radhakrishnan, D. Winston, W. Sun, K. A. Persson, and S. P. Ong, Surface energies of elemental crystals, Scientific Data 3, 10.1038/sdata.2016.80 (2016).
  * Tyson and Miller [1977] W. Tyson and W. Miller, Surface free energies of solid metals: Estimation from liquid surface tension measurements, Surface Science 62, 267 (1977).
  * Boer  _et al._ [1989] F. R. d. Boer, R. Boom, W. C. M. Mattens, A. R. Miedema, and A. K. Niessen,  _Cohesion in metals : transition metal alloys_ (North-Holland Amsterdam, Amsterdam, 1989).
  * Einarsdotter  _et al._ [1997] K. Einarsdotter, B. Sadigh, G. Grimvall, and V. Ozoliņš, Phonon instabilities in fcc and bcc tungsten, Physical Review Letters 79, 2073 (1997).
  * de Jong  _et al._ [2015] M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag, J. J. Plata, C. Toher, S. Curtarolo, G. Ceder, K. A. Persson, and M. Asta, Charting the complete elastic properties of inorganic crystalline compounds, Scientific Data 2, 10.1038/sdata.2015.9 (2015).
  * Barmak  _et al._ [2017] K. Barmak, J. Liu, L. Harlan, P. Xiao, J. Duncan, and G. Henkelman, Transformation of topologically close-packed $\beta$-w to body-centered cubic $\alpha$-w: Comparison of experiments and computations, J. Chem. Phys. 147, 10.1063/1.4995261 (2017).
  * 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, Advanced capabilities for materials modelling with quantum ESPRESSO, J. Phys.: Condens. Matter 29, 465901 (2017).
  * Perdew  _et al._ [1996] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl [1994] P. E. Blöchl, Projector augmented-wave method, Physical Review B 50, 17953 (1994).
  * Marzari  _et al._ [1999] N. Marzari, D. Vanderbilt, A. De Vita, and M. C. Payne, Thermal contraction and disordering of the al(110) surface, Physical Review Letters 82, 3296 (1999).
  * Amodeo and Pizzagalli [2021] J. Amodeo and L. Pizzagalli, Modeling the mechanical properties of nanoparticles: a review, C.R. Phys. 22, 1 (2021).
  * Jana and Caro [2023] R. Jana and M. A. Caro, Searching for iron nanoparticles with a general-purpose gaussian approximation potential, Phys. Rev. B 107, 245421 (2023).
  * Besley  _et al._ [1995] N. A. Besley, R. L. Johnston, A. J. Stace, and J. Uppenbrink, Theoretical study of the structures and stabilities of iron clusters, J. Mol. Struct: THEOCHEM 341, 75 (1995).
  * Swaminarayan  _et al._ [1994] S. Swaminarayan, R. Najafabadi, and D. Srolovitz, Polycrystalline surface properties from spherical crystallites: Ag, au, cu and pt, Surf. Sci. 306, 367 (1994).
  * Abdul-Hafidh [2022] E. H. Abdul-Hafidh, Surface free energy and structural transition of tungsten nanosolid, J. Nanopart. Res. 24, 10.1007/s11051-022-05638-6 (2022).
  * Ozolins [2009] V. Ozolins, First-principles calculations of free energies of unstable phases: The case of fcc w, Phys. Rev. Lett. 102, 065702 (2009).
  * Kittel [1976] C. Kittel,  _Introduction to Solid State Physics_ (John Wiley & Sons, inc., New York, 1976).