\tolerance=10000
\documentstyle[preprint,aps]{revtex}
\begin{document}
\draft

\title{Free-energies of point defects in sodium from first-principles}

\author{Enrico Smargiassi\cite{present} and Paul A. Madden}

\address{Physical Chemistry Laboratory,
University of Oxford, Oxford OX1 3QZ (GB)}

\maketitle

\begin{abstract}
We calculate the formation free-energies of the major point defects in
sodium using a fully {\sl ab-initio}, self-consistent scheme. Our results
are in good agreement with existing theoretical and experimental values,
whenever available. We find that the vacancy has the lowest free-energy of
formation at all temperatures except possibly close to the melting point,
but that vacancies experience an attractive interaction and may stay bound to
each other, forming a divacancy, which might be expected to be stable
for a significant period of time. Contrary
to previous expectations we find that at temperatures close to the melting
point the interstitial may be a significant presence.
\end{abstract}
\pacs{PACS numbers: 61.72.Ji}

\section{Introduction}

The field of defects and self-diffusion in bcc metals has attracted
considerable attention in the past 20 years. This is partly due to the
anomalous behaviour that many of these materials display in the
self-diffusion curve at high temperature\cite{kohler88}. At temperatures
close to the melting point, the Arrhenius plot of the diffusion coefficient
bends upwards, suggesting either a temperature dependence of the diffusion
mechanism, or the simultaneous presence of more than one mechanism. 
  Although
most of the interest has focused on the technologically important
transition metals like Zr or V, the study of simple metals having a bcc
structure is important from a conceptual viewpoint. For example, it is
also believed that the anomalous has a common origin with the martensitic
transition undergone by many of these materials (in itself an intriguing
issue) and in particular that an important contribution to both processes
comes from the ${\bf q}=(2/3,2/3,2/3)$ phonon  which is characteristically
 at very low frequency in the bcc phase.
 The very simplicity of
the electronic structure of simple metals allows one to disentangle
geometrical factors from those depending on the electronic
structure\cite{kohler88}. In spite of a considerable body of theoretical,
numerical and experimental work, many questions remain unanswered.
\par
Molecular-dynamics (MD) studies of defect properties are a major tool for
the study of dynamical properties of matter. In the past, MD studies on
simple metals have been conducted using pair potentials (notably of the
Dagens-Rasolt type\cite{rasolt73,dafano77,vandoan87}) which however are not
completely reliable for defect properties. For example, the activation
energy for vacancy-induced self-diffusion in Na was given\cite{dafano77}
quite correctly as 0.37 eV, but this was a result of a cancellation of
errors between a formation energy too low by 30\% (0.25 eV vs 0.35 eV) and
an activation energy four times too large (0.12 eV vs 0.03 eV). It has been
argued\cite{jacucci81} that the agreement between the calculated activation
energies and the experimental ones (which is found also in Al) is not
accidental, but is a consequence of a systematic cancellation between
many-body forces during the process of creating and subsequently moving the
defect. However, the exactness of the separate formation and migration
energies is of importance. First, the density of defects, a function of the
formation free-energy, is important in itself. Second, an incorrect
migration energy may lead to spurious predictions, and it is difficult to
separate these from the genuine ones. Unexpected phenomena have been
suggested by a few calculations\cite{dafano77,vandoan87} and, in the case
of ref.\cite{dafano77}, not been confirmed experimentally\cite{brunger80}.
\par
  The
atomic environment around a defect is very different from the bulk and thus
a successful description of a defect requires a tailor-made potential. This
new potential however is likely to be inadequate for the description of
bulk properties, like phonons (see the discussion in ref.\cite{jacucci81});
this in turn affects the accuracy of the calculation of the free-energy of
formation of the defect. In other words, in calculations on point defects
the transferability problems that affect classical potentials are
highlighted. A consistent description of all the defect properties seems to
require some sort of explicit treatment of the electronic degrees of
freedom.
\par
Nowadays, {\sl ab-initio} techniques like the Car-Parrinello
method\cite{carpar85} are available, which explicitly and self-consistently
deal with the electronic side of the problem, hence avoiding these
difficulties. Application of these ``exact'' schemes to metals (especially
in the solid phase) on the other hand has been held back by their
extremely large computational requirements. For example, in a
single-particle, Kohn-Sham (KS) representation, reasonable results 
 for simple metals can only
be obtained by using a k-point sampling of the Brillouin zone corresponding
to at least 30 points in the irreducible zone\cite{breier}. A reasonable
description of a point defect however requires typically 50 atoms, which
means that more than 100000 plane waves must be used. Moreover, meaningful
thermal averaging can be done, as we shall see, only on 
 molecular dynamics runs about 5 or 6
ps long. This effectively rules out any finite
temperature simulation of the KS type on the computers now available.
\par
In this paper, we will illustrate the results of the a first-principles MD
study of defect energies and free-energies of defects in Na. We used a new
numerical scheme which gives a first-principles, self-consistent
description of alkali metals at the same level of accuracy, but at a
fraction of the computational cost, of conventional KS schemes. 
 Rather than use a representation of the kinetic energy functional via the introduction of orbitals (as in Kohn-Sham) we use an explicit functional 
 of the density. In the case
of Na we find that the differences between our method, the conventional KS
method and the experiment are very small and probably lie within the basic
uncertainties introduced by the Local Density Approximation plus
pseudopotential scheme. The major advantage of the orbital-free scheme
 is that  we can
easily use cells containing up to 250 atoms (more than 400 for static
calculations) and run dynamics for 100 ps if needed \cite{Foley}.
  This makes our method
ideally suited for defect calculations in alkalis.
\par
The method has been described in detail elsewhere\cite{our}, and so has the
method of calculation of free energies\cite{our1}. In this paper we will
only show and discuss the results for the major point defects: the vacancy,
the divacancy and the interstitial. There are in fact several different
divacancies to be studied; it makes a difference if the two vacant sites
are first neighbours, second neighbours etc.. We will call $D_n$ a
divacancy in which the vacant sites are $n$-th neighbours. We studied the
cases $n=1$ to 4, as they are probably the most relevant in the diffusion
and related processes. Moreover, larger values of $n$ imply a substantial
computational effort (due to the need for larger cells) and would give
results not much different from the case of a couple of isolated vacancies.
As always, the ground-state structure of the interstitial is not a priori
obvious. We performed an explicit study to find the lowest-energy atomic
configurations.
\par
In section II we present the static results regarding formation energies at
the absolute zero and the corresponding electronic and ionic structures for
all the defects. This enables us to assess the validity of our choice by
making systematic comparisons with state-of-the-art KS
calculations\cite{breier}. Section III contains the results for the free
energies and a discussion of the results. Our conclusions are in section
IV.

\section{Energetics of defect formation}

\subsection{Preamble}

We wish first to compare the prediction of our model for the $T=0$ energetics
of defects in Na with those of a recent
accurate study within the KS scheme (ref.~\cite{breier}). In this way we can
check the adequacy of the approximations made in constructing our energy
functional and appropriately choose the convergence parameters. Moreover, as we
are able to use much larger cells than used in the KS calculations \cite{breier}, we may check
for the absence of size effects.
\par
Once  the energy functional has been selected,
one is left with three main approximations
(more details on both the formalism and the parameters used can be found in
ref.\cite{our1}). First, one has to choose a suitable pseudopotential to
describe the electron-ion interaction. We have chosen a Topp-Hopfield
pseudopotential\cite{topp73} which gives a good description of sodium in a
wide variety of circumstances. Second, a plane-wave cutoff must be imposed
in order to have a finite number of electronic degrees of freedom. We have
found that a cutoff corresponding to a plane-wave energy of about 9 Rydberg is
sufficient for all our purposes. Third, the size of the simulation cell is
of importance. An isolated defect will be properly simulated only if the
cell is large enough to make the interaction with the periodic images
negligible. To be sure that this is the case a separate study for each kind
of defect considered is required, especially when ionic relaxation is
allowed. In all of what follows, we consider an energy converged if its
estimated error is less than 0.01 eV, unless otherwise stated. It is doubtful
that a pseudopotential plus local density approximation scheme can give
more accuracy than this.
\par
A formation (free-)energy is calculated as the difference between the energy
of a system with the defect and the one of the perfect crystal, at constant
number of particles. Of course,
the energy functional used and the convergence level
must be the same in both cases. There are however several possible
different choices for the cell volumes. For the perfect crystal, we used
the experimental lattice parameters at all temperatures (except at
$T=100~K$ where we used the same as $T=0$), which is the most common and
sensible choice. Because of the use of the LDA and the choice of pseudopotential
 our scheme overestimates of the lattice parameter of Na
by 2\%; this means that the system is under slight pressure, which
  we verified is always less than 7 Kbar at any temperature.
\par
Once the lattice parameter for the perfect crystal has been chosen, one has
to decide the value appropriate for the defective crystal. A ``constant
lattice parameter'' calculation is one which keeps the cell volume equal to
that of the perfect crystal. The unrelaxed calculations of ref.
\cite{breier} are of this kind. We choose to perform ``constant atomic
volume'' calculations instead, i.e. to scale the volume with the number of
atoms actually present in the cell. One might improve by allowing the
volume of the cell to relax until the pressure is equal to the pressure
experienced by the perfect crystal to which it is compared; this would
allow one to calculate the formation volume of the vacancy. In Ref.
\cite{breier} however it was found that in sodium the change in energy
caused by this relaxation is negligible, and our results confirm this: the
pressure difference between the perfect crystal and the defective crystals
with constant atomic volume never exceeds 0.5 Kbar, suggesting very little
change in the energy upon volume relaxation.

\subsection{Vacancy and divacancy}

For an isolated vacancy with fixed ions, we found that the formation energy
$E_V$ is already converged at a cell size of 16 atoms. Similar conclusions
were reached by Frank {\sl et al.} \cite{breier,frank93} in their studies
of lithium and sodium. The increase in the cell size therefore brings
changes in the total energy only via ionic relaxation. Even so, a 54 atom
cell was found to give convergence for the vacancy and all the divacancies,
as seen in table \ref{tabvac}. We feel that most of the small variations
that the energies show at increasing $N$ are due to numerical effects
rather than physical changes with the cell. The results agree very well
with the available experimental \cite{landolt} and theoretical
\cite{breier} results.
\par
In fig. \ref{figvac} we show the magnitude of the relaxation of the first
shells of neighbours of the vacancy for two cell sizes (54 and 128 atoms
per cell). The relaxation turns out to be symmetric around the vacancy,
except for some of the shells which are incompletely included in the cell.
However, whenever the cell is enlarged so as to include the full shells the
asymmetric part of the relaxation disappears. The asymmetric part of the
relaxation is therefore
a finite-size effect and it has been excluded from fig. \ref{figvac}.
Again, we see that the convergence is good with a 54 atom cell. The results
are also within a few percent of the Kohn-Sham results of
ref.\cite{breier}.
\par
In Fig. (\ref{figdensdiv}) we show the electron density around $D_1$. The
spherical ``holes'' in the electron density are an artifact due to the use
of a pseudopotential to represent the electron-ion interaction. The
bacterium-like hole in the centre is instead real and it corresponds to the
alteration in the electronic density due to the presence of the divacancy.

\subsection{Interstitial}

Close-packed solids like metals usually do not allow the interstitial to
exist in a pure form (by pure we mean an atom sitting in an interstice of
the crystal, while the others sit in a recognisable distortion of the
perfect lattice). Indeed, this is not the case even for several
tetrahedrally coordinated materials, whose structure leaves much more room.
Instead, it usually happens that an atom close to the interstitial is
pushed away from its perfect lattice position. The defect thus is best
described as a pair of atoms occupying a single lattice site, with the
other atoms still quite close to their original sites. This kind of defect
is called a ``dumbbell'' and it might be thought as a molecule embedded in
a crystal. As such, the bond length of the pair may be expected to be of
the order of the bond length of the molecule {\sl in vacuo} (5.8 a.u. for
Na). A dumbbell is characterized by the direction of the bond between the
atoms in question. In a cubic crystal one typically considers the (100),
(110), and (111) dumbbells. However, other kind of interstitial configurations cannot be discarded {\sl a priori}.
\par
We studied the three dumbbells mentioned above, plus two configurations in
which a pure interstitial was either in a tetrahedral or an octahedral
site. In all cases all atoms were allowed to relax using a steepest-descent
procedure until convergence in the energies (an estimated error of less
than $10^{-4}$ hartree) was reached. As already found in other studies
\cite{breier}, we found that the change in energy between the initial,
unrelaxed configuration and the final one is massive (as large as 2 or 3
eV), in fact much greater than the relaxed formation energy itself. This
shows that the ``unrelaxed'' structures of the interstitial are unphysical
and the convergence has to be studied on the relaxed energies directly.
>From table \ref{tabint} we argue that, due to the small differences in
energy between some of the competing structures, a 128 atom cell is
advisable, although a 54 atom cell is probably sufficiently large for many
purposes.
\par
Our for the formation energies  agree well with KS results of ref. \cite{breier}, although ours
look systematically lower. We think that at least part of the difference
may be explained by a different degree of convergence in the ionic
relaxation: as we showed, the ionic relaxation energy is very large and a
slightly imperfect convergence may cause sizeable errors. Indeed, these are
almost to be expected if the relaxation is performed by moving the atoms
``by hand''.
\par
As table \ref{tabint} shows, when a 54 atom cell is used, the (110) and
(111) dumbbells turn out to be energetically almost equivalent. The bond
length between the two atoms is similar too: 5.78 and 5.74 a.u.,
respectively (to be compared with KS results of 5.51 and 5.47 a.u.
\cite{breier}). However when a 128 atom cell is used the (111)
configuration becomes somewhat lower in energy (the bond length increases
to 5.83 a.u., that of the (110) dumbbell remains almost the same). Whether
this is significant or not remains to be seen: a change of about 0.02 eV in
the formation energies as the cell becomes larger may be expected.
Moreover, with such a small energy difference it is almost certain that at
finite temperature the two forms are both present simultaneously. 
\par
In order to check that no other configuration can have a lower energy than
those studied above, we used simulated annealing. We put the extra atom in a
tetrahedral site, heated the system to 1500 K, and cooled it down
gradually ($\Delta T/T \approx 0.01$/step) until the instantaneous
temperature was of the order of 1 K, then it was quenched until converged.
At the end of the run it was found that the atoms forming the dumbbell did
not include the original interstitial, which shows that the system had lost
memory of the original configuration. In the 54 atom cell we found that the
ground-state was a (110) dumbbell (bond length 5.79 a.u.), with a formation
energy almost equal to that given by the relaxation study. Of course the
annealing cannot be expected to discriminate between this structure and the
almost degenerate (111) dumbbell. In a 128 atom cell the ground state
configuration was found to be the (111) dumbbell, again with the same
formation energy of the relaxation study. We conclude that the static
calculations are indeed sufficient to study the ground state of the
interstitial.
\par
The electron density along a (100) slice of the crystal crossing a (111)
dumbbell is shown in Fig. (\ref{figdensint}). The contours corresponding to
regions with the lowest electronic density have been omitted, as they only
appear in the pseudopotential core region (where the calculated
density is unphysical
anyway) and would only confuse the figure. The slice crosses the dumbbell
approximately midway between the two atoms in question, the rightmost atom
being slightly higher than the plane of the figure and the leftmost
slightly lower. The buildup of charge between the atoms is clearly
reminiscent of the chemical bond in a diatomic molecule, however there are
differences. In particular the pockets of electrons on the outside of the
dumbbell appear to be an effect of the crystalline environment. It is
interesting to note that in spite of the importance of the ionic relaxation
in determining the interstitial formation energy, the distortion of the
surrounding lattice is almost undetectable to the naked eye.

\section{Free energies}

\subsection{Preliminaries}

\subsubsection{Defect configurations}

As discussed in \cite{our1}, for a defect which can exist in $m$ equivalent
orientations it is convenient to calculate its free energy when in a specific
form, then add a term which we call ``orientational entropy''. This term
is given by $S^{or} = k_B \ln m$ and all one needs to do is to count the
number of equivalent configurations.
\par
The vacancy is fully specified by the site in which it sits, hence $m=1$
and $S^{or}=0$. The divacancy $D_1$ is specified first by the lattice site
of one of the pair of vacant sites (let us suppose it is $(0,0,0)$); then,
once the first site is fixed there are 8 possibility for finding the
second, as it must lie at a lattice site $a(\pm \frac{1}{2},\pm
\frac{1}{2},\pm \frac{1}{2})$ ($a$ is the lattice parameter). However only
half of these must be counted, due to the indistinguishability of the
vacancies. Therefore for $D_1$, $m=4$ and $S^{or} \simeq 1.4$. For the
divacancy $D_2$, the second vacant site is $a(\pm 1, 0, 0)$ plus all the
permutations of axes, giving $m=6/2=3$ and $S^{or} \simeq 1.1~k_B$. For the
divacancy $D_3$, one has $a(\pm 1,\pm 1, 0)$ (plus permutations), so
$m=12/2=6$ and $S^{or} \simeq 1.8~k_B$. Finally, for $D_4$, one has $a(\pm
\frac{3}{2},\pm \frac{1}{2},\pm \frac{1}{2})$ and $m=24/2=12$ and $S^{or}
\simeq 2.5~k_B$.
\par
In principle, all the possible configurations of the interstitial should be
studied. We have seen from the energetic considerations that only two of
them can be of importance, the (110) and (111) dumbbells. To simplify
matters, we decided to take the (111) dumbbell as representative of both.
In this case, we consider the (110) and (111) dumbbells to be equivalent,
hence $m=6+4=10$ and $S^{or} \simeq 2.3~k_B$.

\subsubsection{$\lambda$ integration}

The technique for calculating the formation free-energy $F^N_f$ of a defect
has been described in detail in ref. \cite{our1}. We summarize here quickly
the relevant formulae for the vacancy and the divacancy and derive those
for the interstitial. For a vacancy or a divacancy in an $N$ site cell:
\begin{equation}
   F^N_f = - \int^1_0  {\rm d}\lambda \, \left\langle \Phi^N (\lambda)
             \right\rangle_\Omega
           + L \Omega_a \int^1_0  {\rm d}\lambda \, P (\lambda) + L \Delta F,
  \label{effelambda}
\end{equation}
where $\Phi^N = \frac{\partial U^N}{\partial \lambda}|_\Omega$, $\Omega_a$
the atomic volume, $P$ is the thermodynamical pressure exerted by the
system, and $L$ is the number of atoms to be decoupled to form the defect
in question. The cell volume $\Omega$ is taken to be a linear function of
$\lambda$. The potential function $U^N=U^N(\lambda)$ is parametrized so
that $U^N(\lambda=1)=$ ionic potential of the perfect crystal, and
$U^N(\lambda=0)=$ is the 
 ionic potential of the crystal with a vacancy (if $L=1$)
or divacancy (if $L=2$) plus the potential energy of
  $L$ atoms decoupled from the others and put in
$L$ harmonic wells $u$. These wells are centered in the vacant sites and
switched off as $\lambda$ approaches 1. The
oscillation frequency of an atom in $u$ is
chosen to be the Debye frequency for sodium: $\omega_D\approx
4.75\,10^{-4}$ a.u.. $\Delta F$ is the difference between the free-energy
per atom in the perfect crystal and the free-energy of an atom in a
harmonic well, at the same temperature: $\Delta F = \frac{1}{N} F^N_C -
\varphi$, where $\varphi$ is the free-energy of an atom acted
upon by $u$. The values of $\Delta F$ at various temperatures have been
calculated previously and published in ref.~\cite{our1}.
\par
The case of the interstitial is analogous. In this case however the defect
is best described as two atoms (the dumbbell) gradually appearing in place
of a regular lattice atom, as $\lambda$ approaches 1. The simplest approach
is then to take the vacancy crystal and not the perfect crystal as the
reference state. In this case, the atoms appear in the site of the vacancy,
and
\begin{eqnarray}
  U^{N+1} (\lambda=1) & = & U^{N+1}_I ,
\end{eqnarray}
and
\begin{eqnarray}
  U^{N+1} (\lambda=0) & = & U^{N-1}_V + 2 u,
\end{eqnarray}
where $U_I$ and $U_V$ refer to the potential energy of the interstitial
crystal and a vacancy crystal respectively. We then have
\begin{equation}
  \int_0^1 \frac{\partial F^{N+1}}{\partial \lambda} {\rm d}\lambda =
    F^{N+1} (\lambda=1) - F^{N+1} (\lambda=0) =
    F^{N+1}_I - F^{N-1}_V - 2 \varphi.
\end{equation}
In this case, the formation free-energy of the interstitial $F^N_{fI}$ is
calculated relative to that of the vacancy $F^N_{fV}$ and Eq.
\ref{effelambda} has to be suitably modified. By definition:-
\begin{eqnarray}
  F^N_{fI} & = & F^{N+1}_I - \frac{N+1}{N} F^N_C
             =   F^{N+1}_I - F^N_C - \frac{F^N_C}{N} \label{fnfi},
\end{eqnarray}
and
\begin{eqnarray}
  F^N_{fV} & = & F^{N-1}_V - \frac{N-1}{N} F^N_C
             =   F^{N-1}_V - F^N_C + \frac{F^N_C}{N}.\label{fnfv}
\end{eqnarray}
Subtracting eq. (\ref{fnfv}) from eq. (\ref{fnfi}) we have, in analogy with
eq. (\ref{effelambda}) \cite{our1},
\begin{equation}
   F^N_{fI} = F^N_{fV} + \int^1_0  {\rm d}\lambda \,
              \left\langle \Phi (\lambda) \right\rangle_\Omega
            - L \Omega_a \int^1_0  {\rm d}\lambda \, P (\lambda)
            - L \Delta F.
  \label{effelambdint}
\end{equation}
\par
In the case of the interstitial it is not obvious {\sl a priori} where the
dumbbell should appear, i.e. where the harmonic wells $u$ have their
centre. Considerations of consistency with the $T=0$ results show that they
must be centered at the place where the dumbbell is found to sit in a
static calculation after full relaxation of the crystal.
\par
A further set of parameters of importance for the free-energy calculations
is then the set of $\lambda$ points used for the calculation of the
integrals above. We checked that two Gauss-Legendre points ($\lambda_1^2 =
0.2113248654$ and $\lambda_2^2 = 1 - \lambda_1^2 = 0.7886751346$, both with
weight 1/2) are sufficient to give convergence for the vacancy and all the
divacancies. The interstitial however causes a much larger ionic relaxation
and this requires the use of three points: $\lambda_1^3 = 0.1127016654$,
$\lambda_3^3 = 1 - \lambda_1^3 = 0.8872983446$, both with weight 5/18, and
$\lambda_2^3 = 0.5 $, with weight 8/18. The integrals over the pressure  
 gave the desired convergence with just one or (for the
interstitial) two $\lambda$ points.
\par
The molecular dynamics 
 run lengths were chosen to be about 6 ps, including thermalization. To
check that this is sufficient for the desired statistical convergence,
after each run we calculated the averages by discarding sections of
different length from the beginning of the run. We found that discarding
the first 2 or 4 ps (i.e. retaining 4 or 2 ps for thermal averaging,
respectively) leads to results very similar to each other.
 However, discarding 1 or 5 ps leads to significant differences.
Allowing 2 ps of thermalization and 2 ps for thermal averaging
therefore seems necessary and sufficient. Test runs as long as 12 ps did not
change this conclusion.

\subsection{Results}

\subsubsection{Vacancy and divacancies}

The free-energy of formation $F_{fV} - T S^{or}_V$ of the vacancy is shown
in Fig. \ref{figfv}, as a function of the temperature. Although from the
slope of the curve it is in principle possible to calculate the formation
entropy, in practice the statistical uncertainty prevents this, even if
more $T$ points had been calculated. All that can be said is that $S_f$ is
of order 1 or 2 $k_B$ at low temperature and probably increases to about 4
or 5 $k_B$ in the vicinity of the melting point. One conspicuous
characteristic of Fig. (\ref{figfv}), repeated in almost all the
calculations, is the increase in the statistical fluctuation of the results
at $T=370~K$, and the bending downwards of the $F(T)$ curve, indicating a
definite increase in entropy. We interpret this as a sign that the melting
is approaching ($370~K$ is the experimental melting temperature of sodium).
Work is in progress to locate the melting point within our model.
\par
Figs \ref{figfdiv} (a,b,c,d) show our results for the various forms of the
divacancy. Again the fluctuations increase sharply at $370~K$. The
formation entropies are 4 -- 6 $k_B$ at low temperatures, and increase
probably to about 8 $k_B$ close to the melting point. The energy
differences between the various forms of the divacancy are small but, in
our opinion, significant. In particular, the conclusion that $D_2$ and
$D_4$ are the forms with lowest free-energy is in agreement with the
experimental findings of Ref. \cite{brunger80}.

\subsubsection{The interaction between vacancies}

We may regard a point defect as a kind of quasiparticle, dressed by a cloud
of virtual phonons (the ionic relaxation), somewhat analogous to a polaron
in an ionic crystal. The defects can interact with each other via the
lattice and with line and surface defects, and can annihilate themselves by
recombining with another defect or with a surface. An interesting point to
study is whether two vacancies repel or attract, or, equivalently, whether
the divacancy is bound or not. In this picture the role of the interaction
potential is played by the difference $F_{int} = F_{fn} - 2 F_{fV}$. The
orientational entropy plays the role of the quasiparticle entropy.
To see this, one may just write the equation for the number
concentration $c_n$ of two $n$-th neighbour vacancies in the case of
$F_{int} = 0$:
\begin{equation}
  c_n = e^{(T S^{or} - F_{fn})/k_B T} = m e^{-2 F_{fV}/k_B T} = m c_V^2,
  \label{vacinteract}
\end{equation}
where $m$ is the number of possible ways in which two $n$-th neighbour
vacancies can be found. Therefore, $m=e^{S^{or}/k_B}$ is the
configuration space available to the motion of the two vacancies. 
\par
The results for $F_{int}$ are shown in Table \ref{tabinteract}. The
comparatively large size of the various uncertainties with respect to these
values make any quantitative conclusion dubious, in particular, the results
at $370~K$ are probably unreliable and are not shown. We believe however
that the qualitative picture is quite clear: at all temperatures there
exists at least one divacancy configuration which is stable with respect to
 dissociation into a pair of separate vacancies. This implies that
vacancies in sodium tend to lump together. The small concentration 
however (about $10^{-5}$ per site at 300 K) means that there is much more
phase space volume available for isolated vacancies than for bound ones and
at equilibrium is likely that the isolated vacancies outnumber the
divacancies.
\par
At large distances the vacancies do not interact, hence
we must have $\lim_{n\to \infty} F_{fn} = 2 F_{fV}$. This limit is already
reached for $n=3$ in the static calculations if the ionic relaxation is not
allowed, in agreement with the findings of the convergence study for the
single vacancy. The difference between $F_3$ and $F_4$ is thus entirely due
to the ionic relaxation, and so it is the interaction between two well
separated vacancies. The data in tables \ref{tabvac} and \ref{tabinteract}
suggest that the interaction between two vacancies is nonmonotonic,
oscillating around zero and becomes negligible only when they are as far
apart as to be sixth or seventh neighbours.

\subsubsection{Interstitial}

Figure \ref{figfi} shows the formation free-energy of the (111) dumbbell
with $S^{or}=k_B \ln 10$. Arguments have been presented 
  in the past (see e.g.
\cite{brunger80} and references therein) in favour of a negative formation
energy for the interstitial, as its presence
  compresses the lattice and  might be
expected to lead to higher phonon frequencies and thus lower entropy. Our
data do not support this view. Indeed, the formation free-energy drops
rapidly with $T$, suggesting a quite large, positive formation entropy,
possibly 6 or 7 $k_B$. At $T=370~K$, $F_{fI}$ is close to zero; although
this is not conclusive, due to the large statistical fluctuations at that
temperature, it clearly suggests that the interstitial may play a
significant role close to the melting point. Indeed, one may suggest a view
of melting in Na, if not an actual mechanism, as the destabilization of the
crystal lattice towards the formation of interstitials.

\section{Conclusions}

We have shown the results of the first fully {\sl ab initio},
self-consistent calculation of defect free-energies in a metal. The method
we used allows one to reach easily the large cell sizes and long simulation
times required by this kind of problems: all the calculations were
performed on a Convex 200 and are confortably within reach of a powerful
workstation.
\par
Our results are in good agreement with previous calculations at $T=0$, and
represent a significant extension of them. In particular, it has been shown
that the role of interstitials at equilibrium cannot be dismissed
as static calculations suggested. Indeed, the increase in the number of
interstitials might be at the root of the increase in the self-diffusion
coefficient of sodium close to melting. To resolve such questions, a dynamical
calculation of the defect-induced self-diffusion for the various mechanisms
should be performed. Work is in progress along these lines.
\par
Experimental evidence of simultaneous contributions from vacancies and
divacancies to the self-diffusion find an interpretation as the presence of
a ``vacancy gas'' in equilibrium with ``vacancy molecules'' which
continuously bind and dissociate after some time. A quantitative theory
however appears difficult to formulate, as the energy differences are tiny
and one would need very accurate results.

\acknowledgments

We thank Professor M. F\"ahnle for providing us with the preprint
\cite{breier} and for further useful information. The work was supported by
SERC and the MOD (UK) under grant GR/H38256. It has also benefited from
further support from SERC under grant no. GR/H10276.

\begin{references}

\bibitem[*]{present} Present address: Laboratoire de Physique,
Ecole Normale Sup\'erieure de Lyon, 46 All\'ee d' Italie,
69364 Lyon C\'edex 07 (F).

\bibitem{kohler88} U. Kohler and C. Herzig, Phil. Mag. A {\bf 58}, 769
(1988).

\bibitem{rasolt73} M. Rasolt and R. Taylor, Jou. of Phys. F {\bf 3}, 67
(1973); C. Dagens, M. Rasolt and R. Taylor, Phys. Rev. B {\bf 11}, 2726
(1975).

\bibitem{dafano77} A. Da~Fano and G. Jacucci, Phys. Rev. Lett. {\bf 39},
950 (1977).

\bibitem{vandoan87} N. Van~Doan and Y. Adda, Phil. Mag. A {\bf 56}, 269
(1987).

\bibitem{jacucci81} G. Jacucci, R. Taylor, A. Tenenbaum, and N. Van~Doan,
Jou. of Phys. F {\bf 11}, 793 (1981).

\bibitem{brunger80} G. Br\"unger, O. Kanert, and D. Wolf, Phys. Rev. B {\bf
22}, 4247 (1980).

\bibitem{carpar85} R. Car and M. Parrinello, Phys. Rev. Lett. {\bf 55},
2471 (1985).

\bibitem{breier} U. Breier, W. Frank, C. Els\"asser, and M. F\"ahnle,
to be published.

\bibitem{Foley} M. Foley, E. Smargiassi and P.A. Madden, Jou. of Phys. --
 Cond. Matter -- to be published.

\bibitem{our} M. Pearson, E. Smargiassi, and P.~A. Madden,
Jou. of Phys. -- Cond. Matter {\bf 5}, 3321 (1993); E. Smargiassi
and P.~A. Madden, Phys. Rev. B {\bf 49}, 5220 (1994).

\bibitem{our1} E. Smargiassi and P.~A. Madden, submitted for publication.

\bibitem{topp73} W.~C. Topp and J.~J. Hopfield, Phys. Rev. B {\bf 7}, 1295
(1973).

\bibitem{frank93} W. Frank, U. Breier, C. Els\"asser, and M. F\"ahnle,
Phys. Rev. B {\bf 48}, 7676 (1993).

\bibitem{landolt} Landolt-B\"ornstein, {\it Numerical Data
and Functional Relationships in Science and Technology} (1991).

\end{references}

\begin{figure}
\caption[]{Relaxation of the ionic shells around the vacancy. Squares: 54
atom cell. Crosses: 128 atom cell. Open crosses: KS results
(ref.~\cite{breier}). The lines are guidelines for the eye.}
\label{figvac}
\end{figure}

\begin{figure}
\caption[]{Electron density around a $D_1$ divacancy, in a 250 atom cell.
The slice is taken along the (110) direction.}
\label{figdensdiv}
\end{figure}

\begin{figure}
\caption[]{Electron density around a (111) dumbbell, in a 432 atom cell
(see text).}
\label{figdensint}
\end{figure}

\begin{figure}
\caption[]{Formation free-energy of the vacancy, in eV (including
$S^{or}$). The error bars show the statistical uncertainty of the results,
estimated by the spread in results obtained by starting the statistical
averaging at different times along the run.}
\label{figfv}
\end{figure}

\begin{figure}
\caption[]{As Fig. \ref{figfv}, for the four kinds of divacancy: (a),
$D_1$; (b), $D_2$; (c), $D_3$; (d), $D_4$.}
\label{figfdiv}
\end{figure}

\begin{figure}
\caption[]{As Fig. \ref{figfv}, for the interstitial.}
\label{figfi}
\end{figure}

\mediumtext
\begin{table}
\caption[]{Formation energy of the vacancy and divacancy at $T=0$,
including ionic relaxation (unrelaxed energies in brackets). For the
vacancy, KS gives 0.36 eV\cite{breier}; the experiments 0.35 eV
\cite{landolt}.}
\begin{tabular}{cccc}
   $N$   &      54     &     128     &     250     \\
\tableline
 vacancy & 0.37 (0.51) & 0.37 (0.52) & 0.35 (0.51) \\
  $D_1$  & 0.81 (1.04) & 0.80 (1.04) & 0.78 (1.03) \\
  $D_2$  & 0.66 (1.00) & 0.66 (1.00) & 0.65 (0.99) \\
  $D_3$  & 0.75 (1.04) & 0.76 (1.03) & 0.74 (1.03) \\
  $D_4$  & 0.69 (1.03) & 0.70 (1.04) & 0.67 (1.03) \\
\end{tabular}
\label{tabvac}
\end{table}

\begin{table}
\caption[]{Formation energy of the interstitial at $T=0$, in eV. Rightmost
column: KS data (ref. \cite{breier}), obtained using $N=54$. Full ionic
relaxation is allowed.}
\begin{tabular}{c|cccc}
   $N$    &  54  & 128  & 250  & 54 (KS) \\
\tableline
  (100)   & 0.67 & 0.66 &      & 0.73 \\
  (110)   & 0.55 & 0.55 &      & 0.59 \\
  (111)   & 0.55 & 0.52 & 0.50 & 0.60 \\
   oct    & 0.67 & 0.66 &      &      \\
   tet    & 0.59 & 0.59 &      &      \\
 annealed & 0.55 & 0.52 &      &      \\
\end{tabular}
\label{tabint}
\end{table}

\begin{table}
\caption[]{Interaction free-energy $F_{int}$ of two vacancies (eV).}
\begin{tabular}{c|cccc}
 $T~(K)$ &   0   &  100  &  200  &  300  \\
 \tableline
  $D_1$  &  0.06 &  0.04 &  0.04 &  0.00 \\
  $D_2$  & -0.08 & -0.10 & -0.14 & -0.11 \\
  $D_3$  &  0.02 &  0.02 &  0.01 & -0.02 \\
  $D_4$  & -0.04 & -0.05 & -0.04 & -0.02 \\
\end{tabular}
\label{tabinteract}
\end{table}

\end{document}
