text:
# Hyperelliptic Integrable Systems on K3 and Rational Surfaces

## 1 Introduction

Some ten years ago, Beauville  presented a construction of algebraically completely integrable Hamiltonian system (ACIHS) associated with a K3 surface. Beauville’s construction is based on Mukai’s work  on symplectic geometry of moduli spaces of sheaves on K3 surfaces. Beauville considered a special case of Mukai’s moduli spaces, and discovered that this moduli space (with a symplectic structure in the sense of Mukai) is an ACIHS, namely, it is a symplectic variety with a Lagrangian fibration whose fibers are Abelian varieties. An interesting byproduct of his result is the fact that this ACIHS has another realization as a symmetric product of the K3 surface.

Beauville’s construction suggests a new approach to the issue of construction of integrable systems, namely, an approach from algebraic surfaces. Hurtubise  indeed found such a general framework. His result covers many classical integrable systems related to loop algebras , Hitchin’s integrable systems of Higgs pairs , etc., as well as Beauville’s integrable systems, on an equal footing. Furthermore, Gorsky et al.  pointed out a close relation between this approach and Sklyanin’s separation of variables .

As for the case of K3 surfaces, however, very few examples of Beauville’s integrable systems seem to be explicitly known. Almost the only example in the literature is the case of quartics in $`^3`$ that Beauville mentioned in his paper to illustrate his construction. (A more detailed description of this example can be found in the paper of Gorsky et al. )

In this paper, we report a few examples of integrable systems related to Beauville’s construction. Two of them are direct application of Beauville’s construction to an elliptic K3 surface and a K3 surface given by a ramified double covering of $`^2`$. Another example is an integrable system associated with a rational elliptic surface, which is also interesting in comparison with elliptic K3 surfaces. A common characteristic of these examples is that they are described by hyperelliptic curves and Jacobi varieties. This enables us to compare these examples with more classical integrable systems that are related to hyperelliptic curves.

## 2 Beauville’s construction

Following Beauville, we assume that the K3 surface $`X`$ in consideration has an projective embedding or finite morphism $`\varphi :X^g`$ by which a $`g`$-dimensional family of genus $`g`$ curves $`C_u=X\varphi ^1(H_u)`$ are cut out by hyperplanes $`H_u`$, $`u(^g)^{}`$. In order to avoid delicate problems, we restrict the hyperplane parameters $`u`$ to a suitable open subset $`𝒰`$ of $`(^g)^{}`$. We thus have a family $`𝒞𝒰`$ of curves fibered over $`𝒰`$.

The phase space of Beauville’s integrable system is the associated relative Picard scheme $`h:𝒫ic^g𝒰`$ of the degree $`g`$ component $`\mathrm{Pic}^g(C_u)`$ of the Picard group of $`C_u`$. This relative Picard scheme, upon suitably compactified, can be identified with a moduli space of simple sheaves on $`X`$, so that, by Mukai’s work, it has a symplectic structure. As Beauville discovered, $`h`$ is a Lagrangian fibration with respect to Mukai’s symplectic structure. Since each fiber $`h^1(u)=\mathrm{Pic}^g(C_u)`$ is obviously an Abelian variety, the total space $`𝒫ic^g`$ of the relative Picard scheme becomes an ACIHS (see the review of Donagi and Markman  for the notion of ACIHS).

Another realization of this integrable system is given by the $`g`$-fold symmetric product $`X^{(g)}=\mathrm{Sym}^g(X)`$ equipped with a symplectic structure induced by the complex symplectic structure of $`X`$. In order to see the relation with the relative Picard scheme, recall that $`\mathrm{Pic}^g(C_u)`$ consists of the linear equivalence $`[D]`$ of an effective divisor $`D=p_1+\mathrm{}+p_g`$ of degree $`g`$. The unordered $`g`$-tuple $`(p_1,\mathrm{},p_g)`$ of points then becomes a point of $`X^{(g)}`$. Conversely, the $`g`$-tuple $`(p_1,\mathrm{},p_g)`$ of points in a general position uniquely determines a hyperplane $`H_u^g`$ (hence a curve $`C_u=XH_u`$) that passes through the $`g`$ points. One can thus define the map

$`\begin{array}{ccc}\hfill X_{}^{(g)}& & 𝒰\hfill \\ \hfill (p_1,\mathrm{},p_g)& & u\hfill \end{array}`$ (3)
from an open subset $`X_{}^{(g)}`$ to $`𝒰`$. This map can be lifted up to an open embedding

$`\begin{array}{ccc}\hfill X_{}^{(g)}& & 𝒫ic^g\hfill \\ \hfill (p_1,\mathrm{},p_g)& & (u,[p_1+\mathrm{}+p_g])\hfill \end{array}`$ (6)
by adding the linear equivalence class $`[p_1+\mathrm{}+p_g]\mathrm{Pic}^g(C_u)`$ as the second data. Beauville observed that this is a symplectic mapping. Thus, in particular, $`X^{(g)}`$ is fibered by the $`g`$-fold symmetric product $`C_u^{(g)}=\mathrm{Sym}^g(C_u)`$ of the curves $`C_u`$, and these fibers are Lagrangian.

## 3 Hyperelliptic integrable system on elliptic K3 surface

We now apply Beauville’s construction to an elliptic K3 surface $`X^1`$. The Weierstrass model of such an algebraic surface can be written, in affine coordinates, as

$`y^2=z^3+f(x)z+g(x),`$ (7)
where $`x`$ is an affine coordinate of $`^1`$ and $`f(x)`$ and $`g(x)`$ are polynomials of degree $`8`$ and $`12`$, respectively. The 2-form

$`\omega ={\displaystyle \frac{dzdx}{y}},`$ (8)
is nowhere-vanishing and holomorphic, so that $`X`$ becomes a complex symplectic surface.

### 3.1 Curves, divisors and symmetric product

We define a five-parameter family of curves $`C_u`$, $`u=(u_1,\mathrm{},u_5)`$, cut out from $`X`$ by the equation

$`z=P(x)={\displaystyle \underset{k=1}{\overset{5}{}}}u_kx^{5k}.`$ (9)
In other words, the curves $`C_u`$ are hyperelliptic curves of genus $`5`$ defined by the equation

$`y^2=P(x)^3+f(x)P(x)+g(x).`$ (10)
Note that the genus and the number of parameters of curves coincide. Furthermore, these curves are hyperplane sections of the following finite morphism $`\varphi :X^5`$:

$`\begin{array}{ccc}\hfill X& & ^5\hfill \\ \hfill (x,y,z)& & (z:x^4:x^3:x^2:x^1:x:1)\hfill \end{array}`$ (13)
Although we have used the affine coordinates for simplicity, it is not difficult to projectivize this definition. The five parameters $`u_1,\mathrm{},u_5`$ thus turn out to be affine coordinates of $`^5`$.

The mapping $`X_{}^{(5)}𝒰`$ takes a particularly simple form in the present setting. Given an unordered $`5`$-tuple $`(p_1,\mathrm{},p_5)`$ of points of $`X`$, let $`(x_j,y_j,z_j)`$ denote the affine coordinates of the $`j`$-th point $`p_j`$. The curve $`C_u`$ passes through the five points if and only if the equations

$`P(x_j)=z_j`$ (14)
are satisfied for $`j=1,\mathrm{},5`$. The polynomial $`P(x)`$ is uniquely determined by these conditions if the $`x_j`$’s are distinct: $`x_jx_k`$ ($`jk`$). In fact, the Lagrange interpolation formula provides an explicit expression of $`P(x)`$:

$`P(x)={\displaystyle \underset{j=1}{\overset{5}{}}}z_j{\displaystyle \underset{kj}{}}{\displaystyle \frac{xx_k}{x_jx_k}}.`$ (15)
This determines $`u_j`$’s as rational functions of $`(x_1,\mathrm{},x_5,z_1,\mathrm{},z_5)`$, e.g.,

$`u_1={\displaystyle \underset{j=1}{\overset{5}{}}}z_j/{\displaystyle \underset{kj}{}}(x_jx_k).`$ (16)
These functions are to be a commuting set of Hamiltonians.

### 3.2 Abel-Jacobi mapping and symplectic form

The symplectic structure of $`X^{(5)}`$ is defined by the 2-form

$`\mathrm{\Omega }={\displaystyle \underset{j=1}{\overset{5}{}}}{\displaystyle \frac{dz_jdx_j}{y_j}},`$ (17)
where $`(x_j,y_j,z_j)`$ denote the affine coordinates of the $`p_j`$’s. According to Beauville, this symplectic structure is mapped to Mukai’s symplectic structure on $`𝒫ic^5`$.

Following Hurtubise, we now slightly modify Beauville’s construction: We shift the divisors $`p_1+\mathrm{}+p_5`$ to $`p_1+\mathrm{}+p_55p_{\mathrm{}}`$, where $`p_{\mathrm{}}`$ is the point $`(x,y,z)=(\mathrm{},\mathrm{},\mathrm{})`$ of $`C_u`$, and consider the relative Jacobian $`𝒥ac𝒰`$ of the Jacobi varieties $`\mathrm{Jac}(C_u)`$ fibered over $`𝒰`$ rather than the relative Picard. The Abel-Jacobi mapping on $`C_u`$ then induces a mapping $`X_{}^{(5)}𝒥ac`$, which plays the role of Beauville’s mapping.

The Abel-Jacobi mapping sends the symmetric product $`C_u^{(5)}`$ onto the Jacobi variety $`\mathrm{Jac}(C_u)=^5/L_u`$ realized as a complex torus with period lattice $`L_u`$:

$`\begin{array}{ccc}\hfill C_u^{(5)}& & \mathrm{Jac}(C_u)=^5/L_u\hfill \\ \hfill (p_1,\mathrm{},p_5)& & (\psi _1,\mathrm{},\psi _5)\hfill \end{array}`$ (20)
In order to define this mapping, one has to choose a basis of holomorphic differentials on $`C_u`$. We take the standard holomorphic differentials

$`\sigma _k={\displaystyle \frac{x^{5k}dx}{y}}(k=1,\mathrm{},5).`$ (21)
The $`\psi _k`$’s are then given by a sum of integrals of the $`\sigma _k`$’s:

$`\psi _k={\displaystyle \underset{j=1}{\overset{5}{}}}{\displaystyle _{(\mathrm{},\mathrm{})}^{(x_j,y_j)}}{\displaystyle \frac{x^{5k}dx}{y}}`$ (22)
Now $`(u_1,\mathrm{},u_5,\psi _1,\mathrm{},\psi _5)`$ give a local coordinate system on the relative Jacobian $`𝒥ac`$. Remarkably, $`\mathrm{\Omega }`$ turns out to take a canonical form in these coordinates:

$`\mathrm{\Omega }={\displaystyle \underset{k=1}{\overset{5}{}}}du_kd\psi _k.`$ (23)
In other words, the “toroidal” coordinates $`\psi _k`$’ are nothing but the canonical conjugate variables of the hyperplane parameters $`u_k`$.

### 3.3 Proof of (23)

Since

$$z_j=P(x_j)=\underset{k=1}{\overset{5}{}}u_kx_j^{5k},$$

$`\mathrm{\Omega }`$ can be written

$`\mathrm{\Omega }={\displaystyle \underset{j,k=1}{\overset{5}{}}}{\displaystyle \frac{du_kx_j^{5k}dx_j}{y_j}}.`$ (24)
It is convenient choose $`(x_1,\mathrm{},x_5,\psi _1,\mathrm{},\psi _5)`$, rather than $`(u_1,\mathrm{},u_5,\psi _1,\mathrm{},\psi _5)`$, as independent variables here. The total differential of $`d\psi _k`$ can be now expressed as

$$d\psi _k=\underset{j=1}{\overset{5}{}}\frac{\psi _k}{x_j}dx_j+\underset{\mathrm{}=1}{\overset{5}{}}\frac{\psi _k}{u_{\mathrm{}}}du_{\mathrm{}},$$

so that

$`{\displaystyle \underset{k=1}{\overset{5}{}}}du_kd\psi _k={\displaystyle \underset{j,k=1}{\overset{5}{}}}{\displaystyle \frac{\psi _k}{x_j}}du_kdx_j+{\displaystyle \underset{k,\mathrm{}=1}{\overset{5}{}}}{\displaystyle \frac{\psi _k}{u_{\mathrm{}}}}du_kdu_{\mathrm{}}.`$ (25)
Differentiating the explicit definition

$$\psi _k=\underset{j=1}{\overset{5}{}}_{\mathrm{}}^{x_j}\frac{x^{5k}dx}{\left(P(x)^3+f(x)P(x)+g(x)\right)^{1/2}}$$

of the Abel-Jacobi mapping gives

$`{\displaystyle \frac{\psi _k}{x_j}}`$ $`=`$ $`{\displaystyle \frac{x_j^{5k}}{y_j}},`$ (26)
$`{\displaystyle \frac{\psi _k}{u_{\mathrm{}}}}`$ $`=`$ $`{\displaystyle \underset{j=1}{\overset{5}{}}}{\displaystyle _{\mathrm{}}^{x_j}}{\displaystyle \frac{x^{5k}x^5\mathrm{}\left(3P(x)^2+f(x)\right)dx}{2\left(P(x)^3+f(x)P(x)+g(x)\right)^{3/2}}}.`$ (27)
In particular,

$`{\displaystyle \frac{\psi _k}{u_{\mathrm{}}}}={\displaystyle \frac{\psi _{\mathrm{}}}{u_k}},`$ (28)
therefore

$`{\displaystyle \underset{k,\mathrm{}=1}{\overset{5}{}}}{\displaystyle \frac{\psi _k}{u_{\mathrm{}}}}du_kdu_{\mathrm{}}=0.`$ (29)
This implies the equality

$`{\displaystyle \underset{k=1}{\overset{5}{}}}du_kd\psi _k={\displaystyle \underset{j,k=1}{\overset{5}{}}}{\displaystyle \frac{x_j^{5k}}{y_j}}du_kdx_j=\mathrm{\Omega },`$ (30)
thus completing the proof of (23).

### 3.4 Remarks

#### 1.

Essentially the same expression of the symplectic form is stated in Hurtubise’s paper  in a more general form. This expression is thus not specific to the case of hyperelliptic curves, but rather universal.

#### 2.

(23) is a concise expression of the fact that the relative Jacobian $`𝒥ac𝒰`$ is an ACIHS. Firstly, the involutivity (Poisson-commutativity) of the “Hamiltonians” $`u_k`$ is obvious from this expression of $`\mathrm{\Omega }`$. Accordingly, the intersection of the level surfaces of the $`u_k`$’s, which are nothing but the Jacobi varieties $`\mathrm{Jac}(C_u)`$, are Lagrangian subvarieties. Lastly, these Hamiltonians generate a commuting set of linear flows on these complex tori. The Hamiltonian flow of $`u_j`$ is indeed given by

$`\psi _k(t)=\psi _k(0)+\delta _{jk}t.`$ (31)
#### 3.

$`\mathrm{\Omega }`$ should be invariant under the shift

$`\psi _k\psi _k+e_k`$ (32)
by any element $`e=(e_1,\mathrm{},e_5)`$ of the period lattice $`L_u`$. This implies the nontrivial equations

$`{\displaystyle \underset{k=1}{\overset{5}{}}}du_kde_k=0`$ (33)
among the components of the vector $`e`$ (which depends on $`u`$). This is an avatar of the “cubic condition” of Donagi and Markman . Actually, we can verify these conditions directly: Note that $`e_k`$ is given by the period integral

$`e_k={\displaystyle _\gamma }{\displaystyle \frac{x^{5k}dx}{y}}`$ (34)
along a cycle $`\gamma `$. From this expression, one can derive the equality

$`{\displaystyle \frac{e_k}{u_{\mathrm{}}}}={\displaystyle \frac{e_{\mathrm{}}}{u_k}}`$ (35)
in the same way as the calculations in the proof of (23).

#### 4.

A one-parameter family of hyperelliptic curves of genus $`5`$ on the elliptic K3 surface appears in the work of Fock et al. on duality of integrable systems . In fact, their hyperelliptic curves are the subfamily of ours with $`u_1=u_2=u_3=u_4=0`$.

## 4 Hyperelliptic integrable system on non-elliptic K3 surface

We here present an example associated with a double covering of $`^2`$. Let $`(x,z)`$ be affine coordinates of $`^2`$ and consider an algebraic surface $`X`$ defined by the equation

$`y^2=f(x,z),`$ (36)
where $`f(x,z)`$ is a polynomial of degree six whose zero-locus $`f(x,z)=0`$ is a nonsingular sextic curve on $`^2`$. This is a K3 surface that appears in Mukai’s work as an interesting example . A nowhere-vanishing holomorphic 2-form is given by

$`\omega ={\displaystyle \frac{dzdx}{y}}.`$ (37)
A two-parameter family of curves $`C_u`$, $`u=(u_1,u_2)`$, can be cut out from $`X`$ by the equation

$`z=P(x)=u_1x+u_2.`$ (38)
In other words, $`C_u`$ is the hyperelliptic curve of genus $`2`$ defined by the equation

$`y^2=f(x,u_1x+u_2).`$ (39)
Note that these curves are nothing but the hyperplane sections associated with the double covering $`\varphi :X^2`$.

The rest of the construction is fully parallel to the case of the elliptic K3 surface. The integrable system is realized on (an open subset of) the $`2`$-fold symmetric product $`X^{(2)}`$ or on the relative Jacobian $`𝒥ac`$. Let $`(p_1,p_2)`$ denote an unordered pair of points of $`X`$, $`(x_j,y_j,z_j)`$ the coordinates of $`p_j`$, and $`\mathrm{\Omega }`$ the symplectic form induced on $`X^{(2)}`$:

$`\mathrm{\Omega }={\displaystyle \underset{j=1,2}{}}{\displaystyle \frac{dz_jdx_j}{y_j}}.`$ (40)
The Abel-Jacobi mapping

$`\begin{array}{ccc}\hfill C_u^{(2)}& & \mathrm{Jac}(C_u)=^2/L_u\hfill \\ \hfill (p_1,p_2)& & (\psi _1,\psi _2)\hfill \end{array}`$ (43)
is now given by

$`\psi _1={\displaystyle \underset{j=1,2}{}}{\displaystyle _{(\mathrm{},\mathrm{})}^{(x_j,y_j)}}{\displaystyle \frac{xdx}{y}},\psi _2={\displaystyle \underset{j=1,2}{}}{\displaystyle _{(\mathrm{},\mathrm{})}^{(x_j,y_j)}}{\displaystyle \frac{dx}{y}}.`$ (44)
The symplectic form $`\mathrm{\Omega }`$ takes the canonical form

$`\mathrm{\Omega }={\displaystyle \underset{k=1,2}{}}du_kd\psi _k`$ (45)
in the coordinates $`(u_1,u_2,\psi _1,\psi _2)`$.

## 5 Hyperelliptic integrable system on rational surface

### 5.1 Integrable system associated with rational elliptic surface

If the polynomials $`f(x)`$ and $`g(x)`$ in the Weierstrass model of elliptic K3 surfaces are replaced by polynomials of degree $`4`$ and $`6`$, resplectively, the outcome is an rational elliptic surface. The 2-form

$`\omega ={\displaystyle \frac{dzdx}{y}}`$ (46)
is nowhere-vanishing and holomorphic in the affine part of the surface, having poles at the compactification divisor at infinity. Let us examine this case in comparison with the case of K3 surfaces.

Simple power counting suggests that the most natural choice of $`P(x)`$ in this case will be a polynomial of degree $`2`$:

$`P(x)=u_0x^2+u_1x+u_2.`$ (47)
The equation $`z=P(x)`$ then defines a hyperelliptic curve of genus $`2`$ in $`X`$:

$`y^2=P(x)^3+f(x)P(x)+g(x).`$ (48)
Ganor  considered the same family of curves in the context of string theory. The foregoing construction, however, does not work literally, firstly because the number of parameters (three) and the genus of the curves (two) do not match.

A correct prescription is to fix $`u_0`$ to a constant $`c`$, and to take the two-parameter family of curves $`C_u`$, $`u=(u_1,u_2)`$, cut out from $`X`$ by the equation

$`z=P(x)=cx^2+u_1x+u_2.`$ (49)
In other words, $`u_0=c`$ should be treated as a Casimir function on a Poisson manifold.

Now we consider $`X^{(2)}`$ (rather than $`X^{(3)}`$ that Ganor argued) with the symplectic form

$`\mathrm{\Omega }={\displaystyle \underset{k=1,2}{}}{\displaystyle \frac{dz_jdx_j}{y_j}},`$ (50)
where $`(x_j,y_j,z_j)`$ are the coordinates of an unordered pair $`(p_1,p_2)`$ of points on $`X`$. The Abel-Jacobi mapping $`(p_1,p_2)(\psi _1,\psi _2)`$ is given by

$`\psi _1={\displaystyle \underset{j=1,2}{}}{\displaystyle _{(\mathrm{},\mathrm{})}^{(x_j,y_j)}}{\displaystyle \frac{xdx}{y}},\psi _2={\displaystyle \underset{j=1,2}{}}{\displaystyle _{(\mathrm{},\mathrm{})}^{(x_j,y_j)}}{\displaystyle \frac{dx}{y}},`$ (51)
and sends $`\mathrm{\Omega }`$ the canonical form:

$`\mathrm{\Omega }={\displaystyle \underset{k=1,2}{}}du_kd\psi _k.`$ (52)
Thus the situation for the rational elliptic surfaces is slightly different from the case of K3 surfaces. Namely, we have to fix the leading coefficient of $`P(x)`$ to a constant $`c`$ so as to match the number of parameters and the genus of curves.

### 5.2 Neumann system and rational surface

It is instructive to compare the foregoing construction with the Neumann system. The Neumann system is a classical integrable system that was solved in the 19th century  using the theory of hyperelliptic integrals developed in those days, and revived by Moser  from a modern point of view.

The hyperelliptic curves of the Neumann system take the special form

$`y^2=P(x)Q(x),`$ (53)
where $`P(x)`$ and $`Q(x)`$ are polynomials of degree $`N`$ and $`N+1`$ of the form

$`P(x)=x^N+{\displaystyle \underset{k=1}{\overset{N}{}}}u_kx^{Nk},Q(x)={\displaystyle \underset{n=1}{\overset{N+1}{}}}(xc_n),`$ (54)
and $`N`$, the genus of the curves, is equal to the degrees of freedom of the Neumann system. The coefficients $`u_k`$ of $`P(x)`$ are integrals of motion in involution. The $`c_n`$’s are non-dynamical structural constants (i.e., Casimir functions) of the system.

One will soon notice that a rational surface is hidden behind these hyperelliptic curves. This rational surface $`X`$ is defined by the equation

$`y^2=zQ(x).`$ (55)
The 2-form

$`\omega ={\displaystyle \frac{dzdx}{y}}`$ (56)
is nowhere-vanishing and holomorphic in the affine part of $`X`$. The curves $`C_u`$ cut out from $`X`$ by the equation $`z=P(x)`$ are exactly the aforementioned hyperelliptic curves. It is not difficult to confirm that the construction based on the symmetric product $`X^{(N)}`$ or the relative Jacobian $`𝒥ac`$ does reproduce the classical method for solving the Neumann system.

It is interesting to note that the leading coefficient of $`P(x)`$, which is now set to $`1`$, is also a kind of structural constant. This constant is related to the radius of an $`N`$-dimensional sphere on which the Neumann system is realized as a mechanical system. If the radius takes a different value, the leading coefficient of $`P(x)`$ accordingly varies as

$$P(x)=cx^N+\underset{k=1}{\overset{N}{}}u_kx^{Nk}.$$

This is the same situation that we have encountered in the case of rational elliptic surfaces.

## 6 Concluding Remarks

We have shown several integrable systems that illustrate Beauville’s construction. The construction is considerably simplified by the fact that the curves arising therein are hyperelliptic. Nevertheless, the most essential part, i.e., the two expressions of $`\mathrm{\Omega }`$ on the symmetric product $`X^{(g)}`$ and on the relative Jacobian $`𝒥ac`$, are rather universal and can be extended to a non-hyperelliptic case. Let us stress that the coordinates $`(x_j,y_j,z_j)`$ of the unordered $`g`$-tuple $`(p_1,\mathrm{},p_g)X^{(g)}`$ correspond to “separated variables” in Sklyanin’s separation of variables.

We have pointed out a delicate difference between the cases of K3 and of rational surfaces. Namely, some of the hyperplane parameters $`u_k`$ in the latter case have to be interpreted as a Casimir rather than an integral of motion. This issue deserves to be pursued further.

The present approach from algebraic surfaces can be a promising alternative to the conventional method based on a Lax representation. We have seen such a possibility in the case of the Neumann system. Let us mention another example, which is related to Witten’s eleven-dimensional interpretation  of the Seiberg-Witten curves for four-dimensional $`N=2`$ supersymmetric QCD.

The Seiberg-Witten curves for $`\mathrm{SU}(N_c)`$ gauge theory coupled to $`N_f`$ ($`N_f2N_c`$) hypermultiplets are hyperelliptic curves of genus $`g=N_c1`$ defined by the equation

$`y^2P(x)y+Q(x)=0.`$ (57)
Here $`P(x)`$ is a polynomial of the form

$`P(x)=x^{N_c}+{\displaystyle \underset{k=2}{\overset{N_c}{}}}u_kx^{N_ck},`$ (58)
and $`Q(x)`$ a polynomial of degree $`N_f`$ of the form

$`Q(x)=\mathrm{\Lambda }^{2N_cN_f}{\displaystyle \underset{\mathrm{}=1}{\overset{N_f}{}}}(x+m_{\mathrm{}}).`$ (59)
The $`u_k`$’s are the “moduli” of the Seiberg-Witten curves, which we interpret as (part of) the parameters of hyperplane sections. If $`N_f=0`$, $`Q(x)`$ reduces to a constant, and the curve is exactly the the spectral curve of the $`N_c`$-periodic Toda chain.

In Witten’s interpretation, these curves are embedded into the affine rational surface $`X`$ defined by the equation

$`yz=Q(x),`$ (60)
A nowhere-vanishing holomorphic 2-form on $`X`$ is given by

$`\omega ={\displaystyle \frac{dydx}{y}}={\displaystyle \frac{dxdz}{z}}.`$ (61)
In fact, this complex surface is a realization of the multi-Taub-NUT space (gravitational instanton) in one of its $`2`$-parameter family of complex structures . As noted by Nakatsu et al.  and de Boer et al. , the Seiberg-Witten curves are cut out from $`X`$ by the equation

$`y+z=P(x).`$ (62)
One can now repeat the foregoing construction with a slightest modification to obtain an integrable system realized on $`X^{(g)}`$ or on $`𝒥ac`$. Presumably, this integrable system will be equivalent to the spin chain that Gorsky et al.  proposed as an integrable system for these Seiberg-Witten curves. Our construction, however, is more direct and seemingly more natural in the context of brane theory .

### Acknowledgements

I am grateful to Yasuhiko Yamada for useful comments. This work is partly supported by the Grant-in-Aid for Scientific Research (No. 12640169), the Ministry of Education, Science and Culture.


format: text

text:
# Electron Cotunneling in a Semiconductor Quantum Dot

## Abstract

We report transport measurements on a semiconductor quantum dot with a small number of confined electrons. In the Coulomb blockade regime, conduction is dominated by cotunneling processes. These can be either elastic or inelastic, depending on whether they leave the dot in its ground state or drive it into an excited state, respectively. We are able to discriminate between these two contributions and show that inelastic events can occur only if the applied bias exceeds the lowest excitation energy. Implications to energy-level spectroscopy are discussed.

PACS numbers: 73.23.Hk, 73.40.Gk, 73.61.Ey

Quantum-dot devices consist of a small electronic island connected by tunnel barriers to source and drain electrodes . Due to on-site Coulomb repulsion, the addition of an electron to the island implies an energy change $`U=e^2/C`$, where $`C`$ is the total capacitance of the island. Hence the number of confined electrons is a well-defined integer, $`N`$, that can be controlled by varying the voltage on a nearby gate electrode. Transport of electrons through the dot is allowed only at the transition points where the $`N`$\- and ($`N+1`$)-states are both energetically accessible. Otherwise, $`N`$ is constant and current is strongly suppressed. This is known as Coulomb blockade . At low temperature, however, higher-order tunneling events can become dominant. These are commonly known as $`cotunneling`$ events since they involve the simultaneous tunneling of two or more electrons . Electron cotunneling has received considerable attention over the last decade. Initially it was recognized as a limitation to the accuracy of single-electron devices. More recently, it has acquired a broader relevance, especially since an increasing activity has been focused on quantum dots with a deliberately enhanced tunnel coupling to their leads. These systems allow the investigation of high-order transport processes and many-body phenomena, such as the Kondo effect . The latter can be regarded as the result of a coherent superposition of higher-order cotunneling events. Here, we will focus on the lowest order contribution to cotunneling.

Previous experiments were performed with metallic islands  or large semiconductor dots , where the energy spectrum is essentially continuous and many levels contribute to cotunneling. Here, we study cotunneling through a small quantum dot where the energy levels are well separated, and where the absolute value of $`N`$ is precisely known. A cotunneling event is called inelastic when it leaves the dot in an excited state. Otherwise it is classified as elastic. We identify two regimes: one consisting of elastic processes only, and one including both elastic and inelastic contributions. We note that the transition between these regimes can be sharper than the characteristic life-time broadening of the dot states. In such a case, the onset of inelastic cotunneling can be exploited to measure the energy spectrum of a quantum dot with improved resolution.

The stability diagram of a generic quantum dot can be obtained by plotting the differential conductance (d$`I`$/d$`V_{sd}`$) as a function of bias, $`V_{sd}`$, and gate voltage, $`V_g`$. Coulomb blockade occurs within the diamond-shaped regions in Fig. 1a. The diamond size is proportional to the addition energy, defined as $`E_{add}(N)\mu _{dot}(N+1)\mu _{dot}(N)`$, where $`\mu _{dot}(N)`$ is the electrochemical potential of an $`N`$-electron dot. Inside the $`N`$-electron diamond, $`\mu _{dot}(N)<\mu _L,\mu _R<\mu _{dot}(N+1)`$, with $`\mu _L`$, $`\mu _R`$ the Fermi energies of the leads. The diamond edges correspond to level alignment: $`\mu _{dot}(N)`$ = $`\mu _L`$ or $`\mu _R`$ (see angled solid lines). This alignment determines the onset for first-order tunneling via the ground state of the dot, leading to a peak in d$`I`$/d$`V_{sd}(V_{sd})`$. The onset for first-order tunneling via the first excited state occurs at a somewhat higher bias (see dot-dashed lines in Fig. 1a, and the corresponding diagrams in Fig. 1b and 1e). These first-order processes have been exploited as a spectroscopic tool on the discrete energy spectrum of dots .

Here, we are interested in second-order tunneling of charge, which becomes more apparent when the tunnel coupling between the dot and the leads is enhanced. We neglect contributions from spin that could give rise to the Kondo effect. Elastic cotunneling is the dominant off-resonance process at low bias. It gives rise to current inside the Coulomb diamond (light-grey region in Fig. 1a). The corresponding two-electron process (Fig. 1c) transfers one electron from the left to the right lead, thereby leaving the dot in the ground state.

For $`e|V_{sd}|\mathrm{\Delta }(N)`$, where $`\mathrm{\Delta }(N)`$ is the lowest on-site excitation energy for a constant $`N`$ , similar two-electron processes can occur which drive the dot into an excited state. For instance, an electron can leave the dot from the ground state to the lowest Fermi sea, while another electron from the highest Fermi sea tunnels into the excited state (see Fig. 1d). Although this type of process is called inelastic , the total electron energy is conserved. The on-site excitation is created at the expense of the energy drop $`eV_{sd}`$. To first approximation, the onset of inelastic cotunneling yields a step in d$`I`$/d$`V_{sd}(V_{sd})`$ . This step occurs when $`e|V_{sd}|=\mathrm{\Delta }(N)`$, which is not or only weakly affected by $`V_g`$ (see also Fig. 1c in Ref. ). As a result, inelastic cotunneling turns on along the vertical (dotted) lines in Fig. 1a. At the edge of the Coulomb diamond the condition for the onset of inelastic cotunneling connects to that for the onset of first-order tunneling via an excited state (dot-dashed lines).

Our device has the external shape of a 0.5-$`\mu `$m-high pillar with a $`0.6\times 0.45`$ $`\mu `$m2 rectangular base (inset to Fig. 2). It is fabricated from an undoped AlGaAs(7 nm)/InGaAs(12 nm)/AlGaAs(7 nm) double barrier heterostructure, sandwiched between n-doped GaAs source and drain electrodes . The quantum dot is formed within the InGaAs layer. The lateral confinement potential is close to that of an ellipse . Its strength is tuned by a negative voltage, $`V_g`$, applied to a metal gate surrounding the pillar. A dc bias voltage, $`V_{sd}`$, applied between source and drain, drives current vertically through the pillar. In addition, we apply a small bias modulation with rms amplitude $`V_{ac}=3`$ $`\mu `$V at 17.7 Hz for lock-in detection. Measurements are carried out in a dilution refrigerator with a base temperature of 15 mK. We find an effective electron temperature $`T_e=25\pm 5`$ mK, due to residual electrical noise.

Figure 2 shows d$`I`$/d$`V_{sd}`$ in grey-scale versus ($`V_{sd}`$,$`V_g`$) at 15 mK. Diamond-shaped regions of low conductivity (light grey) identify the Coulomb blockade regimes for $`N`$ = 1 to 4. The diamonds are delimited by dark narrow lines (d$`I`$/d$`V_{sd}e^2/h`$) corresponding to the onset of first-order tunneling. For $`N=1`$, as well as for $`N=3`$, sub-gap transport is dominated by elastic cotunneling with no evidence for inelastic cotunneling. The differential conductance is uniformly low inside the Coulomb diamond. (Slight modulations are seen due to a weak charging effect in the GaAs pillar above the dot .) This is different for $`N=2`$, where the onset of inelastic cotunneling is clearly observed. As argued before, this onset follows (dotted) lines, nearly parallel to the $`V_g`$ axis . At the diamond edges they connect to (dot-dashed) lines where first-order tunneling via an excited state sets in. Similar considerations apply to $`N=4`$.

The different behavior observed for $`N=even`$ and $`N=odd`$ stems from the fact that inelastic cotunneling occurs only if $`E_{add}(N)>\mathrm{\Delta }(N)`$, as apparent from Fig. 1a. In the case of non-interacting electrons, $`E_{add}(N)=U(N)`$ for $`N=odd`$ and $`E_{add}(N)=U(N)+\mathrm{\Delta }(N)`$ for $`N=even`$, where $`U(N)`$ is the charging energy for $`N`$ electrons; energy levels are spin degenerate and consecutively filled with pairs of electrons. This is a reasonable picture if the level spacing exceeds the exchange interaction energy . In our small quantum dot the first three levels are indeed widely spaced as already discussed in Ref. . For low $`N`$, $`\mathrm{\Delta }(N)`$ exceeds not only the exchange energy but also $`U(N)`$. This implies that for $`N=odd`$, $`\mathrm{\Delta }(N)`$ lies outside the Coulomb diamond (i.e. $`\mathrm{\Delta }(N)>E_{add}(N)`$) and thus inelastic cotunneling is not observed. (Note that cotunneling between spin-degenerate states is an elastic process as initial and final state have the same energy.) For $`N=even`$, $`\mathrm{\Delta }(N)`$ is always smaller than $`E_{add}(N)`$ and inelastic cotunneling can be observed.

We now discuss the difference in life-time broadening between first- and higher-order tunneling. At the onset of first-order tunneling a certain level is aligned to one of the Fermi energies. In this case, an electron can escape from the dot, which leads to a finite life-time broadening of the observed resonance by an amount $`\mathrm{}\mathrm{\Gamma }`$. Here, $`\mathrm{\Gamma }=\mathrm{\Gamma }_L+\mathrm{\Gamma }_R`$, where $`\mathrm{\Gamma }_L`$ and $`\mathrm{\Gamma }_R`$ are the tunneling rates through the left and the right barrier, respectively. (Note that these rates are independent of $`V_{sd}`$, since our bias window ($``$meV) is much smaller than the height of the AlGaAs tunnel barriers ($``$50 meV).

The onset of inelastic cotunneling is also characterized by a width. In the zero-temperature limit, this is determined by the life-time broadening of the excited state. Two types of situations can occur. First, the excited state can be between $`\mu _L`$ and $`\mu _R`$ (see right inset to Fig. 3) so that inelastic cotunneling can be followed by first-order tunneling. Such a decay event leads to a life-time broadening of at least $`\mathrm{}\mathrm{\Gamma }_R\mathrm{}\mathrm{\Gamma }/2`$. Second, the ground and excited state are both well below $`\mu _L`$ and $`\mu _R`$, implying that only higher-order tunneling is allowed (see right inset to Fig. 4). Now, decay from the excited state can only rely on cotunneling. For this higher-order perturbation, the corresponding rate, $`\mathrm{\Gamma }_{co}`$, is much smaller than $`\mathrm{\Gamma }`$, leading to a reduced life-time broadening. To illustrate these arguments, we select different d$`I`$/d$`V_{sd}vsV_{sd}`$ traces and analyse their shape in detail.

Figure 3 shows two traces for $`N=2`$, taken at 15 mK for gate voltages at the horizontal lines in the left inset. The dashed trace has several peaks. The two inner ones, at $`|V_{sd}|1.1`$ mV, correspond to first-order tunneling of the 3rd electron via the 3-electron ground state; i.e. $`\mu _{dot}(3)`$ = $`\mu _L`$ or $`\mu _R`$. The right (left) peak has a full width at half maximum (FWHM) of $``$200 ($``$400) $`\mu `$V. This is somewhat larger than the width, $`\mathrm{}\mathrm{\Gamma }/e150`$ $`\mu `$V, measured in the zero-bias limit. Indeed a finite bias allows non-energy-conserving tunneling events leading to additional broadening. The most likely source for energy relaxation is acoustic-phonon emission . The following pair of peaks, at $`|V_{sd}|2`$ mV, corresponds to the onset of first-order tunneling via the first excited state for $`N=2`$ (see Fig. 1b). Because of the larger bias voltage, these peaks are visibly broader than the inner ones. Additional peak structures occur near the edges of the bias window. The origin of these peaks can not be precisely identified.

The solid trace contains structure from both first- and second-order tunneling. The peaks labeled by open squares arise from first-order tunneling at the edges of the Coulomb diamond (see left inset). Steps, labeled by open circles, identify the onset of inelastic cotunneling and correspond to the open circles in the left inset. Their different heights are probably due to a left-right asymmetry in the tunnel coupling to the leads. Their $`V_{sd}`$-position, which is symmetric around zero bias, provides a direct measure of $`\mathrm{\Delta }(2)`$. The width of these steps is $``$150 $`\mu `$V . Since $`\mathrm{\Delta }(2)U(2)`$, the first excited state lies unavoidably within the bias window when $`|V_{sd}|=\mathrm{\Delta }(2)/e`$ and hence is allowed to decay into the lowest-energy lead (see the right inset to Fig. 3). As argued above, this situation leads to a step-width exceeding $`\mathrm{}\mathrm{\Gamma }_R/e`$, consistent with our finding. Another structure occurs at $`V_{sd}2.6`$ mV and is probably due to the onset of inelastic cotunneling via the second excited state for $`N=2`$. The corresponding line in the stability diagram is hardly visible due to its vicinity to the diamond edge.

To study inelastic cotunneling when both ground and excited state lie well below the Fermi energies of the leads (Fig. 1d) we need $`\mathrm{\Delta }(N)E_{add}(N)`$. To this aim we move to $`N=6`$, since $`\mathrm{\Delta }(6)`$ can be effectively tuned by a magnetic field applied along the vertical axis. We tune the field to 0.35 T, such that $`\mathrm{\Delta }(6)0.1`$ meV, i.e. several times smaller than $`E_{add}(6)`$. From a previous study we know that the ground state is a spin singlet, and the first excited state is a spin triplet . The d$`I`$/d$`V_{sd}vsV_{sd}`$ traces shown in Fig. 4 are taken at two different temperatures, but for the same $`V_g`$, at the horizontal line in the left inset. The solid trace (15 mK) shows a broad minimum around $`V_{sd}=0`$, where transport is dominated by elastic cotunneling via the ground state (see also the light-grey region in the left inset). The differential conductance increases rapidly at the onset of inelastic cotunneling with a step-width of $``$20 $`\mu `$V, i.e. much smaller $`\mathrm{}\mathrm{\Gamma }/e`$. This reduced width stems from the fact that the excited state can not decay directly into the lower energy lead (see right inset). The corresponding life-time broadening, $`\mathrm{}\mathrm{\Gamma }_{co}`$, can be estimated from the cotunneling current, $`I_{co}`$, at $`V_{sd}=\mathrm{\Delta }(6)/e`$. We find $`\mathrm{}\mathrm{\Gamma }_{co}=\mathrm{}I_{co}/e(\mathrm{}/e)_0^{\mathrm{\Delta }(6)/e}`$ d$`I`$/d$`V_{sd}(V_{sd})`$ d$`V_{sd}10`$ $`\mu `$eV, consistent with the observed step-width. At $`T_e=25`$ mK the thermal broadening of the Fermi distribution functions leads to a step-width of $`5.44k_BT_e/e12`$ $`\mu `$eV . Hence life-time broadening has been reduced here to the thermal limit.

The cotunneling onset in Fig. 4 shows a peak structure at low temperature (solid trace) in addition to the expected step structure. This is likely due to Kondo correlations, as discussed in Ref. . On increasing temperature to 200 mK these Kondo correlations are suppressed such that only lowest order cotunneling contributes. This recovers the step structure (dashed trace).

We thank Yu. V. Nazarov, M. R. Wegewijs, M. Eto, K. Maijala, and J. E. Mooij for discussions. We acknowledge financial support from the Specially Promoted Research, Grant-in-Aid for Scientific Research, from the Ministry of Education, Science and Culture in Japan, from the Dutch Organisation for Fundamental Research on Matter (FOM), from the NEDO joint research program (NTDP-98), and from the EU via a TMR network.


text:
# Electrical transport properties of bulk NicFe1-c alloys and related spin-valve systems

## I Introduction

Permalloy is perhaps the most commonly used magnetic material in a variety of devices. This arises primarily from its magnetic properties which combine a reasonably high magnetic moment with a low coercivity. In many of its applications one passes a current through the material, so that its electrical transport properties are also important. There have been several calculations of the transport properties of permalloy, but none has been completely successful. Therefore, one of the outstanding challenges has been an ab initio calculation of the electrical resistivity and anisotropic magnetoresistance across the entire NicFe1-c series. Here we present the first part which addresses the resistivity of bulk permalloy and also the giant magnetoresistance of spin-valve structures containing layers of permalloy. As the lattice structure changes from bcc to fcc at $`c35\%`$ we present results for both lattice structures. We are able to reproduce the overall variation of the resistivity across the series; in particular we find the signature changes in resistivity about $`c35\%`$ as the lattice structure changes from bcc to fcc. Remarkably we find the magnitude of the resistivity is within a factor of two of the experimental values; noteworthy our calculated values are larger, so that when one relaxes the approximations we made one reduces the resistivities and approaches the measured ones.

Here we outline our ab initio calculation and present our results for the resistivity across the entire NicFe1-c series. Then we conclude with results on prototypical spin-valve structures containing permalloy layers, e.g., NicFe1-c(100Å) / Co(6Å) / Cu(9Å) / Co(6Å) / NicFe1-c(5Å). Here we find current in the plane of the layers (CIP) magnetoresistance ratios that compare favorably with those observed in these structures.

## II Method of calculation

In terms of the Kubo-Greenwood approach for disordered layered systems with growth direction along the surface normal ($`z`$-axis) the in-plane conductivity is given by

$$\begin{array}{c}\sigma _{\mu \mu }(n;𝐜;𝒞)=\underset{i,j=1}{\overset{n}{}}\sigma _{\mu \mu }^{ij}(𝐜;𝒞),\mu \{x,y\},\hfill \\ \\ \sigma _{\mu \mu }^{ij}=\frac{\mathrm{}}{\pi N_0\mathrm{\Omega }_{\mathrm{at}}}\mathrm{Tr}J_\mu ^i\mathrm{Im}G^+(ϵ_F)J_\mu ^j\mathrm{Im}G^+(ϵ_F),\hfill \end{array}$$
(1)
where $`n`$ denotes the number of layers considered, $`𝐜=\{c_1,c_2,\mathrm{},c_n\}`$ is a set containing the layer-wise compositions, e.g., the layer-wise concentrations in an inhomogeneously disordered binary alloy system A$`_{c_i}`$B$`_{1c_i}`$ (in a homogeneously disordered alloy $`c_i=c,i`$; whereas in a spin-valve system Am / Bn / Am the $`c_i=1,i\mathrm{A}_m`$, and $`c_i=0,i\mathrm{B}_n`$). In principle the magnetization in each layer of a multilayered structure can be different, and $`𝒞=\{𝐌_1,𝐌_2,\mathrm{},𝐌_n\}`$ denotes a particular magnetic configuration, where the $`𝐌_i`$ refer to the orientations of the magnetizations in the individual layers. $`\sigma _{\mu \mu }^{ij}`$ is the conductivity that describes the current in layer $`i`$ caused by an electric field in layer $`j`$, $`N_0`$ the number of atoms per plane of atoms, $`\mathrm{\Omega }_{\mathrm{at}}`$ is the atomic volume, $`\mathrm{}`$ denotes an average over configurations, $`J_\mu ^i`$ is the $`\mu `$-th component of the current operator referenced to the $`i`$-th plane, and $`G^+`$ is the electron propagator (one-particle Green’s function) from plane $`i`$ to $`j`$ at the Fermi energy $`ϵ_F`$. The resistivity corresponding to the conductivity given by Eq. (1) is then defined by

$$\rho _{\mu \mu }(n;𝐜;𝒞)=1/\sigma _{\mu \mu }(n;𝐜;𝒞),\mu \{x,y\}.$$
(2)
By assuming two different magnetic configurations, $`𝒞`$ and $`𝒞^{}`$, the relative change in the resistivities between them is given by the ratio

$$R_{\mu \mu }(n;𝐜;𝒞;𝒞^{})=\frac{\rho _{\mu \mu }(n;𝐜;𝒞^{})\rho _{\mu \mu }(n;𝐜;𝒞)}{\rho _{\mu \mu }(n;𝐜;𝒞^{})}.$$
(3)
For example in spin-valve structures $`𝒞^{}`$ and $`𝒞`$ refer to “antiparallel” and “parallel” alignments of the magnetizations of the magnetic slabs, and this ratio is referred to as the giant magnetoresistance (GMR) ratio. The ratio in Eq. (3) also serves as a definition for the anisotropic magnetoresistance of homogeneously disordered bulk alloys, $`c_i=c,i`$, with uniform magnetization, $`𝐌_i=𝐌,i`$, with $`𝒞^{}`$ and $`𝒞`$ referring to configurations with the magnetization pointing uniformally parallel and perpendicular to the current axis (however, in this case the ratio is usually taken with respect to an averaged resistivity).

In the present paper all calculations are based on self-consistent effective potentials and effective exchange fields as obtained previously  by using the fully relativistic, spin-polarized, screened Korringa-Kohn-Rostoker method together with the coherent-potential approximation for layered systems. For permalloy a fully relativistic spin-polarized calculation of the resistivity is essential in order to properly account for the spin-orbit coupling of the electron states as well as for the scattering by the disordered constituents of these binary alloys.  The surface-Brillouin-zone integrals needed in the evaluation of the electrical conductivity within the Kubo-Greenwood approach  were obtained by considering 1830 $`𝐤_{}`$-points in the irreducible wedge of the surface Brillouin zone. All scattering channels up to and including $`\mathrm{}_{\mathrm{max}}=2`$ were taken into account. Here we have chosen the magnetization to be uniformly perpendicular to the planes of atoms, i.e., the magnetization points along the $`z`$-axis. With this choice the $`xx`$\- and $`yy`$-components in the above equations are identical, hence we have taken the direction of the current (electric field) to be along the $`x`$-axis so as to simulate the CIP transport geometry (with the magnetization pointing uniformally perpendicular to the current axis).

Consider the following systems:

| (a) |  | $`\text{Ni}_c\text{Fe}_{1c}(100)`$ / | $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_n`$ / | $`\text{Ni}_c\text{Fe}_{1c}(100)`$ |
| --- | --- | --- | --- | --- |
| (b) |  | $`\text{Ni}_c\text{Fe}_{1c}(100)`$ / | $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_n`$ / | Vacuum |
| (c) |  | Vacuum / | $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_n`$ / | Vacuum . |

(4)
All three cases refer to the in-plane conductivity for a film with $`n`$ monolayers of permalloy; what distinguishes them are the boundary conditions.  In case (a) we have outgoing boundary conditions on both sides as there is nothing that separates the film from the substrate it is deposited on nor from the overlayer that caps it, i.e., some electrons will leak out into the semi-infinite leads and their momentum information will be lost hence producing an additional contribution to the resistivity. In (b) we have a free surface on one side which makes it reflecting (the workfunction is finite so that the potential barrier is not infinitely high and there will be a slight leaking out of electrons into the vacuum that has to be taken into account). In the third case (c) we have a perfectly flat free standing film with reflecting boundary conditions on both sides. Only in (c) the calculated in-plane conductivities/resistivities for a perfectly flat film are principally independent of the film thickness $`n`$, at least for $`n>12`$; roughness would act like an outgoing boundary condition. However, for all three cases one obtains the so-called bulk residual resistivity $`\rho ^0(𝐜;𝒞)`$ of substitutionally disordered binary alloys with respect to a particular magnetic configuration by taking the infinite-thickness limit

$$\rho ^0(𝐜;𝒞)=\underset{n\mathrm{}}{lim}\rho _{\mu \mu }(n;𝐜;𝒞).$$
(5)
This infinite-thickness limit ($`n\mathrm{}`$) can be obtained easily since $`n\rho _{\mu \mu }(n;𝐜;𝒞)`$ becomes linear in $`n`$, if $`n`$ is large enough; its slope simply refers to the extrapolated resistivity (see also the discussion of Fig. 2 in the next Section).

In addition, for computational purposes a finite imaginary part $`\delta `$ to the Fermi energy has to be used in the calculation of the conductivity $`\sigma _{\mu \mu }(n;𝐜;𝒞;ϵ_F+i\delta )`$.  To designate resistivities calculated with finite $`\delta `$ we adopt the notation

$$\rho _{\mu \mu }(n;𝐜;𝒞;\delta )=1/\sigma _{\mu \mu }(n;𝐜;𝒞;ϵ_F+i\delta ),\mu \{x,y\}.$$
(6)
The actual resisitivity is arrived at by numerically taking the limit as the imaginary part goes to zero

$$\rho _{\mu \mu }(n;𝐜;𝒞)=\underset{\delta 0}{lim}\rho _{\mu \mu }(n;𝐜;𝒞;\delta ).$$
(7)
If we are interested in the bulk residual resistivity the ($`\delta 0`$) extrapolation is simplified by the fact that the slope of $`n\rho _{\mu \mu }(n;𝐜;𝒞;\delta )`$ (determined at large enough $`n`$) is itself linear in $`\delta `$.

As a check on the above extrapolation procedure we used to arrive at bulk values of the resistivity from those on films of finite thickness (Eq. (5)) we have calculated in-plane resistivities (with the magnetization pointing uniformly along the $`z`$-axis) for 30 monolayers of permalloy with a homogeneous Ni concentration of $`c=0.85`$ for a fixed $`\delta =2`$ mRyd with different boundary conditions: For the reflecting boundary conditions (case (c)) we find a resistivity of $`\rho _{xx}^{\mathrm{ref}}(30;0.85;\widehat{𝐳};2)=20.6`$ $`\mu \mathrm{\Omega }`$cm. A similar calculation with outgoing boundary conditions (case (a)) yields $`\rho _{xx}^{\mathrm{out}}(30;0.85;\widehat{𝐳};2)=31.5`$ $`\mu \mathrm{\Omega }`$cm which by taking $`n\mathrm{}`$ yields $`\rho _{xx}^{\mathrm{out}}(\mathrm{};0.85;\widehat{𝐳};2)=20.6`$ $`\mu \mathrm{\Omega }`$cm; precisely the same value as found for the reflecting boundary conditions (case (c)). This unequivocally demonstrates that the procedure outlined above is able to remove the resistivity that arises from the boundary conditions on finite structures.  In Fig. 1 we show as an illustration the variation of the layer-wise conductivities

$$\sigma _{xx}^i(n;𝐜;𝒞;ϵ_F+i\delta )=\underset{j=1}{\overset{n}{}}\sigma _{xx}^{ij}(n;𝐜;𝒞;ϵ_F+i\delta ),$$
(8)
across the film with reflecting and outgoing boundary conditions. For the outgoing boundary conditions where some electrons leak into the leads all layer-wise conductivity contributions are substantially smaller than for the reflecting boundary conditions, the difference becoming progressively larger when we get closer to the outer edges of the film. When we apply reflecting boundary conditions the layer-wise conductivities are practically constant with small oscillations from the boundary towards the center of the film resembling the well-known Friedel oscillations. When we take the result for $`\delta =2`$ mRyd corresponding to outgoing boundary conditions to $`n\mathrm{}`$ (Eq. (5)), and then numerically take the limit as $`\delta 0`$ (Eq. (7)), we find $`\rho ^0(0.85;\widehat{𝐳})=\rho _{xx}^{\mathrm{out}}(\mathrm{};0.85;\widehat{𝐳};0)=7.1`$ $`\mu \mathrm{\Omega }`$cm which is the calculated resistivity for Ni0.85Fe0.15 to be compared with experiment. It should be noted, however, that the extrapolation procedures give an upper limit for the resistivity since due to computational constraints our calculations can be done only up to $`n45`$ with $`\delta 2`$ mRyd (i.e., the slopes of $`n\rho _{\mu \mu }(n;𝐜;𝒞;\delta )`$ might be slightly too high, and close to $`\delta =0`$ there might be small deviations from linearity). Parenthetically, a similar extrapolation procedure on a film with no disorder, i.e., $`c=1.0`$, yields finite $`\rho _{xx}(n;1.0;𝒞;\delta )`$ in both cases (a) and (c) for finite $`\delta `$; however, the ($`n\mathrm{}`$, $`\delta 0`$) extrapolated value $`\rho _{xx}^{\mathrm{out}}(\mathrm{};1.0;𝒞;0)`$ for outgoing boundary conditions, as well as the ($`\delta 0`$) limit $`\rho _{xx}^{\mathrm{ref}}(n;1.0;𝒞;0)`$ for reflecting boundary conditions, are both zero.

## III Results

### A Bulk alloys

For case (b) of Eq. (4), Eq. (5) is precisely the procedure pursued in experimental studies; see Fig. 6 of McGuire and Potter  where the thickness dependence of the resistivity of Ni0.82Fe0.18 films is shown. It is well known that the film resistivity increases as films become thinner; the resistivity noticeable increases when the film thickness becomes much smaller than the mean free path of the conduction electrons. It is quite interesting to realize that the measured resistivity reaches its asymptotic bulk value only for film thicknesses well above 500Å. In Fig. 2 we display the thickness dependence of our calculated $`\rho _{xx}^{\mathrm{out}}(n;c;\widehat{𝐳};0)`$ for case (a) of Eq. (4) as a function of the number of layers $`n`$ considered in the summation in Eq. (1). As can be seen, for small values of $`n`$ the corresponding changes in $`\rho _{xx}^{\mathrm{out}}(n;c;\widehat{𝐳};0)`$ are large, whereas —just as for the experimental data— for large enough $`n`$ $`\rho _{xx}^{\mathrm{out}}(n;c;\widehat{𝐳};0)`$ approaches towards its asymptotic value, i.e., the functional behavior of $`n\rho _{xx}^{\mathrm{out}}(n;c;\widehat{𝐳};0)`$ becomes linear in $`n`$. Inspecting our ($`n\mathrm{}`$) extrapolation we find (again in accordance with experiment) that above 585Å ($`\widehat{=}`$ 330 monolayers) the deviations from the bulk residual resistivity are smaller than 1 $`\mu \mathrm{\Omega }`$cm, i.e., in order to talk about “bulk” in terms of resistivities one has to consider a tremendously thick film of several thousands Å. It should be noted that in order to arrive at bulk resistivities the same kind of linear extrapolation procedure to very large $`n`$ (Eq. (5)) is done both with experimental data and with our calculated resistivities (Eqs. (1) and (2)).

In Fig. 3 we show our theoretical asymptotic bulk resistivities $`\rho ^0(c;\widehat{𝐳})=\rho _{xx}^{\mathrm{out}}(\mathrm{};c;\widehat{𝐳};0)`$ obtained by performing both limits ($`n\mathrm{}`$ and $`\delta 0`$) for the entire concentration range together with low temperature (4.2 K) experimental data  for Ni concentrations above 70%. The overall trends of measured data can be obtained from an inspection of Fig. 1 of the section on NicFe1-c in Ref.  which shows the temperature dependence of the resistivity for the whole range of concentrations. From this figure one can see that by lowering the temperature the electrical resistivity becomes increasingly sensitive to the structural phase transition from bcc to fcc; at the temperature of liquid nitrogen ($``$195 C, the lowest isothermal curve displayed) in the Fe-rich bcc $`\alpha `$-phase the resistivity reaches a maximum at about 15% of Ni and then slowly decreases up to 30% of Ni, whereas in the Ni-rich fcc $`\gamma `$-phase the resistivity starts to grow below 50% of Ni and seems to diverge near the critical concentration for the structural phase transition. In the $``$195 C curve we also discern a weak shoulder at about 85% Ni; this feature is well resolved in the low temperature (4.2 K) experimental data.  As can be seen from Fig. 3 our theoretical values reproduce very well the vicinity of the structural phase transition from bcc to fcc at about 35% of Ni. There is indeed a maximum in the bcc phase at about 15% of Ni; the onset of the phase transition is clearly visible. Also —just as the experimental data— in the fcc phase the theoretical values show a kind of shoulder at about 85% of Ni.

For this concentration ($`c=0.85`$) we performed an additional calculation using a maximum angular momentum quantum number $`\mathrm{}_{\mathrm{max}}=3`$ (see also Ref. ). As can be seen from Fig. 3, the difference between the $`\mathrm{}_{\mathrm{max}}=3`$ and the $`\mathrm{}_{\mathrm{max}}=2`$ (that is used in all other calculations) results is rather small; however the trend for the inclusion of higher scattering channels is to lower (slightly) the resistivity.

Although the theoretical resistivities agree rather well with the experimental ones, the difference at a particular concentration being of the order of 3 $`\mu \mathrm{\Omega }`$cm, it seems appropriate to comment on the putative sources of the difference. First of all, in the present use of the Kubo-Greenwood equation for layered systems no vertex corrections arising from the configurational average of the product of two Green’s functions are included (see also the discussion in Ref. ). In another calculation  of the resistivity of bulk NicFe1-c the vertex corrections were calculated and found to be quite small over the entire range of concentrations,  i.e., no larger than the differences between the $`\mathrm{}_{\mathrm{max}}=2`$ and $`\mathrm{}_{\mathrm{max}}=3`$ results at $`c=0.85`$ in Fig. 3. Second, the single site approximation to the coherent-potential approximation is used to describe the electronic structure of substitutional alloys. This in turn implies that short range order or concentration fluctuations are excluded. Both the inclusion of vertex corrections and the inclusion of short range order will reduce the present theoretical resistivity values. Furthermore, in magnetic alloys magnetic ordering or clustering can occur; as it implies some ordering, which we have not taken into account, it will again (slightly) lower the resistivity.

In two previous theoretical papers  the validity of the so-called two-current model for transport in ferromagnetic systems with strong spin-dependent disorder was tested; the resistivity and the anisotropic magnetoresistance for NicFe1-c bulk alloys was calculated by solving the Kubo-Greenwood equation using a spin-polarized, relativistic version of the Korringa-Kohn-Rostoker method together with the coherent-potential approximation. As these calculations were done on a fcc lattice for the whole concentration range their results are not applicable to the Fe-rich bcc alloys. Nonetheless a very weak shoulder at about 80% Ni can be discerned from their Fig. 4  —in agreement with the experimental data and our present calculation, see Fig. 3. However, the resistivities they find for 80% Ni, $`\rho 1.8\mu \mathrm{\Omega }`$cm, as well as the maximum at about 40% Ni, $`\rho 3.2\mu \mathrm{\Omega }`$cm, are too low by about a factor of two as compared to available experimental data.

### B Spin-valve structures

The main purpose of our study has been to show how well resistivities for disordered magnetic systems with only two-dimensional translational symmetry (films) can be described, and to confirm that by taking the appropriate limit we derive the bulk resistivities from those calculated for films. Since a large number of spin-valve structures that are used in applications, e.g., the read heads of hard disk drives, contain permalloy layers, we conclude with some results on their transport properties. Including permalloy layers into a multilayer system the electronic structure and the magnetic properties are considerably more complicated than for bulk permalloy; specifically they are sensitive to the thickness and Ni concentration $`c`$ of the permalloy layers, and to the constituency of the neighboring layers. Therefore for concreteness we focus on a permalloy-based system that has been studied, e.g., Si-substrate / Cu-contacts(1500Å) / NicFe1-c(100Å) / Co(6Å) / Cu(9Å) / Co(6Å) / NicFe1-c(5Å) / NiO(250Å),  and we replace it with the reference system $`\text{Ni}_c\text{Fe}_{1c}(100)`$ / $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_{12}`$ Co4 Cun Co4 $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_3`$ / Vacuum. In all cases shown a fcc(100) stacking was chosen, the lattice spacing refers to that of the corresponding fcc bulk NicFe1-c system. In this transcription the rather thick permalloy layer of about 100Å ($`\widehat{=}`$ 56 monolayers) is considered as a substrate, the insulating NiO layer is replaced by a free surface (vacuum), and all other thicknesses are expressed (in monolayers) as close as possible to the measured structure. The 12 layers of NicFe1-c are deliberately chosen so as to guarantee a smooth matching to the substrate via a selfconsistent calculation; this in turn would also allow us to change the left boundary condition such that in principle a system of the type Vacuum / $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_{56}`$ Co4 Cun Co4 $`\left(\text{Ni}_c\text{Fe}_{1c}\right)_3`$ / Vacuum could be described by assuming that the layers in the very interior of the thick NicFe1-c layer exhibit bulk-like properties.

We have considered two series: in the first the Cu-spacer thickness is confined to $`n=5`$ and we vary the Ni concentration $`c`$ while maintaining an fcc lattice throughout; in the second we set $`c=0.85`$ and vary the number of Cu-spacer layers $`n`$. While the structure of NicFe1-c bulk alloys is bcc for $`c<35\%`$ when there are relatively thin layers in a fcc multilayered structure it is appropriate to maintain the fcc lattice structure even for the Fe-rich alloys. In the parallel configuration of the spin valve the magnetizations of all permalloy and Co layers are aligned uniformly perpendicular to the plane of the layers; in the antiparallel configuration they are still normal to the layers but the magnetizations of the two magnetic slabs (separated by the Cu-spacer) point in opposite directions to one another. The conductivity is evaluated at a finite imaginary part to the Fermi energy, $`\delta =2`$ mRyd; this is necessary for otherwise the calculations would take an inordinate time. In principle, in order to obtain the resistivity we have to take the limit as $`\delta 0`$ (see Eq. (7)). However, it is virtually impossible for technical reasons to do this up till now for a whole series of spin-valve structures; therefore we do not show the resistivities we calculated for finite $`\delta `$. But with the assumption that both resistivities for the parallel and the antiparallel configuration of the spin valve scale uniformly with $`\delta `$, and therefore cancel out in the ratio $`R`$ (Eq. (3)), we have used them to estimate the CIP-MR ratios that one can expect from these spin-valve structures.

In Fig. 4a we show the variation of the calculated CIP-MR ratio of the reference system with the Ni concentration $`c`$ for the fixed Cu-spacer thickness $`n=5`$. As can be seen the CIP-MR varies only very little (between about 17.5–21.5%) over the whole range of concentrations considered. In Fig. 4b we show the variation of the CIP-MR ratio for the reference system corresponding to $`c=0.85`$ with respect to the number of Cu-spacer layers $`n`$. By doubling the number of Cu-spacer layers, i.e., by going from $`n=4`$ to $`n=8`$, the CIP-MR is reduced by almost a factor of two. This dramatic decrease can be understood in terms of a simple model: since on increasing the thickness of the nonmagnetic spacer layer $`t_{\mathrm{nm}}`$ the CIP-MR decreases as $`(1/t_{\mathrm{nm}})\mathrm{exp}(t_{\mathrm{nm}}/\lambda _{\mathrm{nm}})`$,  where $`\lambda _{\mathrm{nm}}`$ is the mean free path of the conduction electrons in the nonmagnetic spacer layer, the ratio $`R_{xx}(t_{\mathrm{nm}})/R_{xx}(2t_{\mathrm{nm}})2`$ for $`t_{\mathrm{nm}}\lambda _{\mathrm{nm}}`$ (which is a necessary condition to find at all a CIP-MR and certainly valid for our $`t_{\mathrm{nm}}`$ varying between 7.1 Å and 14.2 Å). Our findings seem to be confirmed by experiments made on similar systems  albeit at room temperature, while our calculated values are for $`T=0`$ K. Depending on the actual thicknesses of: the bottom and top NicFe1-c layers, of the Co slabs and of the Cu spacer, and the pinning and capping layers used, the experimental CIP-MR ratios in permalloy-related spin-valve systems have been found to be in the range between 10 and 20%. Fig. 4 nicely shows that this is about the range that our theoretical calculations predict when we use the resistivities calculated with finite $`\delta `$, and when we make the assumption that for both the parallel and the antiparallel configurations of the spin valve the resistivity scales uniformly with $`\delta `$ and therefore cancels out in the CIP-MR ratio $`R`$.

## IV Summary

By using our fully relativistic, spin-polarized, screened Korringa-Kohn-Rostoker method together with the coherent-potential approximation for layered systems we have calculated the resistivity across the entire NicFe1-c (permalloy) series of bulk alloys within the Kubo-Greenwood approach. We are able to reproduce the overall variation of the resistivity as a function of the Ni concentration; in particular we reproduce the discontinuity in the vicinity of the structural change from bcc to fcc. Also we obtain very reasonable absolute values for the resistivities. We have also estimated the CIP-MR ratios for permalloy-based spin-valve structures and within the limits of our current calculation find ratios in reasonably good agreement with those measured on prototypical spin valves; there is of course the caveat that most of the data on spin valves is quoted at room temperature while our results are for $`T=0`$ K. It should be stressed that our results are based on ab initio calculations with no adjustable parameters. Viewed in this light the calculated resistivities and GMR ratios are noteworthy.

###### Acknowledgements.

This work was supported by the Austrian Science Foundation (T27-TPH, P12146), the Hungarian National Science Foundation (OTKA T030240, T029813), the Defense Advanced Research Projects Agency and the Office of Naval Research (N00014-96-1-1207, MDA 972-99-C-009), NATO (CRG 960340), and the TMR network (EMRX-CT96-0089). We also wish to thank the CNRS IDRIS Computing Center at Orsay for calculations done on their T3E Cray.


text:
# 1 Introduction

## 1 Introduction

A thorough understanding of epistatic interactions of mutations in genomes is becoming more and more crucial to many areas in population genetics and evolutionary biology. Epistasis affects linkage disequilibria (Charlesworth 1976, Barton 1995), robustness to mutations (Lenski *et al.* 1999) or canalization (Nowak *et al.* 1997, Wagner *et al.* 1997, Rice 1998, Ancel & Fontana 2000, for reviews, see Scharloo 1991, Gibson & Wagner 2000), as well as theories on the maintenance of sex (Kondrashov 1982, Kondrashov 1988, West *et al.*1999). The sign of epistatic effects, that is, whether deleterious mutations are reinforcing (synergistic epistasis) or mitigating (antagonistic epistasis) also influences whether or not deleterious mutations can accumulate in the genome via Muller’s ratchet (Muller 1964, Felsenstein 1974). The consensus seems to be that synergistic epistasis can prevent the accumulation of mutations (Crow & Kimura 1978, Kondrashov 1994, but see Butcher 1995 for a dissenting view). On the other hand, the observation of pervasive compensatory mutations (Moore *et al.* 2000), which also render the ratchet powerless, indicates epistasis, but not its sign.

While the genomes of a number of organisms have been examined for signs of epistasis (de Visser *et al.* 1997a, 1997b, de Visser & Hoekstra 1998, Elena & Lenski 1997, Elena 1999), no general trend can be discerned except to say that interactions between mutations are frequent and of both signs, and that weak synergistic epistasis seems to prevail in eukaryotic genomes while viral and prokaryotic genomes show no net tendency in either direction. Experiments to measure epistatic interactions are difficult and usually yield results of weak statistical significance (West *et al.* 1998). Consequently, even epistasis of considerable strength can conceivably be missed in vitro or in vivo. Here, we investigate deleterious mutations in silico, and study a fact which has not received much attention in the population genetics literature. Namely, epistasis is closely related to the geometry of the phenotype space (Rice 1998), which leads to interesting relations between epistasis and the effects of single mutations. Wagner *et al.* (1998) showed in a two-locus, two-allele model that the average effect of a single mutation is correlated to the degree of interaction between loci, and supported their theoretical model with QTL data on body weight in mice. We demonstrate here a similar relation for a measure of epistasis which is more appropriate for high dimensional sequence spaces. We give both theoretical and experimental evidence that the strength of directional epistasis is correlated with the average deleterious effect of a single mutation. As a corollary of this observation, we argue that in situations in which there is a selective pressure to reduce the average deleterious effect, this correlation leads to a reduction of the number of genomes with strong antagonistic effects in a population.

## 2 Neutrality and Epistasis

It is a common observation that the average fitness at a mutational distance $`n`$ from a given reference sequence decays roughly exponentially with $`n`$ (see, e.g., de Visser *et al.* 1997b, Elena & Lenski 1997, Lenski *et al.* 1999, and also our results on RNA sequences below). The simplest explanation for a *perfect* exponential is a multiplicative landscape in which every mutation diminishes the fitness independently by the same factor $`(1s)`$. Here, however, we focus on fitness landscapes with a considerable amount of lethal mutations, in which case a branching process in the high-dimensional sequence space is a more appropriate model: if every viable sequence has a probability of $`1s`$ to remain viable after suffering a single point mutation, then the mean fitness will decay as (neglecting fitness differences in the viable sequences)

$$w(n)=(1s)^ne^{\alpha n},$$
(1)
where we have defined $`\alpha =\mathrm{ln}(1s)`$, and assuming that the fitness of our reference sequence at $`n=0`$ is $`w(0)=1`$. Both explanations for the exponential decay have in common that the effects of subsequent mutations are independent, i.e., there exists no epistasis. If there are some interactions between mutations, we will still observe an exponential decay if antagonistic and synergistic interactions occur in the same proportion. If however there exists a bias towards either antagonistic or synergistic epistatic interactions (directional epistasis), then this bias will naturally appear as deviations from the exponential decay in Eq. (1). While such deviations have previously been indicated by adding a term quadratic in $`n`$ to the exponent of Eq. (1) (Crow 1970, Charlesworth 1990, de Visser *et al.* 1997b, Elena & Lenski 1997), such a parameterization becomes troublesome at larger $`n`$, because $`w(n)`$ could increase beyond the fitness of the reference sequence (for a positive coefficient in the quadratic term). This is avoided by the ansatz (Lenski *et al.* 1999)

$$w(n)=e^{\alpha n^\beta },$$
(2)
where $`\beta =1`$ means there is no bias towards either form of epistatic interactions. A $`\beta >1`$ indicates synergistic mutations prevail (mutations that are on average “worth” more than one independent hit), while $`\beta <1`$ reflects a bias towards antagonistic mutations (mutations whose “damage” is less then one independent mutation).111Note that in the earlier works where a quadratic term was used, the distinction between the different types of directional epistasis depended on whether $`\beta `$ was larger or smaller than 0, rather than 1. Since expression (2) depends only on two parameters, deviations from that form may arise when $`n`$ grows large.

Naively, one might assume that the decay parameter $`\alpha `$ and the epistasis parameter $`\beta `$ are independent. Instead, we shall see that environments with strong selection force a trade-off between $`\alpha `$ and $`\beta `$, so that one can only be optimized at the expense of the other. The reasoning is as follows. In a strongly selective environment mutations can be classified as either neutral or lethal, and $`w(n)`$ can be thought of as the fraction of neutral sequences in genetic space at mutational distance $`n`$. In particular, the *neutrality* $`\nu `$ of a sequence (the number of sequences at Hamming distance 1 with fitness 1) is related to the decay parameter by $`\nu =\mathrm{}(D1)e^\alpha `$, where $`\mathrm{}`$ is the length of the sequence, and $`D`$ is the number of monomers. If all sequences in genetic space have the same $`\nu `$, it follows that $`\beta =1`$. A deviation from $`\beta =1`$ implies that some sequences have more or fewer neutral neighbors than others, giving rise to a correlation between $`\alpha `$ and $`\beta `$. For a viable sequence with lower than average neutrality (higher than average $`\alpha `$), there are comparatively fewer sequences close-by than there are far away, such that this sequence will have a small $`\beta `$. Conversely, a sequence with a high neutrality (small $`\alpha `$) will have comparatively more sequences close by, and $`\beta `$ will be larger. We can make this argument more formal with a simple “conservation law”, which only reflects that the total number of neutral sequences in genetic space is constant. Since for polymers of fixed length $`\mathrm{}`$ made from $`D`$ monomers there are $`\left(\genfrac{}{}{0pt}{}{\mathrm{}}{n}\right)(D1)^n`$ possible $`n`$-mutants, we must have

$$\underset{n=1}{\overset{\mathrm{}}{}}w(n)\left(\genfrac{}{}{0pt}{}{\mathrm{}}{n}\right)(D1)^n=N_\nu ,$$
(3)
where $`N_\nu `$ is the total number of neutral mutants of this wild type. Inserting $`w(n)`$ from Eq. (2) yields an implicit relation between $`\alpha `$ and $`\beta `$. Indeed, for two decay functions of two different reference sequences, with different parameters $`\alpha `$, the only way in which the sum can yield the quantity $`N_\nu `$ (which is independent of the respective reference sequence) is that the two decay functions must have different parameters $`\beta `$ as well. However, the implicit relation depends on the ansatz Eq. (2) being correct for all $`n`$, which is not necessarily the case. Alternatively, we may consider only sequences with up to $`d`$ mutations, in which case we can write

$$\underset{n=1}{\overset{d}{}}w(n)\left(\genfrac{}{}{0pt}{}{\mathrm{}}{n}\right)(D1)^n=N_{\nu ,d},$$
(4)
where $`N_{\nu ,d}`$ is the number of neutral sequences in a sphere of radius $`d`$ around the reference sequence with $`w(0)=1`$. $`N_{\nu ,d}`$ depends on the particular reference sequence chosen, and therefore cannot be regarded a constant. However, for $`d`$ not too small, we may replace it by its average $`N_{\nu ,d}`$ over all viable reference sequences. If we take into account a sufficiently large region of genotype space, we should find roughly the same number of neutral mutants for each viable sequence, since, as $`d`$ approaches $`\mathrm{}`$, the quantity $`N_{\nu ,d}`$ approaches $`N_\nu `$, which in turn is independent of the reference sequence. Hence, Eq. (4) predicts a similar relation between $`\alpha `$ and $`\beta `$, and the two different predictions approach each other as $`dl`$. Predictions based on Eqs. (3) and (4) are used below to compare to our empirical results.

Although the above argument strictly holds only under the assumption that mutations are either neutral or lethal, it is not unreasonable to assume that a similar (possibly weaker) correlation between $`\alpha `$ and $`\beta `$ exists also in more general cases, where slightly deleterious or even advantageous mutations are possible. In that case, under the presence of epistasis, there will still be regions in genotype space in which the number of less-deleterious mutations is higher, and other regions in which it is lower than average. The decay function $`w(n)`$ of a sequence from a region that is rich in non-lethal mutations would have a lower $`\alpha `$, but would be inevitably more synergistic than the decay function of a sequence from a region poor in non-lethal mutations. Our results with digital organisms (see below) support this reasoning.

## 3 Experimental Evidence

Accurate data for the decay parameter $`\alpha `$ and the epistasis parameter $`\beta `$ for biological organisms are rare, which makes our hypothesis difficult to test. A few well-studied systems have emerged which are accessible in silico.

We studied RNA secondary structure prediction using the Vienna RNA package, version 1.3.1., with the default setup (Hofacker *et al.* 1994). We calculated the decay of the average number of neutral folds as a function of the Hamming distance for 100 random RNA sequences of length $`\mathrm{}=76`$. The parameters $`\alpha `$ and $`\beta `$ were determined as follows. We obtained $`\alpha `$ exactly from the fraction of neutral one-mutants. In addition, we sampled the function $`w(n)`$ for Hamming distances up to $`n=8`$, by calculating the structure of up to $`10^6`$ random neighbors of the required Hamming distance. The quantity $`\beta `$ was then determined from a nonlinear fit of $`\alpha n^\beta `$ to the logarithm of $`w(n)`$. A plot of $`\beta `$ versus $`\alpha `$ (Fig. 1) shows a significant correlation, with a correlation coefficient of $`r=0.817`$ ($`p<0.01`$).

According to Eqs. (3) and (4), we can predict the relationship between $`\alpha `$ and $`\beta `$ if we compare the decay functions of sequences that are mutually neutral. For RNA folding, this means we have to determine $`\alpha `$ and $`\beta `$ for a set of sequences that fold into the same structure. We performed experiments with the RNA sequences of length $`\mathrm{}=18`$ used by van Nimwegen *et al.* (1999) (for these sequences we altered the default setup by setting the free energies of dangling ends to zero.) For the particular case that all bonds are of the purine-pyrimidine type (G-C, G-U, A-U), two separate neutral networks (a neutral network is a network of neutral genotypes connected to each other by one-point mutations) were found by van Nimwegen *et al.*, consisting of 51,028 and 5,169 sequences, respectively. For each such set of neutral sequences, Eq. (3) predicts the correlation without free parameter, as long as the number of all neutral sequences $`N_\nu `$ is known. In order to estimate $`N_\nu `$, we generated $`10^8`$ random sequences of length $`\mathrm{}=18`$, of which 10,961 sequences folded correctly. From this, we estimated $`N_\nu =7.5\times 10^6`$ neutral sequences out of the total $`6.9\times 10^{10}`$ sequences of length $`\mathrm{}=18`$. Using this number, Eq. (3) predicts the solid line in Fig. 2, which describes the correlation well. The approach based on averaging over *local* neutral sequences around the reference sequence \[Eq. (4)\] gives rise to the dotted line in Fig. 2 and shows that this too predicts the correlation fairly well (for this second approach, we used $`d=8`$, and from Monte Carlo simulations for 1000 reference sequences, we obtained $`N_{\nu ,8}=2.4\times 10^5`$).

As our second test case, we analyzed the correlation between $`\alpha `$ and $`\beta `$ in digital organisms (Adami 1998, Adami *et al.* 2000). Digital organisms are self-replicating computer programs that mutate and evolve. Their fitness is determined by the ratio of their CPU speed and their gestation time. The latter is given by the number of instructions that have to be executed in order to generate a fully functional offspring. The CPU speed increases when the digital organisms perform logical operations on numbers that they can acquire from their environment. The gestation time, on the other hand, increases either if the organisms employ less efficient mechansims of self-replication, or if they accumulate instructions that are involved in the completion of the above mentioned logical operations. The digital organisms with the highest fitness thus have a very efficient copy mechanism, and perform a large number of logical operations with a comparatively small number of additional instructions. Lenski *et al.* (1999) measured the decay of the mean fitness as a function of the number of mutations accumulated in such digital organisms, and obtained $`\alpha `$ and $`\beta `$ from a fit of Eq. (2) to the measured decay functions. They studied 174 different genomes consisting of two groups of 87 genomes each. The first group of organisms evolved in 87 independent experiments in a complex environment, while the second group was obtained by transfering these organisms to an environment which favoured simple genomes, and allowing the organisms to adapt to this more simple environment.

A statistical analysis revealed a significant correlation between the decay parameter $`\alpha `$ and the parameter of directional epistasis $`\beta `$ for both the complex and the simple digital organisms (Table 1). However, in addition to said correlation, we found a correlation between the $`\alpha `$ and both the genome length $`\mathrm{}`$ and the log fitness for complex and simple organisms, as well as a correlation between $`\beta `$ and $`\mathrm{}`$ in the case of the complex organisms. Hence, for the complex organisms, we cannot rule out the possibility that the correlation between $`\alpha `$ and $`\beta `$ merely reflects an underlying correlation of both quantities with length. In the case of the simple organisms, where we do not see a correlation between $`\beta `$ and $`\mathrm{}`$, we can assume that the correlation between $`\alpha `$ and $`\beta `$ is genuine. To provide further evidence, we examined a reduced data set of all 48 simple organisms with a length between $`\mathrm{}=14`$ and $`\mathrm{}=16`$. These 48 organisms were also of comparable fitness. In that data set, we found an even stronger correlation between $`\alpha `$ and $`\beta `$, while the correlation between either of the two quantities and length or fitness was insignificant (Table 1). It was not possible to study a similar reduced data set for the complex organisms, because the variations in length were too large (the length varied between $`\mathrm{}=20`$ and $`\mathrm{}=314`$ among the 87 genomes). Although the data from the digital organisms is potentially biased because our reference sequences were evolved rather than random as in the RNA case, we believe this bias to be insubstantial in at least the reduced data set, which contains only organisms with very similar phenotypes. Hence, just as for our RNA data, the data from the digital organisms supports our hypothesis of a genuine correlation between $`\alpha `$ and $`\beta `$.

## 4 Adaptation of Epistasis through Correlated Response

The correlation between neutrality and epistasis implies that if one of them is subject to selective pressures, the other will be affected as well due to correlated response. Van Nimwegen *et al.* (1999) have shown that a population evolving on a neutral network reduces its genetic load by moving into the regions of high neutrality in sequence space. In particular, given a random population of molecules on such a network, evolution tends to increase the neutrality in the population, effectively pushing the population into the center of the neutral network. Because of the correlation between neutrality and epistasis, we expect this dynamic to lead to a reduction of antagonistic epistatic effects. To verify this hypothesis, we carried out evolutionary experiments with the RNA sequences of length 18 from the previous section.

We performed one flow-reactor run for each of the two networks found by van Nimwegen *et al.*, starting with an initial population of 1,000 sequences chosen at random from the respective network. We set the replication rates such that sequences folding into the target structure replicated on average once per unit time, while the replication rate of all other sequences was set to $`10^6`$ per unit time. All sequences had a probability of $`\mu =0.5`$ to suffer one random point mutation per replication event. The possibility of several point mutations per replication event was eliminated, to guarantee that the population could not leave the specified neutral network. The epistasis parameter $`\beta `$ was determined for every sequence in the population every two hundred generations, while the population neutrality was monitored constantly. The population neutrality $`\overline{\nu }`$ is the average neutrality of all sequences currently in the population. In Fig. 3, we present the results from the run on the larger of the two networks. The neutrality of the initial population coincides with the network neutrality (the average neutrality of all sequences on the network), which is to be expected for a random initial population. Over the course of evolution, the average population neutrality rose to the predicted equilibrium value (given by the spectral radius of the connectivity matrix, see van Nimwegen *et al.* 1999.) As expected, the average epistasis parameter $`\overline{\beta }`$ increased significantly as well. Results on the second network were qualitatively identical, with $`\beta `$ increasing from 0.78 to around 0.86. Thus, antagonistic epistasis is reduced during adaptation for reduced mutational load on a neutral network.

The above reasoning depends of course on the assumption that a population remains on a single neutral network. In a more realistic scenario, where peaks of different heights are present, the main effect of selection will be to increase the fitness, rather than to increase $`\alpha `$. However, once a local optimum has been reached, we can expect the dynamic described above to take place, as long as there exists some neutrality at the local peak. If there was a correlation between the fitness and $`\alpha `$ or $`\beta `$, then we would see additionally a correlated response in $`\alpha `$ or $`\beta `$ as the fitness is being maximized. Nevertheless, while the correlation between $`\alpha `$ and $`\beta `$ reported here is a general result that follows from geometric constraints on the landscape, no such constraints exist between the fitness and $`\alpha `$ or $`\beta `$. The parameters $`\alpha `$ and $`\beta `$ measure the amount and the distribution of neutral or nearly neutral sequences in the neighborhood of a reference sequence. There is no reason why sequences with high fitness should generally be found in regions with particularly high or low neutrality, or with a particular distribution of the neutral sequences. This is not to say, however, that such a correlation cannot exist in special cases (see, e.g., our data from the complex digital organisms, Table 1, where $`\alpha `$ is correlated with fitness, but $`\beta `$ is not).

## 5 Conclusions

Epistasis plays an important role in evolutionary theory, but remains empirically largely unexplored. Using secondary structure prediction of RNA sequences as well as digital organisms evolving in silico, we have demonstrated a correlation between two important parameters of realistic genetic fitness landscapes: the average deleterious effect of single mutations and the strength of directional epistasis. In conjunction with the results from Wagner *et al.* (1998) for mice, and the two different theoretical explanations (two-locus model in the case of Wagner *et al.*, sequence space based model here), we can expect this correlation to be an ubiquitious phenomenon, present in many natural and artificial fitness landscapes.

The correlation, coupled with the selective pressure which forces random sequences in a neutral network to cluster in the dense areas of the network, leads to a reduction of strong antagonistic epistasis in a population. The nature of this result is purely geometric: as a population tries to reduce the average effect of single mutations, the effect of multiple mutations is inevitably worsened as long as there exist some inhomogeneities in the effect of single mutations across the genotype space. As the result of this geometric constraint, a member of an evolved population will have, on average, a higher $`\beta `$ than a random sample of the fitness landscape would indicate.

It is well known that antagonistic epistasis favors the accumulation of deleterious mutations as well as the operation of Muller’s ratchet. Since in such a situation sexual recombination (within a fixed environment) tends to worsen the loss of information, recombination is unlikely to evolve or be maintained. The mechanism described here may thus provide a path towards an environment more conducive to the evolution of recombination.

## 6 Acknowledgements

We thank Martijn Huynen for providing access to the sequences used in (van Nimwegen *et al.* 1999), and Walter Fontana for providing us with his flow reactor code. We are grateful to Richard Lenski and two anonymous referees for many useful comments and suggestions regarding the manuscript. The genomes of digital organisms used in this study are available at
http:dllab.caltech.edu/pubs/nature99/nature.shtml. This work was supported by the National Science Foundation under contract No. DEB-9981397. Part of this work was carried out at the Jet Propulsion Laboratory, under a contract with the National Aeronautics and Space Administration.