text:
# Errors due to Finite Rise/fall Times of Pulses in Superconducting Charge Qubits

## I Introduction

Building a practical quantum computer with a large number of qubits has recently attracted much attention since the useful quantum algorithms Shor94 ; Grover97  and the quantum error-correction codes Shor95 ; Steane96  were developed. Basic quantum logic gates were implemented on trapped ions Cirac95 ; Monroe95 , QED cavities Turchette95 , and liquid state NMR Chuang98 ; Cory98 ; Jones98 . Solid-state devices have advantages of large scale integration, flexibility in the design, and easy connection to conventional electronic devices. Some of solid-state devices proposed are as follows: electron spins in quantum dots Loss98 , nuclear spins of donor atom in siliconKane98 , and ultrasmall Josephson junctions Shnirman97 ; Averin98 ; Mooij99 ; Makhlin99 ; Shnirman98 ; Makhlin01 .

Two types of superconducting qubits based on ultrasmall Josephson junctions were proposed. One is to use a number of excess Cooper pairs on a superconducting Cooper-pair box, called the superconducting charge qubit Shnirman97 ; Averin98 ; Makhlin99 . And the other is to utilize a single flux quantum of a superconducting loop, called the superconducting flux qubit Mooij99 . Here, we focus on superconducting charge qubits. In theoretical aspects, the measurement of superconducting charge states using a single electron transistor and decoherence due to coupling of qubits with environments were studied Shnirman98 ; Makhlin01 . Also the quantum leakage of superconducting charge qubits was pointed out Fazio99 . In experiments, Nakamura et al. Nakamura99  demonstrated the coherent oscillations of Cooper pairs on a superconducting Cooper-pair box. This corresponds to a rotation of a single qubit about $`x`$ axis. Thus the realization of two-qubit gates is an important step for making a superconducting quantum computer with a medium size.

In this paper, we investigate how the two-qubit gate could be implemented on superconducting charge qubits. The pulse sequence and operation times, and the simulation of two-qubit gates are reported. We examine errors due to finite rise/fall times of pulses in the realization of two-qubit gates.

Our paper is organized as follows. In Sec. II we introduce the Hamiltonians of superconducting charge qubits and of ideal qubits. A numerical method which allows us to simulate a quantum computer is presented. In Sec. III we explicitly show the pulse sequence necessary to implement two-qubit gates with superconducting charge qubits. Also numerical studies of two-qubit gates are presented. In Sec. IV we analyze errors due to finite rise/fall times of pulses on superconducting charge qubits. Finally, in Sec. V we summarize the results.

## II Hamiltonian of superconducting charge qubits

In this paper we consider two qubit systems: ideal qubits and superconducting charge qubits Makhlin99 ; Makhlin01 . By comparing two qubit systems, one can discover differences and similarities between them, which will be helpful to improve the design of superconducting charge qubits.

First consider an ideal model for a quantum computer. The Hamiltonian of the system of $`N`$ ideal qubits reads

$`H_q^{\mathrm{ideal}}(t)`$ $`=`$ $`{\displaystyle \underset{i=1}{\overset{N}{}}}𝐁_i(t)𝝈^{(i)}`$ (1)
$`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i<j}{}}J_{ij}(t)\left[\sigma _x^{(i)}\sigma _x^{(j)}+\sigma _y^{(i)}\sigma _y^{(j)}\right],`$
where $`𝝈^{(i)}=(\sigma _x^{(i)},\sigma _y^{(i)},\sigma _z^{(i)})`$ are Pauli matrices for the $`i`$-th qubit. Here the Zeeman coupling terms $`𝐁_i(t)`$ and the inter qubit couplings $`J_{ij}`$(t) can be turned on and off between zero and finite values in a controlled way.

Let us consider a system of $`N`$ superconducting charge qubits Shnirman97 ; Makhlin99 ; Makhlin01 . Each qubit consists of a single-Cooper-pair box with two ultrasmall Josephson junctions of capacitance $`C_J^0`$ forming a DC-SQUID ring and a gate electrode with capacitance $`C_g`$. The dynamics of a qubit is characterized by relevant energy scales; superconducting gap $`\mathrm{\Delta }`$, charging energy $`E_Ce^2/2(C_g+2C_J^0)`$, Josephson coupling energy $`E_J^0`$, and thermal fluctuation $`k_BT`$. Assume that the system is in the regime of $`E_C\mathrm{\Delta }`$ and $`k_BTE_C`$ in order to suppress quasi-particle tunneling or excitation. Also suppose that the system operates under the conditions $`E_JE_C`$ and $`C_gV_g/(2e)1`$. Then only two charge states $`\{|0,|1\}`$, no excess Cooper pair and one excess Cooper pair on the box, play a role and represent the qubit. A few ways of coupling charge qubits were proposed Shnirman97 ; Averin98 ; Makhlin99 ; Makhlin01 . Here we consider the coupling between charge qubits via the LC resonant circuit, where $`N`$ charge qubits are connected in parallel to a common inductor with inductance $`L`$. The Hamiltonian of the system of $`N`$ superconducting charge qubits is given by

$`H_q(t)`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{i=1}{\overset{N}{}}}\left[E_{Ci}\sigma _z^{(i)}+E_{Ji}(\mathrm{\Phi }_{Xi})\sigma _x^{(i)}\right]`$ (2)
$`{\displaystyle \underset{ij}{}}{\displaystyle \frac{E_{Ji}E_{Jj}}{E_L}}\sigma _y^{(i)}\sigma _y^{(j)},`$
where $`E_{Ci}=4E_C(C_{gi}V_{gi}/e1)`$ is turn on and off by applying the gate voltage $`V_{gi}`$ on the $`i`$-th gate electrode. Here $`E_{Ji}(\mathrm{\Phi }_{Xi})=2E_J^0\mathrm{cos}(\pi \mathrm{\Phi }_{Xi}/\mathrm{\Phi }_0)`$ is the effective Josephson energy of the $`i`$-th qubit, and controlled by the external flux $`\mathrm{\Phi }_{Xi}`$. $`\mathrm{\Phi }_0=\mathrm{}/2e`$ is a flux quantum. Here is $`E_L=\frac{\mathrm{\Phi }_0^2}{\pi ^2L}\left(\frac{2C_J^0}{C_{\mathrm{qb}}}\right)^2`$ with $`C_{\mathrm{qb}}^1=(2C_J^0)^1+C_g^1`$.

The Hamiltonian of Eq. (2) is similar to the ideal one of Eq. (1) except one fact. The system of superconducting charge qubits has only two kinds of independently controllable parameters, $`\{E_{Ci},E_{Ji}\}`$, whereas the ideal model has four kinds of controllable parameters $`\{B_{xi},B_{yi},B_{zi},J_{ij}\}`$. Coupling between two superconducting charge qubits could be realized by turning on both $`E_{Ji}`$ and $`E_{Jj}`$. This causes somewhat disadvantages in operating two-qubit gates, which will be discussed in Sec. IV.

The time-evolution of $`N`$ qubits is governed by a time-dependent Schrödinger equation

$`i\mathrm{}{\displaystyle \frac{}{t}}|\psi =H_q(t)|\psi ,`$ (3)
where $`|\psi =_{m=0}^La_m(t)|m`$ with $`L2^N1`$ and $`a_m𝐂`$. Its computational basis is represented by a tensor product of individual qubits, $`|m=|q_1|q_2\mathrm{}|q_N`$, where integer $`m=2^{N1}q_1+2^{N2}q_2+\mathrm{}+2^0q_N`$ and $`q_i\{0,1\}`$. Multiplying both sides of Eq. (3) by $`n|`$, one gets a set of $`N`$ coupled first-order ordinary differential equations for the function $`a_n(t)`$

$`\dot{a}_n={\displaystyle \frac{i}{\mathrm{}}}{\displaystyle \underset{m=0}{\overset{L}{}}}H_{nm}(t)a_m,`$ (4)
where $`H_{nm}(t)n\left|H_q(t)\right|m`$. We developed a program which allows us to solve the above initial value problem and to simulate a quantum computer. Our code is based on the Runge-Kutta method press92 .

A quantum logic gate can be realized by controlling the time-evolution operator which acts on selected qubits for a fixed period of time. A way to control the time-evolution operator is to turn on or off each terms in the Hamiltonians of Eq. (1) or (2). For an example, a pulse $`P(t)`$, which is applied as $`B_{\alpha i}(t)=B_{\alpha i}^0P(t)`$ with $`\alpha =x,y,z`$, makes it possible to change a term of the Hamiltonian from zero to finite value $`B_{\alpha i}^0`$ for a finite time $`\tau `$. In this paper, we consider a rectangular pulse $`P_{\mathrm{rec}}(t)`$ with width $`\tau =t_bt_a`$ and a unit height, modeled by a kink and anti-kink pair

$`P_{\mathrm{rec}}(t)={\displaystyle \frac{1}{2}}\left[\mathrm{tanh}\left({\displaystyle \frac{tt_a}{ϵ/2}}\right)+\mathrm{tanh}\left({\displaystyle \frac{t_bt}{ϵ/2}}\right)\right],`$ (5)
where the rise (or fall) time is about $`2ϵ`$. That means $`P_{\mathrm{rec}}(t_a+ϵ)\mathrm{tanh}(2)0.964`$ for $`\tau >2ϵ`$. The smaller $`ϵ`$ one takes, the sharper rectangular pulse one obtains. Fig. 1 illustrates a sequence of two rectangular pulses.

## III CNOT gate with superconducting charge qubits

In this section, we explicitly show how the controlled-not (CNOT) gate with superconducting charge qubits could be implemented by applying a sequence of pulses. What makes the CNOT gate so important is that the CNOT gate (or any nontrivial two-qubit gate) and single-qubit gates form a universal set of logic operation. Also a general two-qubit controlled-$`U`$ gate can be built up of two CNOT gates and three single-qubit gates. The CNOT gate acting on two qubits $`i`$ and $`j`$ is represented by the unitary matrix in the basis of $`\{|0_i0_j,|0_i1_j,|1_i0_j,|1_i1_j\}`$

$`U_{\mathrm{CNOT}}^{ij}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 0& 1\\ 0& 0& 1& 0\end{array}\right),`$ (10)
where $`i`$ and $`j`$ are indices for a control bit and a target bit, respectively.

### III.1 CNOT gate with ideal qubits

First consider the implementation of the CNOT gate for an ideal system whose Hamiltonian is given by Eq. (1). The primary two-qubit gate acting on the $`i`$-th and $`j`$-th qubits could be implemented by turning on the coupling $`J_{ij}(t)`$. It can be written in the basis $`\{|0_i0_j,|0_i1_j|1_i0_j,|1_i1_j\}`$ as

$`U_{2\mathrm{b}}^{ij}(\gamma )=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& \mathrm{cos}\gamma & i\mathrm{sin}\gamma & 0\\ 0& i\mathrm{sin}\gamma & \mathrm{cos}\gamma & 0\\ 0& 0& 0& 1\end{array}\right),`$ (15)
where $`\gamma J_{ij}t/\mathrm{}`$. Then, the CNOT gate for an ideal qubit system can be realized by combining single-qubit gates and the primary two-qubit gate $`U_{2\mathrm{b}}^{ij}`$

$`U_{\mathrm{CNOT}}^{ij}`$ $`=`$ $`\mathrm{H}_ie^{i\frac{\pi }{4}}e^{i\frac{\pi }{4}\sigma _x^{(j)}}e^{i\frac{\pi }{4}\sigma _x^{(i)}}`$ (16)
$`\times U_{2\mathrm{b}}^{ij}(\frac{\pi }{4})e^{i\frac{\pi }{2}\sigma _x^{(i)}}U_{2\mathrm{b}}^{ij}(\frac{\pi }{4})\mathrm{H}_i,`$
where $`\mathrm{H}_i`$ is a Hadamard gate acting on the $`i`$-th qubit.

Fig. 2 depicts the time-evolution of qubits under the action of CNOT gate for the ideal model. The lower part of Fig. (2) illustrates the sequence of rectangular pulses which generates the CNOT gate of Eq. (16). Here is taken $`|\psi _{\mathrm{in}}=|3=|11`$, i.e., $`|a_{11}|^2=1`$, as an initial state (denoted by $``$). After implementing the CNOT gate, one obtains $`|a_{10}|^2=1`$, corresponding to state $`|10`$, labeled by $`+`$. In Fig. 2, it takes a long time to implement the single-qubit gate $`e^{i\frac{\pi }{4}\sigma _x^{(j)}}`$ because $`B_{\alpha i}`$’s are taken to be always positive or zero.

### III.2 CNOT gate with superconducting charge qubits

Let us consider the implementation of the CNOT gate with superconducting charge qubits. The coupling between the $`i`$-th and $`j`$-th qubits can be switched on by turning on Josephson couplings, $`E_{Ji}`$ and $`E_{Jj}`$, and by turning off charging energy terms, $`E_{Ci}=E_{Cj}=0`$. Then, the interaction Hamiltonian between two charge qubits, for an example qubits 1 and 2, can be written as

$`H_{\text{ph}}`$ $`=`$ $`{\displaystyle \frac{E_J}{2}}\sigma _x^{(1)}{\displaystyle \frac{E_J}{2}}\sigma _x^{(2)}E_{\mathrm{int}}\sigma _y^{(1)}\sigma _y^{(2)},`$ (17)
where $`E_{\mathrm{int}}E_J^2/E_L`$ with an assumption of $`E_{J1}=E_{J2}E_J`$. For a moment, suppose that switching $`E_{Ji}`$’s on and off could be done instantaneously. As will be discussed later, a finite time in switching on $`E_{Ji}`$’s gives rise to an error in implementing two-qubit gates.

The basic two-qubit gate is given by the time-evolution operator under the Hamiltonian of Eq. (17) for a finite time $`\tau `$

$`U_{\mathrm{ph}}(\tau )=e^{iH_{\mathrm{ph}}\tau /\mathrm{}}.`$ (18)
Let us transform the above time-evolution operator under the unitary operator $`R_y\mathrm{exp}\left[i\frac{\pi }{4}(\sigma _y^{(1)}+\sigma _y^{(2)})\right]`$ as

$`U_{\mathrm{ph}}^{}R_y^{}U_{\mathrm{ph}}R_y=e^{iH_{\mathrm{ph}}^{}\tau /\mathrm{}},`$ (19)
where the Hamiltonian transformed, $`H_{\mathrm{ph}}^{}R_y^{}H_{\mathrm{ph}}R_y`$ reads

$`H_{\mathrm{ph}}^{}`$ $`=`$ $`{\displaystyle \frac{E_J}{2}}\sigma _z^{(1)}{\displaystyle \frac{E_J}{2}}\sigma _z^{(2)}E_{\mathrm{int}}\sigma _y^{(1)}\sigma _y^{(2)}`$ (20a)
$`=`$ $`E_{\mathrm{int}}\left(\begin{array}{cccc}a& 0& 0& 1\\ 0& 0& 1& 0\\ 0& 1& 0& 0\\ 1& 0& 0& a\end{array}\right).`$ (20b)
Here $`aE_J/E_{\mathrm{int}}=E_L/E_J`$. Note that the Hamiltonian of Eq. (20) does not mix the subspace spanned by $`|00`$ and $`|11`$ with that by $`|01`$ and $`|10`$. Then the time-evolution operator of the transformed Hamiltonian $`H_{\mathrm{ph}}^{}`$ becomes

$`U_{\mathrm{ph}}^{}=`$
$`\left(\begin{array}{cccc}\mathrm{cos}\theta +n_z\mathrm{sin}\theta & 0& 0& in_x\mathrm{sin}\theta \\ 0& \mathrm{cos}\varphi & i\mathrm{sin}\varphi & 0\\ 0& i\mathrm{sin}\varphi & \mathrm{cos}\varphi & 0\\ in_x\mathrm{sin}\theta & 0& 0& \mathrm{cos}\theta n_z\mathrm{sin}\theta \end{array}\right),`$ (21)
where the time evolution corresponds to the rotation about the $`y`$-axis by angle $`\varphi E_{\mathrm{int}}\tau /\mathrm{}`$ in the subspace spanned by $`|01`$ and $`|10`$. Also in the subspace spanned by $`|00`$ and $`|11`$, the time evolution gives rise to the rotation about the axis $`(n_x,0,n_z)`$ by angle $`\theta \sqrt{1+a^2}E_{\mathrm{int}}\tau /\mathrm{}`$, where

$`n_z`$ $``$ $`{\displaystyle \frac{a}{\sqrt{a^2+1}}}={\displaystyle \frac{E_L}{\sqrt{E_L^2+E_J^2}}},`$ (22a)
$`n_x`$ $``$ $`{\displaystyle \frac{1}{\sqrt{a^2+1}}}={\displaystyle \frac{E_J}{\sqrt{E_L^2+E_J^2}}}.`$ (22b)
By choosing the appropriate values $`E_L`$, $`E_J`$, and $`\tau `$, one can control the rotation angles $`\varphi `$ and $`\theta `$. Let us consider special angles $`\varphi =\frac{\pi }{4}(2m1)`$ with $`m=1,2,\mathrm{}`$, and $`\theta =n\pi `$ with $`n=0,\pm 1,\pm 2,\mathrm{}`$. Note that the case of even $`n`$ is same to that of odd $`n`$ up to the global phase $`e^{i\pi }`$. These angles can be obtained by taking the evolution-time $`\tau =\frac{\pi }{4}(2m1)\frac{\mathrm{}}{E_{\mathrm{int}}}`$, which is given by the relation of angle $`\varphi `$. Then one has $`\theta =\frac{\pi }{4}(2m1)\sqrt{1+a^2}=\frac{\pi }{4}\sqrt{1+E_L^2/E_J^2}=n\pi `$. Thus we obtain the relation between $`E_L`$ and $`E_J`$

$`{\displaystyle \frac{E_L}{E_J}}=\sqrt{\left({\displaystyle \frac{4n}{2m1}}\right)^21},`$ (23)
where it could be accomplished by tunning $`E_J`$ or $`E_L`$. Note that one can not take $`n=0`$. Also the evolution-time is given by

$`\tau ={\displaystyle \frac{\pi \mathrm{}}{4E_J}}\sqrt{(4n)^2(2m1)^2}.`$ (24)
We have the basic two-qubit gate $`U_{\mathrm{ph}}=R_yU_{\mathrm{ph}}^{}R_y^{}`$ written by

$`U_{\mathrm{ph}}={\displaystyle \frac{1}{2}}\left(\begin{array}{cccc}1+e^{i\varphi }& 0& 0& 1e^{i\varphi }\\ 0& 1+e^{i\varphi }& 1e^{i\varphi }& 0\\ 0& 1e^{i\varphi }& 1+e^{i\varphi }& 0\\ 1e^{i\varphi }& 0& 0& 1+e^{i\varphi }\end{array}\right),`$ (25)
where $`\varphi =\frac{\pi }{4}(2m1)`$ with $`m=1,2,\mathrm{}`$. By combining $`U_{\mathrm{ph}}`$ with single-qubit gates, the controlled phase flip gate $`U_{\mathrm{CPF}}^{ij}`$, operating on the $`i`$-th and $`j`$-th qubits, can be realized by

$`U_{\mathrm{CPF}}^{ij}`$ $`=`$ $`e^{i\varphi \sigma _z^{(i)}}e^{i\varphi \sigma _z^{(j)}}U_{\mathrm{ph}}^{ij}e^{i\pi \sigma _z^{(i)}/2}U_{\mathrm{ph}}^{ij},`$ (26)
where $`\varphi `$ is the value given above. By using the controlled phase flip gate $`U_{\mathrm{CPF}}^{ij}`$ and the Hadamard gate Hj on qubit $`j`$, the CNOT gate can be implemented by

$`U_{\mathrm{CNOT}}^{ij}`$ $`=`$ $`\mathrm{H}_jU_{\mathrm{CPF}}^{ij}\mathrm{H}_j.`$ (27)
Taking $`m=1`$ and $`n=3`$, one gets $`E_L=\sqrt{143}E_J`$. This value satisfies the physical condition $`E_L10E_J`$ Makhlin99 ; Makhlin01 ; Dreher99 . The typical time-scale for operating the basic two-qubit gate is given by $`\tau =\frac{\pi \mathrm{}}{4E_J}\sqrt{143}`$, which is much longer than the operation time of the single-qubit gate, $`\tau _{\mathrm{op}}=\mathrm{}/E_J`$ or $`\mathrm{}/E_{Ci}`$. Fig. 3 shows time-evolution of qubits when $`|\psi _{\mathrm{in}}=|11`$ is taken as an input (denoted by $``$). After operating the CNOT gate on two qubits, one gets the output $`|\psi _{\mathrm{out}}=|10`$, which is labeled by $`\times `$.

## IV Errors due to finite rise/fall times of pulses

Up to now, it was assumed that the effective Josephson coupling $`E_{Ji}(t)`$ could be switched on (or off) instantaneously. This means a pulse applied is a perfect rectangular one and has no rise/fall times. However, in reality it takes finite times to turn on (or off) pulses fully. While $`E_{Ji}(t)`$’s being switched on (or off), the Hamiltonian of Eq. (17) $`H_{\mathrm{ph}}`$ becomes time-dependent and the $`H_{\mathrm{ph}}`$’s at different times do not commute, i.e., $`[H_{\mathrm{ph}}(t),H_{\mathrm{ph}}(t^{})]0`$ for $`tt^{}`$. This causes an error in the realization of two-qubit gates with superconducting charge qubits. If Hamiltonian’s at different times commute, the shape of the pulse applied is of no importance. The ideal qubits and single-qubit rotations of superconducting charge qubits are free from this problem.

The Magnus expansion Slichter92  provides the means of representing the time-evolution operator of a time-dependent Hamiltonian $`H(t)`$ as

$`U(t)=Te^{\frac{i}{\mathrm{}}_0^\tau 𝑑\tau H(\tau )}=e^{\frac{i}{\mathrm{}}\overline{H}\tau },`$ (28a)
where $`T`$ denotes the time-ordering operator. Here the average Hamiltonian is given by $`\overline{H}=\overline{H}^{(0)}+\overline{H}^{(1)}+\mathrm{}`$ where
$`\overline{H}^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{\tau }}{\displaystyle _0^\tau }𝑑t_1H(t_1),`$ (28b)
$`\overline{H}^{(1)}`$ $`=`$ $`{\displaystyle \frac{i}{2\tau \mathrm{}}}{\displaystyle _0^\tau }𝑑t_2{\displaystyle _0^{t_2}}𝑑t_1[H(t_2),H(t_1)].`$ (28c)
If the rise/fall time is short, then the approximation to the first term $`\overline{H}^{(0)}`$ is good. Otherwise the higher order terms become significant. If the pulse for the rise time is modeled by $`P(t)=t/2ϵ`$, the first term $`\overline{H}^{(1)}`$, which arises when the basic two-qubit gate $`U_{\text{ph}}`$ is implemented, can be written by

$`\overline{H}^{(1)}{\displaystyle \frac{E_J^3ϵ}{15E_L\mathrm{}}}\left(\sigma _z^{(1)}\sigma _y^{(2)}+\sigma _y^{(1)}\sigma _z^{(2)}\right).`$ (29)
The error due to the finite rise/fall times of pulses is quantified by the gate fidelity Poyatos97

$`F=\psi _{\mathrm{in}}|U^{}\rho _{\mathrm{out}}U|\psi _{\mathrm{in}}=|\psi _{\mathrm{out}}|\psi _{\mathrm{out}}^ϵ|^2,`$ (30)
where $`U`$ is the unitary operator corresponding to the ideal gate when perfect rectangular pulses are applied. On the other hand, the unitary operator generated by pulses with finite rise/fall times transforms the input state $`|\psi _{\text{in}}`$ into the density operator of the imperfect output state $`\rho _{\mathrm{out}}`$. In our case the gate fidelity is nothing but the square of the overlap between the perfect output state $`|\psi _{\text{out}}=U|\psi _{\text{in}}`$ and the imperfect output state $`\rho _{\text{out}}=|\psi _{\text{out}}^ϵ\psi _{\text{out}}^ϵ|`$. Assuming that Eq. (29) is small, it is straightforward to calculate the gate fidelity

$$F|\psi _{\text{in}}|e^{i\overline{H}^{(1)}\tau /\mathrm{}}|\psi _{\text{in}}|^21ϵ^2\mathrm{\Delta }\eta ^2$$
(31)
where $`\mathrm{\Delta }\eta ^2\psi _{\text{in}}\left|\eta ^2\right|\psi _{\text{in}}(\psi _{\text{in}}\left|\eta \right|\psi _{\text{in}})^2`$ is the dispersion of $`\eta \overline{H}^{(1)}\tau /\mathrm{}ϵ`$. For small rise/fall time, error grows quadratically in rise/fall time $`2ϵ`$.

Fig. 4 shows the numerical study of the error due to finite rise/fall times in implementing the CNOT gate with superconducting charge qubits. For ideal qubits, if finite idle times between successive pulses are applied, the correct CNOT gate is realized. However, for superconducting charge qubits, although a finite idle time reduces the error, there still exits the error due to finite rise/fall times. As shown in Fig. 4 (b), the probability of getting the correct state reduces quadratically in rise/fall time $`2ϵ`$ in agreement with the theoretical prediction of Eq. (31).

In Nakamura et al.’s experiment Nakamura99 , the rise/fall times of the pulse was about $`3040`$ ps at the top of the cryostat. If the effective Josephson coupling is $`E_J50`$ $`\mu `$eV, then the timescale of single-qubit $`x`$ rotation is about $`\mathrm{}/E_J1`$ ps. This means that the rise/fall times of the pulses should be less than 1 ps in order for two-qubit gates to be implemented correctly. It should be noted that this type of errors is caused by the coupling scheme between two qubits. Thus this problems can be solved by improving the design of devices, i.e., by introducing new ways of two-qubit couplings.

## V Summary

We have studied the dynamics of qubits on the action of CNOT gates for ideal and superconducting charge qubits. We have explicitly shown how the CNOT gate could be implemented for superconducting charge qubits. It is found that the error in implementing two-qubit gates with superconducting charge qubits grows quadratically in finite rise/fall times. Thus it is necessary to keep the rise/fall times small or to find new ways of coupling two qubits other than the coupling scheme via the common inductor.

## Acknowledgments

The author would like to thanks Y. Makhlin, A. Shnirman, and G. Schön for helpful discussions. This work was partially supported by Korea Science and Engineering Foundation.


format: text

text:
# Binder cumulants of an urn model and Ising model above critical dimension

## Abstract

Solving numerically master equation for a recently introduced urn model, we show that the fourth- and sixth-order cumulants remain constant along an exactly located line of critical points. Obtained values are in very good agreement with values predicted by Brézin and Zinn-Justin for the Ising model above the critical dimension. At the tricritical point cumulants acquire values which also agree with a suitably extended Brézin and Zinn-Justin approach.

The concept of universality and scale invariance plays a fundamental role in the theory of critical phenomena . It is well known that at criticality the system is characterized by critical exponents. Calculation of these exponents for dimension of the system $`d`$ lower than the so-called critical dimension $`d_c`$ is a highly nontrivial task . On the other hand for $`d>d_c`$ the behaviour of a given system is much simpler and critical exponents take mean-field values which are usually simple fractional numbers.

However, not everything is clearly understood above the critical dimension. One of the examples is the Ising model ($`d_c=4`$) where despite intensive research serious discrepancies between analytical  and numerical  calculations still persist. Of particular interest is the value of the Binder cumulant at the critical point. Several years ago Brézin and Zinn-Justin (BJ) calculated this quantity using field theory methods  and only recently numerical simulations for the $`d=5`$ model are able to confirm it . Some other properties of the Ising model above critical dimension are still poorly explained by existing theories. For example, the theoretically predicted leading corrections to the susceptibility disagree even up the sign with numerical simulations .

In addition to direct simulations of the nearest-neighbour Ising model, there are also some other ways to study the critical point of Ising model above critical dimension. For example, Luijten and Blöte used the model with $`d3`$ but with long-range interactions . Using such an approach they confirmed with good accuracy the BJ predictions for the Binder cumulant.

In the present paper we propose yet another approach to the problem of cumulants above critical dimension. Namely, we calculate fourth- and sixth-order cumulants at the critical point of a recently introduced urn model . Albeit structureless, this model exhibits a mean-field Ising-type symmetry breaking. Along an exactly located critical line, the obtained values are in a very good agreement with values predicted by BJ. Let us notice that our calculations: (i) are not affected by the inaccuracy of the location of the critical point which is a serious problem in the case of the Ising model (ii) are based on the numerical solution of the master equation which offers a much better accuracy than Monte Carlo simulations. Moreover, we calculate these cumulants at the tricritical point and show that the obtained values are also in agreement with suitably extended calculations of BJ. That both the Ising model and the (structureless) urn model have the same cumulants is a manifestation of strong universality above the upper critical dimension: at the critical point not only the lattice structure but also the lattice itself becomes irrelevant. What really matters is the type of symmetry which is broken and since in both cases it is the same $`Z_2`$ symmetry, the equality of cumulants follows.

Our urn model was motivated by recent experiments on the spatial separation of shaken sand . In the present paper we are not concerned with the relation with granular matter and a more detailed justification of rules of the urn model is omitted . The model is defined as follows: $`N`$ particles are distributed between two urns A and B and the number of particles in each urn is denoted as $`M`$ and $`NM`$, respectively. Particles in a given urn (say A) are subject to thermal fluctuations and the temperature $`T`$ of the urn depends on the number of particles in it as:

$$T(x)=T_0+\mathrm{\Delta }(1x),$$
(1)
where $`x`$ is a fraction of a total number of particles in a given urn and $`T_0`$ and $`\mathrm{\Delta }`$ are positive constants. (For urn A and B, $`x=M/N`$ and $`(NM)/N`$, respectively.) Next, we define dynamics of the model :
(i) One of the $`N`$ particles is selected randomly.
(ii) With probability $`\mathrm{exp}[\frac{1}{T(x)}]`$ the selected particle changes urns, where $`x`$ is the fraction of particles in the urn of a selected particle.

To measure the difference in the occupancy of the urns we define

$$ϵ=\frac{2MN}{2N}=\frac{M}{N}\frac{1}{2}.$$
(2)
In the steady state the flux of particles changing their positions from A to B equals to the flux from B to A. Since the selected particles are uncorrelated, the above requirement can be written as:

$$<M>\mathrm{exp}[\frac{1}{T(<M/N>)}]=<NM>\mathrm{exp}[\frac{1}{T(<(NM)/N>)}],$$
(3)
or equivalently

$$(\frac{1}{2}+<ϵ>)\mathrm{exp}[\frac{1}{T(\frac{1}{2}+<ϵ>)}]=(\frac{1}{2}<ϵ>)\mathrm{exp}[\frac{1}{T(\frac{1}{2}<ϵ>)}].$$
(4)
Analysis of eq. (4) shows  that on the ($`\mathrm{\Delta },T_0`$) phase diagram symmetric ($`ϵ=0`$) and asymmetric ($`ϵ0`$) solutions are separated by the critical line which is given by the following equation

$$T_0=\sqrt{\mathrm{\Delta }/2}\mathrm{\Delta }/2,0<\mathrm{\Delta }<\frac{2}{3}.$$
(5)
The critical lines terminates at the tricritical point: $`\mathrm{\Delta }=\frac{2}{3},T_0=\frac{\sqrt{3}1}{3}`$. Let us notice that a random selection of particles implies basically the mean-field nature of this model. Consequently, at the critical point $`\beta =1/2`$ and $`\gamma 1`$ (measured from the divergence of the variance of the order parameter), which are ordinary mean-field exponents. However, the calculation of the dynamical exponent $`z`$ gives $`z=0.50(1)`$  while the mean-field value is 2. We do not have convincing arguments which would explain such a small value of $`z`$. Presumably, this fact might be related with a structureless nature of our model.

Defining $`p(M,t)`$ as the probability that in a given urn (say A) at the time $`t`$ there are $`M`$ particles, the evolution of the model is described by the following master equation

$`p(M,t+1)=`$ $`{\displaystyle \frac{NM+1}{N}}p(M1,t)\omega (NM+1)+{\displaystyle \frac{M+1}{N}}p(M+1,t)\omega (M+1)+`$ (7)
$`p(M,t)\{{\displaystyle \frac{M}{N}}[1\omega (M)]+{\displaystyle \frac{NM}{N}}[1\omega (NM)]\}\mathrm{for}M=1,2\mathrm{}N1`$
$`p(0,t+1)=`$ $`{\displaystyle \frac{1}{N}}p(1,t)\omega (1)+p(0,t)[1\omega (N)],`$ (8)
$`p(N,t+1)=`$ $`{\displaystyle \frac{1}{N}}p(N1,t)\omega (1)+p(N,t)[1\omega (N)],`$ (9)
where $`\omega (M)=\mathrm{exp}[\frac{1}{T(M/N)}]`$. Supplementing the above equations with initial conditions one can easily solve them numerically.

Cumulants that we calculate are defined as

$$x_4=\frac{<ϵ^4>}{<ϵ^2>^2},x_6=\frac{<ϵ^6>}{<ϵ^2>^3}$$
(10)
where

$$<ϵ^n>=\underset{M=0}{\overset{N}{}}(\frac{M}{N}\frac{1}{2})^np(M,\mathrm{})$$
(11)
and the symbol of infinity indicates that we take the long-time (steady-state) solutions of the master equation (9). Calculations are made for $`\mathrm{\Delta }=\frac{1}{8},\frac{1}{4},\frac{1}{2}`$ and $`\frac{2}{3}`$ and for each $`\mathrm{\Delta }`$ the value of $`T_0`$ is calculated from eq. (5). Thus, the last point is the tricritical point and the remaining ones are critical points. Numerical results are presented in Figs. 1-4.

Before discussing our results further, let us briefly describe the BJ approach. To calculate cumulants above the critical dimension they used the Ginzburg-Landau-Wilson model. Then, they calculate the effective action restricting the expansion only to the homogeneous contributions (the lowest-mode approximation). Since at criticality the quadratic (in the order parameter) term vanishes in such an expansion and the leading term is quartic which implies that the probability distribution has the form $`p(x)\mathrm{e}^{x^4}`$, where $`x`$ is a rescaled order parameter. Calculations of moments for such a distribution are then elementary and one obtains

$$x_4=\frac{1}{8\pi ^2}[\mathrm{\Gamma }(\frac{1}{4})]^42.188440\mathrm{},x_6=\frac{3}{8\pi ^2}[\mathrm{\Gamma }(\frac{1}{4})]^46.565319\mathrm{}.$$
(12)
The fact that one can restrict the expansion of the free energy to the lowest order term is by no means obvious . Such a restriction leads to the correct results but only above critical dimension where the model behaves according to the mean-field scenario with fluctuations playing negligible role. For $`d<d_c`$ additional terms in the expansion are also important and cumulants take different value. Numerical confirmation of the above results requires extensive Monte Carlo simulations, and a satisfactory confirmation was obtained only for $`x_4`$ .

Omitting detailed field theory analysis, we can extend the BJ approach to the tricritical point. At such a point also the quartic term vanishes which makes the sixth-order term the leading one and the probability distribution gets the form $`p(x)\mathrm{e}^{x^6}`$. Simple calculations for such a distribution yield

$$x_4=\frac{\mathrm{\Gamma }(\frac{5}{6})\mathrm{\Gamma }(\frac{1}{6})}{\mathrm{\Gamma }(\frac{1}{2})^2}=2,x_6=\frac{\mathrm{\Gamma }(\frac{1}{6})^3}{6\mathrm{\Gamma }(\frac{1}{2})^3}5.162113\mathrm{}.$$
(13)
The BJ results (12)-(13) are indicated by small arrows in Figs. 1-2. Even without any extrapolation one can see, especially for critical points, a good agreement with our results. Data in Figs. 1-2 shows strong finite-size corrections. To have a better estimations of asymptotic values in the limit $`N\mathrm{}`$ we assume finite size corrections of the form

$$x_{4,6}(N)=x_{4,6}(\mathrm{})+AN^\omega .$$
(14)
The least-square fitting of our finite-$`N`$ data to eq. (14) gives $`x_{4,6}(\mathrm{})`$ which agree with BJ values (12)-(13) within the accuracy better than 0.1%. A better estimation of the correction exponent $`\omega `$ is obtained assuming that $`x_{4,6}(\mathrm{})`$ are given by the BJ values. The exponent $`\omega `$ equals then the slope of the date in the logarithmic scale as presented in Figs. 3-4. Our data shows that for the critical(tricritical) point $`\omega =\frac{1}{2}(\frac{1}{3})`$.

Let us notice that leading finite-size corrections to the Binder cumulant in the $`d=5`$ Ising model at the critical point are also of the form $`N^{0.5}`$ (with $`N`$ being the linear system size) . Moreover, for the tricritical point but $`d<d_c`$ the probability distribution is known to exhibit a three-peak structure , which is different than the single-peak form $`p(x)\mathrm{e}^{x^6}`$.

In summary, we calculated fourth- and sixth-order cumulants at the critical and tricritical points in an urn model which undergoes a symmetry breaking transition. Our results confirm that, as predicted by Brézin and Zinn-Justin, the critical probability distributions of the rescaled order parameter has the form $`p(x)\mathrm{e}^{x^4}`$. Similarly, for the tricritical point our results suggest that $`p(x)\mathrm{e}^{x^6}`$.

Although in our opinion convincing, the results are obtained using numerical methods. It would be desirable to have analytical arguments for the generation of such probability distributions. It seems that for the presented urn model this might be easier than for the Ising-type models. Let us notice that for the simplest urn model, which was introduced by Ehrenfest , the steady-state probability distribution can be calculated exactly in the continuum limit of the master equation and the result has the form $`p(x)\mathrm{e}^{x^2}`$, where $`x`$ is now proportional to the difference of occupancy $`ϵ`$. In the Ehrenfest model there is no critical point and we expect that a distribution of the type $`\mathrm{e}^{x^2}`$ might characterize our model but off the critical line (in the symmetric phase). We hope that when suitably extended, an analytic approach to our model might extract critical and tricritical distributions as well. Such an approach is left as a future problem.

###### Acknowledgements.

This work was partially supported by the Swiss National Science Foundation and the project OFES 00-0578 ”COSYC OF SENS”.


text:
# Heavy Majorana neutrinos at a very large electron-proton collider.

## I Introduction

The recent experimental results on solar and atmospheric neutrino properties NIS  confirm that standard neutrinos have small masses and oscillate. However, the nature of those massive neutrinos - if they are Dirac or Majorana particles - is still an open question. This is connected with lepton number conservation and to the general symmetries of some extention of the standard model. The smallness of neutrino masses is generally understood as a consequence of the see-saw mechanism. This brings the question of the possibility of new heavy neutrino states. So far none of these new states was experimentally observed PDG ; ZUB  with masses up to $`M_N`$ 100 GeV. For higher masses there are many suggestions of experimental possibilities in the next high energy hadron-hadron colliders YSP ; ALI , in electron-positron linear accelerators ARS ; DJO ; ACV , in neutrinoless double-beta decay BIL . Another experimental possibility is found on electron-proton colliders with a positron in the final state BUM . This means a very clear experimental signature that violates the electron leptonic number by two units. In this paper we investigate the discovery potential for recently proposed electron- proton colliders for detecting new neutral heavy Majorana leptons as suggested by extentions of the Standard Model. Since 1998-1999 the option of an electron-proton collider for the Very Large Hadron Collider at Fermilab operating with the proton booster has been considered. We study the reaction $`e^{}pe^++`$ jets and present estimates for the signal and Standard Model background.

## II The Very High Energy Electron-Proton Collider

There is a growing work on the options for a future Very Large Hadron Collider (VLHC)SYP . The low magnetic field (approximately 2 T) collider design in a tunnel can also be used for a future high magnetic field (approximately 10 T) collider. These designs could also be used for electron-proton collidersVLH ; DER . In the low-field version electron-proton collisions would be produced by an $`80`$ GeV electron ring and a $`3`$ TeV proton booster ($`\sqrt{s}=980`$ GeV), with a luminosity $``$= 2600 pb-1/yr, while in the high-field option one could reach a center-of-mass energy $`\sqrt{s}=6320`$ GeV and a luminosity $``$ = 1400 pb-1/yr with the collision of a $`200`$ GeV electron beam with a $`50`$ TeV proton beam. These two options are considered in the present work.

## III The model

New heavy leptons are present in many possible extentions of the standard model. In the superstring inspired scenario SIR  we have new doublet vector leptons. In mirror models MIR  there are new fermions with right handed doublets and left handed singlets. Some see-saw extentions of the standard model include new weak-singlet fermions BET .

The general mixing of light and heavy Majorana neutrinos has new mixing parameters as well as new possible CP violating phasesFRI . It is well known that in neutrinoless double beta decay the theoretical estimates depend strongly on mixing angles and CP phases BEL . As a consequence, the comparison between the experimental limits and theory is model dependent. This is also the case for other processes. In this paper we are interested in the study of the discovery potential of new electron-proton colliders for new heavy Majorana neutrinos. From the experimental point of view,instead of start comparisons of very particular models with a great number of parameters and hypothesis one must start looking for the lightest of the new possible Majorana neutrinos, coupled to the electron family. This reduces the number of unknown parameters to two; the heavy neutrino mass and one new mixing angle, that should be determined experimentally.

In this paper we will adopt a general parametrisation that applies to all the above models. We will suppose that mixing angles for heavy-to-light neutrinos and new heavy neutrino masses are independent parameters BUC . The light neutrinos couplings to the neutral Z are given by $`g_{V,A}=g_{V,A}^{SM}sin^2\theta _{\text{mix}}/2`$. Taking the experimental values for $`g_{V,A}`$ and the Standard Model predictions we obtain a small upper bound for $`\theta _{\text{mix}}`$. A recent estimate ACV  gives $`\mathrm{sin}^2\theta _{\text{mix}}<0.0052`$ with $`95\%`$ C.L. This limit value is used throughout this paper for all curves and distributions. The decay of the Majorana neutrino into muons and taus are also possible. In ref. ACV  it was found that for the muon sector the experimental mixing angle limit with $`95\%`$ C.L. is much smaller than for the electron sector.

Most extended models predict new fermions and new gauge vector bosons. Since we presently have no signal for new interactions, we will make the hypothesis that the new heavy neutrino states behave as $`SU_L(2)U_Y(1)`$ basic representations and interact mostly with the standard model gauge bosons. We can resume the new particle interactions in the neutral and charged current lagrangians:

$$_{nc}=\frac{g}{4c_W}sin\theta _{mix}Z_\mu \overline{\psi _N}\gamma ^\mu \left(1\gamma _5\right)\psi _\nu +h.c..$$
(1)
and

$$_{cc}=\frac{g}{2\sqrt{2}}sin\theta _{mix}W_\mu \overline{\psi _N}\gamma ^\mu \left(1\gamma _5\right)\psi _e+h.c.$$
(2)
and similar terms for the other leptonic families. In the Majorana neutrino case, each completely neutral heavy lepton is coupled to all charged leptons but we are considering only the lower mass state.

The decay modes for these leptons, in the Majorana case YSP  must include both signatures $`N\mathrm{}^{}W^\pm `$ and $`N\nu _{\mathrm{}}(\overline{\nu }_{\mathrm{}})Z`$, with $`\mathrm{}=e,\mu ,\tau `$. In this work we have considered the Majorana neutrino decay only to the first family since for the other families we must add more mixing parameters. The Majorana neutrino branching ratios are shown in Table I and explicit formulas for the widths are given in ref. PAN

## IV Results

In Fig. 1 we show the structure of the Feynman diagrams that contribute for the signal. We have included in our results all the Majorana off-shell and on-shell contributions. The calculations for the signal were done using the high energy package CompHep HEP  with CTEQ95 structure functions, after implementation of the extended model. We call attention to the fact that the N width is very narrow for lower heavy neutrino masses and increase for higher masses ( see Table II). Thus propagator effects must be handled with care.

The cross sections for on and off-shell heavy Majorana neutrino production as a function of its mass are displayed in Fig. 2 for $`\sqrt{s}=`$ 980 GeV and Fig. 3 for $`\sqrt{s}=`$ 6320 GeV. In all our results we have done first general detector cuts for the final positron $`E_{e^+}>5`$ GeV and $`|\eta _{e^+}|<3.5`$, where $`\eta _{e^+}`$ is the pseudo rapidity defined relative to the initial proton.

In order to reproduce the signal and background as close as possible to experimental data we have employed the package PYTHIA 6.1 PIT  to generate the quark hadronization and decay from the Majorana neutrino production as well as the corresponding standard model background. The signal events were generated first at parton level with the CompHep package as above mentioned and treated as external process to the PYTHIA program. In this way we have simulated the W, Z, quark and hadron decays. The Standard Model background was estimated leaving open all the possible channels in the PYTHIA program. The positrons come mainly from charm decay. In both cases we did not allow pions to decay.

In the final visible particles we have imposed the following criteria to select events:

a) For $`|\eta _{e^+}|<3.5`$ and $`E_{e^+}>5`$ GeV we must have the number of $`e^+>0`$,

b)For $`|\eta _{h_i}|<3.5`$ and $`E_{h_i}>5`$ GeV we must have the number of hadrons $`>5`$.

With these two conditions, more than 99 $`\%`$ of the background events were eliminated for both energies.

In order to reconstruct the Majorana neutrino invariant mass we have employed the algorithm cluster PYCELL from the PYTHIA package, with 60 azimuthal divisions and 26 divisions in the electron beam direction. For $`\sqrt{s}`$= 980 GeV we have applied the condition $`0.7<\sqrt{(\mathrm{\Delta }\eta )^2+(\mathrm{\Delta }\varphi )^2}`$ while for $`\sqrt{s}`$= 6320 GeV this condition was changed to $`0.3<\sqrt{(\mathrm{\Delta }\eta )^2+(\mathrm{\Delta }\varphi )^2}`$ with $`|\eta |<3.5`$ for both energies.

After these conditions, we have further demanded that:

c) Number of clusters $`>`$ 1,

d) Number of clusters with 1 track = 1,

e) Cluster invariant mass in the range (70, 90) GeV.

This last requirement was done to identify the W coming from the heavy Majorana neutrino decay. After these additional cuts, the background estimated for 650 k events for $`\sqrt{s}1`$ TeV and 2,4 M events for $`\sqrt{s}6`$ TeV was reduced to practically zero. This point was explicitly verified by generating 10 times more events than expected by the luminosity. For the signal, the same cuts imply the results shown in Table II. We also show in the last column of Table II. the significance $`S/\sqrt{S+B}`$ for each case, with ”S” the number of signal events and ”B” the number of background events.

If we choose the clusters obtained with the surviving events then we can reconstruct the Majorana neutrino invariant mass as shown in Figs. 4 and 5.

## V Conclusions

The present work shows that heavy Majorana neutrinos could be produced in the electron-proton option for the VLHC. The very clear signature of a positron in the final state allows to separate efficiently the signal events from the background and to reconstruct the Majorana neutrino invariant mass. In the $`\sqrt{s}1`$ TeV option we can reach Majorana masses up to 300 GeV. For the $`\sqrt{s}6`$ TeV option this value could be as high as 1500 GeV. This is a conservative estimate since we are taking as limiting value $`S/\sqrt{S+B}>5`$.

###### Acknowledgements.

This work was partially supported by the following Brazilian agencies: CNPq and FAPERJ.
Figure Captions

1. General Feynman graphs for single heavy Majorana contribution to $`e^{}pe^++W+q`$.
2. Cross section for single heavy Majorana neutrino production at $`\sqrt{s}=`$ 980 GeV using basic cuts $`E_{e^+}>5`$ GeV and $`|\eta _{e^+}|<3.5`$.
3. Same as figure 2 for $`\sqrt{s}=`$ 6320 GeV.
4. Reconstruction of the invariant Majorana mass distribution for $`\sqrt{s}=`$ 980 GeV after hadronization using cuts and PYCELL algorithm.
5. Same as figure 4 for $`\sqrt{s}=`$ 6320 GeV.


text:
# First-Principles-Based Thermodynamic Description of Solid Copper Using the Tight-Binding Approach

## I INTRODUCTION

Density-functional theory (DFT) first-principles electronic-structure methods describe anomaly-free solids such as elemental copper successfully. They achieve high accuracy for quantities such as bulk properties, surface relaxation and lattice dynamics of the surface, as well as the epitaxial Bain path and elastic constants. DFT methods are routinely used to compute the zero-temperature internal energy, $`\mathrm{\Phi }_0(V)`$, but also can be used to calculate the free energy contributions from the ions, $`F_I(V,T)`$, and the electrons, $`F_E(V,T)`$, resulting in a complete equation of state,

$$F(V,T)=\mathrm{\Phi }_0(V)+F_I(V,T)+F_E(V,T).$$
(1)
However, the required computational effort is expensive, and an alternative efficient evaluation at all volumes and temperatures would be desirable.

In this paper we use the computationally less demanding tight-binding (TB) total energy model in conjunction with well chosen first-principles calculations. In particular, we use the functional fitting forms developed at the U.S. Naval Research Laboratory (NRL) for computing the total energy within the TB formalism, i.e., without an external potential. The model is fit to and accurately reproduces a set of first-principles calculations with a speed-up of many orders of magnitude. In addition, transferability (i.e., a TB parameterization that is accurate for a wide variety of crystal structures and atomic arrangements) has been successfully demonstrated for semiconductors as well as for simple and transition metals. We believe that the TB method can be used as a highly accurate, but computationally more efficient, surrogate for a full first-principles-based approach to calculate the equation of state for solids.

Copper is frequently used as a test material for theoretical methods. In this paper we have (1) developed an improved fit for copper that is accurate for phonons, and (2) used this model to calculate a wide range of temperature- and volume-dependent thermodynamic quantities.

Copper is furthermore widely employed as a pressure standard in high-pressure research. This use is based on correcting $`P(V)`$ data taken along the shock Hugoniot  to room temperature. Such corrections employ model assumptions about the volume dependence of the Grüneisen parameter $`\gamma (V)`$, which is difficult to measure independently. Shock heating increases with pressure, making the corrections more significant at high pressure. It is therefore important to develop theoretical techniques for accurate prediction of $`\gamma `$ for copper at high pressure.

Phonons play a major role in the calculations of thermodynamic quantities, and the TB fits are adjusted to more accurately calculate them. Structures corresponding to high-symmetry phonon modes are shown here to aid in refining the model; the resulting phonon density of states can then be used to determine the free energy and hence all thermodynamic quantities of interest. The precision required to calculate phonon frequencies is an order of magnitude higher than that for the lattice constant or bulk modulus, making this a stringent test for the validity of the tight-binding approach in general and the copper model in particular.

The ion–ion free energy of Eq. 1 is often separated into harmonic and anharmonic parts,

$$F_I(V,T)=F_H(V,T)+F_A(V,T).$$
(2)
Normally, the harmonic component is not a function of volume, but is calculated from the effect of small displacements about the zero-temperature equilibrium lattice. In our calculations, we use the quasi-harmonic approximation, which considers small displacements at any fixed volume (lattice constant) within the harmonic approximation, and hence our phonon frequencies become volume dependent. However, our phonon frequencies are calculated at zero temperature for any given volume, and are not temperature dependent.

The anharmonic part of the free energy involves terms that arise from the potential energy of the lattice when it is expanded beyond the harmonic part to higher than second order. Such terms are needed at high temperatures, when the phonon amplitudes are large, and ultimately lead to melting. They are also needed to explain thermal expansion effects when the harmonic part is based on the equilibrium volume. The quasi-harmonic approximation can handle thermal expansion and the Grüneisen parameter accurately through the volume dependence of the phonons at low temperature. At sufficiently high temperatures, the quasi-harmonic approximation breaks down when the phonon amplitudes become large, and additional anharmonic phonon-phonon corrections are necessary (as indicated in Eq. 2). We have not included these anharmonic types of effects in our calculations. Hence we always set $`F_A(V,T)=0`$, and our calculations will become less reliable at very high temperatures (near melting).

In the following section we introduce the basic ideas of the tight-binding method and the first-principles method used to generate the fitting database, and then describe our TB fitting procedures. In the subsequent section we present calculated results for the thermodynamic properties and compare them with experiment.

## II FITTING THE MODEL

### A Tight-binding electronic structure

The tight-binding approach is essentially a parameterized version of the first-principles calculations and hence is orders of magnitude more computationally efficient. In DFT methods the secular equation,

$$H\psi _{i,v}=ϵ_{i,v}S\psi _{i,v},$$
(3)
is constructed directly from approximate solutions to the full many-body Hamiltonian, and involves a self-consistent potential that is solved iteratively; whereas in the TB approach the elements of the Hamiltonian (and the overlap matrix) themselves have been parameterized. Only two-center terms are considered. For the non-orthogonal tight-binding model described here this requires 73 fitted parameters.

Of those parameters, thirty each are used to describe the inter-site matrix elements of the Hamiltonian and of the overlap matrix. For each combination of symmetries $`(ll^{}m)`$ the form is

$`h_{ll^{}m}(r)`$ $`=\left(a_{ll^{}m}+b_{ll^{}m}r\right)e^{c_{ll^{}m}^2r}f_c(r),`$ (4)
$`s_{ll^{}m}(r)`$ $`=\left(\overline{a}_{ll^{}m}+\overline{b}_{ll^{}m}r\right)e^{\overline{c}_{ll^{}m}^2r}f_c(r),`$ (5)
where $`f_c=1/(1+e^{2(rr_0)})`$ is a multiplicative factor included to ensure a smooth cutoff with increasing distance. In our calculations we have set $`r_0=16.0`$ Bohr radii.

The remaining 13 parameters determine the on-site terms, which allows the parameterization to be applied to structures not included in the fitting database. A measure of the valence electron density,

$$\rho =\underset{ij}{}e^{\lambda ^2r_{ij}}f_c(r_{ij}),$$
(6)
where $`r_{ij}`$ is the interatomic distance, serves to describe the on-site energy,

$$e_\alpha =e_\alpha ^0+e_\alpha ^1\rho ^{2/3}+e_\alpha ^2\rho ^{4/3}+e_\alpha ^3\rho ^2,$$
(7)
for the three orbital types $`\alpha `$, i.e., $`s`$, $`p`$, and $`d`$. These terms are somewhat similar to an embedded-atom-like form in that the energy changes depending on the nearby arrangements of atoms, and may approximately account for self-consistency effects as the atoms move around.

### B Full potential LAPW method

The first-principles quality of the tight-binding model results from fitting to full potential linear augmented plane wave (LAPW) calculations using the reliable WIEN97 program suite. The parameters for the first-principles calculations are listed in Table I.

The LAPW method divides space into spherical regions centered on the atoms and the remaining interstitial region. The radius of the spheres, the muffin-tin radius $`R_m`$, must be chosen such that the spheres do not overlap. The basis functions used to represent the wave function are adapted to the regions: radial solutions to the Schrödinger equation in the spheres, plane waves in the interstitial region. The wave functions then are found iteratively within density-functional theory, constrained to match at the boundaries of the different regions.

### C Initial Fitting Procedure

We first fit the TB method to predict energy differences between the ground-state and non-equilibrium structures. The fitting database included first-principles energies calculated for the cubic structures. In addition to the total energies of these structures, it proved to be crucial to fit the energy bands at high-symmetry points in reciprocal space. By decomposing the electronic wave function in terms of the symmetry character of the eigenvalues the bands are guided to the correct ordering.

The total energies and the band energies can be calculated by starting with a very crude initial tight-binding model that ignores intersite terms; the errors are then minimized utilizing standard nonlinear least squares algorithms.

Figure 1 shows the $`T=0`$ phonon dispersion for fcc copper calculated with the initial model. The long-wavelength modes near $`\mathrm{\Gamma }`$ are well described, the short-wavelength modes near the zone boundary display somewhat high frequencies, in particular the longitudinal modes.

The reasonable agreement for phonons near the zone center $`\mathrm{\Gamma }`$ can be understood by considering the elements of the fitting database. The bulk modulus, i.e., a linear combination of the elastic constants, is implicitly included in the fit. While this does not guarantee accurate elastic constants, i.e., good agreement for the slopes of the dispersion near $`\mathrm{\Gamma }`$, it does set the right scale. Furthermore, the fit includes the bcc structure, which is related to the fcc crystal by a tetragonal strain corresponding to the long-wave-vector limit of the longitudinal mode in the $`[00\xi ]`$ direction. The database lacks any information related to the short-wave-vector modes.

### D Fitting procedure with distorted structures

In order to construct a model with an improved phonon dispersion the database was expanded to include additional information on the phonons, in particular, structures that are snapshots of the crystal deformed by particular phonon modes, i.e., frozen phonons. The undistorted and distorted crystal structures are treated on the same footing in the first-principles calculations and the fitting procedure, implicitly including the differences in energy and hence the frequencies of the phonon modes.

The longitudinal and the transverse mode at the high-symmetry point X ($`𝐪=(0,0,1)`$) were chosen because of the large discrepancy in frequency (see Fig. 1) and because the distorted structures require only a doubling of the unit cell. These distorted structures are considered as additional, distinct structures in the database, to be fit to over a range of volumes.

The initial fit for copper already contains some of the character of distortions related to the elastic constants: the bulk modulus is explicitly included in the energy as a function of volume, and the tetragonal distortion of the fcc crystal is somewhat reflected by fitting to the bcc structure. For completeness, tetragonally- and trigonally-distorted fcc crystals were added to the fit as distinct structures. These additional structures barely influence the model resulting from the fit; however, the fitting process converges much more quickly when they are included.

The cubic structures that were included in the initial fit differ from each other by an energy scale of fractions of electron volts. Phonons require a model tuned to discern energies on a scale that is approximately an order of magnitude smaller. This could be a problem since the minimization procedure tends to ignore small energy differences. For frozen phonons at the zone boundary, where neighboring atoms move against each other, it turns out that amplitudes which are still within the harmonic regime can produce energies that differ from the undistorted structure by fractions of electron volts. The distortions corresponding to elastic constants, however, need to be exaggerated for them to give large enough energy differences. The trigonal distortion used here compresses the base angle from 90 to 75, while the tetragonal distortion changes the $`c/a`$ ratio from unity to 1.9.

Figure 2 shows the energy values in the fitting database alongside those of the initial and improved tight-binding models. The volumes of the first-principles calculations are limited to structures where the muffin-tin radius $`R_m`$ is smaller than the nearest-neighbor distance, particularly for the strongly-distorted fcc structures the choice of $`R_m=2.0`$ a.u. prohibits strong compression. No such limitations exist for the tight-binding approach; the volumes for which the model is appropriate will become clear in the next section.

Figure 3 shows the errors in the improved model’s fit. Compared to the initial fit, errors for the simple, cubic structures remain about the same. The errors for the tetragonally-distorted structures are small around the equilibrium volume (11.93 Å3), but show a tendency to increase as the crystal is compressed. The form of the matrix elements (Eq. 4) cannot be expected to allow a high-quality fit at all volumes; indeed when only a subset of data points are included in the fit the errors show no radical change.

Including the distorted structures in the fit improves the transferability of the model. Figure 4 shows the improved agreement between tight-binding and first-principles energies for the diamond structure, which is not included in the fit. The transferability to a structure of such a different coordination is not guaranteed, and our initial model did not reproduce the diamond energies well, nor did the model of the NRL group.

Figure 5 shows the phonon dispersion calculated with the improved model. Including the distorted fcc structures clearly refines the agreement with the measured values, though the curves do not overlap perfectly: the dispersion of the low-lying transverse modes in the $`[0\xi 1]`$ direction shows a different character, and the high-frequency longitudinal modes remain somewhat large. The discrepancy of the longitudinal frequency at L suggests including this data point in the fit. However, a first-principles, frozen-phonon calculation of this mode shows better agreement with the tight-binding model than with experiment and was therefore not added to the database.

Figure 6 shows the phonon density of states calculated with the improved model. The general shape agrees with the data calculated from the Born-von Kármán force constants fitted to the experimental phonon dispersion along high-symmetry directions. The difference in maximum frequencies and the peak near 7 THz can be attributed to the discrepancy in the dispersion of the longitudinal mode near L in the $`[\xi \xi \xi ]`$ direction. The tight-binding density of states displays more structure around 4 THz, which may be due to modes in low-symmetry directions that are not part of the experimental force-constant model.

The distorted structures added to the fit indeed make for a model that is better suited for phonon calculations. However, while the additional constraints improve the total energies described by the model, the electronic band structure deteriorates. Figure 7 shows the electronic band structure along two sample high-symmetry directions of fcc copper at the experimental volume. While the initial model agrees well with the first-principles band structure, the model improved for thermodynamic quantities loses the good agreement. The resulting electronic density of states, shown in Fig. 8, shows the same discrepancy; however, the density of states at the Fermi energy is quite similar, which is important for the temperature-dependent influence of the electrons (see below). It is possible that a better or more flexible functional form for the distance dependence of the intersite Hamiltonian and overlap matrices are necessary to keep the good transferability and the good agreement with the individual energy bands.

## III CALCULATIONS WITH THE TB MODEL

### A Force Constants

The force constants are calculated from the tight-binding model by the direct-force method, which relies on evaluating the forces on all atoms in a simulation cell in which a reference atom $`(\mathrm{𝟎},i)`$ has been displaced. The large simulation cell consists of primitive cells transposed by vectors $`\mathrm{}`$. Due to periodic boundary conditions on the simulation cell, the force on an atom $`(\mathrm{},j)`$ is in response to the displaced reference atom $`(\mathrm{𝟎},i)`$ as well as its images transposed by vectors $`𝐋`$,

$$𝐅(\mathrm{},j)=\underset{𝐋}{}\underset{¯}{\varphi }(\mathrm{},j;𝐋,i)𝐮(\mathrm{𝟎},i),$$
(8)
i.e., we are actually calculating the cumulant force constant,

$$\underset{¯}{\varphi }^C(\mathrm{},j;\mathrm{𝟎},i)=\underset{𝐋}{}\underset{¯}{\varphi }(\mathrm{},j;𝐋,i).$$
(9)
The components of the cumulant force constant are, in the harmonic regime, given by

$$\varphi _{\alpha \beta }^C(\mathrm{},j;\mathrm{𝟎},i)=\frac{F_\beta (\mathrm{},j)}{u_\alpha (\mathrm{𝟎},i)}\frac{F_\beta (\mathrm{},j)}{u_\alpha (\mathrm{𝟎},i)}.$$
(10)
The calculation of the Hellmann-Feynman forces in the tight-binding approach is achieved by evaluating the analytic derivatives of the Hamiltonian and overlap matrix elements and of the onsite terms.

The frequencies of the phonons with a given wave vector $`𝐪`$ are found by diagonalizing the corresponding dynamical matrix, $`\underset{¯}{D}\left(𝐪\right)`$, which in turn is the Fourier transform of the system’s force constants,

$$D_{\alpha \beta }\left(𝐪\right)=\frac{1}{M}\underset{\mathrm{},j,i}{}\varphi _{\alpha \beta }^C(\mathrm{},j;\mathrm{𝟎},i)e^{i𝐪(𝐫_j𝐫_i+\mathrm{})}.$$
(11)
If the force constants are known for every pair of atoms, the dynamical matrix and hence the frequencies are easily evaluated for any wave vector.

The direct-force method is exact within the quasi-harmonic approximation and the particular model if the forces vanish inside the simulation cell. Computational resources often limit the system size such that the forces do not vanish at the boundary, as in the calculations presented here; the calculations are correct for wave vectors commensurate with the simulation cell and are a good approximation for intermediate values of q.

The evaluation of Eq. 11 involves only the atoms in the simulation cell. In cases where the edge of the simulation cell limits the range of the forces before they actually vanish special care must be taken to ensure that no symmetries are lost. In particular, the inversion symmetry is lost if the reference atom $`(\mathrm{𝟎},i)`$ is not in the center of the simulation cell. Atoms that break the inversion symmetry with respect to the reference atom have to be duplicated (with adjusted weight) and transposed with a basis vector of the simulation cell to reinstate the symmetry.

The cubic symmetry of the fcc crystal allows the calculation of the force constants at a particular volume with a single displacement of the basis atom. Distorted fcc structures no longer have the cubic symmetry, the calculation of the force constants therefore requires the forces to be evaluated for the basis atom displaced in all three Cartesian directions separately. For all calculations the simulation cell contained 108 atoms and a mesh of $`4\times 4\times 4`$ k-points was used.

### B Thermodynamics

As indicated by Eq. 1, the free energy is the internal energy from the tight-binding calculation with entropic terms added from the electrons and the ions. In both terms the relevant physical quantity is the density of states (DOS). The electronic DOS, $`n(E)`$, the occupation of which is given by the Fermi distribution $`f(E,T)=[e^{(EE_f)/(k_BT)}+1]^1`$, determines the electrons’ contribution to the entropy,

$$S_{el}(T)=k_B[f\mathrm{ln}f+(1f)\mathrm{ln}(1f)]n(E)𝑑E.$$
(12)
The phonon DOS, $`g(\omega )`$, contributes through the zero-point energy,

$$U_{zero}=\frac{1}{2}_\mathrm{\Omega }\mathrm{}\omega g(\omega )𝑑\omega ,$$
(13)
as well as the temperature-dependent free energy,

$$F_H(T,V)=k_BT_\mathrm{\Omega }\mathrm{ln}[1e^{\mathrm{}\omega /k_BT}]g(V,\omega )𝑑\omega .$$
(14)
Both DOS are calculated in the tight-binding calculations, the phonon DOS comes from evaluating the dynamical matrix for a fine mesh of wave vectors in the first Brillouin zone, and the electronic DOS results similarly from evaluating the eigenvalues on a fine mesh. The contribution to the free energy from the electrons is on the order of a few percent of that of the phonons, although at low temperatures (where both contributions are very small) and small volumes the percentage rises to about 10%.

Figure 9 shows the resulting free-energy as a function of volume for temperatures between 0 K to 1400 K (at ambient pressure copper melts at 1356 K; melting is an anharmonic effect that lies outside the scope of the quasi-harmonic treatment) in 100 K increments. A comparison with the free energy for the bcc phase shows the fcc structure at lower free energy for all temperatures and volumes, indicating that the model agrees with experiment in that respect.

The free energy as a function of volume and temperature determines the thermal expansion. The temperature-dependent lattice constant derived from the tight-binding model is shown in the inset of Fig. 9 along with the experimental values. As is typical for GGA-calculations, the tight-binding model overestimates the equilibrium volume by 1.4%. The calculated linear expansion coefficient is compared to experimental data in Fig. 10 and shows good agreement, in particular the characteristic temperature, which is determined by the phonon characteristic temperatures (see below).

The shape of the free energy as a function of volume and temperature directly provides the temperature-dependence of the bulk modulus, $`B(T)`$, which is calculated by fitting a second order Birch equation of state.

The bulk modulus is related to two of the elastic constants by $`B(T)=\frac{1}{3}(C_{11}(T)+2C_{12}(T))`$. The temperature-dependence of the other elastic constants, determined by tetragonal ($`C_{11}(T)C_{12}(T)`$) and trigonal strain ($`C_{44}(T)`$), is calculated in two steps. The $`T=0`$ value results from calculating the energy of the appropriately strained crystal and finding the quadratic change. The temperature-dependent value is calculated by following the same procedure but with the free energy of the strained crystal, which is found by calculating the phonon DOS of that structure and then evaluating Eq.14.

Figure 11 shows the calculated temperature-dependence of the elastic constants. The thermal expansion is implicitly included in the calculation of the bulk modulus, a derivative with respect to the volume. For the trigonal and tetragonal strain the thermal expansion is accounted for by finding the equilibrium volume for each temperature and then calculating the effect of the strain on the free energy of that volume. Figure 12 shows the calculated $`T=0`$ elastic constants as a function of volume.

The phonon characteristic temperatures, which are defined as moments of the phonon density of states,

$$\mathrm{ln}(k_B\theta _0)=\mathrm{ln}(\mathrm{}\omega )_{BZ},$$
(15)
$$k_B\theta _1=\frac{4}{3}\mathrm{}\omega _{BZ},$$
(16)
and

$$k_B\theta _2=\left[\frac{5}{3}(\mathrm{}\omega )^2_{BZ}\right]^{1/2},$$
(17)
are shown in Fig. 13. The approximate rule of thumb $`\theta _2\theta _1e^{1/3}\theta _0`$ holds nicely for the calculated values (inset).

At temperatures below the phonon characteristic temperatures individual phonon modes must be considered separately, because they contribute to the crystal’s thermal properties with weights depending on their frequency relative to the temperature. The weight of a mode of branch $`s`$ with wave vector q is determined by the heat capacity for that mode,

$$c_s(𝐪)=\frac{}{T}\frac{\mathrm{}\omega _s(𝐪)}{e^{\beta \mathrm{}\omega _s(𝐪)}1}.$$
(18)
The sum of these individual heat capacities as a function of temperature agrees well with calorimetric data; the comparison is plotted in Fig. 14 in terms of the Debye temperature $`\theta _D`$, which is found such that the Debye model’s heat capacity

$$c_V=9k_B\left(\frac{T}{\theta _D}\right)^3_0^{\theta _D/T}\frac{x^4e^x}{(e^x1)^2}𝑑x$$
(19)
is the same as the heat capacity calculated for the tight-binding model at the same temperature.

The shape of the Debye temperature plotted against temperature remains very similar with compression; the curve itself is shifted upwards with the same volume-dependence as the characteristic phonon temperatures.

The heat capacity of each individual phonon mode, combined with the Grüneisen Parameter of that mode,

$$\gamma _{𝐪,s}=\frac{d\mathrm{ln}\omega _s(𝐪)}{d\mathrm{ln}V},$$
(20)
determines the Grüneisen Parameter,

$$\gamma =\frac{\underset{𝐪,s}{}\gamma _{𝐪,s}c_{v,s}(𝐪)}{_{𝐪,s}c_{v,s}(𝐪)}.$$
(21)
At high temperatures ($`T>\theta _2`$), where all phonon modes contribute equally, $`\gamma \gamma _0=d\mathrm{ln}\theta _0/d\mathrm{ln}\rho `$. At low temperatures only the acoustic phonon modes contribute.

Figure 15 compares the tight-binding results for the Grüneisen parameter with available data. For densities up to near 13 g/cm3 the results roughly agree with the rule of thumb that $`\gamma \rho =`$constant. Our values are slightly below the experimental values, indicating that the phonon frequencies do not increase with compression as rapidly as they should.

Figure 15 also shows the calculated temperature-dependence of the Grüneisen parameter. At low temperatures ($`T\genfrac{}{}{0pt}{}{<}{}`$40 K) the plot shows a fair amount of structure relative to the high-temperature curve. This can be understood from the phonon dispersion shown in Fig. 5, where the lowest branch is in the $`[\xi \xi \xi ]`$ direction and becomes flat around 3 THz, frequencies that become relevant in their contribution to the specific heat at temperatures around a third of their energy, i.e., around 50 K. This branch is the lowest and hence appears first with increasing temperature, furthermore it appears with a lot of weight as there are eight spatial directions corresponding to these modes.

At low temperatures the phonon contribution to the heat capacity is proportional to $`T^3`$ and vanishes more rapidly than the electronic contribution, which is linear in temperature. Figure 16 shows the calculated coefficient of the electronic contribution to the heat capacity,

$$\gamma _{\text{el}}=\frac{\pi ^2}{3}k_B^2n(E_F),$$
(22)
which is proportional to the density of states at the Fermi energy, $`n(E_F)`$. Compression of the crystal reduces $`n(E_F)`$, i.e., $`\gamma _{\text{el}}`$ decreases monotonically.

## IV SUMMARY

The work presented here is aimed at (1) improving the tight-binding fit of copper specifically for the calculation of thermodynamic properties, and (2) investigating the transferability and range of applicability of the improved model.

For the model to be reliable in calculating thermodynamic properties, it must produce a phonon dispersion in good agreement with experiment. The initial model was fit to first-principles calculations of the total energy at a series of different volumes for the cubic crystal structures. The database of first-principles calculations was extended here to include fcc structures distorted to reflect high-symmetry phonon modes and the elastic constants; fitting to the extended database yields the improved model which indeed delivers phonon frequencies significantly closer to the experimental values.

From the phonon density of states the free energy was calculated, in the quasi-harmonic approximation, as a function of volume and temperature. The temperature-dependence of the minimum of the free energy directly yields the thermal expansion and the linear expansion coefficient, both in good agreement with experiment. The elastic constants are somewhat improved over the initial model, though discrepancies with experiment remain evident.

The quantities in the previous paragraph depend on volumes only in the vicinity of the $`T=0`$ equilibrium volume. The volumes used for the cubic and the distorted fcc structures in the fit extend over a wide range; the equilibrium volume is not treated any differently than other values (down to 9.7 Å3, the smallest volume for which distorted structures were fit). This gives some confidence that the model applies to a range beyond the equilibrium volume and its immediate vicinity.

Within the quasi-harmonic approximation the volume dependence of the phonon frequencies gives a non-zero Grüneisen parameter; the results calculated from the TB model roughly agrees with the empirical $`\gamma \rho =\text{constant}`$. The magnitude is somewhat low, i.e., the compression-induced stiffening of the crystal remains somewhat weaker than is experimentally measured.

The compression at which the model clearly fails can be seen from the Grüneisen parameter as well as the volume dependences of the elastic constants, the electronic contribution to the heat capacity, and the characteristic phonon temperatures. All of these entities vary monotonically with compression until the volume reaches approximately 8 Å3, i.e., a density of roughly 13 g/cm3, at which point unphysical behavior appears.

The unphysical behavior points to the limitations of the model. The Hamiltonian and overlap matrix elements are described by a functional form which can at best approximate the actual behavior within a limited range. For an extended range either the functional form must be modified, e.g. by including higher-order terms in Eq. 4, as has been done in a more recent NRL TB copper potential used in Ref. . The need for modification can also be seen in the electronic band structure, which is degraded by the fitting to distorted fcc structures.

## V Acknowledgment

We thank Jon Boettger, Matthias Graf, David Schiferl, and Duane Wallace for helpful and encouraging discussions. This research is supported by the Department of Energy under contract W-7405-ENG-36. All FLAPW calculations were performed using the Wien97 package. Some of the calculations were performed at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098