\documentclass[twocolumn,showpacs,preprintnumbers,amsmath,amssymb]{revtex4}
\usepackage{graphicx}% Include figure files
\usepackage{dcolumn}% Align table columns on decimal point
\usepackage{bm}% bold math

\begin{document}

\title{Order-disorder transition in a Cd-Ca quasicrystal approximant}

\author{M. Widom}
\author{M. Mihalkovic}
 \altaffiliation[Also at ]{Slovakian address}
\affiliation{
Department of Physics\\
Carnegie Mellon University\\
Pittsburgh, PA  15213
}

\date{\today}

\begin{abstract}
Recent experiments discovered an order-disorder transition occuring at
very low temperatures in large unit cell approximants of the stable
Cd-based binary alloy quasicrystals.  The transition is related to
correlations among orientational degrees of freedom whose separations
are 12~\AA.  The orientational degrees of freedom represent a type of
phason disorder within the unit cells of these large approximants.  We
analyze the interactions between the degrees of freedom using {\em
ab-initio} calculations and derive an equivalent antiferromagnetic
Ising model which shows a similar phase transition.
\end{abstract}

\pacs{Valid PACS appear here}
\maketitle

\section{\label{sec:Intro}Introduction}

Thermodynamically stable binary icosahedral quasicrystals occur in the
alloy systems Cd$_{5.7}$Yb and Cd$_{5.7}$Ca~\cite{Tsai,Guo,Jiang}.  In
each case the phase diagram of the binary alloy contains a 1/1 cubic
approximant of structure type Cd$_6$Yb at a close-by
composition~\cite{Pearson,BinaryPD,Gomez}.  This Cd$_6$Yb structure
can be represented in a conventional simple cubic unit cell with a
lattice parameter of $a=15.7$~\AA. However, the Pearson symbol for
Cd$_6$Yb is cI176 indicating that it is cubic, body-centered, with 176
atomic positions per simple cubic unit cell.  Since 176 is not a
multiple of the basic 7-atom stoichiometric unit (6 Cd and 1 Ca, Y or
Yb), there must be partial occupancy. In fact, the unit cell contains
144 Cd atoms, 24 Ca or Yb atoms and 8 vacancies.

The basic structural unit is shown in Fig.~\ref{fig:struct}.
Vacancies concentrate in the innermost cluster, marked $(a)$, in which
4 of a possible 8 sites are occupied.  The 8 sites belong to a cube,
and the 4 occupied sites belong to one out of the two tetrahedral
subsets of the cube.  Each simple cubic cell of Cd$_6$Yb contains two
such basic icosahedral clusters, one centered at the vertex and the
other centered at the body center.  These structures are closely
related to Cd$_6$Y.cI184, differing only in the number of partially
occupied Cd sites, and the nominal orientations of the cluster center
tetrahedra.

\begin{figure}
\includegraphics{struct}
\caption{\label{fig:struct}The basic icosahedral cluster of the
Cd$_6$Y structure type consists of three concentric shells of
icosahedral symmetry (b-d). The internal cluster (a) is a
tetrahedron. (Figure taken from Physics Today~\cite{PhysicsToday})}
\end{figure}

In order for the overall structure to be truly body-centered, it is
necessary that the tetrahedron at each center be identically aligned
or else that all body centers be randomly oriented.  In each case a
body-center translation leaves the structure invariant.  However,
recent experiments indicate order-disorder transitions in which a
supercell ordering develops at low temperatures, breaking the
body-centered cubic symmetry~\cite{Tamura,Tamura2}.

\section{Total energy calculations}

We carry out total energy calculations to determine if the
body-centered cell is energetically favorable.  Because of the low
transition temperature the energy difference will likely be small.
Additionally, any coupling between tetrahedron orientations at cell
vertex and body center must be at a distances of 10-14~\AA.
Interatomic pair potentials might not reliably address this subtle
problem (and they are not available presently in any case).

Instead we carry out {\em ab-initio} calculations using the
pseudopotential-based program VASP~\cite{VASP1,VASP2}.  This approach
uses ultrasoft pseudopotentials~\cite{USPP1,USPP2} or PAW
potentials~\cite{PAW1} to represent the effective interaction of
valence electrons with ionic cores, and solves the many-body quantum
mechanical band structure of these electrons using electronic density
functional theory.  This approach yields reasonably accurate total
energies.  We choose to model Cd-Ca rather than Cd-Yb because the
alkali earth element Ca is easier to treat from first principles than
the rare earth element Yb.

To assess the accuracy, we use VASP to reproduce the sequence
of low-temperature stable phases in the established Ca-Cd binary phase
diagram~\cite{BinaryPD}.  To do this we calculate the cohesive energy
for each known structure, and several hypothetical ones.  Each
structure is fully relaxed in both unit cell parameters and atomic
coordinates.  All energies are converged to an accuracy of 1 meV/atom
or better by increasing the k-point mesh density.  We then subtract
each cohesive energy from the tie-line joining the pure elements in
their ground states.

The result of this calculation is a plot of the enthalpy of formation
(at T=0K) for all known and many hypothetical Ca-Cd structures.
Agreement between our calculation and the established phase diagram
requires that all known low temperature structures lie on the convex
hull of enthalpy versus composition.  Additionally, all hypothetical
structures must lie above the convex hull, as must all known high
temperature, high pressure and metastable phases.

\begin{figure}
\includegraphics*[width = 3 in]{pdiag}
\caption{\label{fig:pdiag}This figure illustrates the enthalpy of
formation of Ca-Cd compounds.  Our notation is heavy circles indicate
known low temperature phases, light circles indicate known high
temperature phases, squares indicate either compound or structure
unknown in Ca-Cd system. Line sigments connect vertices of convex
hull.}
\end{figure}

This is shown in Fig.~\ref{fig:pdiag}.  We label each structure with
its name followed by its Pearson symbol.  Agreement is excellent, with
most low temperature Ca-Cd compounds lying on the convex hull, high
temperature structures slightly above and hypothetical structures at
higher energies yet.  There is one seeming difficulty with CaCd$_2$
which we now address.

The low and high temperature CaCd2 phases are reversed in energy
relative to the experimental report that oI12 is stable at high
temperatures and hP12 at low.  However, the transition between the
hP12 and oI12 variants has not been well established.  Most likely,
according to our findings, the presumed low temperature hP12 phase is
actually a metastable phase, and the nominal high temperature phase is
actually stable all the way to low temperatures.

The established phase diagram also lists the phase Ca$_3$Cd$_{17}$,
which is claimed to exist from T=0K up to melting, at a composition
extremely close to CaCd$_6$.  The structure of this phase is unknown.
Possibly this phase is actually the quasicrystal phase, whose energy
we do not calculate because we do not know its structure.  Note that
if the phase truly is the icosahedral quasicrystal, then either the
quasicrystal or the cubic CaCd$_6$.cI176 would almost surely be a high
temperature phase.  This is because thermodynamic rules governing
alloy phase diagrams dictate that the probability is low for two
distinct phases to coexist over an extended temperature range.

Another candidate for the Ca$_3$Cd$_{17}$ phase (based on atomic size
ratio and chemical similarity) is Sc$_3$Zn$_{17}$.cI160.  We find the
energy of this structure is substantially above the convex hull and
thus unfavorable.  Actually, the Sc$_3$Zn$_{17}$.cI160 structure is
identical to CaCd$_6$.cI176 except that the innermost tetrahedral
shell (shell (a) in Fig.~\ref{fig:struct}) is missing.  We thus
confirm that occupancy of the inner shell is energetically favorable.

\section{Orientational correlations}

The other apparent difficulty is the ``hypothetical'' structure marked
Cd$_6$Ca.cP176 appearing slightly below Cd$_6$Ca.cI176 in energy.
Let's describe these two structures.  They both contain 168 atoms in a
simple cubic unit cell, and they are identical except for the
orientations of the tetrahedron at cluster centers.  Note that, since
the tetrahedron takes two orientations, we can assign an Ising-like
spin variable $\pm$ to each one.  Assign a tetrahedron the $+$ sign if
one of its vertices falls along the direction $(1,1,1)$, and assign it
a $-$ sign if instead one of its vertices falls along the direction
$(-1,-1,-1)$.  In both structures the tetrahedron at the unit cell
vertex takes the $+$ orientation.  In the body-centered cI184
structure the cell center tetrahedron is also in the $+$ orientation,
while in the simple cubic cP176 structure this tetrahedron assumes the
$-$ orientation.

\begin{figure}
\includegraphics*[width = 1.5 in]{cubeA}
\includegraphics*[width = 1.5 in]{cubeB}
\caption{\label{fig:Ising}This figure illustrates the Ising-like configurational
degree of freedom.  Red=occupied by Cd, Green=vacant. (Left) ``+'' configuration,
(Right) ``-'' configuration.}
\end{figure}

To guage the validity of this result it is important to check
convergence in the density of the k-point mesh, the cutoff energy and
the sensitivity to choice of pseudopotential and the density
functional. Table~\ref{tab:cubedata} presents our study.  First we
vary the k-point mesh, from $1\times 1\times 1$ (the $\Gamma$ point)
up to $4\times 4\times 4$ (all Monkhorst-Pack meshes).  The table
presents the convergence of each structure energy separately as well
as the energy difference.  All other computational parameters were
held fixed during these calculations: medium precision (specifies
cutoff energy 168 eV); ultrasoft pseudopotential; Ceperly-Alder LDA;
no atomic relaxation.  In the next series we hold constant the k-point
density (we use $\Gamma$ for speed) and test the convergence in
precision going from low (cutoff 126 eV) to medium (cutoff 168 eV) to
high (cutoff 210 eV).  Next, continuing with medium precision and the
$\Gamma$ point, we compare unrelaxed energies with partial relaxation
(only relaxing clusters (a) and (b) as defined in
Fig.~\ref{fig:struct}) and full relaxation in which all atoms can
move.  The maximum displacement is 0.23~\AA~ for partial relaxation
and 0.30~\AA~ for full relaxation, always concentrated in shell
(b). Specifically, those Cd atoms in shell (b) that adjoin a vacancy
see the largest relaxations.

\begin{table}
\caption{\label{tab:cubedata}Energy convergence studies for CaCd6.cI176
and CaCd6.cP176, and their difference $\Delta=E_{cI}-E_{cP}$. All units
are eV per simple cubic cell.}
\begin{tabular}{l|r|r|r}
\hline
Setting             & E$_{cI}$ & E$_{cP}$ &  $\Delta$ \\
$1\times 1\times 1$ & -309.161 & -309.260 & 0.099 \\
$2\times 2\times 2$ & -309.814 & -309.975 & 0.161 \\
$3\times 3\times 3$ & -309.664 & -309.754 & 0.090 \\
$4\times 4\times 4$ & -309.653 & -309.780 & 0.127 \\
\hline
low               & -278.848 & -278.956 & 0.108 \\
medium            & -309.161 & -309.260 & 0.099 \\
high              & -309.662 & -309.761 & 0.099 \\
\hline
unrelaxed         & -309.161 & -309.260 & 0.099 \\
partial           & -311.902 & -311.995 & 0.093 \\
full              & -312.243 & -312.317 & 0.074 \\
\end{tabular}
\end{table}

Given the Ising-like $Z_2$ symmetry of the order parameter (one of two
orientations) it is appropriate to model the energy using an
Ising-model Hamiltonian. Including sufficiently far-neighbor
interactions we can surely capture the energetics accurately.
However, we have only a single energy difference to work with here, so
this is sufficient to extract only one coupling.  Assume this is the
nearest-meighbor coupling, along the cube body diagonal, and call it $J_1$.

Each BCC lattice site has 8 nearest neighbors, each of which reverses
sign when going from cI to cP.  There are two lattice sites per simple
cubic cell, but we must avoid over counting the bonds, since each bond
is shared by two lattice sites.  Hence we conclude that
$\Delta=16J_1$, or $J_1\approx +0.004$ eV (using the fully relaxed
value).  Because the value of $J_1$ is positive the interaction is
antiferromagnetic.

To determine the transition temperature for this system we wrote a
simple Monte Carlo program to simulate the BCC Ising antiferromagnet.
Actually, by simply reversing the sign convention for spins at
body-center sites, the BCC antiferromagnet can be seen to be
equivalent to the BCC ferromagnet.  According to our simulations, the
transition temperature should be around T=350 K.

Because this temperature is well above the reported transition
temperature of T=100K, we investigated the role of further neighbor
interactions.  If further-neighbor couplings have appropriate signs
the magnetism can become frustrated, lowering the transition
temperature and leading to spatial modulation of the low temperature
structure.

The next nearest neighbor lies along the cube edge, and we will call
this coupling $J_2$.  In order to extract values of $J_2$ we need to
study a larger cell, so we doubled the cell along the x-axis, and
considered the configurations denoted
${\sigma_1^v\sigma_1^c\sigma_2^v\sigma_2^c}$ in which the cube vertex
of the first cell has spin $\sigma_1^v$, etc.  Owing to the large
number of atoms present we report here only the results of $\Gamma$
point calculations.

\begin{table}
\caption{\label{tab:longdata}Energies of Ising decorations in double-length
structures. Units are eV per doubled cell.}
\begin{tabular}{l|r|r|r}
Config &  E$_0$   &  dE     & H \\
\hline
+~+~+~+   & -617.958 &  0      & $C+16J_1+12J_2$ \\
+~+~+~-   & -618.031 & -0.072  & $C+8J_2$ \\
+~+~-~-   & -617.986 & -0.028  & $C+4J_2$ \\
+~-~+~-   & -618.178 & -0.219  & $C+-16J_1+12J_2$ \\
\hline
\end{tabular}
\end{table}

Fitting these energies to a two-coupling Ising model yields values of
$J_1=+0.0068$ and $J_2=-0.0102$ eV.  Because $J_2$ is negative, the
next-nearest interaction proves ferromagnetic, which will increase the
transition temperature and also will not lead to superlattice
ordering.  In other words, the Ising model presented so far is rather
inconsistent with the experimental findings.

\section{Discussion}

Given that our results disagree with experiment in both superlattice
ordering and transition temperature it is clear that additional study
is needed.  We do not believe the basic methodology is flawed.  Rather
the most likely explanation is deficiency of the structural model
being used.  We restricted our attention to tetrahedron orientations
belonging to subsets of a basic cube, in order to have a simple
Ising-model description.  While the cubic sets of possible tetrahedron
oreintations seem appropriate given the overall cubic symmetry of the
cI176 structure, we must recognize that the second (b) shell is just a
slightly distorted dodecahedron, so there are many more plausible
tetrahedron orientations than the two we employ.

A recent study of the $M$Cd$_6$ structure family by Gomez and
Lidin~\cite{Gomez} finds the electron density distribution of the
central tetrahedra in CaCd$_6$ actually favor cube {\em edges} rather
than vertices.  They introduce additional partially occupied sites,
proposing the new Pearson type CaCd$_6$.cI232.  We have carried out
VASP calculations for these new structures and confirm they are
energetically advantageous over the original cI176/cP176
structures. If an order-disorder transition exists in this new
structure, it will be at a lower temperature than found in our initial
study because: (1) the energy differences among orientations are
generally lower; (2) the order parameter has a higher symmetry than
the Ising spins so there is a higher orientational entropy.

\begin{acknowledgments}
We wish to acknowledge useful discussions with R. Tamura and
C.L. Henley. This work was supported in part by NSF grant DMR-0111198
\end{acknowledgments}

\bibliography{CdCa}

\end{document}

