text:
# 3-Local Hamiltonian is QMA-complete

## 1 Introduction

Complexity theory is one of the cornerstones of theoretical computer science, formalizing the notion of an efficient algorithm (see, e.g., ). With the advent of quantum computing a plethora of new complexity classes have entered the field. One of the major challenges for theoretical computer science is to understand their structure and the interrelation between classical and quantum classes.

A seminal result in classical complexity theory is the celebrated Cook-Levin theorem which states that SAT is NP-complete. Namely, we are given a set of clauses (disjunctions) over a set of $`n`$ variables and asked whether there exists an assignment to the variables that satisfies all clauses. Moreover, the 3SAT problem in which each clause contains at most three literals is also NP-complete. It turns out that the 2SAT problem (where each clause contains at most two literals) can be solved in polynomial time (actually, there is even a linear time algorithm). However, the MAX2SAT problem in which we are given an extra number $`d`$ and asked whether there exists an assignment that satisfies at least $`d`$ clauses is still NP-complete.

In this paper we will be interested in the quantum analogues of the above results. For a good introduction the reader is referred to a recent survey by Aharonov and Naveh  and to a book by Kitaev, Shen and Vyalyi . Kitaev defined a quantum analogue of the classical class NP and named it BQNP. Strictly speaking, this class is the quantum analogue of MA, the probabilistic version of NP, and hence we will call it QMA (as was done in ).

QMA is naturally defined as a class of promise problems: A promise problem $`L`$ is a pair ($`L_{yes},L_{no}`$) of disjoint sets of strings corresponding to “Yes” and “No” instances of the problem. The problem is to determine, given a string $`xL_{yes}L_{no}`$, whether $`xL_{yes}`$ or $`xL_{no}`$. Let $``$ be the Hilbert space of a qubit.

###### Definition 1.1 (QMA)

Fix $`\epsilon =\epsilon (|x|)`$ such that $`2^{\mathrm{\Omega }(|x|)}\epsilon \frac{1}{3}`$. Then, a promise problem $`LQMA`$ if there exists a quantum polynomial time verifier $`V`$ and a polynomial $`p`$ such that:

* $`xL_{yes}|\xi ^{p(|x|)}Pr\left(V(|x,|\xi )=1\right)1\epsilon `$
* $`xL_{no}|\xi ^{p(|x|)}Pr\left(V(|x,|\xi )=1\right)\epsilon `$

where $`Pr\left(V(|x,|\xi )=1\right)`$ denotes the probability that $`V`$ outputs $`1`$ given $`|x`$ and $`|\xi `$.

By using amplification methods, it was shown in  that for any choice of $`\epsilon `$ in the above range the resulting classes are equivalent. In this paper we will assume that $`\epsilon `$ is $`2^{\mathrm{\Omega }(|x|)}`$.

We would also like to find an analogue of the SAT problem. One natural choice is the Local Hamiltonian problem. As we will see later, this problem is indeed a complete problem for QMA:

###### Definition 1.2

We say that an operator $`H:^n^n`$ on $`n`$ qubits is a $`k`$-local Hamiltonian if $`H`$ is expressible as $`H=_{j=1}^rH_j`$ where each term is a Hermitian operator acting on at most $`k`$ qubits.

###### Definition 1.3

The (promise) problem $`k`$-local Hamiltonian is defined as follows:

* Input: A $`k`$-local Hamiltonian on $`n`$-qubits $`H=_{j=1}^rH_j`$ with $`r=poly(n)`$. Each $`H_j`$ has a bounded operator norm $`H_jpoly(n)`$ and its entries are specified by $`poly(n)`$ bits. In addition, we are given two numbers $`a`$ and $`b`$ (with $`poly(n)`$ precision) such that $`ba>1/poly(n)`$. We are promised that the smallest eigenvalue of $`H`$ is either at most $`a`$ or larger than $`b`$.
* Output:

$`1`$ if $`H`$ has an eigenvalue not exceeding $`a`$,

$`0`$ if all eigenvalues of $`H`$ are larger than $`b`$.

We note that the original definition required that $`0H_j1`$ (i.e., that both $`H_j`$ and $`IH_j`$ are nonnegative, meaning that they only have nonnegative eigenvalues). However, it is easy to see that the two definitions are equivalent: given $`H_j`$’s such that $`H_jpoly(n)`$ for each $`j`$, normalize $`a`$, $`b`$ and all the $`H_j`$’s by a factor of $`1/poly(n)`$ such that $`H_j\frac{1}{2}`$. Then, add half the identity to each $`H_j`$ (such that $`0H_j1`$) and $`\frac{r}{2}`$ to $`a`$ and $`b`$ where $`r`$ is the number of terms in $`H`$.

It can be seen that the $`k`$-local Hamiltonian problem is NP-hard for all $`k2`$. This was recently shown by Wocjan and Beth  (see also ). One possible proof is to show that for any $`k2`$ the problem is at least as hard as MAX-$`k`$-SAT. The idea is to represent the $`n`$ variables by $`n`$ qubits and represent each clause by a Hamiltonian. Each Hamiltonian acts on the $`k`$ variables that appear in its clause. It ‘penalizes’ the assignment which violates the clause by increasing its eigenvalue. Therefore, the lowest eigenvalue of the sum of the Hamiltonians corresponds to the maximum number of clauses that can be satisfied simultaneously.

However, the only known QMA-completeness result was due to Kitaev which showed that the 5-local Hamiltonian problem is QMA-complete . An interesting open question which was already mentioned in  is whether the locality $`5`$ is optimal. Given that classically MAX$`2`$SAT is NP-complete we might hope to reduce the locality of the Hamiltonians. Our main theorem is the following:

###### Theorem 1.4

The problem 3-local Hamiltonian is QMA-complete.

We note that the 1-local Hamiltonian problem can be solved in polynomial time by a classical algorithm and is therefore unlikely to be QMA-complete. We leave the case of 2-local Hamiltonian as an open problem. Finally, we mention that using the methods of  one can show that the 2-local Hamiltonian problem is QMA-complete if we allow higher dimensional systems instead of qubits.

## 2 Kitaev’s Construction

In this section we will recall Kitaev’s proof that $`O(\mathrm{log}n)`$-local Hamiltonian is QMA-complete (his proof that $`5`$-local Hamiltonian is QMA-complete follows by a simple modification and we will mention it later). The proof begins by showing that $`k`$-Local Hamiltonian is indeed in QMA for any $`k=O(\mathrm{log}n)`$:

###### Lemma 2.1 ()

The $`k`$-local Hamiltonian problem is in QMA for any $`k=O(\mathrm{log}n)`$.

Then, it is enough to show that any problem $`L`$ in $`QMA`$ can be reduced to $`O(\mathrm{log}n)`$-local Hamiltonian. Let $`U_x=V(|x,)=U_T\mathrm{}U_1`$ be a quantum circuit of size $`T=poly(|x|)`$ operating on $`N=poly(|x|)`$ qubits. Notice that the input $`xL`$ is encoded into the circuit. We assume without loss of generality that $`TN`$ and that each gate $`U_i`$ operates on two qubits. Moreover, we assume that initially, the first $`m=p(|x|)`$ qubits contain the proof and the remaining ancillary $`Nm`$ qubits are zero (see Definition 1.1). Finally, we assume that the output of the circuit is written in the first computation qubit (i.e., it is 1 if the circuit accepts). The Hamiltonian $`H`$ that is constructed operates on a space of $`n=N+\mathrm{log}(T+1)`$ qubits. The first $`N`$ qubits represent the computation and the last $`\mathrm{log}(T+1)`$ qubits represent the possible values $`0,\mathrm{},T`$ for the clock. The Hamiltonian is constructed of three terms,

$$H=H_{in}+H_{out}+H_{prop}.$$
(1)
The terms are given by

$`H_{in}`$ $`=`$ $`{\displaystyle \underset{i=m+1}{\overset{N}{}}}|1_i1|_i|00|`$
$`H_{out}`$ $`=`$ $`|0_10|_1|TT|`$
$`H_{prop}`$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{T}{}}}H_{prop,t}`$ (2)
and

$$H_{prop,t}=\frac{1}{2}(I|tt|+I|t1t1|U_t|tt1|U_t^{}|t1t|)$$
(3)
for $`1tT`$ where $`|\alpha _i\alpha |_i`$ is the projection on the subspace in which the $`i`$’th qubit is $`|\alpha `$. It is understood that the first part of each tensor product acts on the space of the $`N`$ computation qubits and the second part acts on the clock qubits. $`U_t`$ and $`U_t^{}`$ in $`H_{prop,t}`$ act on the same computational qubits as $`U_t`$ does when it is employed in the verifier’s circuit $`U_x`$. Intuitively, each Hamiltonian ‘checks’ a certain property by increasing the eigenvalue if the property doesn’t hold: The Hamiltonian $`H_{in}`$ checks that the input of the circuit is correct (i.e., none of the last $`Nm`$ computation qubits is $`1`$), $`H_{out}`$ checks that the output bit indicates acceptance and $`H_{prop}`$ checks that the propagation is according to the circuit. Notice that these Hamiltonians are $`O(\mathrm{log}n)`$-local since there are $`\mathrm{log}(T+1)=O(\mathrm{log}n)`$ clock qubits. The proof is completed by the following lemmas and recalling that $`\epsilon `$ is chosen to be $`2^{\mathrm{\Omega }(|x|)}`$ so that $`\frac{c}{T^3}\frac{\epsilon }{T+1}>1/poly(n)`$:

###### Lemma 2.2 ()

Assume that the circuit $`U_x`$ accepts with probability more than $`1\epsilon `$ on some input $`|\xi ,0`$. Then the Hamiltonian $`H`$ has an eigenvalue smaller than $`\frac{\epsilon }{T+1}`$.

###### Lemma 2.3 ()

Assume that the circuit $`U_x`$ accepts with probability less than $`\epsilon `$ on all inputs $`|\xi ,0`$. Then all the eigenvalues of $`H`$ are larger than $`\frac{c}{T^3}`$ for some constant $`c`$.

Although the proof of this lemma will not be used in this paper, we sketch it here for completeness:

* We write $`H=H^{}+H_{prop}`$ where $`H^{}`$ denotes $`H_{in}+H_{out}`$. We start by noticing that both $`H^{}`$ and $`H_{prop}`$ are non-negative Hamiltonians. We can lower bound the smallest non-zero eigenvalue of $`H^{}`$ by $`1`$ since it is the sum of commuting projections. It can also be shown that the smallest non-zero eigenvalue of $`H_{prop}`$ is at least $`\mathrm{\Omega }(1/T^2)`$. This, however, is not enough to prove the lemma since, for example, the null-spaces of $`H^{}`$ and $`H_{prop}`$ might have a non-trivial intersection (i.e., there exists a non-zero vector in their intersection).

The next step is to show that since the circuit $`U_x`$ accepts with small probability, the angle between the null-spaces of $`H^{}`$ and $`H_{prop}`$ is not too small (in particular, this implies that the intersection of the two null-spaces is trivial). More specifically, we define the angle $`\theta `$ between the null-spaces of $`H^{}`$ and $`H_{prop}`$ by

$$\mathrm{cos}\theta =\mathrm{max}|\eta _1|\eta _2|$$

where the maximum is taken over all $`\eta _1`$ in the null-space of $`H^{}`$ and $`\eta _2`$ in the null-space of $`H_{prop}`$. Then, one can prove that $`\mathrm{sin}^2\theta \mathrm{\Omega }(1/T)`$. Finally, it can be shown that the smallest eigenvalue of $`H=H^{}+H_{prop}`$ can be lower bounded by the smallest eigenvalue among the non-zero eigenvalues of $`H^{}`$ and $`H_{prop}`$ times $`2\mathrm{sin}^2\frac{\theta }{2}`$. Hence, we get the lower bound

$$\mathrm{\Omega }(1/T^2)2\mathrm{sin}^2\frac{\theta }{2}$$

which is at least $`\frac{c}{T^3}`$ for some constant $`c>0`$.

## 3 The Construction

The result of the previous section can be improved to $`5`$-local Hamiltonian by using a unary representation for the clock and noting that three clock qubits are enough to identify the current time step (and since two computation qubits are also required, we get 5-local Hamiltonians). In addition, one has to add a Hamiltonian that penalizes clock qubits which are ‘illegal’, i.e., that do not represent a legal unary encoding. For more detail, see . In this section, we show how to use the result of the previous section to obtain the $`3`$-local Hamiltonian result. Our construction follows the ideas of Kitaev’s 5-local proof. The main difference is that our Hamiltonians use only one clock qubit instead of three. This requires another modification, namely, the penalty for illegal clock representations has to be considerably higher.

According to Lemma 2.1, $`3`$-local Hamiltonian is in QMA. Hence, it is enough to show that any problem $`L`$ in QMA can be reduced to the 3-local Hamiltonian problem. We are given a circuit $`U_x=U_T\mathrm{}U_1`$ as in the previous section. We construct a Hamiltonian $`H`$ that operates on a space of $`N+T`$ qubits. The first $`N`$ qubits represent the computation and the last $`T`$ qubits represent the clock. The Hamiltonian is constructed of four terms,

$$H=H_{in}+H_{out}+H_{prop}+H_{clock}.$$
(4)
The first three terms check that the input of the circuit is correct, that the output bit indicates acceptance and that the propagation is according to the circuit. As before, tensor products separate the computation qubits from the clock qubits:

$`H_{in}`$ $`=`$ $`{\displaystyle \underset{i=m+1}{\overset{N}{}}}|1_i1|_i|0_10|_1`$
$`H_{out}`$ $`=`$ $`|0_10|_1|1_T1|_T`$
$`H_{prop}`$ $`=`$ $`{\displaystyle \underset{t=1}{\overset{T}{}}}H_{prop,t}`$ (5)
$`H_{prop,t}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|10_{t,t+1}10|_{t,t+1}+I|10_{t1,t}10|_{t1,t}U_t|1_t0|_tU_t^{}|0_t1|_t)`$
for $`2tT1`$ and

$`H_{prop,1}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|10_{1,2}10|_{1,2}+I|0_10|_1U_1|1_10|_1U_1^{}|0_11|_1)`$
$`H_{prop,T}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|1_T1|_T+I|10_{T1,T}10|_{T1,T}U_T|1_T0|_TU_T^{}|0_T1|_T).`$
For any $`0tT`$, let $`|\widehat{t}`$ denote the state

$$|\underset{t}{\underset{}{1\mathrm{}1}}\underset{Tt}{\underset{}{0\mathrm{}0}}.$$

These are the legal unary representations. The last term is chosen to give a high penalty to states which do not contain a legal unary representation in the clock qubits:

$$H_{clock}=T^{12}\underset{1i<jT}{}|01_{ij}01|_{ij}$$
(6)
We denote the sum $`H_{in}+H_{prop}+H_{out}`$ of the computation related Hamiltonians by $`H_{comp}`$. Note that $`H`$ is a sum of 3-local Hamiltonians of bounded norm which can be specified by a polynomial number of bits, as required by Definition 1.3. We note that some of the terms in $`H_{prop}`$ are negative, but this is allowed by Definition 1.3.

###### Lemma 3.1 (Completeness)

Assume that the circuit $`U_x`$ accepts with probability more than $`1\epsilon `$ on some input $`|\xi ,0`$. Then $`H`$ has an eigenvalue smaller than $`\frac{\epsilon }{T+1}`$.

* Consider the vector

$$|\eta \stackrel{def}{=}\frac{1}{\sqrt{T+1}}\underset{t=0}{\overset{T}{}}U_t\mathrm{}U_1|\xi ,0|\widehat{t}.$$
(7)
Then,

$`\eta |H|\eta =\eta |H_{in}|\eta +\eta |H_{prop}|\eta +\eta |H_{clock}|\eta +\eta |H_{out}|\eta `$ (8)
and it is easy to see that the first three terms are zero. Moreover, since $`U_x`$ accepts with probability higher than $`1\epsilon `$,

$$\eta |H_{out}|\eta <\frac{\epsilon }{T+1}.$$
(9)

###### Lemma 3.2 (Soundness)

Assume that the circuit $`U_x`$ accepts with probability less than $`\epsilon `$ on all inputs $`|\xi ,0`$. Then all the eigenvalues of $`H`$ are larger than $`\frac{c}{T^3}`$ for some constant $`c`$.

* Let $`_{legal}`$ denote the subspace spanned by states whose clock qubits represent a unary encoding. The orthogonal space is denoted by $`_{illegal}`$. We will use a simple upper bound on the operator norm of $`H_{comp}`$ given by

$$H_{comp}H_{in}+H_{out}+\underset{t=0}{\overset{T}{}}H_{prop,t}N+1+2T4T.$$
(10)

We will show that for any unit vector $`|\eta `$, $`\eta |H|\eta \frac{c}{T^3}`$. Write $`|\eta =\alpha _1|\eta _1+\alpha _2|\eta _2`$ with $`|\eta _1_{legal}`$, $`|\eta _2_{illegal}`$, $`|\eta _1=|\eta _2=1`$ and $`\alpha _1,\alpha _2[0,1]`$ with $`\alpha _1^2+\alpha _2^2=1`$. If $`\alpha _2\frac{1}{T^5}`$ then

$$\eta |H|\eta \eta |H_{clock}|\eta H_{comp}\alpha _2^2T^{12}4T>1.$$
(11)
It remains to consider the case $`\alpha _2<\frac{1}{T^5}`$. Noting that $`H_{clock}0`$ we get:

$`\eta |H|\eta =\eta |H_{clock}|\eta +\eta |H_{comp}|\eta \eta |H_{comp}|\eta =`$
$`\alpha _1^2\eta _1|H_{comp}|\eta _1+2\alpha _1\alpha _2\text{Re}(\eta _1|H_{comp}|\eta _2)+\alpha _2^2\eta _2|H_{comp}|\eta _2=`$
$`\eta _1|H_{comp}|\eta _1\alpha _2^2\eta _1|H_{comp}|\eta _1+2\alpha _1\alpha _2\text{Re}(\eta _1|H_{comp}|\eta _2)+\alpha _2^2\eta _2|H_{comp}|\eta _2`$
$`\eta _1|H_{comp}|\eta _1{\displaystyle \frac{1}{T^{10}}}H_{comp}{\displaystyle \frac{2}{T^5}}H_{comp}{\displaystyle \frac{1}{T^{10}}}H_{comp}`$
$`\eta _1|H_{comp}|\eta _1{\displaystyle \frac{8}{T^9}}{\displaystyle \frac{8}{T^4}}>\eta _1|H_{comp}|\eta _1{\displaystyle \frac{9}{T^4}},`$ (12)
where we used the bound on the operator norm $`H_{comp}`$. Therefore, it is enough to show that for any $`\eta _{legal}`$, $`\eta |H_{comp}|\eta \frac{c}{T^3}`$. We will show that by using Lemma 2.3:

$`\eta |H_{comp}|\eta =\eta |\mathrm{\Pi }H_{comp}\mathrm{\Pi }|\eta =`$
$`\eta |\mathrm{\Pi }H_{in}\mathrm{\Pi }|\eta +\eta |\mathrm{\Pi }H_{out}\mathrm{\Pi }|\eta +{\displaystyle \underset{t=1}{\overset{T}{}}}\eta |\mathrm{\Pi }H_{prop,t}\mathrm{\Pi }|\eta `$ (13)
where $`\mathrm{\Pi }`$ is the projection on the subspace $`_{legal}`$. We compute $`\mathrm{\Pi }H_{comp}\mathrm{\Pi }`$:

$`\mathrm{\Pi }H_{in}\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \underset{i=m+1}{\overset{N}{}}}|1_i1|_i|\widehat{0}\widehat{0}|`$
$`\mathrm{\Pi }H_{out}\mathrm{\Pi }`$ $`=`$ $`|0_10|_1|\widehat{T}\widehat{T}|`$
$`\mathrm{\Pi }H_{prop,t}\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|\widehat{t}\widehat{t}|+I|\widehat{t1}\widehat{t1}|U_t|\widehat{t}\widehat{t1}|U_t^{}|\widehat{t1}\widehat{t}|)`$ (14)
for $`2tT1`$ and

$`\mathrm{\Pi }H_{prop,1}\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|\widehat{1}\widehat{1}|+I|\widehat{0}\widehat{0}|U_1|\widehat{1}\widehat{0}|U_1^{}|\widehat{0}\widehat{1}|)`$ (15)
$`\mathrm{\Pi }H_{prop,T}\mathrm{\Pi }`$ $`=`$ $`{\displaystyle \frac{1}{2}}(I|\widehat{T}\widehat{T}|+I|\widehat{T1}\widehat{T1}|U_T|\widehat{T}\widehat{T1}|U_T^{}|\widehat{T1}\widehat{T}|).`$
The Hamiltonian $`\mathrm{\Pi }H_{comp}\mathrm{\Pi }`$ acts on the Hilbert space $`_{legal}`$ whose dimension is $`2^N(T+1)`$. The Hamiltonian presented in Section 2 acts on a Hilbert space of the same dimension. In fact, notice that the two Hamiltonians are equivalent up to a renaming of the basis elements. Therefore, Lemma 2.3 implies that for any $`\eta _{legal}`$, $`\eta |H_{comp}|\eta \frac{c}{T^3}`$ which completes the proof.

## Acknowledgments

We wish to thank Dorit Aharonov and Fréderic Magniez for useful discussions. JK’s effort is sponsored by the Defense Advanced Research Projects Agency (DARPA) and Air Force Laboratory, Air Force Materiel Command, USAF, under agreement number F30602-01-2-0524. OR’s research is supported by NSF grant CCR-9987845.


format: text

text:
# Untitled Document

New Zero-Resistance States in Heterostructures: Anderson-Brinkman-Phillips Charge-Density Waves and New States Keshav N. Shrivastava School of Physics, University of Hyderabad,
Hyderabad 500046, India

Recently observed zero-resistance states are interpreted by Anderson-Brinkman and Phillips as charge density waves. We point out the existence of charge-density waves, superconducting states and new states including a state of zero charge and finite spin.

Corresponding author: keshav@mailaps.org
Fax: +91-40-2301 0145.Phone: 2301 0811.

1. Introduction

Recently zero-resistivity has been observed at small magnetic fields when the heterostructures are irradiated by microwaves. Phillips has suggested that the observation may be explained by sliding charge density waves which can give very low resistivity in the presence of a small magnetic field. Anderson and Brinkman have considered two possible mechanisms. According to one of the mechanisms, there are both positive as well as negative contributions to the resistivity so that the resultant value comes near zero. The negative resistance may imply charge-density fluctuations which may be observable in the power spectrum. In the past whenever a near zero resistance was found, it was assigned to charge-density waves and in a few cases, superconductivity was identified by using Meissner effect. When both the experiments, i.e., the zero-resistivity and negative susceptibility were successful, it was thought to be superconductivity with spin-singlet states. In recent years, a complex gap of a superconductor has been discussed. However, in the case of $`{}_{}{}^{3}He`$, singlets, triplets with $`\pm `$ 1 component, and triplets with 0 as well as $`\pm `$ 1 components have been found to phase separate.

Anderson-Brinkman as well as Phillips agree that there are some kind of charge density waves. We have been working on such problems for some time so that we pray that our interpretation may also be considered. We are not against the charge-density waves which is one of the components in our solution. Our solution follows.

2. Description

When we try to calculate the effective magnetic moment we find several solutions. One of the solutions has zero g-value so that the frequency of this mode is zero. This solution is therefore similar to the charge-density waves. When spin is reversed, the zero energy solution changes to a finite energy solution. If we can have one particle with spin up and zero energy and the other particle with spin down and finite energy, we can have a new type of superconductivity so that the resistivity is zero due to zero orbital angular momentum and zero spin. The orbital angular momentum of each of the particles is zero, $`l=0`$, and spin for one is $`down`$ while for the other is $`up`$ so that there is a zero momentum, spin-singlet type state which is superconducting. Thus for $`l=0`$, only these two solutions are found, namely, the charge density waves or the superconductivity. These states are distinguished by the sign of susceptibility. The charge-density waves are paramagnetic whereras superconductivity is diamagnetic.

Let us now imagine that l becomes finite. Here the contribution of Anderson model becomes very important. It may mean that going from l=0 to l=1 is an Anderson type transition. Once l=1, we have an interesting situation in which one spin configuration such as -1/2 gives an effective charge of 1/3 and +1/2 gives a charge of 2/3. In fact, varying the values of l reproduces “Stormer’s” values exactly. It is a wonder that this theory gives all of the fractional charges exactly as found in the experimental data on the quantum Hall effect. For very large values of the orbital angular momentum, even the 1/2 filled Landau level is produced and there is a Bose-Einstein condensation at this point. The width of the plateau is determined by the ordinary electron-phonon type interaction which can be reduced to an activation energy type process.

In this way, we have many solutions, one of which is a charge density, the other is superconductivity and there are others with effective fractional charges. This point of view is slightly different from Anderson-Brinkman-Phillips interpretation, yet it has all of the features of the experimental data. We have a particle with zero charge and spin 1/2 so that its frequency is also zero which is a characteristic of the charge-density wave, in agreement with Anderson, Brinkman and Phillips.

3. Conclusions.

Some new states occur. These states are superconducting but as the magnetic field is varied, the spin flips and the state is destroyed and then another state is created. One of the solutions has a charge density wave with zero frequency but the resistivity results are dominated by non-charge density waves. Therefore, the results of refs.1 and 2 have new results.

4. References

1. R. G. Mani, et al, Nature 419, 646(2002)
2. M. A. Zudov, et al, Phys. Rev. Lett. 90, 046807(2003).
3. P. W. Anderson and W. F. Brinkman, cond-mat/0302129.
4. J. C. Phillips, cond-mat/0212416.
5. K.N. Shrivastava, Introduction to quantum Hall effect,
Nova Science Pub. Inc., N. Y. (2002).
6. K. N. Shrivastava, cond-mat/0212552.
7. H. L. Stormer, Rev. Mod. Phys. 71, 875 (1999).

Note: Ref.5 is available from:
Nova Science Publishers, Inc.,
400 Oser Avenue, Suite 1600,
Hauppauge, N. Y.. 11788-3619,
Tel.(631)-231-7269, Fax: (631)-231-8175,
ISBN 1-59033-419-1 US$`\$69`$.
E-mail: novascience@Earthlink.net


text:
# Experimental implementation of an adiabatic quantum optimization algorithm

## Abstract

We report the realization of a nuclear magnetic resonance computer with three quantum bits that simulates an adiabatic quantum optimization algorithm. Adiabatic quantum algorithms offer new insight into how quantum resources can be used to solve hard problems. This experiment uses a particularly well suited three quantum bit molecule and was made possible by introducing a technique that encodes general instances of the given optimization problem into an easily applicable Hamiltonian. Our results indicate an optimal run time of the adiabatic algorithm that agrees well with the prediction of a simple decoherence model.

Since the discovery of Shor’sShor (1994) and Grover’s Grover (1997) algorithms, the quest of finding new quantum algorithms proved a formidable challenge. Recently however, a novel algorithm was proposed, using adiabatic evolution Farhi et al. ; Hogg (2000). Despite the uncertainty in its scaling behavior, this algorithm remains a remarkable discovery because it offers new insights into the potential usefulness of quantum resources for computational tasks.

Experimental realizations of quantum algorithms in the past demonstrated Grover’s search algorithm Jones et al. (1998); Vandersypen et al. (2000a), the Deutsch-Jozsa algorithm Chuang et al. (1998); Marx et al. (2000); Jones and Mosca (1998), order-finding Vandersypen et al. (2000b), and Shor’s algorithm Vandersypen et al. (2001). Recently, Hogg’s algorithm was implemented using only one computational step Peng et al. (2002), however a demonstration of an adiabatic quantum algorithm thus far has remained beyond reach.

Here, we provide the first experimental implementation of an adiabatic quantum optimization algorithm using three qubits and nuclear magnetic resonance (NMR) techniques Gershenfeld and Chuang (1997); Cory et al. (1997). NMR techniques are especially attractive because several tens of qubits may be accessible, which is precisely the range that could be crucial in determining the scaling behavior of adiabatic quantum algorithms Farhi and Goldstone (2001). Compared to earlier implementations of search problems Jones et al. (1998); Vandersypen et al. (2000a), this experiment is a full implementation of a true optimization problem, which does not require a black box function or ancilla bits.

This experiment was made possible by overcoming two experimental challenges. First, an adiabatic evolution requires a smoothly varying Hamiltonian over time, but the terms of the available Hamiltonian in our system cannot be smoothly varied and may even have fixed values. We developed a method to approximately smoothly vary a Hamiltonian despite the given restrictions by extending NMR average Hamiltonian techniques Rhim et al. (1970). Second, general instances of the optimization algorithm may require the application of Hamiltonians that are not easily accessible. We developed methods to implement general instances of a well known classical NP-complete optimization problem given a fixed natural system Hamiltonian.

We provide a concrete procedure detailing these methods. We then apply the results to our optimization problem which is known as Maximum Cut or maxcut Garey et al. (1976). Our experimental results indicate there exists an optimal total running time which can be predicted using a decoherence model that is based on independent stochastic relaxation of the spins.

The evolution of a quantum state during an adiabatic quantum algorithm is determined by a slowly varying, time-dependent Hamiltonian. Suppose we are given some time-dependent Hamiltonian $`H(t)`$ where $`0tT`$, and at $`t=0`$ we start in the ground state of $`H(0)`$. By varying $`H(t)`$ slowly, the quantum system remains in the ground state of $`H(t)`$ for all $`0tT`$ provided the lowest two energy eigenvalues of $`H(t)`$ are never degenerate Messiah (1976). Now suppose we can encode an optimization problem into $`H(T)`$. Then the state of the quantum system at time $`t=T`$ represents the solution to the optimization problem Farhi et al. . The total run-time $`T`$ of the adiabatic algorithm scales as $`g_{min}^2`$ where $`g_{min}`$ is the minimum separation between the lowest two energy eigenvalues of $`H(t)`$ Farhi et al. ; van Dam et al. . It is the scaling behavior of $`g_{min}`$ that will ultimately determine the success of adiabatic quantum algorithms. Classical simulations of this scaling behavior are hard due to the exponentially growing size of Hilbert space. In contrast, sufficiently large quantum computers could simulate this behavior efficiently.

Smoothly varying some time-dependent Hamiltonian appears straightforward but it contrasts with the traditional picture of discrete unitary operations including fault tolerant quantum circuit constructions Aharonov and Ben-Or (1997). Fortunately, we can approximate a smoothly varying Hamiltonian using methods of quantum simulations Trotter (1958) and recast adiabatic evolution in terms of unitary operations.

Discretizing a continuous Hamiltonian is a straightforward process and changes the run time $`T`$ of the adiabatic algorithm only polynomially van Dam et al. . For simplicity, let the discrete time Hamiltonian $`H[m]`$ be a linear interpolation from some beginning Hamiltonian $`H[0]=H_b`$ to some final problem Hamiltonian $`H[M]=H_p`$ such that $`H[m]=(m/M)H_p+(1m/M)H_b`$. The unitary evolution of the discrete algorithm can be written as:

$$U=\underset{m}{}U_m=\underset{m}{}\mathrm{e}^{\mathrm{i}((1m/M)H_b+(m/M)H_p)\mathrm{\Delta }t}$$
(1)
where $`\mathrm{\Delta }t=T/(M+1)`$, and $`M+1`$ is the total number of discretization steps. The adiabatic limit is achieved when both $`T,M\mathrm{}`$ and $`\mathrm{\Delta }t0`$.

Full control over the strength of $`H_b`$ and $`H_p`$ is needed to implement Eq.(1). However, this may not necessarily be a realistic experimental assumption. We will next show how the discrete time adiabatic algorithm can still be implemented when $`H_b`$ and $`H_p`$ cannot both be applied simultaneously and when they are both fixed in strength.

When both $`H_b`$ and $`H_p`$ are fixed, we can approximate $`U_m`$ to second order by using the Trotter formula $`\mathrm{exp}((A+B)\mathrm{\Delta }t)=\mathrm{exp}(A\mathrm{\Delta }t/2)\mathrm{exp}(B\mathrm{\Delta }t)\mathrm{exp}(A\mathrm{\Delta }t/2)+𝒪(\mathrm{\Delta }t^2)`$ Trotter (1958). Higher order approximations can be constructed if more accuracy is required.

Now suppose $`H_b`$ and $`H_p`$ are both constant. Since any unitary matrix is generated by an action $`\mathrm{i}H\mathrm{\Delta }t`$, we can increase the effect of a constant Hamiltonian $`H`$ by lengthening the time $`\mathrm{\Delta }t`$. Thus, we can implicitly increase the strength of $`H_b`$ and $`H_p`$ even when they are constant by simply increasing the time during which they are applied.

This technique also allows cases when the accessible Hamiltonians are not of the required strength, for example when we are given $`H_b^{}=gH_b`$ and $`H_p^{}=hH_p`$ but still wish to implement $`H_b`$ and $`H_p`$. Using all of the described techniques, we can now write $`U_m`$ as:

$$U_m\mathrm{e}^{\mathrm{i}H_b^{}[(1m/M)\mathrm{\Delta }t/2g]}\mathrm{e}^{\mathrm{i}H_p^{}[(m/M)\mathrm{\Delta }t/h]}$$
(2)
where $`AB=ABA`$. Each discretization step is of length $`(1m/M)\mathrm{\Delta }t/g+(m/M)\mathrm{\Delta }t/h`$, which is not constant when $`gh`$. As an illustration consider Fig. 1a.

In this experiment we choose $`\mathrm{\Delta }t=T/(M+1)`$ to be constant as we vary the number of discretization steps $`M+1`$. This way, the total run time $`T`$ increases with $`M+1`$, allowing us to test the behavior of the algorithm when approaching one of the conditions for the adiabatic limit. Even when the discrete approximation is not close to the adiabatic limit, the implemented algorithm can often find solutions using relatively few steps but lacks the guaranteed performance of the adiabatic theorem Hogg .

Adiabatic evolution has been proposed to solve general optimization problems, including NP-complete ones. In this general setting, the algorithm can depend on the existence of a black box function or the usage of large amounts of workspace. Our goal here is to optimize a hard natural problem in a way that avoids these difficulties. We will first describe which problem we chose and later on explain why it does not require ancilla qubits.

We found the maxcut problem to be a well-suited problem to demonstrate an adiabatic quantum algorithm because it allows a variety of interesting test cases. It also has applications in the study of spin glasses Barahona (1985) and VLSI design Cho et al. (1998), among others. The decision variant of the maxcut problem is part of the core NP-complete problems Garey et al. (1976), and even the approximation within a factor of 1.0624 of the perfect solution is NP-complete Ausiello et al. (1999).

The maxcut problem can be understood as follows. A cut is defined as the partitioning of an undirected $`n`$-node graph with edge weights into two sets. We define the payoff as the sum of weights of edges crossing the cut. The maximum cut is a cut that maximizes this payoff. By assigning either $`s_i=0`$ or $`s_i=1`$ to each node $`i`$, depending on its location with respect to the cut, the maxcut problem can be restated as finding the $`n`$ bit number $`s`$ that maximizes the payoff. An extension of the maxcut problem is to let the nodes themselves carry weights, which can be regarded as the nodes having a preference on their location. As an illustration consider a graph with three nodes as drawn in Fig. 1b.

The payoff as a function of the cut defined by $`s`$ is given by

$$P(s)=\underset{i}{}w_is_i+\underset{i,j}{}s_i(1s_j)w_{ij}$$
(3)
where $`w_{ij}`$ are the edge weights, $`w_i`$ denotes the preference of the nodes to be on the $`1`$ side of the cut, and $`s_i`$ is the value of the $`i`$-th bit of $`s`$, for $`0s2^n1`$.

The smallest meaningful test case of the maxcut problem requires $`3`$ nodes and admits a variety of interesting cases by varying $`w_i`$ and $`w_{ij}`$. We aimed at two goals when choosing a representative set of weights. First, we wanted the minimum energy gap $`g_{min}`$ to be smaller than the one for a $`3`$-qubit adiabatic Grover search. Second, we wanted a resulting energy landscape with both a global and local maximum such that a greedy classical search would incorrectly find the local maximum half the time gre . These goals are met by the choice $`w_1=w_2=w_3=2`$, $`w_{12}=2`$, $`w_{13}=1`$, $`w_{23}=3`$. The payoff function for this set of weights is $`P(s)=[\mathrm{0\; 6\; 7\; 7\; 5}\mathbf{9}\mathrm{8\; 6}]`$ where $`s`$ = \[000 001 010 011 100 101 110 111\]. The global maximum lies at $`s=101`$ so the answer on the quantum computer following measurement should be $`|101`$, and not at the local maximum $`s=110`$.

In the quantum setting, this payoff function $`P(s)`$ can be encoded into the Hamiltonian $`H_p`$ by rewriting Eq.(3) using Pauli matrices:

$$H_p=\underset{i}{}w_i(I\sigma _{zi})/2+\underset{i<j}{}w_{ij}(I\sigma _{zi}\sigma _{zj})/2$$
(4)
where $`I`$ is the $`2^n`$x$`2^n`$ identity matrix and $`\sigma _{zi}`$ is the Pauli $`Z`$ matrix on spin $`i`$. The identity matrices in the equation above only lead to an overall phase which cannot be observed, and hence they can be ignored. The diagonal values of Eq.(4) are equal to $`P(s)`$. Because of the direct encoding of $`P(s)`$ into $`H_p`$ no black box function or ancilla qubits are required, which makes this a full implementation of an optimization problem.

Similar to Eq.(4), the natural Hamiltonian of $`n`$ weakly coupled spin-1/2 nuclei subject to a static magnetic field $`B_0`$ is well approximated by Ernst (1994)

$$=\underset{i}{}\omega _i\sigma _{zi}/2+\underset{i<j}{}\pi J_{ij}\sigma _{zi}\sigma _{zj}/2+_{env}$$
(5)
where the first term represents the Larmor precession of each spin $`i`$ about $`B_0`$, and $`\omega _i`$ is its Larmor frequency. The second term describes the scalar spin-spin coupling of strength $`J_{ij}`$ between spins $`i`$ and $`j`$. The last term represents coupling to the environment, causing decoherence. Note the resemblances between $``$ and $`H_p`$.

Despite the similarities, the spin-spin couplings of Eq.(5) are generally different from a randomly chosen set of weights. Therefore, we require a procedure to turn the fixed $`J_{ij}`$ into any specified weights $`w_{ij}`$. This is achieved using refocusing schemes that are typically used to turn on only one of the couplings while turning all others off Ernst (1994).

We have modified a refocusing scheme to effectively change the couplings to any arbitrary value. Consider the pulse sequence drawn in Fig. 2. Based on this scheme, we can derive the under-constrained system $`(\alpha +\beta \gamma \delta )J_{12}=w_{12}`$, $`(\alpha \beta \gamma +\delta )J_{13}=w_{13}`$, and $`(\alpha \beta +\gamma \delta )J_{23}=w_{23}`$, which can be solved for positive $`\alpha `$, $`\beta `$, $`\gamma `$, and $`\delta `$ such that $`J_{ij}w_{ij}`$.

The single weights $`w_i`$ are implemented by introducing a reference frame for each spin $`i`$ which rotates about $`B_0`$ at frequency $`(\omega _iw_i)/2`$. In order to apply the single qubit rotations of our refocusing scheme on resonance, we apply the reference frequency shift only during the delay segment $`\alpha `$, which we can always choose to be a positive value. Thus, $`H_p`$ is implemented by applying the refocusing scheme from Fig. 2 while going off-resonance during the delay segment $`\alpha `$.

A full implementation of an adiabatic algorithm also requires a proper choice of $`H_b`$. We choose $`H_b=_i\sigma _{xi}`$ for several reasons. First, its highest two excited states are non-degenerate. Second, it can be easily generated using single qubit rotations, and third, its highest excited state is created from a pure state with all qubits in the $`|0`$ state by applying a Hadamard gate on all qubits (we require the initial state to be the highest excited state of $`H_b`$ because we are optimizing for the maximum value of $`H_p`$).

The full adiabatic quantum algorithm is now implemented by first creating the highest excited state of $`H_b`$. We then apply $`M+1`$ unitary matrices as given by Eq.(2) and illustrated by Fig.1a. Accordingly, from slice to slice, we decrease the time during which $`H_b`$ is active while increasing the time during which $`H_p`$ is active. Finally, we measure the quantum system and read-out the answer.

We selected 13C-labeled CHFBr2 for our experiments. The Hamiltonian of the 1H-19F-13C system is of the form of Eq.(5) with measured couplings $`J_{HC}=224`$ Hz, $`J_{HF}=50`$ Hz, and $`J_{FC}=311`$ Hz. The interaction with the Br nuclei is averaged out, contributing only to $`_{env}`$. Experiments were carried out at MIT using an 11.7 Tesla Oxford Instruments magnet and a Varian Unity Inova spectrometer with a triple resonance (H-F-X) probe from Nalorac.

The experiments were performed at room temperature at which the thermal equilibrium state is highly mixed and cannot be turned into the required initial state by just unitary transforms. We thus first created an approximate effective pure state as in ref. Vandersypen et al., 2000a by summing over three temporal labeling experiments.

In our experiments, we actually implemented $`0.5H_p`$ and $`0.5887H_b`$ instead of $`H_p`$ and $`H_b`$. This ensures that the error due to the 2nd order Trotter approximation is bounded by $`\sqrt{_i|ϵ_i|^2}<0.0356`$ where $`ϵ_i`$ is the contribution of the $`i`$-th undesired Pauli matrix. We also choose $`g`$ so the applied RF field does not heat the sample, and $`gh`$ so $`J_{ij}`$ can be ignored when applying $`H_b`$. All of these choices result in a total experimental time that is within the shortest $`T_2`$ decoherence time Vandersypen et al. (2000a). We reconstructed the traceless deviation density matrices upon completion of the experiments using quantum state tomography Vandersypen et al. (2000a).

We executed this algorithm for several $`M`$ (with $`w_i`$ and $`w_{ij}`$ as listed above Eq.(4)). Since we chose $`\mathrm{\Delta }t`$ to be constant, this meant increasing the run-time $`T`$ of the algorithm. The reconstructed deviation density matrices are shown in Fig. 3. The plots clearly display the expected pure state $`|101`$. The local maximum at $`s=110`$ has a decreasingly small probability of being measured for increasing $`M`$. Simulations using Eq.(2) show that this optimization algorithm performs better for increasing $`M`$. We wanted to verify whether this is indeed true experimentally.

For this purpose, we estimate the error of our obtained deviation density matrices compared with the ideal case of $`M=\mathrm{}`$. Fig. 4a plots the trace distance as a function of $`M`$, using the same arbitrary scale as in Fig. 3. From the plot, we observe there exists an optimal run-time of the algorithm, corresponding to $`0.226`$ seconds in our experiment. This optimal run time is in good agreement with the prediction of a previously developed simple decoherence model Vandersypen et al. (2001). Predicting the impact of decoherence has already provided invaluable insight into estimating errors in previous experiments Vandersypen et al. (2001), and we believe continued effort towards understanding decoherence will greatly benefit experimental investigations of quantum systems.

In conclusion, we have provided the first experimental demonstration of an adiabatic quantum optimization algorithm. We show a concrete procedure turning a continuous time adiabatic quantum algorithm into a discrete time version, even when certain restrictions apply to the accessible Hamiltonians. Our results indicate that there exists an optimal run-time of the algorithm which can be roughly predicted using a simple decoherence model. We believe this implementation opens the door to a variety of interesting experimental demonstrations and investigations of adiabatic quantum algorithms.

###### Acknowledgements.

We wish to thank A. Childs, A. Landahl, and E. Farhi for useful discussions. This work was supported by the NSF grant CCR-0122419, the DARPA QuIST program, the HP/MSRI postdoc fellowship.


text:
# Bonding in MgSi and AlMgSi Compounds Relevant to AlMgSi Alloys

## I Introduction

Precipitation or age hardened alloys are today one of the most important alloy types in industry. In the AlMgSi alloy system, Mg-Si and Al-Mg-Si precipitates formed during specific heat treatments give rise to a very significant increase in strength. The precipitation sequence is generally accepted to be:

$$SSSS\text{Mg/Si Clusters}GPZ\beta ^{\prime \prime }\beta ^{}\beta ,$$
(1)
where GPZ refers to a Guinier Preston Zone and SSSS refers to a Super Saturated Solid Solution. Very little is known about the early stages of the precipitation process. However, it is believed that when the SSSS is heat treated Mg and Si atoms quickly diffuse substitutionally to form small clusers due to the large amount of quenched-in vacancies Porter and Easterling (1997) (a large part of the vacancies move to interfaces like surfaces and grain boundaries in the later stages of the heat treatment). Although the details are hard to investigate experimentally, several studies of such clustering have been carried out using atom probe microscopy Murayama and Hono (1999); Edwards et al. (1998). The first phase which can be resolved using high resolution electron microscopy (HREM) is the GPZ. From this, a model of its crystal structure has recently been proposed Marioara et al. (2001). The structure of the proceeding phase, $`\beta ^{\prime \prime }`$, was also solved using electron microscopy techniquesZandbergen et al. (1997), a result which has later been supported by *ab initio* calculationsDerlet et al. (2002).

Earlier it was believed that the GPZ, $`\beta ^{\prime \prime }`$ and $`\beta `$ phases all had the stoichiometry of the $`\beta `$ phase, $`\mathrm{Mg}_2\mathrm{Si}`$, and that alloys should be optimized accordingly. Thus, the terminal $`\beta `$ phase was of primary importance. Later however, it had been confirmed that the $`\beta ^{\prime \prime }`$ gives a greater contribution to the hardness due to its semi-coherent interface with the aluminium matrix and needle-like shape, which is more effective for dislocation pinning Takeda et al. (1998). The $`\beta ^{\prime \prime }`$ stoichiometry has been shown to be $`\mathrm{Mg}_5\mathrm{Si}_6`$Zandbergen et al. (1997).

Recently, several additional phases have been identified experimentallyMatsuda et al. (2000); Andersen et al. (2002a, b), giving the extended precipitation sequence:

$`SSSSMg/Si\text{ Clusters}GPZ\beta ^{\prime \prime }`$ (2)
$`(\beta ^{}+U2+U1+B^{})\beta .`$
In ref. Matsuda et al., 2000 the phases $`U1`$, $`U2`$ and $`B^{}`$ are referred to as type *A*, *B* and *C* respectively. These three phases, in addition to $`\beta ^{}`$, are often grouped together since little is known about their interdependence. However it is believed that the peak of concentration with respect to time for each of these structures follows the ordering given by the above precipitation sequence Matsuda et al. (2000). They all form relatively late in the precipitation sequence usually at temperatures in the range 200-300 $`{}_{}{}^{}\mathrm{C}`$, and in Si-rich alloy compositions.

It has been considered a general rule of thumb that successful aluminium precipitation hardening alloys contain secondary and ternary alloying elements which are larger and smaller than aluminium Ryum (1986). Now, the concept of atom size in this context must be based on the type of bonding involved. One can have either ionic, metallic or covalent radii for the constituent elements giving dramatically different values for atomic size Kittel (1996). It is clear that when studying the electronic density of compound structures this type of concept may lead to an oversimplification. It it therefore interesting to carry out a theoretical study of the electronic structure and bonding characteristics of the relevant phases with respect to their relative stability. This is the purpose of the present work, where we employ full potential *ab initio* methods based on Density Functional Theory (DFT) to investigate the bonding within the above mentioned structures. In section II we describe the Linear Augmented Plane Wave + local orbitals approach, (L)APW+(lo), used for all calculations. Section III describes the results obtained for the three models of precipitate phases containing only Mg and Si ($`\beta ^{\prime \prime }`$, $`\beta ^{}`$, and $`\beta `$), and section IV deals with the phases containing Al, Mg and Si ($`U1`$ and $`U2`$).

## II Method

The *ab initio* calculations were performed using WIEN2k, a program package implementing the full potential (L)APW+(lo) DFT methodBlaha et al. (1999). The Augmented Plane Wave approach (APW)+(lo) methodSjostedt et al. (2000) differs from the LAPW method in the linearization of the APW’s. In slater’s original APW method Slater (1964) the unit cell is partitioned into non-overlapping atomic spheres and an interstitial region. The basis functions are plane waves for the interstitial region and radial wave functions within the atomic spheres. In the LAPW method the basis functions inside the spheres are linearized with respect to the energy $`E_l`$:

$$\psi _{𝒌_𝒏}(𝒓)=\underset{lm}{}(A_{lm,𝒌_𝒏}u_l(r,E_l)+B_{lm,𝒌_𝒏}\dot{u}_l(r,E_l))Y_{lm}(\widehat{𝒓})$$
(3)
where $`u_l(r,E_l)`$ are radial functions and $`Y_{lm}(𝒓)`$ are spherical harmonics. $`\dot{u}_l(r,E_l)`$ is the energy derivative of $`u_l(r,E_l)`$. $`A_{lm,𝒌_𝒏}`$ and $`B_{lm,𝒌_𝒏}`$ are determined by matching the above basis to the value and derivative of the plane waves for each k-vector at the sphere boundary. An alternative way for doing the linearization is the APW+lo method. One starts with the original APW’s and adds local orbitals to obtain the variational flexibility in the basis functions.

$`\psi _{𝒌_𝒏}(𝒓)={\displaystyle \underset{lm}{}}(A_{lm,𝒌_𝒏}u_l(r,E_l)+\varphi _{lo})Y_{lm}(\widehat{𝒓})`$ (4)
$`\varphi _{lo}=B_{lm}u_l(r,E_l)+C_{lm}\dot{u}_l(r,E_l)`$ (5)
At first sight this looks very similar to the LAPW basis. However, the $`B_{lm}`$ and $`C_{lm}`$ are no longer dependent on the wave-vector and are determined by the requirement that the local orbital is zero at the sphere boundary and normalized. The great advantage of this scheme is that the calculations converge to results almost identical with those of the LAPW method but for dimensioning parameters which effectively leads to a smaller basis set Madsen et al. (2001). The WIEN2k code uses a mixed APW+lo/LAPW basis set, exploiting the advantages of both methods.

To make the results for the different structures comparable, we used the same set of APW+(lo) parameters for all calculations: $`R_{mt}=2.1`$ Bohr, $`R_{mt}K_{max}=7`$ and $`G_{max}=14`$ $`\mathrm{Ry}^{1/2}`$. Here, $`R_{mt}`$ is the muffin tin radius, $`K_{max}`$ is the plane wave cut-off and $`G_{max}`$ is the maximum Fourier component of the electron density. For all calculations we used the modified tetrahedron methodBlochl et al. (1994) for Brillouin zone integrations. All k-point meshes were checked for convergence. Thus in general, the highly symmetric structures (like bulk Al, Si and Mg) with few symmetrically inequivalent atoms in the unit cell, require a denser k-point mesh than the precipitate phases with larger unit cells and a lower symmetry. For the exchange-correlation potential we used the Generalized Gradient Approximation (GGA) of Perdew et. al Perdew et al. (1996)

## III $`\mathrm{MgSi}`$ Phases

### III.1 $`\beta `$

The $`\beta `$ phase (fluorite $`\mathrm{Mg}_2\mathrm{Si}`$) is the terminal equilibrium structure of the precipitation sequence. It has a fcc primitive unit cell (space group $`\mathrm{Fm}\overline{3}\mathrm{m}`$ (225)), with an experimental lattice parameter a=6.39 Å Marioara et al. (2001). It forms precipitates of a plate-like or cubic shape up to 20 $`\mu \mathrm{m}`$ in diameter. Its interface with the Al matrix is fully incoherent Marioara et al. (2001). The conventional unit cell, containing 8 Mg atoms and 4 Si atoms, is shown in fig. 1. As can be seen, each of the Si atoms in the structure has 8 Mg nearest neighbours, giving each Mg atom 4 Si nearest neighbours at the same distance. The Si atoms are arranged as an fcc lattice interpenetrated by a sc Mg lattice.

Several *ab initio* studies of the bonding in Fluorite $`\mathrm{Mg}_2\mathrm{Si}`$ have been carried out in the past Baranek et al. (1997, 1999); Wood and Zunger (1986). However, to our knowledge, none of these involved the LAPW or (L)APW+(lo) DFT method. For our calculations we use the experimentally observed lattice constant for this phase in the Al matrix: 6.39 Å. Performing a volume relaxation, this value differed by only 0.25 % from the calculated optimized value 6.37 Å, and by 0.6 % from the experimental value for bulk $`\mathrm{Mg}_2\mathrm{Si}`$, 6.338 ÅAnastassakis and Hawranek (1972). The calculated bulk modulus, derived from a second order Birch fit, was 54.3 GPa, compared to 59 GPa from experiment Anastassakis and Hawranek (1972). Using the 6.39 Å lattice constant the nearest neighbour (nn) Mg-Si distance is 2.77 Å, the Si-Si nn distance is 4.51 Å and the Mg-Mg nn distance is 3.29 Å.

Past linear-combination-of-atomic-orbitals Hartree Fock and DFT pseudopotential calculations Baranek et al. (1997); Wood and Zunger (1986) indicate considerable charge transfer from the electropostive Mg atoms to the electronegative Si atoms resulting in a partly ionic Si-Mg bond for the $`\beta `$ phase. This result is also supported by the present work. Fig. 2 displays the bonding charge density for the $`(\overline{1}01)`$ plane of fig. 1. We define the bonding charge density as the difference between the converged valence charge density from DFT and the charge density derived from the isolated neutral atoms. Thus the bonding charge density indicates the charge transfer resulting from the converged electronic structure. In fig. 2, the thick lines represent positive contours corresponding to charge transfer to the region, while the thinner lines represent negative contours corresponding to charge transfer away from the region. The concentration of thick lines around the Si atoms thus displays a net build up of charge around the Si atoms counterbalanced by a net reduction of charge around the central Mg atoms, indicating that some ionicity is indeed at play.

Fig. 3 displays the corresponding total and partial electronic density of states (DOS) for the $`\beta `$ phase. Here the partial DOS represents the “Muffin tin decomposed” DOS for both the Mg and Si atoms, in which the occupied states are projected onto the muffin tin eigenvectors (see section II) of the particular atom. In this way, some information can be gained on the local DOS at each atom, however such a procedure is not complete since the interstitial region cannot be locally resolved. On the other hand, the total DOS represents the DOS derived from all the“Muffin tin decomposed” DOS and the interstitial DOS of the entire computational cell. Note that the calculation is based on the primitive unit cell, containing 1 symmetrically inequivalent Si atom and 2 symmetrically inequivalent Mg atoms, thus the partial Mg DOS contains the contribution for 2 Mg atoms, and the Si DOS the contribution from 1 Si atom.

A general feature of the total and partial DOS in fig. 3 is the two broad bands below the Fermi energy (indicated by the horizontal line at 0.26 Ry) separated by an energy of approximately 0.4 Rydbergs. That the occupancy is dominated by Si states, where the partial Mg DOS represents the two inequivalent Mg atoms, is a further indication of the presence of ionicity whereby the Mg has simply donated electrons to the Si. In addition this phase has a band-gap. The magnitude of the gap at the $`\mathrm{\Gamma }`$-point is 0.13 Ry (1.77 eV), which is somewhat smaller than the experimental value of 2.27 eV for a lattice constant of 6.338 Å Madelung (1983).

Figs. 4a and b now display the $`s`$, $`p`$ and $`d`$ decomposed partial DOSs for both Si and the two inequivalent Mg atoms. In this case, the total DOS represents the total muffin tin decomposed DOS for the particular atom, that is, the sum of the $`s`$, $`p`$ and $`d`$ partial DOSs. Inspection of fig. 4a reveals that the lower band consists of $`s`$ states where as the upper band consists of predominantly $`p`$ states with a little $`s`$ state character. For the case of the Mg, fig. 4b, the DOS below the Fermi level largely mimics that of Si apart from both occupied bands being of mixed $`s`$ and $`p`$ character.

Insight into the general features of the DOSs shown in fig. 4 can be understood from the perspective of beginning with pure fcc Si with a lattice constant of 6.39Å, giving a rather large nn Si-Si separation of 4.52Å. For such a system there can be no strong hybridization between the atomic valence $`s`$ and $`p`$ states, maintaining a gap of approximately 0.3 Ry in the corresponding DOS (not shown). The $`s`$ and $`p`$ band centers differ by about 0.4 Ry, which is not so dissimilar from the 0.49 Ry difference between the isolated atom 3s and 3p valence states. With the addition of the interpenetrating cubic array of Mg atoms, both bands broaden reducing the gap to about 0.2 Ryd. Thus in fig. 4a, the dominant s character of the filled lower band and the dominant $`p`$ character of the upper filled band arises from the strong onsite Si orthogonality requirement, whereas in fig. 4 the mixed $`s`$ and $`p`$ character of the corresponding Mg bands arises from the Mg-Si $`s`$ and $`p`$ matrix interaction elements. We note that the corresponding heights of the Mg total DOS are significantly less than the Si DOS, following approximately the square-root dependence of a metallic DOS. Thus, through charge transfer to the Si, the Mg plays the role of “strengthening” the Si backbone lattice, providing an explanation of the origin of the $`\beta `$ phase band gap through the pulling down of the excited $`s`$, $`p`$ and $`d`$ states of atomic Si.

### III.2 $`\beta `$

The $`\beta ^{\prime \prime }`$-phase has a base-centred-monoclinic conventional unit cell (space group $`\mathrm{C2}_\mathrm{m}`$) with experimental lattice parameters $`a=15.16`$ Å, $`b=6.74`$ Å, $`c=4.05`$ Å, and $`\gamma _{ab}`$=$`105.3^{}`$ Zandbergen et al. (1997). Fig. 5 displays the conventional cell for two viewing orientations. This phase forms precipitates of a needle-like shape, typically with a thickness of 30-50 Å and a length of 300-400 Å. The needle length, which is along the $`<001>`$ direction of the $`\beta ^{\prime \prime }`$-phase unit cell, runs parallel to the $`<001>`$ direction of the Aluminium matrix. For the present work we use the experimentally derived lattice parameters and the relaxed atomic positions of the inequivalent Si and Mg atoms by minimizing the forces by a Newton-Coates procedure with respect to the $`\mathrm{C2}/\mathrm{m}`$ space group. The resulting inequivalent atomic coordinates are listed in table 1. The given space group of this phase was first solved by high resolution electron microscopy techniques Zandbergen et al. (1997) and later supported by FLAPW DFT calculations Derlet et al. (2002).

Fig. 6 displays the bonding charge density contour plot for the (002) plane. A dominant feature is the concentration of charge between the Si nearest neighbours, and to a lesser extent between the Si and Mg nearest neighbours, indicating that covalency is at play in this system Derlet et al. (2002). Such charge transfer to the bonding regions originates from the core regions of both the Si and Mg atoms, in addition to the homogeneous interstitial region between the Mg atoms. The depletion of charge from the Mg is comparable to that seen in the $`\beta `$-phase, indicating that for this system both ionicity and covalency is present in the bonding. The nearest neighbour distance range in this structure is 2.39-2.53 Å for Si-Si and 2.61-2.84 Å for Si-Mg. This is not dissimilar to the nearest neighbour distance of 2.33 Å in covalently bonded diamond cubic Si and 2.76 Å in the covalent/ionic equilibrium $`\beta `$ phase. Moreover, one of the inequivalent Si atoms is sp3 (tetrahedrally) coordinated and another is sp2 coordinated Derlet et al. (2002).

Fig. 7 displays the Total DOS (including the muffin tin and interstitial regions) for the $`\beta `$ phase. Unlike the $`\beta `$ phase, there is no band gap. However, depressions can be identified at -0.1 and 0.4 Rydbergs reminiscent of the gaps seen in fig. 3. By inspection of the partial DOS for each inequivalent atom we find that for the Si atoms there is a dominant $`s`$ character for the lower energy states ($`<0.2`$ Ry), whereas, at about 0.2 Ry, the $`p`$ character dominates. In between these regimes, a mixture of $`p`$ and $`s`$ character exists, indicating strong hybridization. Indeed fig. 8a displays the partial DOS for Si2, which from past work is locally sp3 coordinated Derlet et al. (2002), showing approximate correspondence with the diamond cubic Si DOS Liu et al. (1988). In particular the dominant $`s`$ character of the L$`_2^{}`$, the $`s`$ and $`p`$ character of the L1 and the dominant $`p`$ character of X4 the high symmetry points. Alternatively, fig. 8b shows the partial DOS for Mg3, which from charge density profile analysis Derlet et al. (2002) is in a homogeneous metallic environment, displaying an approximate shifted metallic square-root behaviour.

## IV $`\mathrm{AlMgSi}`$ Phases

### IV.1 $`U1`$

The $`U1`$ conventional/primitive unit cell is trigonal (space group $`P_{\overline{3}m1}`$) with experimentally derived lattice parameters $`a=b=4.05`$ Å, and $`c=6.74`$ Å. It contains 1 Mg, 2 Al, and 2 Si atoms giving the formula $`\mathrm{MgAl}_2\mathrm{Si}_2`$. The unit cell and atomic positions are shown in table 2 and fig. 9. Experiments show that the $`<001>`$ direction of the $`U1`$ unit cell runs parallel to the $`<310>`$ direction of the fcc Al matrix Andersen et al. (2002a). The $`U1`$ phase usually forms rod-like precipitates with a length of 50-500 nm and a width of approximately 50 nm. Performing a volume relaxation, we found the calculated unit cell equilibrium volume to differ by only 1% from the experimentally derived unit cell volume. From a second order Birch fit we obtained a bulk modulus of 71 GPa.

The $`U1`$ phase can be categorized as belonging to a class of structures given the name $`\mathrm{CaAl}_2\mathrm{Si}_2`$-type Zintl compounds Cheng et al. (1986). In Zintl compounds, the electropositive elements are thought of as merely electron donors to the electronegative elements which thereby are able to fulfill the octet rule. For $`\mathrm{MgAl}_2\mathrm{Si}_2`$ this implies that each Mg atom donates two electrons to $`\mathrm{Al}_2\mathrm{Si}_2`$ for each unit cell. This gives $`\mathrm{Al}_2\mathrm{Si}_2^2`$, or two units of $`\mathrm{AlSi}^{}`$ which is isoelectronic to AlN, fulfilling the octet rule. The resulting layered structure is sketched in fig. 10. $`\mathrm{AlSi}^{}`$ constitutes a double layer of tightly bound “chair-like” six-membered rings separated by layers of hexagonal $`\mathrm{Mg}^{2+}`$. In the Al-Si network, each Si atom bonds to 4 Al atoms forming an umbrella-like structure, while each Al atom forms the more common tetrahedral structure with 4 nearest neighbour Si atoms The Al-Si bond length within each AlSi layer is 2.48 Å, while the length of the Al-Si bonds connecting the two $`\mathrm{AlSi}^{}`$ layers is 2.62 Å. The length of the Mg-Si bond is 2.86 Å (see below and fig. 11), somewhat larger than the length of the partly ionic bond in the $`\beta `$ equilibrium structure, but similar to Mg-Si bonds found in the $`\mathrm{Mg}_5\mathrm{Si}_6`$ $`\beta ^{\prime \prime }`$ phase Derlet et al. (2002).

The calculated bonding charge density for the (110) plane is shown in fig 11. There is a clear concentration of charge between the Al and Si atoms indicating strong Al-Si bonding. There is also a small build-up of charge between the Si and corner Mg atoms which constitutes a coupling between Mg2+ and AlSi- layers.

Evidence of the layered structure can also found in the DOS. Figs. 12 and 13 show the muffin tin decomposed total and partial DOS. The total DOS below the fermi level is characterized by two broad bands, separated by a 0.1 Ry band gap (fig. 12a). It is also evident that the major contribution to the total DOS stems from the Al and Si states. This is an indication that the Mg donates charge to the Al-Si network and is not strongly involved in the electronic bonding. Furthermore, the Si states give a relatively larger contribution to the total DOS than the Al states. From the angular decomposed Si DOS (fig. 12b) one can see that the contribution to the lower band comes mainly from Si $`s`$ states, while the higher band consists mainly of $`p`$ states with a modest mixing in of $`s`$-states. This can be seen as an indication of a partial hybridization of the Si states resulting in the Si bonding environment being partly covalent. This hybridization can also be identified in the Al DOS (fig. 13a). Compared to a metallic DOS Kittel (1996), there is a peak concentration of $`s`$ states with mixed in $`p`$ states slightly below 0 Ry resembling the $`\mathrm{L}_1`$ lobe of a covalent DOS. However for both Al and Mg (fig. 13b), the overall shape of the total atomic DOS seems roughly to follow the square-root relationship characteristic of a metallic DOS.

The origin of the band gap at the centre of the occupied states can be understood from the same line of reasoning as for the $`\beta `$ phase. In the $`U1`$ phase, the Si atoms are arranged as an hcp lattice with lattice parameters $`a=4.05`$ Å and $`c=6.74`$ Å. If we once again consider only the Si system (with a nearest neighbour distance of 4.05 Å), then little hybridization of the atomic Si $`s`$ and $`p`$ valence states can be expected and a similar DOS to the $`\beta `$ is seen. That is, the occupied valence states will be separated into two relatively narrow bands. By introducing the Al and Mg, the Al forms a tightly bound bonding network with the silicon via the donation of charge from the Mg atoms. This coupling results in a broadening (with respect to the artificial Si sub-structure) of the two occupied bands and in the removal of the band gap seen in $`\beta `$. This has its origins in the stronger mixed $`s`$ and $`p`$ character of the second occupied band (centred at $`0.2`$ Ry) of Si (fig. 12b), when compared to the strong $`p`$ character of the corresponding band for $`\beta `$ (fig. 4) — indicating strong hybridization between the Si $`p`$ states and the Al $`s`$ and $`p`$ states.

### IV.2 $`U2`$

The $`U2`$ primitive/conventional unit cell is orthorhombic (space group $`P_{nma}`$) with experimentally derived lattice parameters $`a=6.75\AA `$, $`b=4.05\AA `$, and $`c=7.94\AA `$. It contains 4 Mg, 4 Al, and 4 Si atoms giving the formula $`\mathrm{Mg}_4\mathrm{Al}_4\mathrm{Si}_4`$ (fig. 14). The coordinates of the symmetrically inequivalent atoms of the unit cell are shown in table 3. The morphology and size of the precipitate formed by the $`U2`$ phase is usually the same as for the $`U1`$ phase. The $`<100>`$ direction of the unit cell is oriented parallel to the $`<310>`$ direction of the fcc Al matrixAndersen et al. (2002b). By performing a volume relaxation of the unit cell, we get an equilibrium volume differing by less than 1% from the experimentally derived unit cell volume. A second order Birch fit gives a bulk modulus of 69.1 GPa, very close to the value obtained for the $`U1`$ phase.

The $`U2`$ phase is similar to the TiNiSi structure type Nuspl et al. (1996); Landrum et al. (1998). As for the $`U1`$ phase, the bonding in the $`U2`$ phase is characterized by the electropositive Mg atoms donating charge to the electronegative Al and Si atoms forming a tightly bound bonding network (fig. 15). However, unlike for the $`U1`$ phase, the electropostitive Mg atoms retain a large amount of charge and (using a simple chemical picture) there is not enough transfer for the electronegative elements to fulfill the octet rule Landrum et al. (1998). The Al-Si bond lengths are all $`2.59`$ to $`2.61\AA `$, while the bond lengths of the Mg-Si neighbour coupling is $`2.78`$ and $`2.86\AA `$.

Fig. 16 displays the bonding charge density maps for the (040), (020) and (0$`\overline{4}`$0) planes of the $`U2`$ unit cell. As for the $`U1`$ phase, there is a concentration of charge between Al and Si nearest neighbours making up the AlSi bonding network. This charge concentration appears to have more of a covalent, directional character, than the $`U1`$ Al-Si bond. Furthermore, one can identify a more pronounced charge concentration in between Mg-Si neighbours than for the $`U1`$ phase. This fits well with the general assumption for this structure class — that a large amount of Mg charge is retained. As for the U1 phase there is no indication of charge build-up between Al and Mg neigbours, indicating that the AlSi network is mainly coupled to the Mg atoms via the Si atoms.

This picture is further justified by the total and partial DOS shown in figs. 17 and 18. Comparing the DOS of the $`U2`$ phase with the DOS of the $`U1`$ phase one can identify several common features, supporting the proposition that the electronic bonding picture is very similar for the two phases. Note that the total DOS is the sum of 1 Mg, 1 Al, and 1 Si atom, as opposed to 1 Mg, 2 Al, and 2 Si atom for the $`U1`$ phase (fig. 12a). From the atom decomposed DOS of the $`U2`$ phase (fig. 17a) one can see that the Al and Si states are dominating with respect to the Mg states, and that the Si states give a slightly larger contribution than the Al states. This is consistent with the donation of charge from Mg to Al+Si atoms, as in the case of $`U1`$. However, the Mg DOS gives a larger contribution to the total DOS than for the $`U1`$ phase. This is an indication that there is less Mg to Al and Si charge transfer, than for the $`U1`$ Zintl-type phase. We remind the reader that it is difficult to quantify the charge transfer because the (L)APW+(lo) method does not use an atomic/localised basis set. Furthermore, the magnitude of the band gap separating the two occupied bands is slightly smaller for the $`U2`$ phase than for the $`U1`$ phase (0.075 Ry vs. 0.1 Ry for the $`U1`$ phase). Consequently, the region of hybridization around 0 Ry is increased for the Al and Si DOS (fig. 18 a and b) consistent with the pronounced covalent character of the Al-Si bonds found in the bonding charge density.

As for the $`U1`$ and $`\beta `$ phase, the origin of the band gap in the occupied states can be explained by the separation of the Si atoms with respect to each other, which for this structure forms a slightly distorted h.c.p structure (not shown here) with lattice parameters comparable to those of the $`U1`$ Si hcp sub-lattice.

## V Discussion and Concluding remarks

We have made a comparative study of the electronic structure of the $`\beta `$, $`\beta ^{\prime \prime }`$, $`U1`$, and $`U2`$ phases in the AlMgSi alloy precipitation sequence. The bonding in the $`\beta `$ phase is characterized by the covalent bonding between Si-Si nearest neighbour pairs and ionic/covalent bonding between Si-Mg nearest neighbour pairs. For the $`\beta `$ phase the bonding is dominated by the partly ionic Mg-Si bond. By calculating the heat of formation for the various precipitate phases further insight can be gained in the observed precipitation sequence. Presently we employ the formula:

$$\mathrm{\Delta }H=E_{\text{AlMgSi}}x_{\text{Mg}}E_{\text{Mg}}x_{\text{Si}}E_{\text{Si}}x_{\text{Al}}E_{\text{Al}},$$
(6)
where $`\mathrm{\Delta }H`$ is the formation energy/enthalpy per atom, $`E_{\text{AlMgSi}}`$ is the energy of the given AlMgSi compound, and $`E_{\text{Mg}}`$, $`E_{\text{Si}}`$, and $`E_{\text{Al}}`$ are the equilibrium ground state (zero temperature) energies per atom at of hcp Mg, dc Si, and fcc Al respectively. $`x_i`$ is the relative content of element $`i`$ in the compound. Results for the various phases in the precipitation sequence are shown in table 4 and figure 19. For the $`\beta `$ and $`\beta ^{\prime \prime }`$ phases we used the experimentally observed lattice parameters. For the $`\beta ^{\prime \prime }`$ we used the relaxed coordinates found from force minimization (table 1). We have also included the corresponding energy from a past model for the $`\beta ^{}`$ phase proposed by Matsuda et al.Matsuda et al. (2000).

It is difficult to draw any quantitative conclusions based on the relative formation energy of these phases in Aluminium, at finite temperatures, without taking into account the effects of the interface and entropy. The calculated energies are nonetheless very reasonable. As can be seen, only 3 phases ($`\beta `$, $`U1`$ and $`U2`$), give a negative formation energy. The $`U1`$ and $`U2`$ are lower in energy than $`\beta ^{\prime \prime }`$, but higher than the $`\beta `$ phase. $`\beta `$ shows the most negative energy, which is line with expectations, since it is the terminal equilibrium structure of the precipitation sequence. The energy of $`\beta ^{\prime \prime }`$ is slightly positive, but not markedly different from $`\beta `$, $`U1`$ and $`U2`$. $`\beta ^{\prime \prime }`$, which can be derived from an fcc super-cell, forms early in the precipitation sequence at temperatures above 150 C, but gradually disappears after further heat treatment Marioara et al. (2001).

We note the unusually high heat of formation for Matsuda’s $`\beta ^{}`$ model which has a hexagonal unit cell with space group $`\mathrm{P}\overline{6}2\mathrm{m}`$ with experimentally derived lattice parameters $`a=b=4.05\AA `$, $`c=6.74\AA `$ and the stoichiometry $`\mathrm{Mg}_2\mathrm{Si}`$. For the calculated formation energy, we used the experimentally derived lattice parameters and published coordinates. Upon volume relaxation, we however found that the structure reduced in volume by 26% corresponding to a heat of formation of 25 mRy. Furthermore, from a second order Birch fit, we find that for the experimental lattice constant, the Matsuda’s $`\beta ^{}`$ phase is under a negative pressure of approximately 6.9 GPa. Since the experimentally derived unit cell parameters are regarded as being extremely accurate it becomes difficult to see how such as phase might exist within the Al matrix.

To this date $`U1`$ and $`U2`$ are the only phases in the precipitation sequence of the AlMgSi alloy system containing Aluminium in addition to Magnesium and Silicon. The crystallographic structure type of both phases belong to structure families containing an abundance of documented compounds. A characteristic for both phases is the formation of a tightly bound network of Al and Si atoms made possible by at charge transfer from the electropositive Mg atoms. Considering the wide range of structures belonging to this family it is possible that these these type of structures containing other hardening elements than Mg and Si would form stable precipitates in Aluminium alloys. For example, it is quite remarkable that the $`U1`$ phase with Magnesium replaced by the rare earth element Europium produces a $`\mathrm{CaAl}_2\mathrm{Si}_2`$-type structure with lattice parameters and coordinates differing from the $`U1`$-phase by only few percent Schobinger-Papamantellos and Hulliger (1989). Suggested further work include the calculations on these two phases with iso-valent species substituted for the Si and Mg sites. In the future these calculations could also be used together with interfacial energies and entropy calculations to estimate the solid solubility of these phases in Aluminium.

Furthermore, the rule of thumb that successful aluminium precipitation hardened alloys have secondary and ternary elements that are larger and smaller than aluminium is not sufficiently accurate when applied to the stability of bulk phases in the precipitation sequence. The atomic radii of Al, Si, and Mg vary with the bonding environment. For example, in the equilibrium $`\beta `$ phase the atomic size of Si is increased and that of Mg correspondingly decreased in forming a dominating ionic/covalent Mg-Si bond. Also, using the appropriate values for ionic radii Kittel (1996), the Mg-Si bond length in the $`\beta `$ and $`\beta ^{\prime \prime }`$ phases are predicted to be considerably smaller than what the present calculations show. A more precise formulation would involve the ability of the secondary(A) and ternary(B) elements to form stable bonds (A-B, A-A, and B-B) with bond lengths which are comparable that of the Al-Al bond.

###### Acknowledgements.

This study has been performed as part of the research programs “KMB Heat Treatment Fundamentals” and “SUP Micro- and Nanostructure Based Materials Development”, supported by the Norwegian Research Council and collaborators from the industry. Parts of this work has been supported by the Norwegian Research Council through CPU time on the NOTUR Origin 3800.


text:
# A First Look at White Dwarf - M Dwarf Pairs in the Sloan Digital Sky Survey (SDSS)

## 1 Introduction

The formation and evolution of a low mass star in a binary system is a common phenomenon, and can lead to interesting results such as cataclysmic variables (Warner 1995), which are the most numerous UV and X-ray sources in the Galaxy, and Type Ia supernovae (Langer et al 2000), which are standard candles for cosmology. While the low mass secondaries play a pivotal role in determining the evolution of the system, it is not well understood how the binary environment affects the evolution and properties of the secondary.

Studies of pre-cataclysmic variable stars (or Pre-Cataclysmic Binaries – PCBs) (Bond et al 1985; Hillwig et al 2000) envision the following scenario for their evolution. The more massive primary in a wide binary star system evolves through the giant and planetary nebula phases, leading to significant mass loss and a “common envelope” encompassing both primary and secondary. A large amount of the system’s orbital angular momentum is lost during this phase, leading to a much closer binary orbit. Tidal interactions in the binary system will have led to synchronous rotation of the secondary for separations closer than 5-7 Rg, where Rg is the radius of the primary during its giant phase (Soker 2002), corresponding to an orbital period of $``$50 days (assuming an M5 secondary). The remaining binary may undergo further orbital decay via magnetic braking or gravitational radiation (see Warner 1995). Eventually, systems with very short periods will interact and become cataclysmic variables (CVs). CVs therefore represent only the short-period subset of all these binaries. In this paper we add to the growing database of PCBs, focusing on faint WD/M pairs. We hope to eventually build a large sample of PCBs which represents various phases of binary evolution, in order to test the validity of the above scenario, and the efficiency of the various physical mechanisms for angular momentum loss.

Note that we define the WD, the original high-mass component, as the primary in these systems. This contrasts with previous surveys, many of which defined the primary as the brighter star in a given filter, whether it be the WD or the M dwarf.

The Sloan Digital Sky Survey (SDSS; York et al 2000) database, with its accurate photometry in five broad bands covering the entire optical spectrum, is an outstanding resource for finding these pre-interacting systems. We have identified 109 WD/M binaries with g $`<`$ 20, using data available through 11/14/2001. Our efficiency at finding these objects using our current photometric targeting algorithm is $``$ 60 %, and our sample is augmented by serendipitous observations obtained through other target selection pipelines (e.g. quasars). Forty three of our systems are contained in the June 2001 SDSS Early Data Release (Stoughton et al 2002).

Currently known are 42 PCBs (Hillwig et al 2000), and roughly 800 widely separated, non-interacting WD/M binaries (Silvestri, Oswalt & Hawley 2002 and references therein). We expect to accumulate a sample of 700-1000 systems by the end of the survey, including $``$ 50 PCBs – doubling the current samples of these objects.

For each of our systems, we analyzed the SDSS discovery spectrum to determine the temperature and distance of the white dwarf primary and the spectral type and magnetic activity level of the low mass secondary. We are particularly interested in studying the dependence of the magnetic activity level on rotation rate, which increases in close binary systems where the secondary is in synchronous rotation. The magnetic activity of the secondary has been widely cited as a primary driver for angular momentum loss as stars evolve toward the period gap (Patterson 1984). We also investigated the influence of irradiation from the hot white dwarf on the activity level of the secondary.

To further compare activity levels across spectral type and rotation rate, we need to determine which pairs are in close binaries (thus exhibiting synchronous rotation). Population models (de Kool 1992; Warner 1995 and references therein) indicate that about 10% of the WD/M pairs should be close enough to enter the common envelope stage, with periods $``$ 4 hrs. We targeted the most active systems for followup observations at the 3.5m ARC telescope at Apache Point Observatory (APO), to search for radial velocity variations, and obtain orbital periods. For a 2 M WD progenitor, a secondary closer than $``$0.25 AU (P $``$ 50 days) will rotate synchronously (Soker 1996). The radial velocity of an M5 dwarf (M $``$ 0.2 M) in such a binary is $``$40 km s-1. With our radial velocity errors of $`\pm `$10 km s-1, we can therefore detect both synchronously and non-synchronously rotating binaries.

In Section 2 we describe our technique for selecting WD/M pairs from the SDSS photometric database. Section 3 presents our methods for deriving system parameters from the SDSS spectra. In Section 4 we discuss the results of this analysis, including our ongoing followup at APO. Section 5 summarizes our findings.

## 2 SDSS and APO Data

### 2.1 Photometric Selection of WD/M Pairs with SDSS

Our goal is to find the evolutionary precursors of CVs, systems on the verge of interacting, or only weakly interacting. For these systems we expect no significant flux from an accretion disk, and the composite spectrum of the candidate pair will be simply the superposition of the spectra of each component. Therefore we wish to isolate objects with blue spectral energy distributions (SEDs) at short optical wavelengths (hot WDs), and red SEDs at long wavelengths (M dwarfs). To accomplish this, we use the Sloan filter system (Fukugita et al 1996) to devise a photometric targeting algorithm.

Through a trial and error procedure during the commissioning year of the SDSS, we optimized a set of color selection criteria in the 4 SDSS colors: u - g $`<`$ 0.45 (blue), g - r $`<`$ 0.7 (blue), r - i $`>`$ 0.3 (red), and i - z $`>`$ 0.4 (red). We found that a practical magnitude cutoff of g $`<`$ 20 was required to obtain targeted SDSS spectra of our candidates with a sufficient signal to noise ratio for subsequent analysis.

Figure 1 shows the color distribution of our full sample of 109 systems. About half of the sample was found by visual inspection of SDSS spectra targeted for other purposes (e.g. galaxies, quasars), and many of these objects lie outside of our color cuts. However, we are unable to broaden the color selection criteria further, as too many interlopers from the stellar locus would be targeted, lowering our efficiency.

When our color criteria are applied to all SDSS objects for which spectra have been taken, $``$ 60% of these are WD/M pairs. We have repeat observations of 16 systems, one of which was observed 3 times. Note that our sample is not at all complete, as we take an average of one spectrum per SDSS field, which is 7 square degrees. The completeness of and biases in our sample will be addressed in an upcoming paper based on a larger SDSS WD/M sample.

Table 1 lists the SDSS photometry for the full sample, and Figure 2 illustrates the wide variety of systems we observed. Panels a) and b) are classic WD/M pairs, with SEDs that are blue at short wavelengths and red at long wavelengths. In some cases, the binary spectrum is dominated by the hot WD primary (panel c) or the cool M dwarf secondary (panels d and f), while panel e) shows nearly equal contributions over the whole spectrum. Panels a) and c) show very active systems with strong H$`\alpha `$ emission and emission cores inside the pressure-broadened WD Balmer absorption lines. The WD primaries in the systems from panels d) and f) are similar, however, d) is active while f) is not.

### 2.2 APO Data

We have observed 31 of our 109 pairs with the ARC 3.5-m telescope at APO, using the high resolution grating on the Double Imaging Spectrograph (DIS), which produces spectra in both the blue and red portions of the optical spectrum simultaneously. The dispersion is about 1.6Å pixel-1 in the blue, and 1.1 Å pixel-1 in the red.

Table 3 shows information for the 12 systems with APO spectra on at least two different nights. We tracked 5 objects for several hours, to look for variations on $``$ 4 hour timescales.

## 3 Analysis of SDSS Data

The temperature of the WD primary in each spectrum was determined by fitting the spectrum to hydrogen WD models (Hubeny & Lanz, 1995). Only one object, SDSS144258.48$`+`$001031.5, showed helium absorption lines. We fit the spectra between 3800 and 5000 Å, a region with very little contamination from the M dwarf secondary, to models with temperatures between 8,000K and 98,000K. The models are normalized to the spectrum at 4200Å, which is a fairly flat region of continuum between the strong H$`\delta `$ and H$`\gamma `$ absorption lines. The best fit was found using a simple chi-squared minimization technique. The errors are, conservatively, $`\pm `$ 2,000K for temperatures below 20,000K, and $`\pm `$ 4,000K above 20,000K. Once the best fit temperature was determined, the distance was found from the normalization factor that matches the model flux to the spectrum, assuming a WD radius of 8 $`\times `$ 108 cm (MWD = 0.6 M). We then subtracted the best-fit WD model from the binary spectrum to extract the spectrum of the low mass secondary.

The spectral types of the secondaries were determined from these extracted spectra, using primarily the TiO5 spectral index defined in Reid, Hawley and Gizis (1995). This index measures the strong TiO band near 7050 Å, and is well correlated with spectral type over the range K7 - M6. There are no secondaries with later types in our sample, which is expected given the extreme faintness of these very late M dwarfs (precluding their selection in our sample).

The distances obtained from the WD analysis were checked using the spectroscopic parallaxes of the M dwarf secondaries. We obtained absolute magnitudes for the M dwarfs as a function of spectral type using relations from Hawley et al (2002). The distances from the two methods generally agreed within 20%. However, the spectroscopic parallax of late-type stars using SDSS colors is under revision (West et al in preparation). We therefore use the distances derived from the WD temperature fitting.

The activity level of an M star is best measured by the ratio of its H$`\alpha `$ luminosity to its bolometric luminosity, L/Lbol. The H$`\alpha `$ equivalent width (EW) of each star was measured interactively using the IRAF66affiliation: IRAF is distributed by the National Optical Astronomy Observatories, which are operated by the Association of Universities for Research in Astronomy, Inc., under cooperative agreement with the National Science Foundation. “splot” routine, as continuum placement is difficult in these subtracted spectra, making automatic routines unreliable. The H$`\alpha `$ luminosities of our M dwarfs were calculated from the equivalent widths using the average H$`\alpha `$ continuum fluxes as a function of spectral type from the PMSU sample (Hawley, Gizis & Reid 1996). The bolometric luminosities were calculated from the spectral type averages given in Reid and Hawley 2000 (New Light on Dark Stars). Table 1 contains all of the WD and M dwarf parameters determined from our analysis. Table 2 bins the data by the spectral type of the secondary.

## 4 Results

Figure 3 shows the temperature distribution of the white dwarfs in our sample. The average temperature is $``$16,000K, slightly higher than the average for field WDs, which is $``$12,000K. The highest temperature seen is $``$44,000K $`\pm `$4,000K. Thus, using the cooling curves from Fontaine et al (2001) and assuming a mass of 0.6 M, the youngest (hottest) WD has a post main-sequence age of less than 0.1 Gyr. The coolest model was 8,000K, implying a post main-sequence age between 1 and 2 Gyr.

The average M spectral type of the low mass secondaries is $``$M4.5. The distribution of spectral types of the sample is shown in Figure 4. The rise toward spectral type M4-5 is due to the mass function of the local neighborhood, while the decrease toward later spectral types is due to the lower detection efficiency for fainter objects.

The distribution of derived distances is shown in Figure 5, assuming radii of 8 $`\times `$ 108 cm for all WDs, which should be a valid assumption for at least $``$ 80 % of WDs (Bergeron et al 1992). All but one of the systems lie within 200 pc, indicating that we are sampling the local disk population.

Roughly 60% of the systems are magnetically active (see Table 2). This is a higher fraction for these early-mid M dwarfs than that measured by the Palomar/MSU (PMSU) Survey of nearby stars (Hawley, Gizis & Reid 1996). As noted in the PMSU Survey and for later types by Gizis et al (2000), the fraction of stars that are active increases with later spectral type, until virtually every M7 star is active. Figure 6 shows that a higher fraction of early-type M dwarf components (M1 - M4) in WD/M pairs are active than their field counterparts. The level of activity as measured by L/Lbol (Figure 7) is similar for active M dwarfs both in the field and in our binaries, but there is a possible suggestion that the M2 - M4 dwarfs in the WD/M pairs are somewhat more active. However, our sample is relatively small, with few points earlier than M2.5 or later than M6, so these results must be confirmed by further studies with a larger sample.

Surprisingly, Figure 8 shows only a very weak correlation between M dwarf activity and WD temperature. A stronger correlation would imply that irradiation by the WD was enhancing the H$`\alpha `$ emission from the M dwarf.

The source of magnetic activity in late type stars is not well understood, and information on the rotation dependence of the activity would provide a constraint on existing dynamo models. Close binaries, with large radial velocity variations, should have strong enough tidal effects to induce the M dwarf secondary to rotate synchronously with the white dwarf primary, thereby significantly increasing its rotation rate. If the M dwarf dynamo is rotation dependent, the activity level should thus increase in these close binary systems.

For the 12 objects with at least two followup APO spectra, we can calculate a lower limit for the radial velocity variation. Figure 9 compares their activity levels with Vmin (i.e. their minimum $`v\mathrm{sin}i`$). No correlation is immediately apparent, but again more data are needed to provide a significant test. One object out of the five that we tracked for 3 continuous hours, SDSS112909.50$`+`$663704.4, was found to have a period of $``$4.1 hours (see Figure 10). No clear periodicity was found in the four other systems.

## 5 Conclusions

We have compiled a sample of 109 WD/M binaries from the first year of the Sloan Digital Sky Survey. We expect the sample to contain 700-1000 systems by the end of the five year survey. These systems, which are the precursors of cataclysmic variable stars and type Ia supernovae, span a wide range in WD temperature, M dwarf spectral type and chromospheric activity level.

The average temperature of our sample is $``$16,000K, considerably hotter than the average temperature for field WDs. The mean spectral class of the low-mass secondaries is $``$M4-M5. The fraction of active early-type M dwarfs in binaries appears higher than the active fraction among field M dwarfs (Figure 6). Their magnetic activity levels may also be elevated relative to those seen in the field M dwarfs (Figure 7). More data are required to confirm and quantify these suggestive results.

Since a higher orbital velocity implies a closer orbit, in which tidal interactions spin up the M dwarf secondary and possibly increase the efficiency of a rotation-dependent dynamo, we have begun to obtain orbital velocities for the brightest, most active systems in our sample.

Population models (de Kool 1992) predict that $``$ 10% of WD/M binaries should be pre-interacting with orbital periods of $``$ 4 hours. We have discovered one such system from the five which were followed for 3 or more hours. Eleven other systems show some evidence for velocity variations during scattered observations. We see no correlation between the minimum orbital velocity and the M dwarf activity level. However, much more data are needed to reach any firm conclusions. Our sample size is continuously increasing, and further analysis and followup observations are in progress. We encourage others to pursue velocities of the identified systems.

Funding for the creation and distribution of the SDSS Archive has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Aeronautics and Space Administration, the National Science Foundation, the U.S. Department of Energy, the Japanese Monbukagakusho, and the Max Planck Society. The SDSS Web site is http://www.sdss.org/.

The SDSS is managed by the Astrophysical Research Consortium (ARC) for the Participating Institutions. The Participating Institutions are The University of Chicago, Fermilab, the Institute for Advanced Study, the Japan Participation Group, The Johns Hopkins University, Los Alamos National Laboratory, the Max-Planck-Institute for Astronomy (MPIA), the Max-Planck-Institute for Astrophysics (MPA), New Mexico State University, University of Pittsburgh, Princeton University, the United States Naval Observatory, and the University of Washington.