text:
# Data Collection for the Sloan Digital Sky Survey — A Network-Flow Heuristic

## 1 The Sloan Digital Sky Survey

> “The Sloan Digital Sky Survey \[is\] a joint project of the Astrophysical Research Consortium. … The goal of the project, which is scheduled to begin in 1997 and take five years, is to make a much better map of the universe than is currently available. The volume of the universe to be surveyed will be 100 times larger than the volume of previous surveys. The number of galaxies with known distances is expected to increase by a factor of 100 to 1,000,000 galaxies and the number of quasars to increase to 100,000.
>
> “The Sloan Foundation … has contributed $8 million to the $18 million capital costs of the project. …
>
> “In order to do the survey, ARC is designing and building a special purpose 2.5 meter (100-inch) telescope at its Apache Point Observatory. …
>
> “\[The Sky Survey will proceed in two phases. In the first phase, a two-dimensional map of the sky will be made. For the second phase, the\] million brightest stars and the one hundred thousand brightest quasars will be selected for spectroscopic analysis from the two-dimensional map…
>
> \[29, Savani, 1994\]

To gather the spectroscopic data in the second phase, the telescope will be pointed repeatedly at the sky to take a series of “snapshots”. Each snapshot will capture data for up to $`660`$ galaxies and quasars in the circular portion of the sky visible through the telescope. For each captured galaxy, light from that galaxy will enter the telescope and travel through an optical fiber to a spectral analyzer. The optical fibers (one for each galaxy) will be held in place by a “plug plate” drilled to hold the up to $`660`$ fibers, each aligned to accept the light of its respective galaxy .

### 1.1 A Capacitated Covering Problem.

The second phase of the survey is expected to cost on the order of $4-5 million. This cost will depend primarily on the number of snapshots taken. This paper concerns the following problem: given the “2-dimensional” locations of the desired galaxies, determine a minimum-size set of snapshots that capture them. Formally:

> Euclidean Capacitated Covering by Disks (ECCD)
>
> Given a collection of points on the unit sphere, a radius $`r`$, and a capacity $`c`$, find a small set of discs of radius $`r`$ (located on the sphere) such that each given point can be assigned to a disc containing it, with no disc being assigned more than $`c`$ points.

The sphere corresponds to the view-sphere centered at the telescope. The points correspond to the images of the galaxies projected on the view-sphere. Each disc represents one snapshot to be taken through the telescope; the points assigned to that disc correspond to those galaxies for which data will be collected in that snapshot. The capacity $`c`$ is the maximum number of galaxies for which spectral data can be gathered in a single snapshot (due to limitations in packing the optical fibers).

The ECCD problem is NP-hard . The instances we need to solve will have hundreds of thousands of points. Luckily, as Figure 1 illustrates, the instances we need to solve are nicely structured.

In this paper we describe a heuristic algorithm for the problem. The algorithm is effective for instances arising in the Survey and will be used for it. The basic idea behind the algorithm is to start with a near-uniform cover of the sphere by discs  and then to iteratively improve the cover. The key observation is that a given cover can be improved by first solving a relaxation of the problem in which the “point-in-disc” constraints are replaced by penalties for assigning points to discs not containing them, and then moving the discs to minimize the cost of the assignment found. The relaxed problem reduces to the minimum-cost flow problem. In our tests, the algorithm runs in nearly linear time and finds covers that are roughly 20% better than comparable near-uniform covers.

## 2 Related Work

The NP-completeness of the variant when the points lie in the plane was proven by Megiddo and Supowit . The proof adapts easily to our problem. The NP-completeness of the planar problem when the discs are required to be centered on the given points was proven by Marchetti-Spaccamela . When the covering regions are rings, instead of discs, Maass  showed the problem NP-complete even if the points all lie on a single line.

Papadimitriou  (improving results by Fisher and Hochbaum ), considered the related $`p`$-medians problem in the plane, which is that of covering the given points with $`p`$ discs (of arbitrary radii, but centered at $`p`$ of the given points) so as to minimize the sum of the disc radii. He showed the problem to be NP-complete and presented average-case analyses of several algorithms. One of the heuristics is a uniform (“honeycomb”) covering of the points by discs, which he shows gives a near-optimal solution with high probability when $`p`$ is $`\omega (\mathrm{log}n)`$ and $`o(n/\mathrm{log}n)`$ and the points are randomly distributed in the unit square.

The problem can be modeled as a capacitated set-cover problem. The well-known greedy algorithm of Johnson  and Lovász , as modified for the capacitated case by Bar-Ilan, Kortsarz, and Peleg , would yield a $`\mathrm{ln}n`$-approximate solution, where $`n`$ is the number of galaxies. This algorithm is not good enough in practice. For this particular set-cover problem the dual of the set system has bounded VC-dimension; in this case an improved approximation algorithm is known for the uncapacitated case , but, judging from a few small trials, this algorithm does not appear to take sufficient advantage of the structure of our problem instances to perform well in practice.

Numerous generalizations of our problem have been considered under various names, including “(un)capacitated facility (or plant) location,” “$`p`$-centers”, and “minimax facility location”. These problems have been studied under various metrics and also in general graphs. In general, polynomial-time exact algorithms are known only when the number of covering regions (in our case, discs) is small (e.g., ) or when the underlying metric space (or network) is tree-like (e.g., ). Generally, these algorithms are for uncapacitated problems.

There is a large literature on these problems in Operations Research. Relevant books include . Much of this research has concentrated on adapting integer-programming techniques to fairly general formulations of the problem. For example, recent works on the Capacitated Facility Location Problem (a generalization of our problem to arbitrary networks) include . Quoting from the conclusion of “Approximate Solutions to Large Scale Capacitated Facility Location Problems” (1990) :

> The problem of locating facilities has inspired a rich body of literature which spans well over two decades. Numerous algorithms have been devised and successfully applied to problems with as many as 200 customers and 100 facilities. The computational experience on larger problems, however, has been virtually non-existent… In the work leading to this paper, the objective was to develop a heuristic algorithm that can be used to generate effective solutions for large scale facility locations problems. The computational results obtained so far seem to indicate that this requirement can be met for problems with as many as 1000 customers and 100 facilities.

## 3 The Algorithm

The instances arising in the Sky Survey exhibit particular structure. Within any given region, the galaxies are distributed densely throughout the region, somewhat uniformly but with clustering tendencies and variation in density. The density of the galaxies means that virtually the entire region must be covered by discs. The variation in density means that more discs must be concentrated within densely populated regions. As a reference point, consider the sparsest possible covering of the area by discs (resembling a “honeycomb”). This cover provides roughly the right total capacity and does well in sparse areas, but in dense areas does not provide sufficient capacity. Any good solution will have to maintain a honeycomb-like structure in sparse areas while bunching discs more densely in dense areas.

The outer loop of the algorithm does a binary search for the smallest value of a density parameter $`\delta `$ that leads to success in the inner loop. The inner loop begins with a near-uniform cover of normalized density $`1+\delta `$ and iteratively improves it (see Figure 1 for “before” and “after” covers). Each iteration of the loop perturbs the discs, as described below, in an attempt to improve the cover (Figure 2 shows the results of such a series of improvement steps). If the desired coverage is obtained, the inner loop stops (successfully). If the perturbation ceases to improve the cover, the inner loop stops (unsuccessfully).

Next we describe how the algorithm perturbs a given cover in order to improve it. We start with the observation that for a given set of discs (with known locations), the problem of finding the maximum number of galaxies that can be assigned reduces to a generalized maximum bipartite matching problem in a graph $`G=(U,V,E)`$, where the vertices in $`U`$ correspond to galaxies, the vertices in $`V`$ correspond to the discs, and edge $`(u,v)`$ is present if $`u`$’s galaxy is in $`v`$’s disc. A maximum legal assignment of galaxies to discs then corresponds to a maximum size set $`S`$ of edges such that each $`uU`$ is incident to at most one edge in $`S`$ while each $`vV`$ is incident to at most $`c`$ edges in $`S`$.

Since the latter problem reduces in a standard way to the maximum flow problem , which can be efficiently solved, it follows that for a given set of discs, one can efficiently find a maximum legal assignment of galaxies to discs.

Of course, the maximum legal assignment may still leave many galaxies unassigned, even though many discs are not used to capacity. In this case, how can discs be moved to improve the coverage? Consider the following relaxation of the problem:

> Relaxed Problem
>
> Given a set of discs, a set of galaxies, and a capacity $`c`$, find a minimum-penalty assignment of the galaxies to discs such that no disc is assigned more than $`c`$ galaxies.

Here a galaxy can be assigned to a disc not containing it, but there is a penalty for doing so that encourages assignments of galaxies to nearby discs (details of the penalty function are in § 3.2).

The relaxed problem can be solved efficiently (even for arbitrary penalties) by reducing it to the assignment problem or to minimum-cost maximum flow. We reduced it to the latter, more general, problem in anticipation of having to incorporate more complex constraints on the assignment (that no sufficiently close pairs of galaxies should be assigned to the same disc) at a later point. As described below, even the more general problem can be solved quickly enough for our purposes.

A solution to the relaxed problem will assign all galaxies to discs, but a given disc may be assigned galaxies outside of it. The advantage of the relaxed problem is that a solution to it can give information about how to improve a given set of discs. The intuition is that if excess demand (i.e. a high density of galaxies relative to discs) exists in one area, and excess capacity exists in another, then a disc between the two areas will tend to be assigned galaxies that are outside of the disc and that lie towards the area of excess demand. Figure 3 illustrates this.

Once a minimum-penalty solution to the relaxed problem has been found, the algorithm moves the discs to minimize the cost of the particular assignment of galaxies to discs specified by the minimum-penalty solution. This problem can be solved independently for each disc. For a given disc, for a fixed set of galaxies assigned to it, the sum of the penalties for those assignments is a function $`f(x,y)`$ of the coordinates $`(x,y)`$ of the center of the disc. As long as the penalty function is convex and reasonably smooth, $`f`$ will be also. Starting with the current location $`(x_0,y_0)`$ of the disc, a simple gradient-based method (described in § 3.2) is used to find $`(x,y)`$ maximizing $`f(x,y)`$.

This gives us the essentials of the inner loop of the algorithm. It starts with a near-uniform cover of some specified (normalized) density $`1+\delta `$. It improves the cover by finding a minimum-penalty assignment of the galaxies to the discs and then moving the discs to their optimal locations given that assignment. It continues, alternately improving the assignment and then moving the discs, until the net penalty ceases to decrease appreciably. At the end of the inner loop, the algorithm finds a legal (not relaxed) assignment of galaxies to discs maximizing the number of assigned galaxies. Figure 2 shows a sequence of covers generated by a single run of the inner loop.

The outer loop performs a binary search to find the smallest $`d`$ that causes the inner loop to successfully cover the galaxies. The presentation here is a slight simplification of the actual algorithm, in that the actual algorithm uses a “polishing” heuristic before terminating the inner loop, and a heuristic is applied to reduce size of the network-flow problem before solving it. These heuristics and other details about convergence of the inner and outer loops, and starting conditions for the outer loop, are described in § 3.2.

### 3.1 Example Run of Inner Loop.

The sample instance in Figure 1 contains 12642 points — a random $`10\%`$ of the points in a subregion of the sky previously scanned. The size of this subregion is about $`10\%`$ of that of the region that will be mapped by the Survey. A uniform cover of 218 discs of capacity 60 (total capacity 13080) allows 81% of the galaxies to be assigned. After 16 iterations of the inner loop of the algorithm, the improved cover captures 97.8% of the points. Figure 1 shows the initial near-uniform cover and the final cover; Figure 2 shows a composite of the successive covers. Section 4 describes comprehensive testing of quality of solutions given by the algorithm and its running time.

### 3.2 Implementation Details.

For the initial near-uniform covers, we use Hardin, Sloane, and Smith’s catalogue of packings of points on the sphere . These packings give covers of the entire sphere, but we need a cover of only a (usually rectangular) subregion of the sky. To prune a “global” cover $`C`$ the algorithm first finds a maximum legal assignment of galaxies to discs in $`C`$, then discards all discs having at most a few assigned galaxies. (The cutoff for discarding a disc is chosen so that the resulting number of discs is as desired.)

The inner loop of the algorithm is implemented in C++ using LEDA  for basic data structures. We use a scaling algorithm by Andrew Goldberg to solve the minimum-cost flow problems . We use TCL for the outer loop of the algorithm and to collect performance data.

#### Penalty function:

The penalty for assigning a galaxy to a disc whose center is distance $`d`$ away is proportional to

$$p(d)=\{\begin{array}{cc}d^2r^2\hfill & \text{if }dr\hfill \\ 100(d^2r^2)\hfill & \text{if }dr\text{.}\hfill \end{array}$$

Recall that $`r`$ is the disc radius. When solving the relaxed problem, the algorithm first rounds the penalties. Rounding so that there are few distinct penalties allows a heuristic reduction in the size of the resulting flow problem. (This heuristic is discussed further below.) Figure 5 shows plots of $`p`$ and the rounded penalties. The rounding is chosen to preserve the distance between the galaxy and the edge of the disc within roughly a factor of 2. The edge of the disc is important because the penalty function is least smooth for points near the edge. The factor of 2 is somewhat arbitrary, it was chosen to balance between the advantages of rounding and the resulting loss of accuracy. After rounding, only 14 or so distinct penalties (each an integer power of 2) arise.

#### Reducing the size of the flow problem:

We expected the bottleneck in the algorithm to be solving the minimum-cost flow problems. To minimize this time, the algorithm uses a heuristic to reduce the minimum-cost flow problem to a smaller, approximately equivalent, problem. This is the “approximation” of the flow problem mentioned in the high-level description of the algorithm. First, the algorithm only considers assigning each galaxy to discs whose centers are within a distance of 2 disc radii, and of these at most the three closest discs. Second, it rounds the penalties as described above to reduce the number of distinct penalties. Finally, instead of having vertices for individual galaxies, it has vertices for equivalence classes of galaxies, where two galaxies are equivalent if they have the same assignable discs with the same rounded penalties. With these heuristics, even for very dense sets of galaxies, the number of equivalence classes will be proportional to the number of discs as long as each disc intersects $`O(1)`$ other discs. This is true in our case.

The precaution of using equivalence classes turned out to be unnecessary for two reasons. First, the average number of galaxies per equivalence class was typically no more than $`3`$. More fundamentally, solving the flow problems was not in fact a substantial bottleneck (see the data in § 4.2 and the subsequent discussion).111It is conceivable that the rounding of the penalties decreased the time used by the minimum-cost flow algorithm, as the latter works by scaling.

#### Constructing the flow problem:

The algorithm stores all the discs in a two-dimensional array so that discs near any given point can be found rapidly. To construct the flow problem, the algorithm iterates through the galaxies. For each galaxy, it finds the discs whose centers are within 2 disc radii. It selects the three nearest of these discs and computes the rounded distances to each. These discs and their rounded distances determine the equivalence class of the galaxy. The equivalence class is found (or created if necessary). From the equivalence classes, the flow network is constructed.

So that the equivalence classes can be found quickly, each equivalence class is stored in a hash table maintained at its nearest disc. The hash table for a disc $`D`$ contains those equivalence classes whose nearest disc is $`D`$. This method preserves locality of reference. In an earlier implementation, a single large hash table held all the equivalence classes. For large problems, this table was too large to fit in main memory. This slowed the algorithm by a factor of roughly 50.

#### Moving the discs:

After the minimum-penalty relaxed assignment is found, recall that each disc is moved individually to minimize the penalty associated with that disc. The “simple gradient-descent method” used to do this is as follows. To minimize $`f(x,y)`$, starting at a point $`(x_0,y_0)`$, compute the gradient (direction of maximum rate of increase), then move $`(x,y)`$ in steps of $`\alpha `$ (approximately $`16/1000`$ disc radii, chosen to balance speed and accuracy) in the direction opposite the gradient until such steps ceased to decrease the value of $`f(x,y)`$. Recompute the gradient at the new location and repeat the process with $`\alpha `$ halved. Continue in this fashion, halving $`\alpha `$ each time, until $`\alpha `$ is decreased to approximately $`2/1000`$ disc radii.

#### Convergence and “Polishing”:

The outermost loop of the algorithm does a binary search on the size of the uniform starting cover. Within this loop, the inner loop iteratively improves the given cover.

We describe convergence of the inner loop first. Recall that the inner loop starts with a given cover and improves it until the desired number of galaxies are legally covered, or until “convergence” occurs. Convergence is determined as follows: after each iteration, if the gap between the actual number of galaxies covered and the desired number did not decrease by at least 5%, then the algorithm considers the process “stuck”. At this point it changes the basic improvement step (this is the “polishing” heuristic mentioned in the high-level descriptions of the algorithm) as follows: it solves the relaxed problem as if the disc radius were 2% smaller. It continues with this heuristic until it also becomes stuck. Every time the process becomes stuck, the algorithm alternates between the standard improvement step and the modified one. If the process is ever stuck for at least two sequential rounds, it is considered to have converged. The purpose of the polishing heuristic is that in the original relaxed problem, a disc may be assigned galaxies that are just barely outside of it at little penalty. These galaxies cannot be legally assigned, yet may “hold” discs in place in the subsequent disc-moving step. “Shrinking” the effective radius of the disc for a few rounds encourages these galaxies to be assigned elsewhere.

Next we describe initial conditions and the convergence criterion for the outer loop. The outer loop maintains a lower bound $`L`$ and an upper bound $`U`$ on the minimum sufficient cover size. It also maintains covers $`C_L`$ and $`C_U`$ obtained by starting with a uniform cover of size $`L`$ or $`U`$ (respectively) and applying the basic algorithm to improve the cover until the desired coverage is obtained or convergence occurs. Initially $`L`$ and $`U`$ are taken to be $`1.05`$ and $`1.15`$, respectively, times the number of galaxies divided by the capacity per disc. The binary search maintains the invariant that $`C_L`$ and $`C_U`$ are, respectively, insufficient and sufficient to achieve the desired coverage. If this invariant does not hold initially, $`L`$ and/or $`U`$ are adjusted in increments of 5% to achieve the invariant. The algorithm halts the binary search as soon as the following condition ceases to be met: $`C_U`$ has more than one more disc than $`C_L`$, $`C_U`$ is at least 0.5% bigger than $`C_L`$, and $`C_U`$ legally covers at least 0.5% more galaxies than $`C_L`$. Once the search halts, the algorithm returns $`C_U`$.

## 4 Performance of the Algorithm

We tested the running time and the quality of the solutions found by the algorithm on sample instances. In this section we describe the results.

The Survey will map roughly 25% of the sky — the region having right ascension zero through $`360`$ degrees and declination $`30`$ degrees through $`90`$ degrees. Roughly one million galaxies will be mapped. Because the two phases of the Survey will be pipelined (the second will be started before the first is done), the second phase will be done in pieces.

We generated the problem instances from data from a region of the sky that had been previously scanned for a different purpose. We selected 6 subregions, and for these subregions we generated 4 problem instances by randomly sampling 30, 50, 70 or 100% of the galaxies. This gave us 24 sample problems. We took the disc radius to be 1.5 arc-seconds and the capacity to be 600 times 0.3, 0.5, 0.7, or 1 corresponding to the sampling percentage above. (The base capacity is 600 instead of 660 because approximately 60 points in each disc will be reserved for quasars not in the sample.) The largest region has an area roughly 4% of the entire sky. For each subregion, the right ascension and declination ranges and the number of galaxies are shown in Figure 6.

### 4.1 Quality of solutions.

Figure 6 illustrates the quality of the solutions returned by the algorithm on the 24 problem instances. The figure plots the size of the cover needed to assign 98% of the galaxies in each region, normalized by dividing by the number of discs needed just to provide enough capacity to hold 98% of the galaxies. The plot shows the same information for covering by near-uniform covers. The algorithm (very roughly) requires 5% to 15% extra capacity, whereas using uniform covers requires 25% to 35% extra capacity.

### 4.2 Running time.

Plots of the time per galaxy to solve each problem instance as a function of the number of galaxies appear in Figure 7. This net time includes all of the iterations needed to find the final cover for the given problem instance, including the binary search “outer loop”.

The three main components of the running time are the time building the graphs (including finding the equivalence classes of galaxies), the time solving the flow problems, and the time moving the discs. These plots show that the net running time is on the order of $`0.1`$ cpu seconds per galaxy ($`850000`$ galaxies per day), with the three main components each taking a substantial fraction of the time. These tests were carried out on a Silicon Graphics machine with 6 150 MHZ processors, a 16 Kbyte data cache, a 1 Mbyte secondary cache, and 256 Mbytes of main memory.

Most of the variation in the time per galaxy is due to the number of iterations, which varied from 30 to 100 per problem instance, and which increases (in our implementation) with the problem size. Figure 7 plots the average time per galaxy per iteration versus the problem size. Each iteration represents the solution of one relaxed problem and one perturbation of one set of discs. The time per iteration grows linearly with the number of galaxies. This is as expected, except for the surprising speed of the flow computations.

We note that the binary search is fairly naive, given that in principle a fairly precise guess about the correct size of the starting cover could be made. Similarly, we feel that a more careful and less conservative estimate of convergence, possibly interleaving the two loops in some fashion, might be warranted. These improvements might reduce the total number of iterations substantially.

#### Time to solve flow problems.

The heuristics for keeping the flow problems small appear to be effective. Figure 8 plots the average number of edges per galaxy in each flow problem as a function of the sampling density of the instance.

Figure 8 plots the average time per edge to solve the individual flow problems. The time appears to grow only near-linearly with the number of edges. Better than worst-case behavior on certain classes of problems is not uncommon; further the flow problems arising here are not particularly hard ones. See  for computational studies related to this issue.

## 5 Retrospective

The Euclidean capacitated covering problem arising here is very natural. Looking in the Operations Research literature, we found numerous capacitated covering algorithms based on integer linear programming, but these were not fast enough problems of the desired size. The Computer Science literature had a number of efficient approximation algorithms for covering that had provable worst-case performance guarantees, yet these algorithms would not produce good enough solutions in practice.

Nonetheless, our final solution rests on theoretical foundations. Our algorithm works in the spirit of Lagrangian relaxation. We decompose the problem into two parts: finding the cover $`C`$ and finding the assignment $`A`$. We relax the constraints on $`A`$ by replacing the “disc-containment” constraint by a penalty function. Then, for any given cover, finding the minimum-penalty assignment is a tractable problem. Likewise, for any given assignment, finding the minimum-penalty cover is tractable. Thus, the relaxation yields a scheme that iteratively reduces the minimum penalty and so drives the pair $`(C,A)`$ closer to feasibility. Lagrangian relaxation is a common technique in both the Operations Research  and the Computer Science  literature.

Finding the decomposition required an understanding of network flow theory. The ability to solve large problems hinges on a fast network flow algorithm. Classic augmenting paths algorithms are far too slow for our purpose. Goldberg’s algorithm incorporates both recent research within the worst-case model and heuristics discovered by empirical studies (in the spirit of ).

Also useful were Hardin, Sloane, and Smith’s sphere covers . These enabled us to start with better uniform covers than we might have otherwise. Finally, in prototyping and testing ideas, it helped to have a pre-existing library of relevant high-level data types and algorithms. For this we used LEDA .

Worst-case analysis did side-track us slightly. Although worst-case analysis suggested that network flow would be the bottleneck for large problems, it was not at all. Ironically, as described in § 3.2, our first attempt to keep the flow problems small by using equivalence classes backfired: our original implementation used a single hashing data structure to hold all the equivalence classes; although the standard worst-case model suggests hashing is quite fast, its incautious use slowed the solutions of large problems by a factor of 50 due to the lack of locality of reference.

In conclusion, our experience suggests that a successful approach rested on theoretical understanding, but required that it be creatively adapted to take advantage of the particular structure of our problem instances.

## 6 Acknowledgements

Thanks to Ken Steiglitz for introducing two of the coauthors and to an anonymous referee for helpful suggestions.


format: text

text:
# 1 Introduction

## 1 Introduction

As witnessed by the introduction of a session on the so–called left–handed materials in this conference, materials with negative phase velocity have attracted much attention during the last two years. Over three decades ago, Veselago  suggested many unusual properties of materials with negative real relative permittivity and negative real relative permeability at a certain frequency, including inverse refraction, negative radiation pressure, inverse Doppler effect. However, his considerations were completely speculative in view of the lack of a material, or even a nonhomogeneous composite medium, with a relative permittivity having a negative real part and a very small imaginary part. A breakthrough was achieved by Smith et al. , who, developing some earlier ideas by Pendry et al. , presented evidence for a weakly dissipative composite medium displaying negative values for the real parts of its effective permittivity and effective permeability. Their so–called meta–material consists of various inclusions of conducting rings and wires embedded within printed circuit boards. Other types of nanostructural combinations with similar response properties can also be devised .

Experimental results published last year by Shelby et al.  on the transmission of a $`10`$–GHz beam through a wedge provided impetus to the electromagnetics community for a discussion on the concept of a negative index of refraction. Doubts have emerged on the homogeneity and the isotropy of the composite materials used by Shelby et al. as well as on the adequacy of their measurement setup . Despite those doubts, the crucial flipping of the transmission pattern about the normal to the exit surface, when a teflon wedge was replaced by a wedge made of the ring–wire material, appears to be unexplainable in any way other than by resorting to the essence of Veselago’s suggestion. In view of the considerable literature accumulated during the past few months and the explosive nature of the current scientific scene , we take this opportunity to present our thoughts on a variety of related issues.

## 2 What’s in a name?

The emergence of a clear terminology is often a difficult process with regards to scientific findings relating to novel effects, something that is also apparent in the present instance. The first label for the candidate materials is left–handed materials . But chiral materials are important subjects of electromagnetics research and the terms left–handedness and right–handedness have been applied to the molecular structure of such materials for well over a century . The continued use of the term left–handed materials (LHMs) for achiral materials will thus confuse the crucial issues .

The term backward (BW) medium has been proffered by Lindell and colleagues . This term presumes the a priori definitions of forward and backward directions. Whatever be the merits of this term for planewave propagation, it would founder for problems involving nonplanar interfaces.

Ziolkowski and Heyman  recently provided the most extensive theoretical and numerical analysis of the negative index of refraction to date. They introduced the technical term double negative (DNG) medium to indicate that the real parts of both permittivity and permeability are negative. While sensible enough, such nomenclature conceals the importance of dissipative effects.

After a careful study of the relevant constitutive parameters, we have come to the conclusion that the term negative phase–velocity (NPV) medium is unambiguous and covers all possible situations that we could think of. It also provides a contrast to the emerging negative group–velocity (NGV) mediums, reports on which are now emerging with regularity .

## 3 The condition for NPV

Consider an isotropic dielectric–magnetic medium with relative permittivity $`ϵ_r=ϵ_r^{}+iϵ_r^{\prime \prime }`$ and relative permeability $`\mu _r=\mu _r^{}+i\mu _r^{\prime \prime }`$. Dissipation is reflected in the imaginary parts $`ϵ_r^{\prime \prime }`$ and $`\mu _r^{\prime \prime }`$, whilst causality dictates that $`\mu _r^{\prime \prime }>0`$ and $`ϵ_r^{\prime \prime }>0`$, so that $`ϵ_r`$ and $`\mu _r`$ lie in the upper half of the complex plane. The phase velocity is opposite to the direction of power flow, whenever the inequality

$$\left[+\left(ϵ_{r}^{}{}_{}{}^{2}+ϵ_{r}^{\prime \prime }{}_{}{}^{2}\right)^{1/2}ϵ_r^{}\right]\left[+\left(\mu _{r}^{}{}_{}{}^{2}+\mu _{r}^{\prime \prime }{}_{}{}^{2}\right)^{1/2}\mu _r^{}\right]>ϵ_r^{\prime \prime }\mu _r^{\prime \prime }$$
(1)
holds.  Clearly, the simultaneous satisfaction of both $`ϵ_r^{}<0`$ and $`\mu _r^{}<0`$ is a sufficient, but not necessary, requirement for the phase velocity to be negative. This result has been illustrated by frequency–domain  as well as time–domain  calculations of planewave reflection at the planar interface of free space and a NPV medium with Lorentzian characteristics.

A plane electromagnetic wave polarized parallel to the $`x`$ axis, and propagating along the $`z`$ axis in a medium characterized by $`ϵ_r`$ and $`\mu _r`$, is described by

$$𝐄(z)=A\mathrm{exp}(ik_0n_\pm z)𝐮_x,$$
(2)
where $`n_\pm =\pm \sqrt{ϵ_r\mu _r}`$. The choice of the sign of the refractive index is mandated by the direction of power flow. If the criterion (1) is satisfied, then $`n_+`$ (resp. $`n_{}`$) applies for power flow along the $`+z`$ (resp. $`z`$) axis; accordingly, $`\mathrm{Re}\left[n_+\right]<0`$ (resp. $`\mathrm{Re}\left[n_{}\right]>0`$), where $`\mathrm{Re}\left[\text{ }\text{ }\text{}\text{ }\right]`$ denotes the real part. Thus, the phase velocity is oriented parallel to the $`z`$ (resp. $`+z`$) axis. Some confusion in the literature  emerges from the claimed inadmissability of either $`n_+`$ or $`n_{}`$ on grounds other than the dictates of power flow.

## 4 Perfect lenses

Pendry  presented the possibility of fabricating a perfect (i.e., aberration–free) lens from a material with $`ϵ_r=\mu _r=1`$, which attracted enormous attention from such luminaries as science reporters attached to various newspapers. Attention came from researchers as well , . In particular, Ziolkowski and Heyman  concluded from extensive two–dimensional simulations that the condition $`ϵ_r=\mu _r=1`$ cannot be met by realistic meta–materials, even in some narrow frequency range.

Aberrations due to dissipation inside the desired NPV medium (with $`ϵ_r=\mu _r=1`$) would prove to be a stumbling block in fabricating the desired perfect lenses . However, the possible use of active (i.e., non–passive) elements in meta–materials may provide some relief from the glorious tyranny of the principle of conservation of energy . Chromatic aberrations due to non–fulfilment of the required conditions outside some narrow frequency range will also be important .

## 5 Distributed Bragg reflectors

The Bragg regime of a multilayer distributed Bragg reflector (DBR) would undergo a blue–shift, if a conventional positive phase–velocity (PPV) constituent were to be replaced by its NPV counterpart . An underlying cause may be the reversal of phase of the reflected and the transmitted plane waves, when a PPV medium is replaced by a NPV medium . Anyhow, multilayer DBRs with NPV constituents could be useful in wavelength regimes that are inaccessible with DBRs employing PPV constituents exclusively.

## 6 Homogenization

The meta–materials wherein the phase velocity can be directed opposite to the power flow are composite materials comprising inclusions of various kinds dispersed in some host medium. At sufficiently low frequencies, a homogeneous medium can be prescribed as effectively equivalent to a particulate composite medium for certain purposes –W̧00.

Incorporation of NPV mediums in such well–known homogenization approaches as the ones named after Maxwell Garnett and Bruggeman  is mathematically trivial, and we forecast many theoretical publications thereon. Whether or not the theoretical predictions of those approaches will be physically realized is another matter. We do, however, note that surprising results can emerge from the consideration of NPV mediums in homogenization approaches. For instances ,

* the Bruggeman approach forecasts that a certain mixture of a NPV medium with its prosaic PPV counterpart, both impedance–matched, can function as the medium that Pendry deems desirable for fabricating perfect lenses ; and
* the Maxwell Garnett approach predicts that composite mediums with zero permittivity and zero permeability can be made as electrically small NPV inclusions dispersed randomly but homogeneously in a PPV host medium.

## 7 Complex mediums

Oppositely directed phase velocity and power flow are distinct possibilities in complex mediums, as discussed by Lindell et al. , and could yield interesting phenomenons. As an example, the circular Bragg phenomenon will be reversed in ferrocholesteric materials with negative real permittivities and permeabilities . Even in isotropic chiral mediums, negative $`ϵ_r^{}`$ and $`\mu _r^{}`$ will lead to a reversed circular dichroism .

To conclude, this brief overview of various emerging issues relating to NPV mediums is expected to stimulate different lines of thought among the participants of Complex Mediums III. In view of claims and counter–claims launched respectively by Shelby et al.  and the detractors of their experimental results, the scientific situation is very presently volatile. Though the commonplace Lorentz model does allow the possibility of isotropic, homogeneous, dielectric–magnetic materials exhibiting negative phase velocity , our brief overview definitely does not contain the last word on the topic of artificial materials acting similarly.


text:
# Deprojection of luminosity functions of galaxies in the Coma cluster

## 1 Introduction

In the past large field studies of the photometric properties of galaxies in the Coma cluster have been based on photographic plates (e.g. Godwin et al. godwin (1983); Lugger lugger (1989)). More recent studies have been made with CCD detectors which provide higher photometric accuracy. However, a limitation of the use of CCDs has been the relatively small field of view. Consequently, these studies have been limited to very small areas, ranging from $`50`$ arcmin2 (e.g. Bernstein et al. bernstein\_etal (1995); Mobasher & Trentham mobasher\_trentham (1998)) to 1500 arcmin2 (e.g. Lobo et al. lobo\_etal (1997)), mainly focused on the cluster center. These small fields can not be used to study the dependence of the luminosity function (LF) on the projected distance from the cluster center. Furthermore, lacking the photometric data at larger distances from the core it is not possible to correct a central LF for projection effects. Accordingly, the published LFs of galaxies in the Coma cluster are all projected and are, therefore, luminosity distributions rather than actual LFs. It is important to test how well these correspond to the actual deprojected LFs.

Recently, Valotto et al. (valotto\_etal (2001)) used numerical simulations to analyse systematic effects in the determination of the galaxy luminosity function in clusters. They find that projection effects conspire to produce artificially steep faint end slopes. Beijersbergen et al. (beijersbergen\_etal (2002), hereafter BHDH) presented statistical $`U`$, $`B`$ and $`r`$ band projected LFs for galaxies in the Coma cluster in five annuli out to a projected radius of 1.4 degrees for the $`B`$ and $`r`$ bands and 0.74 degrees for the $`U`$ band. In total an area of 5.2 deg2 in $`B`$ and $`r`$ and 1.3 deg2 in $`U`$ has been surveyed. This data offers the interesting possibility to deproject actually observed LFs and to obtain true 3-d LFs.

In this *Research Note* we use the information of all annuli to deproject the 2-d LFs. Our study extends the work in BHDH by inferring for the first time deprojected LFs as a function of projected distance from the cluster center. In addition, we study how the shapes of the LFs are affected by projection by measuring the faint end slopes as a function of distance to the cluster center and compare this with results from numerical simulations. Throughout this *Research Note* the projected LFs are taken from BHDH and are called 2-d LFs. Deprojected LFs are named 3-d LFs.

## 2 Deprojection of luminosity functions

Extensive literature exists on the uncovering of three-dimensional structure from observations of projected quantities (see e.g. Zaroubi et al. zaroubi\_etal (1998) and references therein). Usually methods have been described which were designed for a specific type of observation, in particular X-ray measurements. Here we use a very simple, analytic geometric deprojection scheme (Fabian et al. fabian (1981)). This scheme has mainly been used for the analysis of X-ray observations, but in principle it is suitable for any type of observations. We use this projection scheme to deproject the observed 2-d luminosity distributions of the galaxies in the Coma cluster to real, 3-d, luminosity functions. This scheme consists of assuming that the galaxies are symmetrically distributed in a sphere. This sphere is divided up in a number of concentric shells, in which the galaxy density per magnitude is assumed to be constant. The value of the 3-d luminosity function in a particular shell can then be found by subtracting the contributions from the other shells due to projection. In Fig. 1 we show an outline of this model.

The number of projected galaxies per magnitude $`m`$, Ni,m, in annulus $`i`$, and the galaxy density per magnitude, $`\varphi _m`$, are related as:

$`N_{i,m}=\varphi _{i,m}F(R_i,R_{i1},R_i)+{\displaystyle \underset{j=i+1}{\overset{5}{}}}\varphi _{j,m}`$ (1)
$`\left\{[F(R_j,R_{i1},R_i)F(R_{j1},R_{i1},R_i)]\right\}`$
where Ri are the radii of the annuli and the function F is defined by:

$$F(\alpha ,\beta ,\gamma )\frac{4}{3}\pi \alpha ^3\left[(1\beta ^2/\alpha ^2)^{3/2}(1\gamma ^2/\alpha ^2)^{3/2}\right]$$
(2)
With the available 2-d information we can solve for all $`\varphi _{i,m}`$, i.e. the true 3-d LFs. The $`N_{i,m}`$ are dependent on the areas of the annuli which are accurately known. However, the outermost annuli are not entirely filled by data: annulus V is filled for $`80`$% in $`B`$ and $`r`$ and annulus IV is filled for $`45`$% in $`U`$. We omitted annulus IV for the $`U`$ band and corrected $`N_{5,m}`$ for the $`B`$ and $`r`$ bands by assuming a constant projected density throughout this annulus. The function F is solely dependent on geometry and has values of order $`10^3`$ to $`10^1`$. Figs. 2 and 3 show both the 2-d (open symbols) and 3-d (solid symbols) LFs for all bands and annuli (all magnitudes are corrected for atmospheric extinction). The error bars of the 2-d LFs are not shown for clarity, but are somewhat smaller than the 3-d error bars. The 2-d errors are estimates of the uncertainty in the subtraction of the number of field galaxies. These estimates are based on very large control fields and therefore very accurate (see BHDH for an extensive discussion). The 3-d errors have been calculated using Eq. 1. The corresponding errors are therefore dependent on each other and fully determined by the uncertainty in the subtraction of the number of field galaxies. The vertical scales represent volume densities and the 2-d LFs are scaled as to match the bright ends of the 3-d LFs. This simplifies the comparison of their shapes. In annulus V the galaxy density is much lower than in the center, but still higher than in the field. That is to say, the derived 3-d $`B`$ and $`r`$ band LFs still suffer from (minor) projection effects due to cluster galaxies at radii even larger than 1.7$`h_{100}^1`$ Mpc. The contamination by projection for the 3-d $`U`$ band LFs is larger since there we miss photometric information for radii larger than 0.5$`h_{100}^1`$ Mpc.

## 3 Discussion

In the literature LFs of galaxies in clusters are usually given in units of number of galaxies per magnitude per area. What actually is measured is a luminosity distribution, denoted here as 2-d LF. A luminosity *function* of a sample of galaxies has units of galaxies per magnitude per volume. For observations of the entire projected volume of a cluster the two have the same shape since for all luminosities the same volume is sampled. BHDH studied changes in the LF as a function of distance to the cluster center. However, they did not correct for the fact that fore- and background galaxies, residing in the outer parts of the cluster, are projected onto the central cluster regions. The question is how much their results have been influenced by this effect. Below we study how the shapes of the luminosity distributions given in BHDH are affected by projection effects.

In the outermost annulus the projected area is converted to volume through a factor depending on the adopted geometry (radii), i.e. for this annulus the shapes of the 2-d and 3-d LFs remain identical. From all LFs of the inner annuli the contributions to the galaxy counts from surrounding shells is subtracted using the model described in Sect. 2. The impact on the resulting LFs can be seen in Figs. 2 and 3. When the 2-d and 3-d LFs are compared binwise the differences, generally, are small. Yet, the error bars become larger and it should be kept in mind that the scale is logarithmic. Our result can be explained by the small values of the F integrals in Eq. 1 combined with the density profile of the cluster. In all expressions we find combinations of relatively high values of the integrals times low densities, and vice versa, resulting in a minor net correction term. However, this correction term is different for each magnitude bin causing the shapes of corrected and uncorrected LFs to be different. Below, we quantify the differences by studying the change in the shapes of the 3-d LFs as a function of position in the cluster. The results are compared to the exact same analysis for the 2-d LFs.

### 3.1 Dependence of luminosity functions on radial distance from the cluster center

BHDH found that the 2-d central LFs could be well represented by Schechter functions. In Fig. 2 the best-fitting Schechter (schechter (1976)) functions of the 3-d LFs are shown as solid lines. In Table 1 we compare the best-fitting Schechter parameters, $`\alpha `$ and $`M^{}`$, for the 2-d and 3-d case. For the $`r`$ and $`B`$ bands the Schechter fits are again reasonable representations of the data, though the best-fitting parameters have changed. The $`U`$ band fit remains a poor representation of the data.

One of the main results of BHDH was the finding of a significant steepening of the faint end slopes with increasing cluster radius. A straight line fit ($`sr+b`$), where $`r`$ is the clustercentric radius, to the weighted $`U`$, $`B`$ and $`r`$ band slopes at each radius gives $`s=0.11\pm 0.02`$ Mpc-1 for the 2-d LFs. In order to compare this with the situation where all 2-d LFs have been corrected for projection we fitted power laws ($`b10^{aM}`$) to the faint ends of the 3-d LFs in Figs. 2 and 3. The infalling group around NGC 4839 (Neumann et al. neumann\_etal (2001)) does not have a significant influence on any of the derived LFs. The LFs were fitted for $`M_\mathrm{r}>20`$, $`M_\mathrm{B}>20`$ and $`M_\mathrm{U}>19`$. We omitted the slope of the $`B`$ band LF for area IV since there the LF is not well constrained. The results are shown in Fig. 4. A straight line fit to the weighted slopes of the 3-d LFs at each radius gives $`s=0.07\pm 0.04`$ Mpc-1 and is indicated by the dashed line. Although this is a somewhat weaker trend than derived from the 2-d results, it is different at $`1\sigma `$ only and the trend is still clearly there.

Valotto et al. (valotto\_etal (2001)) used a deep mock catalogue to analyse systematic effects in the determination of the galaxy luminosity function in clusters. They find that projection effects produce artificially steep faint end slopes. The cluster core suffers most from objects that are projected onto it. The LFs of area I shows similar ($`r`$ and $`B`$ band) or shallower ($`U`$ band) faint end slopes when they are corrected for projection. We conclude that the impact of projection effects on the shapes of the 2-d LFs in Coma is not as severe as their simulations suggest. However, the 2-d LFs derived for relatively small areas suffer from projection effects and should therefore be properly corrected before one can reliably interpret their shapes.

## 4 Summary and Conclusions

In this *Research Note* we have used a simple analytic model to deproject 2-d LFs of cluster galaxies. We have applied this to the LFs of the galaxies in the Coma cluster measured by BHDH. Under the assumptions of a spherically symmetric galaxy distribution and a constant galaxy density in each shell the 2-d LFs have been corrected for projection effects. The main results are summarized by Figs. 2 and 3. Projection effects affect the shapes of the 2-d LFs and in general make the faint end slopes artificially steeper, although the differences are small. The alterations are not significant enough as to arrive at different conclusions compared to those inferred from 2-d LFs. Furthermore, the previously derived steepening of the faint end slopes with increasing cluster radius is only slightly weakened when projection effects are taken into account. Overall, we have demonstrated by using a simple model that before interpreting shapes of 2-d LFs one needs to correct for projection.

###### Acknowledgements.

We thank the referee for constructive comments.


text:
# Instability of the Gravitational 𝑁-Body Problem in the Large-𝑁 Limit

## 1 Introduction

Miller (1964) first noted the remarkable sensitivity of the gravitational $`N`$-body problem to small changes in the initial conditions. Errors or perturbations in the coordinates or velocities of one or more stars grow roughly exponentially, with an $`e`$-folding time that is of order the crossing time. The implication, verified in a number of subsequent studies (Lecar, 1968; Hayli, 1970), is that $`N`$-body integrations are not reproducible over time scales that exceed a few crossing times. The instability is somewhat reduced when the Newtonian force law is modified by a cutoff (Standish, 1968), indicating that it is driven by close encounters.

Of interest is the behavior of Miller’s instability in the limit of large $`N`$. It is commonly assumed that the $`N`$-body equations of motion go over to the collisionless Boltzmann equation (CBE) as $`N\mathrm{}`$ (Binney & Tremaine, 1987). This would imply, for instance, that a particle trajectory which is integrable in a smooth potential should exhibit increasingly regular behavior as the number of point masses used to represent the smooth potential increases. On the other hand, if the rate of growth of small perturbations remains substantial even for large $`N`$, there would be an important sense in which the CBE does not correctly describe the behavior of $`N`$-body systems.

In fact there are indications that the growth rate of Miller’s instability remains constant or even increases with $`N`$ (Kandrup & Smith, 1991; Heggie, 1991; Goodman Heggie & Hut, 1993), although this result is uncertain since published numerical experiments have been limited to $`N\begin{array}{c}<\hfill \\ \hfill \end{array}10^3`$. Here we describe the application of a new, “systolic” $`N`$-body algorithm to this problem, which allows us to treat systems with $`N`$ as large as $`10^5`$. We observe for the first time a clear $`N`$-dependence of the instability: the growth rate is found to increase approximately as $`\mathrm{ln}N`$. Our methods and results are described in §2 and §3, and the implications for galactic dynamics are discussed in §4.

## 2 Method

Following Miller (1971), we integrated the coupled $`N`$-body and variational equations:

$`\ddot{𝐱}_i`$ $`=`$ $`Gm{\displaystyle \underset{j=1}{\overset{N}{}}}{\displaystyle \frac{𝐱_i𝐱_𝐣}{|𝐱_i𝐱_j|^3}}`$ (1a)
$`\ddot{𝐗}_i`$ $`=`$ $`Gm{\displaystyle \underset{j=1}{\overset{N}{}}}\left[\left(𝐗_i𝐗_j\right)3{\displaystyle \frac{\left(𝐗_i𝐗_j\right)\left(𝐱_i𝐱_j\right)}{|𝐱_i𝐱_j|^2}}\left(𝐱_i𝐱_j\right)\right]\times {\displaystyle \frac{1}{|𝐱_i𝐱_j|^3}}.`$ (1b)
Here $`𝐱_i`$ are the configuration-space coordinates of the $`i`$th particle and $`𝐗_i`$ are the components of its variational vector. The masses $`m`$ are assumed equal. The variational equations represent the time development of the infinitesimal distance between two neighboring $`N`$-body systems.

Equations (1a) and (1b) were integrated using the systolic $`N`$-body algorithm described by Dorband, Hemsendorf & Merritt (2002). This algorithm implements the fourth-order Hermite integration as described by Makino & Aarseth (1992). We adopted their formula for computing the time step of particle $`i`$,

$$\mathrm{\Delta }t_i=\sqrt{\eta \frac{|𝐚_i(t_1)||𝐚_i^{(2)}(t_1)|+|\dot{𝐚}_i(t_1)|^2}{|\dot{𝐚}_i(t_1)||𝐚_i^{(3)}(t_1)|+|𝐚_i^{(2)}(t_1)|^2}}.$$
(2)
Here $`𝐚`$ is the acceleration $`\ddot{𝐱}`$, superscripts denote the order of the time derivative, $`t_1`$ is the system time, and $`\eta `$ is a dimensionless constant; we set $`\eta =0.02`$. The same time step was used to integrate both the $`N`$-body and variational equations. The systolic algorithm distributes the particles equally among $`p`$ processors and computes forces by systematically shifting the particle coordinates between processors in a ring. A single processor was used for small particle numbers while 64 processors were used for the largest-$`N`$ integrations (Table 1). The multi-processor integrations were carried out using the Cray T3E at the Höchstleistungsrechenzentrum in Stuttgart.

Initial conditions were generated randomly from the isotropic Plummer model, whose density and potential satisfy

$$\rho (r)=\frac{3GM}{4\pi }\frac{b^2}{\left(r^2+b^2\right)^{5/2}},\mathrm{\Phi }(r)=\frac{GM}{\sqrt{r^2+b^2}}.$$
(3)
We adopted standard $`N`$-body units (Heggie & Mathieu, 1986) such that $`G=M=1`$, $`E=1/4`$, giving a scale factor $`b=3\pi /16`$. We defined the crossing time $`t_{cr}`$ as $`R/V`$, with $`RGM^2/2E`$ and $`V^22E/M`$; in these units, $`t_{cr}=2\sqrt{2}`$. The variational vectors $`𝐗`$ and $`𝐕`$ were assigned an initial amplitude of $`10^{30}`$ for each particle with randomly chosen directions.

The parameters of the integrations are listed in Table 1. Each integration was continued until a time of 20, or roughly $`7`$ crossing times, with the exception of the largest run, $`N=131072`$ which was terminated at $`t=12`$. This time is short enough that two-body relaxation should not be important except perhaps for the smallest $`N`$, and long enough to show a clearly exponential growth of the solutions of the variational equations.

We focussed on the amplitudes of the Cartesian components of the variational vectors, $`𝐗_i=(X_i,Y_i,Z_i)`$, since the variational velocities tend to exhibit spikes in their time dependence (Miller, 1964). We examined three choices for the amplitude $`\mathrm{\Delta }`$ of the separation:

1. $`\mathrm{\Delta }_m`$, the maximum over $`i`$ of the $`|𝐗_i|`$.

2. $`\mathrm{\Delta }_a`$, the arithmetic mean of the $`|𝐗_i|`$.

3. $`\mathrm{\Delta }_g`$, the geometric mean of the $`|𝐗_i|`$.

The instability growth rate $`\mu _e`$ was defined as

$$\mu _e=\frac{\mathrm{ln}\mathrm{\Delta }(t_2)\mathrm{ln}\mathrm{\Delta }(t_1)}{t_2t_1}.$$
(4)
Except in the case of small $`N`$, $`N\begin{array}{c}<\hfill \\ \hfill \end{array}1024`$, growth in the variation was found to be nearly exponential and hence the computed values of $`\mu _e`$ depended only weakly on $`t_2`$ and $`t_1`$. We chose $`t_1=1`$ and $`t_2=20`$, except for the $`N=131072`$ run for which $`t_2=12`$.

## 3 Results

The three amplitudes $`\mathrm{\Delta }`$ defined above were found to be very similar in most of the integrations, as shown in Figure 1 for the 10 integrations with $`N=32768`$. In what follows we adopt $`\mathrm{\Delta }\mathrm{\Delta }_g`$, the geometric mean.

Figure 2a shows $`t_e\mu _e^1`$ for each of the integrations. The mean value of the $`e`$-folding time and its uncertainty are plotted in Figure 2b; the latter was defined as the standard error of the mean, or $`n^{1/2}`$ times the standard deviation, with $`n`$ the number of distinct $`N`$-body integrations. For the first time, a clear $`N`$-dependence can be seen, in the sense that the average value of $`t_e`$ declines with increasing $`N`$: large-$`N`$ systems are more unstable than small-$`N`$ systems.

A number of predictions have been made for the large-$`N`$ dependence of the $`e`$-folding time. Gurzadyan & Savvidy (1986) estimated $`t_eN^{1/3}t_{cr}`$ based on a geometrical approach. This $`N`$ dependence is clearly inconsistent with Figure 2. Gurzadyan & Savvidy’s approach was criticized already by Heggie (1991) and Goodman, Heggie & Hut (1993) due to its improper treatment of close encounters. The latter authors argued that $`t_e/t_{cr}1/\mathrm{ln}N`$ or $`1/\mathrm{ln}(\mathrm{ln}N)`$; the weaker dependence would hold only after a time long enough that the perturbation from one star was able to propagate through the system.

We tested these predictions against the $`N`$-body data. Figure 3 shows fits of two functional forms to the mean growth rates:

$`\mu _e`$ $`=`$ $`a+b\mathrm{ln}N,`$ (5a)
$`\mu _e`$ $`=`$ $`c+d\mathrm{ln}(\mathrm{ln}N).`$ (5b)
Since any such relation is expected to be valid only in the limit of large $`N`$, we restricted the fits to $`N1024`$. The best-fit parameters are given in Table 2. We used a standard least-squares routine that accounts for errors in the dependent variable ($`\mu _e`$); in the case of the data pointa with $`N=65536`$ and 131072, for which there was only one integration, the uncertainty in $`\mu _e`$ was assumed to be the same as in the integration with $`N=32768`$. These two data points were omitted when computing the values of $`\stackrel{~}{\chi }^2`$ given in Table 2.

While $`\mu _et_{cr}N^0`$ dependence can clearly be ruled out, we can not distinguish between a $`\mathrm{ln}N`$ and a $`\mathrm{ln}(\mathrm{ln}N)`$ dependence. Expressed in terms of $`t_e`$ and $`t_{cr}`$, the best-fit relations are

$`t_e/t_{cr}`$ $``$ $`{\displaystyle \frac{0.99}{\mathrm{ln}\left(1.09\times 10^3N\right)}},`$ (6a)
$`t_e/t_{cr}`$ $``$ $`{\displaystyle \frac{0.112}{\mathrm{ln}\left(0.70\mathrm{ln}N\right)}}.`$ (6b)
Distinguishing between these two functional forms, or other similar ones, would clearly be very difficult; even for $`N=10^6`$, the two relations predict values of $`t_e`$ that differ only by $`5\%`$.

## 4 Discussion

Our results confirm that the gravitational $`N`$-body problem is inherently chaotic and furthermore that the degree of chaos, as measured by the rate of divergence of nearby trajectories, increases with increasing $`N`$. We have directly established a characteristic $`e`$-folding time of $`t_{cr}/20`$ for systems with $`N10^5`$; if the weak $`N`$-dependence found here can be extrapolated, the characteristic time for systems containing $`10^{12}`$ particles would be only slightly smaller, $`t_{cr}/30`$. Hence trajectories in stellar and galactic systems diverge on a time scale that is generically much shorter than the crossing time.

If the collisionless Boltzmann equation (CBE) is a valid representation of large-$`N`$ stellar systems, it should be possible to show that the $`N`$-body trajectories go over, in the limit of large $`N`$, to the orbits in the corresponding smoothed-out potential. The generic instability of the $`N`$-body problem precludes this, since the characteristics of the CBE can not be identified with the $`N`$-body orbits for times longer than $`t_{cr}`$. At the same time, much experience with $`N`$-body integrations demonstrates that in many ways the behavior of large-$`N`$ systems matches expectations derived from the CBE (e.g. Aarseth & Lecar 1975).

A likely resolution of this seeming paradox is that the macroscopic, or finite-amplitude, behavior of trajectories is not well predicted by the rate of growth of small perturbations. For instance, orbits integrated in “frozen” $`N`$-body potentials behave more and more like their smooth-potential counterparts as $`N`$ is increased, even though their Liapunov exponents remain large (Valluri & Merritt, 2000; Kandrup & Sideris, 2001). The initial growth of perturbations is exponential but it saturates, at an amplitude that varies inversely with $`N`$ (Valluri & Merritt, 2000; Hut & Heggie, 2001). Thus the CBE may be a good predictor of the macroscopic dynamics of large-$`N`$ systems even if it does not reproduce the small-scale chaos inherently associated with $`N`$-body dynamics.

This work was supported by NSF grant 00-71099 and by NASA grants NAG5-6037 and NAG5-9046 to DM. We thank S. Aarseth, J. Goodman, D. Heggie and H. Kandrup for useful discussions. We are grateful to the Höchstleistungsrechenzentrum in Stuttgart for their generous allocation of computer time.


text:
# First-order transitions for 𝑛-vector models in two and more dimensions; rigorous proof.

Abstract: We prove that various SO(n)-invariant $`n`$-vector models with interactions which have a deep and narrow enough minimum have a first-order transition in the temperature. The result holds in dimension two or more, and is independent on the nature of the low-temperature phase.

Recently Blöte, Guo and Hilhorst , extending earlier work by Domany, Schick and Swendsen  on 2-dimensional classical XY-models, performed a numerical study of 2-dimensional $`n`$-vector models with non-linear interactions. For sufficiently strong values of the non-linearity, they found the presence of a first-order transition in temperature. In  a heuristic explanation of this first-order behavior, based on a similarity with the high-$`q`$ Potts model, was suggested, explaining the numerical results. A further confirmation of this transition was found by Caracciolo and Pellisetto , who considered the $`n\mathrm{}`$ (spherical limit) of the model, and found the same first-order transition.

On the other hand, various studies, mostly based on Renormalization-Group analyses or Kosterlitz-Thouless type arguments based on the picture of binding-unbinding of vortices, have contested this first-order behaviour (e.g. ).

Here we settle the issue by presenting a rigorous proof of the existence of this first-order transition.

It may seem somewhat surprising that two-dimensional $`n`$-vector models, whose magnetization by the Mermin-Wagner theorem  is always zero, can have such a phase transition. The reason is that the transition we are talking about here is manifested by the long-range order in higher-order correlation functions. Such transitions were discovered by one of us some time ago, see . But the results of  were related to the fact that there the symmetry group was the (disconnected) group $`O\left(2\right),`$ and at the transition point the discrete symmetry $`_2`$ is broken, while the connected part – $`SO\left(2\right)`$ – of the symmetry persists. The nature of the transition we study here, however, is not connected to any type of symmetry breaking, and as such is much closer to the high-$`q`$ $`q`$-state Potts model, or the model studied in , where a first order transition between a low-energy and a high-entropy phase occurs. Thus we confirm the original intuition of .

For XY-spins in two dimensions there can be a low-temperature phase with slow polynomial decay of correlations, while the majority belief in the field, despite the work of Patrascioiu and Seiler , is that for $`n>2`$, the $`n`$-vector models at low finite temperatures have exponentially decaying correlations, just as at high temperatures. Our result unfortunately does not say anything about this question.

Our proof is directly inspired by the existing proofs for low-energy–high-entropy phase transitions, and is indeed an adaptation of those. We employ the method of Reflection Positivity  (RP). For simplicity we write the proof for 2-dimensional XY-spins, the extension to the general case is immediate.

We remark that also the generalization to higher dimensions is immediate. Thus the low-temperature phase can be either magnetized, Kosterlitz-Thouless-like (not magnetized with slow correlation decay), or possibly non-magnetized with exponentially decaying correlations. As such our result contradicts strong “universality” claims, stating that universality classes of interactions exist, all elements of which have the same kind of phase transition between a high-temperature phase and a low-temperature phase, and which are only determined by the dimension of the system, the symmetry of the interaction and whether the interaction is short-range or long-range.

We will consider the Hamiltonian given by

$$H=J\underset{i,j^2}{}\left(\frac{1+\mathrm{cos}\left(\varphi _i\varphi _j\right)}{2}\right)^p.$$
(1)
To formulate our result we have to introduce for every n.n. bond $`b=(i,j)`$ the following bond observables:

$$P_b^<(\varphi _i,\varphi _j)=\{\begin{array}{cc}1\hfill & \text{ if }\left|\varphi _i\varphi _j\right|<\epsilon /2,\hfill \\ 0\hfill & \text{ if }\left|\varphi _i\varphi _j\right|>\epsilon /2,\hfill \end{array}$$
(2)
which project on the ordered bond configurations, and $`P_b^>(\varphi _i,\varphi _j)=1P_b^<(\varphi _i,\varphi _j).`$ Our main result is contained in the following

###### Theorem 1

Suppose the parameter $`p`$ is large enough. Then there exists a transition temperature $`\beta _c=\beta _c(J,p),`$ such that there are two different Gibbs states, $`^<`$ and $`^>,`$ at $`\beta =\beta _c,`$ corresponding to the Hamiltonian (1). For some specific choice of $`\epsilon =\epsilon \left(p\right)`$ in (2), we have for the “ordered” state $`^<`$ that

$$P_b^<^<>\kappa \left(p\right),$$

while in the “disordered” phase $`^>`$

$$P_b^>^>>\kappa \left(p\right),$$

for each bond $`b,`$ with $`\kappa \left(p\right)1`$ as $`p\mathrm{}.`$

Before analysing the model (1), we present an even simpler toy model, which already displays the mechanism, and which is even closer to the Potts model. The single spin space is the circle, $`S^1,`$ the free measure is the Lebesgue measure, normalized such that $`S^1`$ has measure one, so $`S^1=[\frac{1}{2},\frac{1}{2}]`$. ( One can take here any sphere $`𝕊^n`$ instead.) The toy Hamiltonian is

$$H=J\underset{i,j^2}{}U(\varphi _i,\varphi _j)$$
(3)
The nearest neighbour interaction $`U(\varphi _1,\varphi _2)=U\left(\left|\varphi _1\varphi _2\right|\right)`$ is given by

$$U\left(\varphi \right)=\{\begin{array}{cc}1\hfill & \text{ if }\left|\varphi \right|\frac{\epsilon }{2},\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$

Here $`\epsilon `$ plays a similar role to $`\frac{1}{q}`$ in the $`q`$-state Potts model. The Hamiltonian is RP under reflections in coordinate planes. For the case $`^2`$ one has also RP under reflections in lines at 45 degrees, passing through the lattice sites. That case is the easiest.

One has to show that:

* $`P_b^>_\beta `$ is small for large $`\beta ;`$ (the ordered, typical low-temperature phase bonds)
* $`P_b^<_\beta `$ is small for small $`\beta ;`$ (the disordered, typical high-temperature phase bonds)
* $`P_b^{}^<P_{b^{\prime \prime }}^>_\beta `$ is small for all $`\beta ,`$

provided $`\epsilon `$ is small enough. Here $`_\beta `$ is the state with periodic b.c. in the box $`\mathrm{\Lambda }`$ of size $`L,`$ and $`b^{},b^{\prime \prime }`$ are two orthogonal bonds sharing the same site. The estimates have to be uniform in $`L,`$ for $`L`$ large. The first two are straightforward application of RP and the chess-board estimate. So let us get the last one. By the chess-board estimate,

$$P_b^{}^<P_{b^{\prime \prime }}^>_\beta P^\mathrm{\#}_\beta ^{1/\left|\mathrm{\Lambda }\right|},$$
(4)
where the observable $`P^\mathrm{\#}`$ is the indicator

$$P^\mathrm{\#}=\underset{bE_{01}}{}P_b^<\underset{bE_{23}}{}P_b^>.$$

Here $`E_{01},E_{23}`$ is the partition of all the bonds in $`\mathrm{\Lambda }`$ into two halves; $`E_{01}`$ consists of all bonds $`(x,x+e_1)`$ and $`(x,x+e_2),`$ for which $`x^1+x^2=0`$ or $`1\mathrm{mod}\mathrm{\hspace{0.17em}4},`$ while $`E_{23}`$ is the other half; $`e_1`$ and $`e_2`$ are the two coordinate vectors.

To proceed with the estimate (4) we need the estimate on the partition function. We have

$$Z_\mathrm{\Lambda }(\beta ,\epsilon )\left(\frac{\epsilon }{2}\right)^{\left|\mathrm{\Lambda }\right|}e^{2\beta \left|\mathrm{\Lambda }\right|}+\left(14\epsilon \right)^{\left|\mathrm{\Lambda }\right|/2}.$$
(5)
(The first summand is obtained by integrating over all configurations $`\varphi ,`$ such that $`\left|\varphi _x\right|\frac{\epsilon }{4}`$ for all $`x\mathrm{\Lambda }.`$ For the second one we take all configurations $`\varphi `$ which are arbitrary on the even sublattice and which satisfy $`\left|\varphi _x\varphi _y\right|>\epsilon /2`$ for every pair of n.n.; for every $`y`$ on the odd sublattice that leaves the spins to be free in a set of measure $`14\epsilon .`$) Solving

$$\left(\frac{\epsilon }{2}\right)^{\left|\mathrm{\Lambda }\right|}e^{2\beta _0\left|\mathrm{\Lambda }\right|}=\left(14\epsilon \right)^{\left|\mathrm{\Lambda }\right|/2}$$

for $`\beta _0`$ we find

$$e^{4\beta _0}=\left(14\epsilon \right)\left(\frac{\epsilon }{2}\right)^2,$$
(6)
so for $`\beta \beta _0`$ the first term in (5) dominates, while for $`\beta \beta _0`$ the second term dominates. Similarly, the partition function $`Z_\mathrm{\Lambda }^\mathrm{\#}(\beta ,\epsilon ),`$ taken over all configurations $`\varphi `$ with $`P^\mathrm{\#}\left(\varphi \right)=1`$ satisfies

$$Z_\mathrm{\Lambda }^\mathrm{\#}(\beta ,\epsilon )e^{\beta \left|\mathrm{\Lambda }\right|}\epsilon ^{\frac{3}{4}\left|\mathrm{\Lambda }\right|+O\left(\sqrt{\left|\mathrm{\Lambda }\right|}\right)}.$$
(7)
If $`\beta \beta _0,`$ we write, using (6):

$$P_b^{}^<P_{b^{\prime \prime }}^>_\beta \frac{e^\beta \epsilon ^{\frac{3}{4}}}{\frac{\epsilon }{2}e^{2\beta }}=2\frac{1}{\epsilon ^{\frac{1}{4}}e^\beta }2\frac{1}{\left[\epsilon \left(14\epsilon \right)\left(\frac{\epsilon }{2}\right)^2\right]^{\frac{1}{4}}}C\epsilon ^{1/4}.$$

If $`\beta \beta _0,`$ we similarly have

$$P_b^{}^<P_{b^{\prime \prime }}^>_\beta \frac{e^\beta \epsilon ^{\frac{3}{4}}}{\left(14\epsilon \right)^{1/2}}\frac{\left[\left(14\epsilon \right)\left(\frac{\epsilon }{2}\right)^2\epsilon ^3\right]^{\frac{1}{4}}}{\left(14\epsilon \right)^{1/2}}C^{}\epsilon ^{1/4}.$$

So we are done.

For the non-linear models, we employ the fact that for small difference angles $`\mathrm{cos}(\varphi _i\varphi _j)`$ is approximately $`1O((\varphi _i\varphi _j)^2)`$ and furthermore that $`lim_p\mathrm{}(1\frac{1}{p})^p=e^1`$. This suggests to choose $`ϵ(p)=\frac{1}{\sqrt{p}}`$. Because the separation between ordered and disordered bonds is somewhat arbitrary, to obtain an inequality similar to (5), we make a slightly different choice. We consider a bond $`(i,j)`$ disordered if $`|\varphi _i\varphi _j|\frac{C}{\sqrt{p}}`$ for some large $`C`$. So first we choose a sufficiently large constant $`C`$. For the estimate of the ordered partition function we only integrate over the much smaller intervals of “strongly ordered” configurations: $`|\varphi _i|\frac{C^1}{\sqrt{p}}`$ to obtain a lower bound:

$$Z_\mathrm{\Lambda }(\beta ,p)\left(\frac{1}{2C\sqrt{p}}\right)^{\left|\mathrm{\Lambda }\right|}e^{[2\beta (1O(\frac{1}{C}))]\left|\mathrm{\Lambda }\right|}+\left(1\frac{4C}{\sqrt{p}}\right)^{\left|\mathrm{\Lambda }\right|/2}.$$
(8)
This makes use of the fact that the “strongly ordered” bonds all have energy almost equal to $`J`$, whereas the disordered partition function is bounded by that of the toy model, but with $`\epsilon `$ replaced by $`\frac{C}{\sqrt{p}}`$.

For the estimate which shows that ordered and disordered bonds tend not to neighbor each other, we obtain:

$$Z_\mathrm{\Lambda }^\mathrm{\#}(\beta ,p)e^{\beta \left|\mathrm{\Lambda }\right|}(\frac{C}{\sqrt{p}})^{(\frac{3}{4}+O(e^C))\left|\mathrm{\Lambda }\right|+O\left(\sqrt{\left|\mathrm{\Lambda }\right|}\right)}.$$
(9)
The rest of the argument is essentially unchanged. We first chose C big enough (such that $`\frac{1}{C}`$ is small wrt to $`1`$), and we can still choose $`p`$ big enough such that $`p`$ is large wrt to $`C^2`$, which finishes the argument.


text:
# Spin-dependent properties of a two-dimensional electron gas with ferromagnetic gates

## Abstract

A theoretical prediction of the spin-dependent electron self-energy and in-plane transport of a two-dimensional electron gas in proximity with a ferromagnetic gate is presented. The application of the predicted spin-dependent properties is illustrated by the proposal of a device configuration with two neighboring ferromagnetic gates which produces a magnetoresistance effect on the channel current generated by nonmagnetic source and drain contacts. Specific results are shown for a silicon inversion layer with iron gates. The gate leakage current is found to be beneficial to the spin effects.

The generation, control, and measurement of electron spin polarization are the essential ingredients of spintronics. The celebrated giant magnetoresistance effect in metallic multilayers has lead to the realization of remarkable non-volatile magnetic random access memories. A great deal of interest is presently focused on semiconductor spintronics because of the potentially appealing nonlinear functionalities offered by semiconductors review . One particularly relevant system is that with planar confined carriers whose density can be controlled by a gate. Datta and Das Datta  proposed a spintronics device consisting of a gated two-dimensional electron gas (2DEG) with magnetic source and drain contacts. The spin valve behavior of the Datta and Das proposal depends on the orientation of the magnetization in the source and drain contacts and is modulated by the Rashba spin-orbit effect rashba  induced by the gate. Much recent progress in spin injection review  and control of spin-orbit interaction nitta  has been made towards the goal of an eventual demonstration of the Datta and Das device.

In this letter, we investigate the spin-dependence of the energies and lifetimes of a two-dimensional electron gas where the ferromagnetic elements are the gates, while the source and drain are nonmagnetic contacts. The spin-dependence of the electronic transport arises out of the proximity to the ferromagnetic gates, in stark contrast to the injection of spin-polarized electrons from the ferromagnetic source. Recent experiments FMimprinting ; FMcoherence  have shown that large spontaneous spin polarization can be achieved in an optically excited semiconductor layer in contact with a ferromagnet . This phenomenon is explained in terms of spin-dependent reflection at the semiconductor/ferromagnet interface imprinting . For quantum confined semiconductor electrons, the proximity effect with the ferromagnet is equivalent to multiple reflections. A required condition is that the interface barrier is sufficiently thin to allow for sizeable wave-function coupling. This is analogous to the proximity effect of a metal such as Pd next to Fe ivan . Our scheme avoids the restrictions of the Rashba effect, namely narrow gap semiconductors (for large spin-orbit coupling) and one-dimensional channels (to maximize the favorable directions of spin torque). There is a wider range of appropriate systems for our scheme, such as the inversion layer in Metal-Oxide-Semiconductor (MOS) structures, the 2DEG in a modulation doped III-V quantum well, and the InAs accumulation or inversion layers. The paper is structured as follows. We shall discuss (i) the spin-dependent properties of a 2DEG with a single ferromagnetic gate, (ii) a device proposal of a spin valve with two neighboring gates, (iii) model results for a silicon inversion layer with iron gates.

(i) 2DEG with a single ferromagnetic gate
The tunneling coupling of the electron wave function in the 2DEG to the exchange-split band states in the metal produces a spin-dependent broadening and shift of the quantized electron energy in the 2DEG. Namely, we have a complex energy $`E_\pm ^\mathrm{s}=\mathrm{Re}(E_\pm ^\mathrm{s})i\mathrm{}/(2\tau _\pm )`$ where $`E_+^\mathrm{s}`$ ($`E_{}^\mathrm{s}`$) refers to the majority (minority) spin channel. Note that $`\tau _\pm `$ describes the escape time of the quasi-bound semiconductor electrons into the ferromagnet band continuum. Due to the spin-dependent self-energy, the in-plane transport of the 2DEG is described by the spin-dependent conductance $`\sigma _\pm `$ whose expression is

$$\frac{1}{\sigma _\pm }=\frac{m_x}{n_\pm e^2}\left(\frac{1}{\tau _0}+\frac{1}{\tau _\pm }\right),$$
(1)
where $`\tau _0`$ is the ordinary spin-independent Drude scattering time, $`m_x`$ the in-plane effective mass and $`n_\pm `$ the spin-dependent density of the 2DEG due to the spin splitting $`\mathrm{Re}(E_+^\mathrm{s}E_{}^\mathrm{s})`$. However, this is not the only aspect of the problem, because a fraction of the channel current is leaked into the gate (implying that drain and source currents are different). The spin-dependent leakage probabilities $`\alpha _\pm `$ for non-equilibrium electrons injected into the channel are governed by the transit times $`\tau _{\mathrm{tr}}^\pm =L_\mathrm{G}/v_x^\pm `$ where $`L_\mathrm{G}`$ is the length of the gate (along the $`x`$-direction) and $`v_x^\pm `$ are the spin-dependent drift velocities induced by the source-drain bias . Namely,

$$\alpha _\pm =_0^{\tau _{\mathrm{tr}}^\pm }𝑑t\frac{1}{\tau _\pm }e^{t/\tau _\pm }=1\mathrm{exp}\left(\frac{L_\mathrm{G}}{v_x^\pm \tau _\pm }\right).$$
(2)
In the linear regime, the drift velocity is $`v_x^\pm =\mu _\pm _x`$ where $`_x`$ is the in-plane field due to the source-drain bias and $`\mu _\pm =\sigma _\pm /(n_\pm e)`$ are the spin-dependent mobilities. When the lateral extension of the 2DEG channel is shorter than the spin mean-free path, the majority and minority spin channels can be considered decoupled in the active region between the source and drain contacts. Hence, the spin-dependent drain current will be $`j_\pm ^\mathrm{D}\sigma _\pm (1\alpha _\pm )=\stackrel{~}{\sigma }_\pm `$, i.e., the channel conductance $`\sigma _\pm `$ is renormalized by the leakage probability. Note that the shorter the $`\tau _\pm `$, the smaller are both $`\sigma _\pm `$ and $`(1\alpha _\pm )`$. Hence, the leakage current, which ordinarily reduces the channel current, actually enhances the difference between the renormalized resistances $`1/\stackrel{~}{\sigma }_\pm `$ of the two spin channels.

(ii) Spin valve with two ferromagnetic gates
The spin-dependent in-plane transport due to the proximity of a ferromagnetic gate can be conveniently exploited to design a spin valve (see Fig. 1b). If two neighboring ferromagnetic gates are used, the electron current detected at the drain contact will depend on the relative orientation of the magnetizations of the two gates, giving rise to the magnetoresistance ratio $`=(j_{++}^\mathrm{D}j_+^\mathrm{D})/j_{++}^\mathrm{D}`$ where $`j_{++}^\mathrm{D}`$ and $`j_+^\mathrm{D}`$ are the total drain currents for parallel and antiparallel magnetizations respectively. For simplificity, we consider the case where the two ferromagnetic gates are identical and the effect of the channel region between the two gates on the current is negligible. The spin-dependent channel current density, $`j_\pm (x)=\sigma _\pm (x)_x\varphi _\pm (x)`$ is driven by the electrochemical potential $`\varphi _\pm (x)`$ and is constrained by the continuity equation $`_xj_\pm =j_\pm /[v_x^\pm (x)\tau _\pm (x)]`$ accounting for the leakage of current into the gates. Combining the two equations leads to a second order differential equation for $`\varphi _\pm (x)`$. The boundary conditions are (a) $`\varphi _+=\varphi _{}=\varphi _{source}`$ ($`\varphi _{drain}`$) at the edge of the gate near the nonmagnetic source (drain) contact, (b) $`\varphi _\pm (x)`$ and $`j_\pm (x)`$ vary continuously along $`x`$. Under the first gate, $`\sigma _\pm (x)=\sigma _\pm `$. Under the second gate, $`\sigma _\pm (x)=\sigma _\pm `$ in parallel configuration, but $`\sigma _\pm (x)=\sigma _{}`$ in the antiparallel case. We have solved exactly this equation both for the parallel and antiparallel configuration. The analytical expression for the magnetoresistance ratio $``$ is quite involved and will be presented in a longer publication. However, we have verified that a good approximation is given by

$$\left[\frac{\sigma _+(1\alpha _+)\sigma _{}(1\alpha _{})}{\sigma _+(1\alpha _+)+\sigma _{}(1\alpha _{})}\right]^2.$$
(3)
This formula is obtained by a circuit diagram of the two spin channels in parallel and by properly combining in series the leakage-renormalized conductances $`\stackrel{~}{\sigma }_\pm =\sigma _\pm (1\alpha _\pm )`$ for the two gates.

(iii) Model results for a silicon inversion layer
Here we present reasonable estimates for a $`\mathrm{Si}/\mathrm{SiO}_2/\mathrm{Fe}`$ heterostructure. We consider a free-electron Hamiltonian with two exchange split spin bands in the ferromagnet slon , an insulator barrier and the electrostatic confinement potential for the inversion layer. We will focus on the low temperature regime in which only the first quantized subband in the inversion layer is filled and the confining potential can be approximated with a triangular shape (Fig. 1a). As suggested by the recent literature on tunneling Si2 ; Si3 , the conservation of the in-plane wave-vector is relaxed and, hence, a one-dimensional calculation can be carried out. Within this simplified model, the spin-dependent wave-function $`\mathrm{\Psi }_\pm (z)`$ of the quasi-bound electron in the inversion layer can be written in terms of special Airy functions $`\mathrm{Ai}`$ and $`\mathrm{Bi}`$. For $`z>0`$ (the semiconductor region),

$$\mathrm{\Psi }_\pm (z)=A_\pm \mathrm{Ai}\left[k_s(z\overline{z}_\sigma )\right],$$
(4)
where $`k_\mathrm{s}^3=2m_\mathrm{s}e_\mathrm{s}/\mathrm{}^2`$ depends on the electric field $`_\mathrm{s}`$ in the inversion layer and the semiconductor effective mass $`m_\mathrm{s}`$. The position $`\overline{z}_\pm `$ represents the turning point. For $`d<z<0`$ (the insulator barrier region),

$$\mathrm{\Psi }_\pm (z)=\left\{B_\pm \mathrm{Ai}\left[k_\mathrm{b}(z\widehat{z}_\pm )\right]+C_\pm \mathrm{Bi}\left[k_\mathrm{b}(z\widehat{z}_\pm )\right]\right\},$$
(5)
where $`k_\mathrm{b}^3=2m_\mathrm{b}e_\mathrm{b}/\mathrm{}^2`$ and $`_\mathrm{b}=(ϵ_\mathrm{s}/ϵ_\mathrm{b})_\mathrm{s}`$ is the electric field in the barrier region ($`ϵ_\mathrm{b}`$ and $`ϵ_\mathrm{s}`$ being the dielectric constants) and $`\widehat{z}_\pm =(U_\mathrm{b}E_\pm ^\mathrm{s})/(e_\mathrm{b})`$ is the extrapolated turning point. The quantity $`E_\pm ^\mathrm{s}`$ is the quantization energy of the inversion layer subband and has to be calculated self-consistently. Finally, for $`z<d`$,

$$\mathrm{\Psi }_\pm (z)=D_\pm \mathrm{exp}\left[i\left(k_\pm ^{\mathrm{fm}}+i\frac{1}{l_\pm ^{\mathrm{fm}}}\right)z\right],$$
(6)
where $`k_\pm ^{\mathrm{fm}}=\sqrt{(k_\pm ^\mathrm{F})^2+2m_{\mathrm{fm}}(e_\mathrm{b}d+E_\pm ^\mathrm{s})/\mathrm{}^2}`$ in terms of $`k_\pm ^\mathrm{F}`$, the spin-dependent Fermi wave-vector and $`m_{\mathrm{fm}}`$, the ferromagnet band mass. The term $`l_\pm ^{\mathrm{fm}}`$ is the mean-free path in the ferromagnet. The coefficients $`A_\pm `$, $`B_\pm `$, $`C_\pm `$, $`D_\pm `$ are determined by the usual boundary conditions at the interfaces $`z=0`$ and $`z=d`$. The resultant secular equation for $`k\overline{z}_\pm `$ has a complex solution, yielding the spin-dependent complex  energy $`E_\pm ^\mathrm{s}=k_\mathrm{s}\overline{z}_\pm (\mathrm{}^2k_\mathrm{s}^2/2m_\mathrm{z})`$ in the inversion layer.

Fig. 2 depicts the majority and minority spin life times $`\tau _+^\mathrm{s}`$ (solid lines) and $`\tau _{}^\mathrm{s}`$ (dashed lines) as functions of the oxide thickness. The ultrathin range of the oxide thickness and the oxide electron effective mass are based on the recent measurements of electron tunneling across the silicon-based MOS junction Situnneling ; Si2 ; Si3 . The spin-dependent life time decreases exponentially with decreasing barrier thickness. The ratio of the times for two spins, $`\tau _{}^\mathrm{s}/\tau _+^\mathrm{s}`$, is dependent on the parameters of the ferromagnetic junction and can exceed a factor of 2. An increase of the electric field in the oxide (via the gate voltage) dramatically shorten the lifetimes, approaching the picosecond range. The oxide electric fields $`_\mathrm{b}`$ from the top to the bottom curves are 2, 7, and 12 $`\times 10^6`$ V/cm (the highest value approaches the dielectric breakdown). These results show a significant spin-dependence of the channel mobility for ultrathin oxide thickness. The spin splitting (not shown) is comparable to the energy broadening $`(\mathrm{}/2\tau _+^\mathrm{s})+(\mathrm{}/2\tau _{}^\mathrm{s})`$ and is about 0.3 meV when the lifetime is of the order 1 ps. Thus, for this specific system, the spin splitting can be disregarded because it is much smaller than the Fermi energy of the 2DEG. The leakage probability $`\alpha _\pm `$ depends also on the drift velocity according to Eq. (2). For a gate length $`L_\mathrm{G}=0.1\mu `$m, in-plane field $`_x=0.2\mathrm{V}/\mu `$m and Drude scattering time $`\tau _0=1`$ ps, the transit time in the gated region is $`\tau _{\mathrm{tr}}0.5`$ ps. This means that when the spin-dependent lifetimes are $`\tau _+=1`$ ps and $`\tau _{}=2`$ ps, the leakage probabilities are $`\alpha _+40\%`$ and $`\alpha _{}20\%`$. With these parameters, we get a magnetoresistance ratio $`10\%`$. A decrease of the drift velocity (through a decrease of the source-drain bias) will make the leakage of the channel current more severe (thus reducing the drain current), although it will increase $``$ according to Eq. (3).

In summary, we have predicted the spin-dependent properties of a 2DEG system due to the proximity of ferromagnetic gates, showing the potential spintronics applications of this kind of configuration.

###### Acknowledgements.

This work is supported by DARPA/ONR N0014-99-1-1096 and NSF DMR 0099572. CC is also grateful to the Swiss National Science Foundation for additional support. JPM acknowledges support by the California Institute for Telecommunications and Information Technology. We thank Edward Yu, Daniel Schaadt, and Ivan Schuller for helpful discussions.