text:
# Random Geometric Graphs

## 1 Introduction

The interest in complex networks has exploded over the last five years  where data on very large networks like the WWW , collaborations in the scientific community , transportation , movie actor collaborations  etc. have become accessible.

Random graphs are often used to model complex networks . Ever since Erdös and Rényi’s groundbraking work more than forty years ago , intense theoretical research on random graphs has been taking place . In contrast to random graphs the interactions between the sites in a lattice are usually between nearest neighbours, reflecting a myopic world. Lattices are therefore often said to be at the other end of the spectrum of network models .

Properties of real networks like robustness , growth , and topology have attracted much attention, primarily from physicists. It has been consistently shown that many of the networks possess small world characteristics . Like random graphs, small world networks are characterized by short average distances between any two sites, and by a high degree of localness, much like in lattices. However, individually, random graphs and lattice models in their pure forms are poor models of many real world networks. One could argue that high-dimensional lattices have the necessary high clustering and low average path length, though this has not been explored much . In the current paper we provide results on high-dimensional systems.

A random geometric graph (RGG) is a random graph with a metric. It is constructed by assigning each vertex random coordinates in a $`d`$-dimensional box of volume 1, i.e. each coordinate is drawn from a uniform distribution on the unit interval. RGGs have been used sporadically in real networks modeling  and extensively in continuum percolation , but almost exclusively in two and three dimensions. Although RGGs are the continuum version of lattices, they deserve some attention of their own, since percolating continuum systems display behaviour that lattices are incapable of . In addition, the connectivity in RGGs can be increased in a more natural way than by adding new bonds randomly in lattices.

Recently, continuum percolation has been used in the study of the stretched exponential decay of the correlation function in random walks on fractals and the conjectured relation to relaxation in complex systems . However, continuous systems in general and RGGs in particular are relevant whenever we need a multi-dimensional system with a metric, as for example when modeling the spread of diseases .

In this paper we study RGGs in arbitrary dimensions. In low dimensions the systems are dominated by local interactions. For higher dimensions RGGs are usually believed to approach standard random graphs, which we show is true only in some respects. We focus on ‘phase transitions’  at the percolation threshold by looking at the size of the largest cluster, and we determine how the value of the critical parameter in RGGs approaches that of random graphs as the dimension increases. We also extract the distribution of cluster sizes in the critical region. Furthermore, an expression for the cluster coefficient, a quantity that has attracted much interest in network theory recently, is derived. Results relevant for graph bi-partitioning are established. Finally, we discuss how to implement random geometric graphs efficiently.

The layout of this paper is as follows. In Section 2 and 3 we describe random graphs and random geometric graphs, respectively. In Section 4 we present our results, and Section 5 contains the details regarding the implementation. Finally, in Section 6 we sum up.

## 2 Random Graphs

Random graphs consist of $`N`$ vertices (points/sites) and $`K`$ edges (lines) where each possible edge is present with probability $`p`$, i.e. $`K=pN(N1)/2`$.111From here on we consider $`NN1`$ in accordance with the literature , since we are only investigating large systems. To keep the discussion independent of the system size $`N`$, graphs are often characterized by the connectivity (degree) $`\alpha =2K/N=pN`$, i.e. the average number of connections per vertex, instead of $`K`$ or $`p`$. As the connectivity increases clusters of vertices appear, where a cluster consists of all vertices linked together by edges, directly or indirectly.

The size of the largest cluster in the macroscopic limit $`N\mathrm{}`$ can be calculated analytically . It is $`NG(\alpha )`$, where

$$G(\alpha )=1\frac{1}{\alpha }\underset{n=1}{\overset{\mathrm{}}{}}\frac{n^{n1}}{n!}(\alpha e^\alpha )^n.$$
(1)
By the use of

$$y=\underset{n=1}{\overset{\mathrm{}}{}}\frac{n^{n1}}{n!}x^nx=ye^y$$
(2)
we can invert Eq. (1), getting

$$\alpha (G)=\frac{1}{G}\mathrm{log}(1G).$$
(3)
from which it is trivial to show that $`\alpha _c=1`$. With Eq. (3) it is an easy task to plot the fraction of vertices in the largest cluster—the giant component—as done in Fig. 1, where we see the prototype of a phase transition in combinatorial problems.

In random graphs the probability distribution of edges $`p_k`$ is binomial

$$p_k=(\begin{array}{c}N\\ k\end{array})p^k(1p)^{Nk}\frac{\alpha ^ke^\alpha }{k!},$$
(4)
where the approximation resulting in the Poisson distribution is valid for large systems sizes $`N`$, which is exactly the limit in which we are interested. The critical connectivity $`\alpha _c`$ for graphs with arbitrary random degree distribution $`p_k`$ has recently been derived by other techniques than those originally leading to Eq. (1. Unfortunately, we cannot use these results in connection with random geometric graphs, as will become clear in the next section.

## 3 Random Geometric Graphs

A $`d`$-dimensional random geometric graph (RGG) is a graph where each of the $`N`$ vertices is assigned random coordinates in the box $`[0,1]^d`$, and only points ‘close’ to each other are connected by an edge. The degree distribution of a RGG with average connectivity $`\alpha `$ is therefore given by Eq. (4) as well. However, a RGG is a special kind of random graph with properties not captured by the theoretical tools mentioned above. For one thing, the probability that three vertices are cyclically connected is different in random graphs and RGGs, regardless of the degree distribution of the random graph.

RGGs are sometimes named spatial graphs . Fig. 2 illustrates a RGG in $`2D`$. As in lattices, different boundary conditions can be applied. We will see that toroidal (continuous) boundary conditions make a vital difference compared to having open boundary conditions.

The volume of a $`d`$-dimensional (hyper)sphere with radius $`r`$ is

$$V_{sphere}=\frac{\pi ^{d/2}r^d}{\mathrm{\Gamma }(\frac{d+2}{2})},$$
(5)
where $`\mathrm{\Gamma }(x)`$ is the gamma function. This volume is needed in order to find the edges in RGGs.

To ‘visualise’ a RGG in general, one can think of a box filled with small spheres with radius $`r`$ and volume $`V`$ given by Eq. (5), where points are connected by an edge only if the distance between their centers is $`<2r`$, i.e. if the spheres overlap. Since the total volume of our box is 1, the probability that two arbitrarily chosen vertices are connected is equal to the volume of a sphere with radius $`R=2r`$. In continuum percolation theory this volume is denoted the *excluded volume* $`V_{ex}`$, where $`V_{ex}=2^dV`$ in a RGG. The excluded volume is the basic quantity of interest because it is directly related to the connectivity

$$\alpha =Np=NV_{ex},$$
(6)
from which it is clear why the connectivity is frequently called the total excluded volume of the system. Eqs. (5) and (6) give us

$$R=\frac{1}{\sqrt{\pi }}\left[\frac{\alpha }{N}\mathrm{\Gamma }\left(\frac{d+2}{2}\right)\right]^{\frac{1}{d}}.$$
(7)
Fig. 3 shows the radius $`R`$ of the excluded volume as a function of $`N/\alpha =1/p=1/V_{ex}`$. $`R`$ decreases monotonically: for a given connectivity $`\alpha `$ the spheres have to become smaller when more vertices are added to the graph.

Eq. (7) provides us with the required relation between $`\alpha `$ and $`R`$ when creating a RGG. The distance between every pair of vertices must be calculated, and an edge is added if the distance is less than $`R`$. Thus, it seems unavoidable to have a runtime of $`O(N^2)`$ making it unfeasible to investigate as large systems as with random graphs—see Fig. 1—where the number of calculations for a given $`\alpha `$ needed to create all the edges is $`O(N)`$. To overcome this obstacle we have designed a data structure which is described in Section 5, with a runtime of $`O(N^\beta )`$ where $`\beta 1.3`$. This allows us to study RGGs with up to $`N=4^{11}>410^6`$ vertices, which is more than an order of magnitude larger than usually accomplished .

## 4 Results

In our simulations of RGGs we define $`\alpha _c`$ to be the lowest connectivity at which the fraction of vertices in the largest cluster is $`>0`$ in the macroscopic limit. We make the bold claim that the systems we are able to analyse consist of enough points to make the critical connectivity almost as sharply defined as in Fig. 1. However, our main purpose is not to derive high precision percolation thresholds. Instead, we are more interested in the critical connectivity as a function of the dimension of the RGGs.

In this paper we express our threshold values in terms of $`\alpha `$. Other popular choices are the fractional volume $`s`$ occupied by the spheres  or the density $`N`$ of spheres. The relation between these parameters at the percolation threshold is

$$\alpha _c=N_cV_{ex}=2^d\mathrm{ln}(1s_c)$$
(8)
(see e.g.  for a derivation). Usually, in continuum percolation the volume $`V`$ of each sphere is fixed while $`N`$ is the independent variable in a system of size $`[0,L]^d`$. The approach of measuring $`N_c`$ or $`s_c`$ for various values of $`L`$ has been used in both two  and three  dimensions, i.e. for discs and spheres, where the critical values are determined by the use of finite size scaling. This procedure resembles site percolation in lattices. From the previous sections it is clear that we take a route closer to bond percolation in lattices by fixing $`L=1`$ while tuning $`\alpha `$ for different values of $`N`$. In Section 5 we describe how this has been carried out in practice.

### The Size of the Largest Cluster

Let $`G_d(\alpha )`$ denote the fraction of vertices in the largest cluster in $`d`$ dimensions. Since a RGG in the limit of infinite dimension is often assumed equivalent to a random graph, we expect that Eq. (3) provides us with an expression for $`G_{\mathrm{}}(\alpha )`$. But what does $`G_d(\alpha )`$ look like for finite $`d`$? And what is the behaviour of $`\alpha _c(d)`$? How does it approach $`\alpha _c(\mathrm{})`$ as $`d`$ increases? These are the questions addressed in this and the following section.

Figs. 4 and 5 illustrate the average size of the largest cluster in RGGs in 2, 3, 4, and 5 dimensions with and without toroidal boundary conditions. The curves correspond to $`N=4^k`$ vertices with $`k=5,6,\mathrm{},11`$, where the larger systems display the sharpest transitions. The legend in Fig. 5 applies to all diagrams in Figs. 4 and 5. In these 8 diagrams each curve is based on 300 data points. In other words, $`G_d(\alpha )`$ is calculated in intervals of $`\mathrm{\Delta }\alpha =0.005`$ resulting in the smooth lines in the figures. For every data set we have averaged over enough runs for error bars to be completely negligible.

Since continuous boundary conditions mean addition of extra edges, the size of the largest component $`G(\alpha )`$ obviously grows faster in Fig. 5 than in Fig. 4, especially in the smaller systems. These relatively few extra edges make a decisive difference, connecting vertices not already in the same cluster. Since toroidal systems are models of bulk systems, $`G`$ is much less $`N`$-dependent in that case. However ‘unphysical’ RGGs with open boundaries may seem, they are the most popular RGG version in the literature. Consequently, we consider them alongside the continuous case.

From Figs. 4 and 5 we see that the continuous boundary conditions make the transition where $`G>0`$ more abrupt, but that an estimation of $`\alpha _c`$ does not depend much on the boundary conditions if only we base our judgment on large enough systems. This is confirmed in the inset of Fig. 4, where $`\alpha _c=4.53`$ is obtained by finite size scaling, i.e. plotting $`G(x)`$ where $`x=N^{1/\nu }(\alpha \alpha _c)`$. However, it is clearly easier to make precise estimates of the critical connectivity with than without continuous boundary conditions. We note in passing that the exponent $`\nu =3`$ is equal to the value of $`\nu `$ found in random graphs .

### The Critical Connectivity

With numerically obtained knowledge of $`G(\alpha )`$, it is possible to extract $`\alpha _c`$. The procedure is simple. By inspection of Fig. 5 we can estimate $`\alpha _c`$ for $`d5`$. To obtain further data points we have run our algorithm on RGGs with $`N=4^{10}`$ for systems of larger dimensions as well. Though this results in increased runtime per graph, the results get more homogeneous and fewer runs are needed in order to get a decent estimate of $`G_d(\alpha )`$. Our findings presented in Table 1 and Fig. 6 strongly suggest that

$$\alpha _c(d)=\alpha _c(\mathrm{})+Ad^\gamma ,$$
(9)
where $`\alpha _c(\mathrm{})=1`$, $`\gamma =1.74(2)`$ and $`A=11.78(5)`$. As expected, Eq. (9) predicts that $`\alpha _c(\mathrm{})`$ is equal to $`\alpha _c`$ in random graphs, confirming that RGGs and random graphs become more and more similar as $`d`$ increases. However, when we derive the cluster coefficient, we will see that this is not true in all respects.

Finally, we note that our findings are in accordance with the most precise estimates that we know of: $`\alpha _c=4.51223(5)`$  and $`\alpha _c=2.734(6)`$  in $`2D`$ and $`3D`$, respectively, obtained by the use of finite size scaling. For $`d>3`$ we have not been able to find any estimates of $`\alpha _c`$ to compare with .

### The Distribution of Cluster Sizes

Having examined the size of the largest cluster and the critical connectivity, we now look at the distribution of cluster sizes in RGGs.

The inset illustrates the scale free power-law distribution at $`\alpha =2.6`$. Right below $`\alpha _c`$, clusters of all sizes can be encountered. The small hump at large cluster sizes is always present because the clusters cannot contain more than all of the vertices. The clusters pile up when their size approaches this boundary, in this case a cluster size of $`1000`$, just below the inevitable cut-off.

Our simulations show that for $`\alpha `$ significantly below $`\alpha _c`$ the distribution is approximately exponential. As the connectivity increases the distribution becomes power-law-like. As $`\alpha `$ is further increased the distribution is separated in two parts; there are no clusters of medium size, only the largest macroscopic cluster and a few small ones around it. We have observed this overall behaviour in all our tests of the distribution of cluster sizes in various dimensions.

Fig. 7 shows our data in $`3D`$. For $`\alpha =2.1`$ ($``$) the data points lie on an almost straight line indicating an exponential distribution. Increasing the connectivity to $`\alpha =2.4`$ ($`\mathrm{}`$) results in a broader distribution that is no longer exponential. Right at the critical connectivity ($``$) the distribution flattens out. Clusters of all sizes are observed. Right above $`\alpha _c`$ ($``$) two separate regions begin to materialise. Already at $`\alpha =3.3`$ ($``$) the largest cluster makes it highly unlikely that a cluster of medium size can be present as well. The distribution is cut in two.

### The Cluster Coefficient

In network theory the cluster coefficient $`C`$ is an often calculated quantity , which is defined in the following way. Let the vertices $`i`$ and $`j`$ be connected directly to a common vertex $`k`$. $`C`$ is then the probability that vertex $`i`$ and vertex $`j`$ are directly connected as well. From this we see that the cluster coefficient is a measure of the ‘cliquishness’ of the graph. In this section we derive $`C=C_d`$ analytically in arbitrary dimensions $`d`$, showing that $`C_d`$ decreases in an exponential fashion.

To determine $`C_d`$ we make use of the concept of the excluded volume $`V_{ex}`$. If we again use the vertices $`i`$, $`j`$, and $`k`$, then $`i`$ and $`j`$ must both be within the excluded volume of $`k`$. Put differently, the probability that $`i`$ and $`j`$ are connected is equal to the probability that two randomly chosen points in a sphere of volume $`V_{ex}`$ and radius $`R`$ is less than a distance $`R`$ apart. In other words, given the coordinates of vertex $`i`$ the probability that there is an edge between $`i`$ and $`j`$ is equal to the fraction of the excluded volume of vertex $`i`$ that lies inside the excluded volume of $`k`$. By averaging over all points in $`V_{ex}`$ we get the cluster coefficient $`C_d`$.

The task of calculating $`C_d`$ is considerably simplified by the spherical symmetry of the problem. The fractional volume ‘overlap’ $`\rho _d`$ of two spheres only depends on the distance $`r`$ between the centers and not on any angular parts, i.e. $`\rho _d=\rho _d(r)`$. In general, the cluster coefficient can therefore be written as

$$C_d=\frac{1}{V_{ex}}_{V_{ex}}\rho _d(r)\text{d}V.$$
(10)
In Appendix A we derive that

$$C_d=\{\begin{array}{cc}1H_d(1)\hfill & \text{even }d\hfill \\ & \\ \frac{3}{2}H_d(\frac{1}{2})\hfill & \text{odd }d\hfill \end{array}$$
(11)
where

$$H_d(x)=\frac{1}{\sqrt{\pi }}\underset{i=x}{\overset{d/2}{}}\frac{\mathrm{\Gamma }(i)}{\mathrm{\Gamma }(i+\frac{1}{2})}\left(\frac{3}{4}\right)^{i+\frac{1}{2}}.$$
(12)
When $`d`$ is large Eq. (11) reduces to (see Appendix A)

$$C_d3\sqrt{\frac{2}{\pi d}}\left(\frac{3}{4}\right)^{\frac{d+1}{2}}.$$
(13)
The cluster coefficient is plotted in Fig. 8 ($``$) together with the asymptotic solution in Eq. (13) (full line).

Eq. (11) shows that the cluster coefficient is a purely geometric quantity depending only on the dimension $`d`$; neither the connectivity $`\alpha `$ nor the system size $`N`$ are present. In random graphs $`C=\alpha /N`$, since there is per definition no correlation between edges. So, in contrast to what is usually believed, RGGs are *not* identical to random graphs when $`d\mathrm{}`$.

In higher dimensions, the cluster coefficient in RGGs becomes exceedingly small. This peculiar fact can be explained by noting that the distribution of distances between two connected vertices gets more and more peaked at the maximal distance $`R`$ as $`d`$ increases. This implies that if the vertices $`i`$ and $`j`$ are both connected to vertex $`k`$ in a high-dimensional space, then it is highly unlikely that $`i`$ and $`j`$ are directly connected by an edge as well. Only in low dimensions are RGGs dominated by small loops. On the contrary, the way that a standard random graph is designed implies a cluster coefficient which can only be interpreted statistically, and not geometrically. Despite the fact that $`\alpha _c=1`$ in both random graphs and RGGs of infinite dimensionality, they do not have the same topology.

### Graph Bi-partitioning

Random geometric graphs are useful outside network modeling and percolation theory as well. In this section we look at RGGs in relation to graph bi-partitioning, a well known problem in combinatorial optimization.

The NP-hard problem of partitioning a graph with $`N`$ vertices in two subsets with $`N/2`$ vertices each, in such a way that the cutsize $`E`$, i.e. the number of edges between vertices in different subsets, is minimized, is called the graph bi-partitioning (GBP) problem. Fig. 2 illustrates a bi-partitioned RGG, where $`N/2`$ of the points are marked by squares, the other half being dots.

The GBP problem of RGGs with open boundary conditions has been tested by various heuristics . In this section we use our numerical findings to establish the critical connectivity in relation to GBP. Additionally, for $`\alpha >\alpha _c^{GBP}`$ we argue that the cutsize $`E`$ depends on $`N`$ and $`\alpha `$ in a simple way.

In GBP the connectivity is critical when $`G=1/2`$. As soon as the largest cluster contains more than half of the vertices, it becomes impossible to bipartition the graph without violating any edges. For random graphs Eq. (3) immediately gives us $`\alpha _c^{GBP}=2\mathrm{ln}21.386`$.

In RGGs $`\alpha _c^{GBP}(d)`$ can be extracted in the same way as $`\alpha _c`$ was in Section 4. Our numerical findings in RGGs with continuous boundary conditions are presented in Table 2. We stress that the results are valid only for large $`N`$, as a closer look at Fig. 5 reveals. In $`2D`$ the average fraction of vertices in the largest cluster is independent of $`N`$ only for $`\alpha >\alpha _c^{GBP}`$. This means that if one looks at GBP in $`2D`$ with $`N=1000`$, one cannot use the value of $`\alpha _c^{GBP}`$ in Table 2. In higher dimensions the interval around $`\alpha _c`$ where $`G_d(\alpha )`$ is size-dependent gets smaller and does not play a role in relation to GBP.

With open boundary conditions the picture is messy, as Fig. 4 shows. In this case $`G(\alpha )`$ is highly $`N`$-dependent, and it is not possible to speak of a critical connectivity $`\alpha _c^{GBP}`$ without specifying $`N`$. This is true despite the fact that $`G(\alpha )`$ is an averaged quantity, i.e. for small $`N`$ will a fraction of the graphs contain a cluster with more than $`N/2`$ vertices even when $`\alpha <\alpha _c^{GBP}`$. Fig. 4 clearly shows that $`\alpha _c^{GBP}`$ is a decreasing function of $`N`$ for $`d>2`$. In $`2D`$ however, all curves cross at almost the same (pivotal) point, and it is reasonable to speak of $`\alpha _c^{GBP}`$ without specifying $`N`$. As the inset in Fig. 4 shows this would lead to an estimate of $`\alpha _c^{GBP}=4.53(1)`$, close to $`\alpha _c^{GBP}`$ in RGGs with toroidal boundary conditions.

The size of the largest cluster near $`\alpha _c`$ grows so rapidly in $`2D`$ that $`\alpha _c=\alpha _c^{GBP}`$ cannot be ruled out on the basis of our numerical data. This is true with both open and continuous boundary conditions. However, as this would imply that the phase transition is of 1st order in $`2D`$ only, we believe that the two critical connectivities are close but not identical.

When bi-partitioning a RGG, it is obvious that the ‘area of contact’  between the two subsets in the optimal configuration must be close to a minimum. In $`2D`$ this means that the best achievable partition must be close to simply cutting the graph in two at the coordinate values $`x_1=1/2`$ or $`x_2=1/2`$. This observation is especially relevant for large connectivities where the cutsize is, fluctuations neglected, proportional to the length of the dividing line. All this tentatively indicates how the cutsize $`E`$ in GBP behaves as a function of $`N`$ and $`\alpha `$ by looking at RGGs partitioned at $`x_i=1/2`$, where $`1id`$. As we are about to argue, we expect a scaling relation like

$$E_dN^{1/\nu }\alpha ^\beta (d),$$
(14)
where the exponents $`\nu `$ and $`\beta `$ only depend on the dimension of the RGG.

The exponents in Eq. (14) can be determined in the following way. Given the radius $`R`$ of the excluded volume of each vertex, the cutsize must be proportional to $`NR`$, since only vertices with $`1/2R<x_i<1/2`$ contribute to the cutsize (to avoid counting the violated edges twice we only look at the vertices at one side of the partitioning plane at $`x_i=1/2`$), times the average number of violated edges per vertex in this region, which is proportional to $`NR^d`$. In other words,

$$E_dN^2R^{d+1}.$$
(15)
If instead of $`R`$ we want to express the result in terms of $`\alpha (d)NR^d`$, we get

$$1/\nu =1\frac{1}{d},\beta =1+\frac{1}{d}.$$
(16)
Since $`EN^2`$ in Eq. (15), the relation $`1/\nu +\beta =2`$ holds in arbitrary dimensions.

Now, it is obvious that the scaling Ansatz is reasonable only for $`\alpha >\alpha _c^{GBP}`$. As Fig. 2 illustrates, the optimal partition at $`\alpha \alpha _c^{GBP}`$ is highly complex and not at all close to a straight line. If we incorporate that $`E=0`$ for $`\alpha <\alpha _c^{GBP}`$ and replace Eq. (14) with

$$E_dN^{1/\nu }(\alpha (d)\alpha _c^{GBP})^\beta ,$$
(17)
we do not expect Eq. (16) to hold if we focus only on a region near the critical connectivity. By the use of extremal optimization, a heuristic that works particularly well near phase transitions in hard combinatorial problems, Boettcher and Percus  have found $`\alpha _c^{GBP}4.1`$, $`1/\nu 0.6`$ and $`\beta 1.4`$ in $`2D`$ for $`4<\alpha <6`$, not far off our estimates in Eq. (16) valid for large connectivities. Note that the low estimate of $`\alpha _c^{GBP}`$ is expected; the algorithm does not always find the best partition, and some graphs with $`\alpha <\alpha _c`$ does have $`E>0`$.

## 5 Implementation

The implementation is of major importance when studying random geometric graphs, since a straightforward check of all possible edges between the $`N`$ points will result in unfeasible runtimes $`O(N^2)`$. We now outline how our program works and describe how to avoid runtimes $`O(N^2)`$.

The main idea is to divide and conquer. Partition the $`d`$-dimensional box in smaller subboxes and determine which subbox each vertex belongs to. Given the connectivity and thereby the radius $`R`$ of the excluded volume, for each vertex we then only have to look for potential edges to vertices in the subboxes adjacent to the subbox where the vertex itself is located. This leads to a huge reduction in the number of comparisons. And this just gets better when $`N`$ increases, resulting in a decrease in $`R`$ as we saw in Fig. 3. By partitioning the box further as $`N`$ increases we avoid a linear increase in the number of comparisons per vertex, which would lead to the undesirable $`O(N^2)`$ growth.

The algorithm used when looking at RGGs is simple. It works like this:

1. Generate $`d`$ coordinates for each vertex.
2. Partition the space in small subboxes.
3. Find the edges.
4. Calculate the relevant quantities ($`G`$, cluster sizes etc.) as $`\alpha `$ increases.

Obviously, a trade-off in Step 2 is involved when choosing the number of small boxes.

Being the most time consuming part of the algorithm, Step 3 is the main contributor when deciding how the runtime depends on $`N`$. The runtimes for most of our runs are shown in Fig. 9. We see that the runtime is $`O(N^\beta )`$, where $`\beta 1.3`$, resulting in ‘feasible’ runtimes for graphs with up to $`N>410^6`$. Note that the runtime of the much simpler algorithm used on random graphs also grows like a power-law with $`\beta =1.15`$, even though the number of operations is clearly $`O(N)`$. In fact, the number of comparisons with potential neighbours per vertex is very nearly constant in our implementation, i.e. the total number of neighbour tests is $`O(N)`$ in RGGs as well. Of course, this is only possible if the number of subboxes also increases with $`N`$. Managing the partitioning part of the algorithm adds to the runtime. To sum up, the power-law increase in the runtime illustrated in Fig. 9 for both random graphs and RGGs is probably mainly due to cache misses. The slightly higher values of $`\beta `$ in the RGGs stems from the additional time used when partitioning the $`d`$-dimensional box into smaller boxes.

Step 4 is worth a comment. When running the algorithm we are interested in information at certain values of $`\alpha `$. Instead of generating a new graph for every data point needed, we first set up the graph with the minimal connectivity we want to look at. This is easily accomplished with our algorithm. Given an $`\alpha `$-window $`[\alpha _{min},\alpha _{max}]`$ in which we want to examine the graph, we find all the edges belonging to the graph when $`\alpha =\alpha _{max}`$, but we only add the edges corresponding to $`\alpha =\alpha _{min}`$. The rest of the edges, those who are to be added when $`\alpha `$ is gradually increased to $`\alpha _{max}`$, are stored in a priority queue. It is then a simple task to increase $`\alpha `$ as one wishes. As mentioned earlier, in Figs. 4 and 5 each curve is based upon 300 data points, i.e. $`\mathrm{\Delta }\alpha =0.005`$.

The source code, written in C, is available upon request. For a more accurate and technical discussion of fast algorithms in relation to RGGs, see e.g. .

## 6 Summary

In this paper we have illustrated the usefulness of random geometric graphs in network theory and how to implement them efficiently. Several properties of random geometric graphs in the vicinity of the critical connectivity $`\alpha _c`$ have been analysed. We have determined the size of the largest cluster numerically and shown that $`\alpha _c(d)`$ approaches $`\alpha _c(\mathrm{})=1`$ found in random graphs in a power-law fashion. We have verified that the distribution of cluster sizes is cut in two just when the connectivity becomes larger than $`\alpha _c`$. Interestingly, the derivation of the cluster coefficient shows that, even in the limit of infinite dimensionality $`d`$, random geometric graphs are not identical to random graphs.

Random geometric graphs share properties with both lattice models and standard random graphs. Random geometric graphs allow us to work with random graphs with a local structure. In addition, it is straightforward to add ‘long’ edges if one wishes to simulate, e.g., a small world network. With all this in mind, we hope this paper will make random geometric graphs more widely used in network theory.

## Appendix A Derivation of $`C_d`$

In order to determine the cluster coefficient for arbitrary $`d`$, one must find the fractional overlap $`\rho _d`$. Since $`\rho _d`$ has no angular dependence, Eq. (10) reduces to

$$C_d=\frac{d}{R^d}_0^R\rho _d(r)r^{d1}\text{d}r.$$
(18)
Since $`\rho _1=1\frac{r}{2R}`$, $`C_1=\frac{3}{4}`$. From Fig. 10 we see that in $`2D`$ the overlapping area—the area circumscribed by the fat lines—is $`2(AB)`$, where $`A`$ is the area of the part of the circle swept out by the angle $`\theta =2\mathrm{arccos}(r/2R)`$ between the two dashed lines originating from the center of the lowest circle, and $`B`$ is the area of the dashed triangle. Now, $`A=\frac{1}{2}\theta R^2`$ and $`B=R^2\mathrm{cos}(\theta /2)\mathrm{sin}(\theta /2)=\frac{1}{2}R^2\mathrm{sin}\theta `$. The area of the overlap is then $`R^2(\theta \mathrm{sin}\theta )`$, so $`\rho _2=\frac{1}{\pi }(\theta \mathrm{sin}\theta )`$ and $`C_2=1\frac{3\sqrt{3}}{4\pi }.`$

For $`d3`$, the use of cylindrical coordinates and the relation

$$2\pi \underset{i=2}{\overset{n1}{}}_0^\pi \mathrm{sin}^{ni}\theta _i\text{d}\theta _i=\frac{n\pi ^{n/2}}{\mathrm{\Gamma }(\frac{n+2}{2})}$$
(19)
results in

$$\rho _d(r)=\frac{2}{\sqrt{\pi }}\frac{\mathrm{\Gamma }(\frac{d+2}{2})}{\mathrm{\Gamma }(\frac{d+1}{2})}_0^{\mathrm{arccos}(\frac{r}{2R})}\mathrm{sin}^d\theta \text{d}\theta .$$
(20)
By reversing the integration in $`C_d`$ we get

$$C_d=\frac{3}{\sqrt{\pi }}\frac{\mathrm{\Gamma }(\frac{d+2}{2})}{\mathrm{\Gamma }(\frac{d+1}{2})}_0^{\frac{\pi }{3}}\mathrm{sin}^d\theta \text{d}\theta ,$$
(21)
which can be solved by integration by parts. The use of the duplicate formula for the Gamma function then finally leads to Eq. (11).

For large $`d`$, the ratio of the Gamma functions in Eq. (21) is given by Stirling’s approximation. By putting $`x=\mathrm{cos}\theta 1/2`$, the cluster coefficient can therefore be written as

$$C_d\sqrt{\frac{6d}{\pi }}\left(\frac{3}{4}\right)^{\frac{d}{2}}_0^{\frac{1}{2}}\mathrm{exp}[\frac{d1}{2}\mathrm{ln}f(x)]\text{d}x,$$
(22)
where $`f(x)=1\frac{4x}{3}(1+x)`$. Since the contributions to the integral for large $`d`$ are significant only when $`x0`$, $`\mathrm{ln}f`$ can be expanded to 1st order and Eq. (13) is recovered.

Acknowledgements

We thank J. Neil Bearden, Stefan Boettcher and Paolo Sibani for helpful comments. We are especially indebted to Allon Percus for very valuable discussions and correspondence. This work was financially supported by Statens Naturvidenskabelige Forskningsråd.


format: text

text:
# Interaction Of Electrons With Spin Waves In The Bulk And In Multilayers

(Physics Department, Carnegie Mellon University, Pittsburgh, PA 15213)

The exchange interaction between electrons and magnetic spins is considerably enhanced near interfaces, in magnetic multilayers. As a result, a dc current can be used to generate spin oscillations. We review theory and experimental evidence. The s-d exchange interaction causes a rapid precession of itinerant conduction-electron spins $`\stackrel{}{s}`$ around the localized spins $`\stackrel{}{S}`$ of magnetic electrons. This $`\stackrel{}{s}`$ precession has been observed directly (Weber et al., 2001) with electron beams through Fe, Co, and Ni films. Because of it, the time-averaged interaction torque between $`\stackrel{}{s}`$ and $`\stackrel{}{S}`$ vanishes. Thus, electrons do not interact at all with long-wavelength spin waves, in the bulk. An interface between a magnetic layer and a spacer causes a local coherence between the precession phases of different electrons, in a region within 10 nm from the interface (Slonczewski, 1996; Berger, 1996). Also, a second magnetic layer with pinned $`\stackrel{}{S}`$ is used to ”prepare” $`\stackrel{}{s}`$ in a specific direction, before the electrons cross over to the active magnetic layer. The current-induced drive torque of $`\stackrel{}{s}`$ on $`\stackrel{}{S}`$ in the active layer may be calculated from the spin current (Slonczewski) or from the spin imbalance $`\mathrm{\Delta }\mu `$ (Berger). This torque is equivalent to negative Gilbert damping, leading to instability of $`\stackrel{}{S}`$. Spin current and $`\mathrm{\Delta }\mu `$ are proportional to each other (Berger, 2001) and can arise fom Fermi-surface translation, as well as from expansion/contraction. In fields $`\stackrel{}{H}`$ normal to layers, the critical current $`I_c`$ for $`\stackrel{}{S}`$ instability is predicted to be proportional to the ferromagnetic-resonance frequency $`\omega `$ (consistent with the Tsoi et al. 1998, 2000, experiments). However, for in-plane $`\stackrel{}{H}`$, due to elliptic $`\stackrel{}{S}`$ precession, $`I_c`$ is not proportional to $`\omega `$, but linear in H (Katine et al. 2000, experiments). Apart from the current-induced drive torque, an extra Gilbert damping is predicted near the interface even at zero current (Berger, 1996). It has been observed by ferromagnetic resonance (Urban et al., to be published).

I. INTRODUCTION AND INTERACTIONS IN THE BULK

Recently, it has been suggested that the exchange interaction between electrons and magnetic spins is considerably enhanced near interfaces, in magnetic multilayers. As a result, a dc current can be used to generate spin oscillations, or even to switch the spin direction from up to down. Using classical spin models and simple mathematics, we start by describing the situation in bulk samples and single thin films. The effect of interfaces will be discussed in Sections II and III. Then, some of the existing experimental evidence will be covered. The analogy with a semiconductor laser is explained in Section VIII.

In metallic ferromagnets, the s-d exchange interaction $`2J_{sd}\stackrel{}{s}\stackrel{}{S}(\stackrel{}{r})`$ couples 1 the localized spin $`\stackrel{}{S}(\stackrel{}{r})`$ of 3d magnetic electrons to the itinerant spins $`\stackrel{}{s}`$ of 4s conduction electrons. This interaction can be approximately represented by a classical field $`\stackrel{}{H}_{sd}=J_{sd}<\stackrel{}{S}(\stackrel{}{r})>/\mu _B`$ acting on $`\stackrel{}{s}`$. Here, $`\mu _B`$ is a Bohr magneton, and we have $`J_{sd}0.1eV`$ and $`\mu _0H_{sd}1000T`$.

We assume that a spin wave is present, i.e., that $`\stackrel{}{S}`$ is precessing around the field $`\stackrel{}{H}`$ (Fig. 1a), which includes the external field as well as demagnetizing and anisotropy fields. In the 1950s and 1960s, theoretical investigations of electron scattering by spin waves were often based on the first Born approximation. In that approximation, s-d exchange is assumed to be a small perturbation. In classical language, this means that $`\stackrel{}{s}`$ remains nearly fixed and nearly parallel or antiparallel to $`\stackrel{}{H}`$ (Fig. 1a), as $`\stackrel{}{S}`$ precesses around $`\stackrel{}{H}`$.

However, since $`H_{sd}H`$, s-d exchange is clearly too large for the Born approximation to be valid. A more correct picture 2 (Fig. 1b) has $`\stackrel{}{s}`$ precessing rapidly around the instantaneous $`\stackrel{}{S}`$ direction, as $`\stackrel{}{S}`$ itself precesses more slowly around $`\stackrel{}{H}`$. Depending on initial conditions, in many cases $`\stackrel{}{s}`$ follows closely the instantaneous $`\stackrel{}{S}`$ direction, namely, $`\stackrel{}{s}`$ is nearly parallel to $`\stackrel{}{S}`$ for a spin-up electron (Fig. 1b) and antiparallel for spin-down. This could be called the strong-coupling (or adiabatic) approximation. Because the angle between $`\stackrel{}{s}`$ and $`\stackrel{}{S}`$ is now smaller (Fig. 1b), the instantaneous value of the exchange torque between them is much reduced. Also, because of the rapid $`\stackrel{}{s}`$ precession around $`\stackrel{}{S}`$, that torque changes sign frequently. Finally, different conduction electrons precess out of phase with each other in a bulk sample, so that their torques tend to cancel. We conclude 2 that electrons do not interact practically at all with spin waves through ordinary s-d exchange, at least for long-wavelength spin waves!

This interesting conclusion must be modified 2 if the electron collision rate with solutes or phonons is high and comparable to the $`\stackrel{}{s}`$ precession frequency $`\omega _s`$. Also, if the spin-wave wavelength is short and comparable to $`v_F/\omega _s`$, where $`v_F`$ is the electron Fermi velocity; indeed, a current-induced exchange torque leading to a spin-wave instability has been derived 3 quantum-mechanically for finite spin-wave wavelength. The present paper will emphasize the limit of long wavelength (i.e., uniform $`\stackrel{}{S}`$ precession) and long electron collision time.

The precession of $`\stackrel{}{s}`$ around $`\stackrel{}{S}`$ has been observed directly 4 by passing a polarized electron beam through a film of Fe, Co or Ni of thickness 0-10 nm (Fig. 1c). The magnetic film is supported by a gold foil of thickness 20 nm. The electron beam is generated by a gun having a GaAs photocathode. The electron spin $`\stackrel{}{s}`$ is originally perpendicular to the fixed magnetic spins $`\stackrel{}{S}`$ of the uniformly magnetized film. After passing through the film, $`\stackrel{}{s}`$ is measured, and is found to be rotated by an angle $`ϵ`$ around $`\stackrel{}{S}`$. The value of $`ϵ`$ is proportional to magnetic-film thickness. In cobalt, for example, the proportionality constant is $`19^{}/nm`$. In these experiments, the incident electron energy exceeds the Fermi level in the film by only $`47eV`$, and its exact value is found to have little effect on $`ϵ`$.

II. EFFECT OF REFERENCE LAYER AND INTERFACE

S-d exchange torques exerted by conduction electrons on precessing magnetic spins $`\stackrel{}{S}_2`$ in layer $`F_2`$ (Fig. 2a) can be increased 5,6 over their bulk value, by introducing another magnetic layer $`F_1`$ separated from $`F_2`$ by a nonmagnetic spacer N. The magnetic spins $`\stackrel{}{S}_1`$ in $`F_1`$ are assumed to have a fixed direction parallel to the precession axis of $`\stackrel{}{S}_2`$. If the $`F_1`$ thickness is at least a local spin-diffusion length, this layer is thick enough to control the ratio of spin-up current to spin-down current in the multilayer, and to give to this ratio a value as different from one as possible (as needed 6 for current-driven experiments). As an example, the spin-diffusion length is 44 nm in Co nanowires 7 at 77 K. A second role of layer $`F_1`$ is to relax the spin $`\stackrel{}{s}`$ of electrons in $`F_1`$ to a direction parallel (or antiparallel) to $`\stackrel{}{S}_1`$, under the combined influence of s-d exchange, spin-orbit interaction and random collisions. Thus, $`F_1`$ prepares $`\stackrel{}{s}`$ in a fixed direction, and could be called the “reference” layer.

Some electrons originating in $`F_1`$ cross over to $`F_2`$ (Fig. 2a). By solving the Schroedinger equation 5,6,8, while assuming static and uniform $`\stackrel{}{S}_2`$, one finds that the $`\stackrel{}{s}`$ direction changes somewhat as the electron crosses the interfaces between these layers, but remains roughly parallel or antiparallel to $`\stackrel{}{S}_1`$. Thus, the angle between this $`\stackrel{}{s}`$ and $`\stackrel{}{S}_2`$ in $`F_2`$ is large, and comparable to the $`\stackrel{}{S}_2`$ precession-cone angle $`\theta `$ (Fig. 2a). In turn, this implies a large instantaneous torque between $`\stackrel{}{s}`$ and $`\stackrel{}{S}_2`$, as in the case of the Born approximation of Fig. 1a. But, as the electron propagates inside $`F_2`$, away from the $`NF_2`$ interface, $`\stackrel{}{s}`$ precesses rapidly around $`\stackrel{}{S}_2`$ (Fig. 2b). As a result, the total torque exerted on all magnetic spins in $`F_2`$, and obtained by summation over all $`\stackrel{}{s}`$ locations in $`F_2`$, is not in the tangential direction shown as $`\stackrel{}{\tau }_n`$ in Fig. 2a, as one might naively think. Rather, it is 5,6 in a direction $`\stackrel{}{\tau }_d`$ radial to the $`\stackrel{}{S}_2`$ precession orbit; one exception being the case of a very thin layer $`F_2`$ discussed in Section VI.

Electrons with opposite $`\stackrel{}{s}`$ produce opposite torques. Hence, after summation over all electrons, the total torque $`\stackrel{}{\tau }_d`$ is found to be zero in equilibrium. An electric current I normal to layers is needed to populate certain electron states near the Fermi level and to depopulate others, resulting in a torque proportional to I. It is also approximately proportional 5,6 to $`sin\theta `$ for small $`\theta `$, and may therefore be written 5 in the convenient form

$$\stackrel{}{\tau }_d=g(\theta )(\mathrm{}/2eS_1S_2^2)I\stackrel{}{S}_2\times (\stackrel{}{S}_1\times \stackrel{}{S}_2),$$
(1)
where $`g(\theta )`$ is a slowly varying function of $`\theta `$. Note that the direction of a Gilbert or Landau-Lifshitz damping torque is also radial 9. However, it is directed toward the precession axis of $`\stackrel{}{S}_2`$, while the current-induced drive torque $`\stackrel{}{\tau }_d`$ is directed outwards provided I flows (Fig. 2a) from $`F_1`$ toward $`F_2`$, i.e., $`I>0`$.

For two electrons moving in different directions inside $`F_2`$ (Fig. 2b), the phases of $`\stackrel{}{s}`$ precession are equal at first, but become different at a given distance x from the $`NF_2`$ interface, as x increases from zero. Because of boundary conditions at the interface, the phase for an electron moving nearly parallel to the x axis changes with x at a slower rate than the phase for an electron going in an oblique direction, assuming a spherical Fermi surface. Because of this dephasing effect, the torques exerted on $`\stackrel{}{S}_2`$ by different electrons tend to cancel 5 for x larger than a certain characteristic length $`L_0=\pi /(k_{}k_{})`$, where $`k_{},k_{}`$ are the spin-up and spin-down wavenumbers in $`F_2`$ (Fig. 2b). From the data of Ref. 4, and assuming a reasonable free-electron Fermi velocity of $`1\times 10^6m/s`$, we derive $`L_09.5nm`$ and $`H_{sd}1900T`$ for cobalt. The length $`L_0`$ plays in the present problem the same role as the transverse spin relaxation time $`T_2`$ in nuclear magnetic resonance.

Since regions of $`F_2`$ with $`x>L_0`$ do not contribute to the total torque on $`F_2`$ (except 8 for some periodic fluctuations), the current-induced torque per unit volume of $`F_2`$ is inversely proportional 5,6,8 to the $`F_2`$ thickness $`L_2^x`$ if $`L_2^x>L_0`$. Electron scattering by solutes or phonons also reduces 6 the contribution of regions with large x to the torque, and has been observed 4 to change $`\stackrel{}{s}`$ in the transmitted beam to a direction making a smaller angle with $`\stackrel{}{S}_2`$.

III. SPIN CURRENT AND SPIN IMBALANCE

In order to calculate the total drive torque $`\stackrel{}{\tau }_d`$ on layer $`F_2`$ (Section II), one can introduce 5 the concept of angular-momentum current $`\stackrel{}{P}(x)`$, often called spin current. It is related to $`\stackrel{}{\tau }_d`$ by $`\stackrel{}{\tau }_d=\stackrel{}{P}(x=0)\stackrel{}{P}(x=L_2^x)`$. In practice, the second term is of no importance if $`L_2^xL_0`$, for the same reasons (Section II) that regions with $`x>L_0`$ do not contribute to the total torque. This approach has the advantage that the details of the exchange interaction between $`\stackrel{}{s}`$ and $`\stackrel{}{S}_2`$ do not appear explicitly in the calculations. Also, it applies even to very large $`\theta `$ values.

As an alternative approach, one can introduce 6 the difference $`\mathrm{\Delta }\mu =\mu _{}\mu `$ between the Fermi levels $`\mu _{},\mu _{}`$ of spin-up and spin-down electrons. In Fig. 2c, we show the occupation numbers $`f_{},f_{}`$ versus energy $`ϵ`$ for spin-up and spin-down Fermi distributions at $`T0`$, if $`\mathrm{\Delta }\mu <0`$. Quantum theory 9 shows that the energy of a spin wave in $`F_2`$ consists of an integer number $`n_m`$ of energy quanta $`\mathrm{}\omega `$ (magnons), where $`\omega `$ is the $`\stackrel{}{S}_2`$ precession frequency. Note that $`n_m`$ is directly related to the average $`\theta `$ angle, assumed to be small. The simplest collision process (Fig 3a), caused by s-d exchange and conserving angular momentum, is 10 such that an electron flips its spin $`\stackrel{}{s}`$ from down to up while one magnon of energy $`\mathrm{}\omega `$ and wavevector $`\stackrel{}{q}`$ is created. We show a complicated interaction vertex at the center of Fig. 3a, as a reminder that the simple Born approximation (Section I) is not valid for this process. Energy conservation gives $`ϵ_{}ϵ_{}=\mathrm{}\omega `$, where $`ϵ_{},ϵ_{}`$ are the final and initial electron energies in this so-called “spin-flip process”, also represented by an oblique arrow on Fig. 2c.

Because of the Pauli exclusion principle, such processes must start from a filled state and end in an empty state. Fig. 2c shows that the only suitable initial states for these processes are those located in a region of energy thickness $`|\mathrm{\Delta }\mu |\mathrm{}\omega =(\mathrm{\Delta }\mu +\mathrm{}\omega )`$ near the top of the spin-down Fermi sea. Thus, the number of possible initial states and, in turn, the resulting drive torque are proportional to $`\mathrm{\Delta }\mu +\mathrm{}\omega `$:

$$\stackrel{}{\tau }_d=\alpha _s(\theta )(\mathrm{\Delta }\mu +\mathrm{}\omega )\frac{\stackrel{}{S}_2\times (\stackrel{}{S}_1\times \stackrel{}{S}_2)}{\mathrm{}S_1S_2}.$$
(2)
The term $`\mathrm{\Delta }\mu `$ is proportional to the current I, and thus is consistent with Eq.(1). But the new $`\mathrm{}\omega `$ term represents 6 a surface damping torque independent of I, with $`\alpha _s(\theta )`$ as the corresponding dimensionless Gilbert parameter. It is found that $`\alpha _s`$ is inversely proportional to $`L_2^x`$. This surface damping is a consequence of the Pauli exclusion principle, and is not obtained if the spins $`\stackrel{}{S}_2`$ are assumed classical and static as in the Slonczewski theory 5. Although the present argument leading to Eq. (2) is based on a Fermi distribution with sharp Fermi level at $`T0`$ and with $`\mathrm{\Delta }\mu +\mathrm{}\omega <0`$ (Fig. 2c), Eq. (2) turns out to be valid 6 even at high temperature $`k_BT\mathrm{}\omega ,\mathrm{\Delta }\mu `$ where the Fermi level is very blurred, and even when $`\mathrm{\Delta }\mu +\mathrm{}\omega >0`$.

The current I flowing from $`F_1`$ to $`F_2`$ (Fig. 2a) may act 5,6 through a translation of the Fermi surfaces in momentum space. In ferromagnets, the spin-up and spin-down currents, as well as the corresponding Fermi-surface displacements (Fig. 3b), are different, leading to non-zero $`\stackrel{}{P}`$ and $`\mathrm{\Delta }\mu `$ at a given point of the Fermi surface, with $`|\mathrm{\Delta }\mu |10^4eV`$. We saw in Section II that the important electrons are those originating in $`F_1`$ and moving toward $`F_2`$. They are located on the $`k_x>0`$ half of the Fermi surface, and indicated by a filling of dots on Fig. 3b. By summing over all such $`k_x>0`$ states, one can calculate the total spin current $`\mathrm{\Sigma }\stackrel{}{P}(x=0)`$ and the average $`\mathrm{\Delta }\mu `$. Actually, only one component of $`\stackrel{}{P}`$, here called $`P_{x2}`$, is important. It is found 11 that spin current and $`\mathrm{\Delta }\mu `$ remain in a constant ratio as the current I is varied, or as the ratio of spin-up to spin-down conductivities is varied:

$$\frac{\mathrm{\Sigma }P_{x2}(x=0)}{\mathrm{\Delta }\mu sin\theta }=const..$$
(3)
This shows that the Slonczewski theory 5, based on $`\mathrm{\Sigma }P_{x2}`$, and the Berger theory 6, based on $`\mathrm{\Delta }\mu `$, are largely equivalent. The most significant difference between them is that, in their original form, the first one uses phenomenological bulk damping only, and the second surface damping only. This difference has observable consequences, which will be described in Section V and Fig. 4b.

There exists another kind of spin imbalance, associated with isotropic expansion or contraction of the Fermi surfaces (Fig. 3c), which goes back to Aronov and to Johnson and Silsbee 12. It is the same at all points of the Fermi surface. We denote it by $`\overline{\mathrm{\Delta }\mu }`$, and it can reach $`10^3eV`$. Like $`\mathrm{\Delta }\mu `$, it is related linearly to I, but the relation is more nonlocal. It exists only within a spin-diffusion length from interfaces. As before, only electrons on the $`k_x>0`$ half of the Fermi surface are really active in creating torques, and they are shown by a filling of dots in Fig. 3c. Again, the total spin-current component and average spin imbalance, now denoted by $`\mathrm{\Sigma }\overline{P}_{x2}`$ and $`\overline{\mathrm{\Delta }\mu }`$, can be calculated 11. Again, we find $`\mathrm{\Sigma }\overline{P}_{x2}/\overline{\mathrm{\Delta }\mu }=const.`$, and the constant has nearly the same value as before.

Numerical estimates 11, based on the assumption of one-dimensional current flow, suggest that $`\mathrm{\Sigma }\overline{P}_{x2}`$ and $`\overline{\mathrm{\Delta }\mu }`$ for Fermi-surface expansion/contraction are each 10 or 100 times larger than the corresponding quantities for translation. Thus, the expansion/contraction mechanism is dominant if the flow is one-dimensional, i.e., if current leads attached to the multilayer sample do not flare out to a larger diameter too close to the multilayer itself.

IV. CRITICAL CURRENT VERSUS FIELD

In Section II, we saw that the current-induced drive torque $`\stackrel{}{\tau }_d`$ on layer $`F_2`$ (Eq. 1) is directed outwards, i.e., tends to increase $`\theta `$, if $`I>0`$. On the other hand, bulk Gilbert damping, characterized by the Gilbert parameter $`\alpha _b`$, gives an inward torque:

$$\stackrel{}{\tau }_b=\alpha _b\frac{\stackrel{}{S}_2}{S_2}\times \frac{d\stackrel{}{S}_2}{dt}=\alpha _b\omega \frac{\stackrel{}{S}_2\times (\stackrel{}{S}_1\times \stackrel{}{S}_2)}{S_1S_2}.$$
(4)
Strictly speaking, the second equality only holds for circular precession in a field $`\stackrel{}{H}`$ normal to layers. The physical origin of this bulk damping will be discussed briefly in Section VII.

At small I or $`|\mathrm{\Delta }\mu |`$, damping dominates, so that $`\theta `$ actually decreases as time passes (Fig. 3d). However, if I is positive and exceeds a certain critical value $`I_c`$, the drive torque overcomes damping, and $`\theta `$ increases (Fig. 3e). Depending on the function $`g(\theta )`$ in Eq. (1), $`\theta `$ may stabilize at a constant value, or may go all the way to $`180^{}`$ and stay there. Detailed numerical as well as analytical calculations of dc current-induced precession and switching have been done 13 for various thin-film geometries and anisotropy cases. And even when the precessing spin distribution in $`F_2`$ becomes non-uniform 14.

The value of $`I_c`$ is obtained by writing $`\stackrel{}{\tau }_d+\stackrel{}{\tau }_b=0`$, where we take $`\stackrel{}{\tau }_d`$ from Eq. (2) rather than Eq. (1) to include surface damping. Using also Eq. (4), and writing $`\alpha _s=C_s/L_2^x`$ and $`\mathrm{\Delta }\mu =KI`$, where $`C_s`$ and K are positive constants, we obtain

$$I_c=\frac{1}{K}(1+\frac{L_2^x\alpha _b}{C_s})\mathrm{}\omega .$$
(5)
In the case of $`\stackrel{}{H}`$ normal to layers, with circular $`\stackrel{}{S}_2`$ precession, we have 9 a frequency $`\omega =\gamma (H_zM_s)`$ for the uniform mode. Here, $`M_s`$ is the saturation magnetization in SI units. There is an additional $`4\pi `$ factor in front of $`M_s`$ if cgs units are used. Eq. (5) shows that $`I_c`$ itself is proportional to $`\omega `$:

$$I_c\omega =\gamma (H_zM_s).$$
(6)
Here, $`\gamma `$ is the gyromagnetic ratio. Tsoi et al. 15,18 measured $`I_c`$ on a multilayer consisting of many alternating Co and Cu layers, normal to the field. A sharp silver point (Fig. 4a) contacting the sample was used to achieve large local current densities. Consistent with Eq. (6), they found $`I_c`$ to increase linearly with increasing $`H_z`$, though not quite with the predicted vertical intercept.

The case 9 of in-plane field $`\stackrel{}{H}`$ with $`\omega =\gamma (H_z(H_z+M_s))^{1/2}`$ is more common, as it does not require large fields to saturate the sample, but more complicated. Since $`\stackrel{}{S}_2`$ precession is now elliptical (Fig. 3f), the bulk damping torque $`\stackrel{}{\tau }_b`$ (Eq. 4) and the two parts of $`\stackrel{}{\tau }_d`$ vary through a precession cycle, in different manners. As a result, $`I_c`$ is not simply proportional to frequency $`\omega `$, but rather 16 linear in $`H_z`$:

$$I_cH_z+(M_s/2).$$
(7)
Again, an extra $`4\pi `$ factor should be put in front of $`M_s`$ if cgs units are used. The prediction of Eq. (7) is more or less consistent with $`I_c`$ measurements by Katine et al. 16 and Albert et al. 17 on Co/Cu/Co multilayers in in-plane fields, where $`I_c`$ is found to increase linearly with increasing $`H_z`$, with a positive vertical intercept.

Tsoi et al. 18 have succeeded in measuring the frequency $`\omega `$ of current-induced spin oscillations. At fixed $`H_z`$, they find that $`\omega `$ increases with increasing I above $`I_c`$. This is not surprising: Since $`I_c`$ increases with increasing $`\omega `$ (Eq. (6)), the only spin wave excited at $`II_c`$ is the one with lowest $`\omega `$, i.e., the uniform mode with zero wavenumber q (Fig. 3g). At $`I>I_c`$, other modes with larger $`\omega `$ and q may appear. Slonczewski 19 has even proposed that the excited modes at $`II_c`$ already have $`q>0`$, due to the radiation of spin waves away from the point contact; this is rather consistent with the measured $`\omega `$ at $`I=I_c`$ in Tsoi’s Fig. 3.

To measure $`\omega `$, Tsoi et al. adapted a technique well known to radio engineers under the name of “synchronous rectification” or of “direct conversion”. The sample is exposed to a microwave of known and adjustable frequency $`\omega _{ex}`$ (Fig. 4a). While varying the dc current I through the point contact, they measured the dc ohmic voltage V between point and sample. Due to nonlinear effects, a glitch appears 18 in dV/dI whenever $`\omega _{ex}`$ is equal to the frequency $`\omega `$ of current-induced spin precession.

V. CRITICAL CURRENT VERSUS LAYER THICKNESS

According to Eq.(5), $`I_c`$ is a linear function of $`F_2`$ thickness $`L_2^x`$, with positive intercept on the $`I_c`$ axis. This prediction is shown as the upper solid line in Fig. 4b. Actually, it should hold for any direction of $`\stackrel{}{H}`$.

$`I_c`$ has been measured versus $`L_2^x`$ by Albert et al. 20, for in-plane field. They found $`I_c`$ to vary linearly with $`L_2^x`$, but with no appreciable vertical intercept. If we were to ignore the surface-damping term $`\mathrm{}\omega `$ in Eq.(2), the constant first term would disappear from Eq. (5), bringing that equation into better agreement with the experimental findings above. This is shown as the lower solid line in Fig. 4b.

On the other hand, Urban et al. 21 have recently measured the Gilbert damping constant versus precessing-layer thickness in very thin Fe/Au/Fe/Au multilayers with single-crystal layers, using a transmission-resonance technique. They find a considerable contribution inversely proportional to thickness, in other words, a surface damping. Conceivably, the discrepancy between the two experiments may come from a non-magnetic “dead” atomic layer 4 reducing the effective $`L_2^x`$.

VI. NEXI TORQUES IN VERY THIN LAYERS.

We saw in Section II that the total current-induced torque on the spins $`\stackrel{}{S}_2`$ of layer $`F_2`$ is usually in the radial direction $`\stackrel{}{\tau }_d`$ of the $`\stackrel{}{S}_2`$ precession orbit of Fig. 2a. However, for very small $`F_2`$ thickness $`L_2^x<L_010nm`$, the spin $`\stackrel{}{s}`$ of an electron does not have enough time to precess appreciably around $`\stackrel{}{S}_2`$ while crossing $`F_2`$ at the Fermi velocity (Fig 2b). Rather, $`\stackrel{}{s}`$ remains largely in its initial direction parallel or antiparallel to $`\stackrel{}{S}_1`$. Therefore, the total torque is in the tangential direction $`\stackrel{}{\tau }_n`$ (Fig. 2a), like the torque of a current-induced exchange field of fixed direction parallel or antiparallel to $`\stackrel{}{S}_1`$. This effect was first proposed by Heide et al. 22, who called it non-equilibrium exchange interaction (NEXI), and expected it to exist at all layer thicknesses.

We have calculated the magnitude of $`\stackrel{}{\tau }_n`$ as a function of $`L_2^x`$, by the same methods 8 as for $`\stackrel{}{\tau }_d`$. The normalized torque per unit volume of $`F_2`$ is plotted versus $`k_NL_2^x`$ in Fig. 5, for $`k_{}/k_N=1.5`$ and $`k_{}/k_N=1.0`$. Here, $`k_N`$ is the Fermi wavenumber in layer N. The corresponding plot for $`\stackrel{}{\tau }_d`$ was in Fig. 2 of Ref. 8. For $`L_2^x<L_0`$, $`\stackrel{}{\tau }_n`$ is found to be comparable to $`\stackrel{}{\tau }_d`$. However, as expected, $`\stackrel{}{\tau }_n`$ vanishes at larger $`L_2^x`$ values, except for periodic oscillations around zero. The period of these oscillations is approximately $`2L_0=2\pi /(k_{}k_{})`$ on the $`L_2^x`$ scale.

VII. ANISOTROPIC S-D EXCHANGE IN THE BULK

We saw in Section I that ordinary s-d exchange is not active in the bulk. Anisotropic s-d exchange is 23 a combination of spin-orbit interaction and exchange, and does not conserve total angular momentum. It is probably responsible for most of the Gilbert damping (Eq. (4)) observed in bulk and single-film magnetic metals. This theoretical expectation is confirmed by the observation 24 of an increase of damping in Ni and Co at low temperature, of the kind predicted 23 on the basis of anisotropic s-d exchange.

VIII. ANALOGY WITH LASER DIODE

Because of a rather complete analogy with the semiconductor laser, we called 6 SWASER (Spin Wave Amplification by Stimulated Emission of Radiation) the current-induced spin instability described in the present paper. Like the rate of creation of photons in a laser, the rate of creation $`dn_m/dt`$ for magnons is 6 proportional to the number $`n_m`$ of magnons itself, and this is called stimulated emission. By opposition, spontaneous emission is independent of $`n_m`$, and is probably negligible 25 for magnons whenever $`\theta 0.1^{}`$. The energy gap $`\gamma H_z`$ in the spin-wave spectrum plays in our theory a role like that of the semiconductor band gap in a laser. Also, the spin imbalance $`\mathrm{\Delta }\mu `$ is the analog of the difference $`\mu _v\mu _c`$ between the Fermi levels 26 of the valence and conduction bands. The spin-wave damping proportional to $`\mathrm{}\omega `$ (see Section III) also exists for light waves in a laser 26. As a result, our condition 6 $`\mathrm{\Delta }\mu +\mathrm{}\omega =0`$ for the critical current has 26 a correspondent $`\mu _v\mu _c+\mathrm{}\omega =0`$ for the laser threshold. Finally, the $`NF_2`$ interface acts as a source of momentum 8 to bridge the momentum gap $`k_{}k_{}`$ between spin-up and spin-down Fermi surfaces. Similarly, solute atoms and lattice pores are used 27 to bridge the momentum gap between valence and conduction bands in lasers made of an indirect-gap semiconductor such as Si. One difference, however, is that spin waves are far more nonlinear than light waves.

The SWASER can be used 17 to reliably switch small magnetic elements, leading to possible applications for data storage. Also, it constitutes a very compact, tunable, microwave oscillator.

REFERENCES

1. C. Herring, in $`\underset{¯}{Magnetism}`$, edited by G.T. Rado and H. Suhl (Academic, New York, 1966), Vol. 4, p. 162.

2. L.L. Hirst, Phys. Rev. $`\underset{¯}{141}`$, 503 (1966); J. I. Kaplan, Phys. Rev. $`\underset{¯}{143}`$, 351 (1966); C. Herring, Phys. Rev. $`\underset{¯}{85}`$, 1003 (1952); $`\underset{¯}{87}`$, 60 (1952). Here, instead of the s-d exchange model, often all electrons are assumed itinerant, but a spin precession still happens in the common exchange field.

3. Ya.B. Bazaliy, B.A. Jones and S.C. Zhang, Phys. Rev. B $`\underset{¯}{57}`$, R3213 (1998).

4. W. Weber, S. Riesen and H.C. Siegmann, Science $`\underset{¯}{291}`$, 1015 (2001).

5. J.C. Slonczewski, J. Magn. Magn. Mater. $`\underset{¯}{159}`$, L1 (1996).

6. L. Berger, Phys. Rev. B $`\underset{¯}{54}`$, 9353 (1996).

7. L. Piraux, S. Dubois and A. Fert, J. Magn. Magn. Mater. $`\underset{¯}{159}`$, L287 (1996).

8. L. Berger, J. Appl. Phys. $`\underset{¯}{81}`$, 4880 (1997).

9. A.H. Morrish, $`\underset{¯}{ThePhysicalPrinciplesOfMagnetism}`$ (Wiley, New York, 1965). See p. 298 for magnons, p. 549 for Gilbert damping, and p. 545 for the frequency of the uniform mode.

10. A.H. Mitchell, Phys. Rev. $`\underset{¯}{105}`$, 1439 (1957).

11. L. Berger, J. Appl. Phys. $`\underset{¯}{89}`$, 5521 (2001).

12. A.G. Aronov, JETP Letters, $`\underset{¯}{24}`$, 32 (1976); M. Johnson and R.H. Silsbee, Phys. Rev. Lett. $`\underset{¯}{55}`$, 1790 (1985).

13. J.Z. Sun, Phys. Rev. B $`\underset{¯}{62}`$, 570 (2000); Ya.B. Bazaliy, B.A. Jones and S.-C. Zhang, J. Appl. Phys. $`\underset{¯}{89}`$, 6793 (2001); R.H. Koch and C. Heide, paper DF-01, Joint MMM-Intermag Conference, San Antonio, Texas, January 2001.

14. J. Miltat, G. Albuquerque, A. Thiaville and C. Vouille, J. Appl. Phys. $`\underset{¯}{89}`$, 6982 (2001).

15. M. Tsoi, A.G.M. Jansen, J. Bass, W.-C. Chiang, M. Seck, V. Tsoi and P. Wyder, Phys. Rev. Lett. $`\underset{¯}{80}`$, 4281 (1998); $`\underset{¯}{81}`$, 493(E) (1998).

16. J.A. Katine, F.J. Albert, R.A. Buhrman, E.B. Myers and D.C. Ralph, Phys. Rev. Lett. $`\underset{¯}{84}`$, 3149 (2000);

17. F.J. Albert, J.A. Katine, R.A. Buhrman and D.C. Ralph, Appl. Phys. Lett. $`\underset{¯}{77}`$, 3809 (2000).

18. M. Tsoi, A.G.M. Jansen, J. Bass, W. C. Chiang, V. Tsoi and P. Wyder, Nature $`\underset{¯}{406}`$, 46 (2000).

19. J.C. Slonczewski, J. Magn. Magn. Mater. $`\underset{¯}{195}`$, L261 (1999).

20. F.J. Albert, R.A. Buhrman, J. Katine and D.C. Ralph, paper ED-14, Joint MMM-Intermag Conference, San Antonio, Texas, January 2001.

21. R. Urban, G. Woltersdorf and B. Heinrich, submitted to Phys. Rev. Lett..

22. C. Heide, P.E. Zilberman and R.J. Elliott, Phys. Rev. B, $`\underset{¯}{63}`$, 64424 (2001).

23. V. Kambersky, Can. J. Phys. $`\underset{¯}{48}`$, 2906 (1970); V. Korenman and R.E. Prange, Phys. Rev. B $`\underset{¯}{6}`$, 2769 (1972); L. Berger, J. Phys. Chem. Phys. Solids $`\underset{¯}{38}`$, 1321 (1977).

24. S.M. Bhagat and P. Lubitz, Phys. Rev. B $`\underset{¯}{10}`$, 179 (1974), see Fig. 7; B. Heinrich, D.J. Meredith and J.F. Cochran, J. Appl. Phys. $`\underset{¯}{50}`$, 7726 (1979).

25. L. Berger, Phys. Rev. B $`\underset{¯}{59}`$, 11465 (1999).

26. M.G.A. Bernard and G. Duraffourg, Phys. Stat. Solidi, $`\underset{¯}{1}`$, 699 (1961).

27. C.B. Duke and N. Holonyak, Physics Today, December 1973, p. 23; R.T. Collins, P.M. Fauchet and M.A. Tischler, Physics Today, January 1997, p.24.

FIGURE CAPTIONS

FIG. 1. a) Spin configuration assumed in the Born approximation, with $`\stackrel{}{s}`$ nearly static and antiparallel to $`\stackrel{}{H}`$. b) More correct picture with $`\stackrel{}{s}`$ precessing rapidly around $`\stackrel{}{S}_2`$, and following it closely. c) Observation of $`\stackrel{}{s}`$ precession by an angle $`ϵ`$ around a fixed $`\stackrel{}{S}`$, in Fe, Co, Ni films (Ref. 4).

FIG. 2. a) Multilayer with magnetic layers $`F_1,F_2`$ and nonmagnetic $`N,N_2`$. The oblique arrow represents an electron crossing from $`F_1`$ to $`F_2`$. b) Two electrons propagating inside $`F_2`$ in different directions have different $`\stackrel{}{s}`$ precession phases around $`\stackrel{}{S}_2`$ at a given $`x>0`$. c) Fermi distributions $`f_{}(ϵ),f_{}(ϵ)`$ for spin-up and spin-down electrons of energy $`ϵ`$, with different Fermi levels $`\mu _{},\mu _{}`$, at $`T0`$. The oblique arrow represents a spin-flip process where a magnon of energy $`\mathrm{}\omega `$ is emitted.

FIG. 3. a) Spin-flip process where a magnon of energy $`\mathrm{}\omega `$ and wavevector $`\stackrel{}{q}`$ is created. b) Current-induced translations of spin-up and spin-down Fermi surfaces. c) Current-induced expansion or contraction of spin-up and spin-down Fermi surfaces. d) Precession orbit of spin $`\stackrel{}{S}_2`$ if $`I<I_c`$. e) Orbit if $`I>I_c`$. f) Orbit if $`I>I_c`$ and $`\stackrel{}{H}`$ is in-plane. g) Dispersion relation $`\omega (q)`$ for spin waves in a field $`H_z`$. Horizontal lines show the maximum $`\omega `$ value of current-induced spin waves.

FIG. 4. a) Multilayer with many magnetic layers and a point contact, used by Tsoi et al. (Refs. 15 and 18). b) The upper solid line represents the predicted dependence of $`I_c`$ on $`F_2`$ thickness $`L_2^x`$ if both surface and bulk dampings are introduced. The lower solid line applies if only bulk damping is introduced. The crosses represent schematically the experimental data of Albert et al. (Ref. 20).

FIG. 5. Normalized tangential torque component $`\tau _n`$ per unit volume of $`F_2`$, plotted versus normalized layer thickness $`L_2^x`$, for fixed I, $`\theta `$ and $`\omega `$. The point corresponding to $`L_2^x=L_0`$ is marked on the horizontal scale.


text:
# Possible Recovery of SN 1961V In Hubble Space Telescope Archival Images1footnote 11footnote 1Based on observations made with the NASA/ESA Hubble Space Telescope, obtained from the data archive of the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555.

## 1. Introduction

Supernova 1961V in NGC 1058 is an unusual object among the nearly 2100 historical supernovae (SNe), and its nature is unclear. Discovered in 1961 July by P. Wild, SN 1961V had perhaps the most bizarre light curve ever recorded for a SN, with pre-maximum and post-maximum plateaus of $`0.5`$ yr, a sustained three-year plateau at $`m_{\mathrm{pg}}18.5`$ mag, and a subsequent decline to $`m_{\mathrm{pg}}21`$–22 mag (see Humphreys & Davidson 1994; Humphreys, Davidson, & Smith 1999). Spectroscopically, it is one of very few members of the archaic Zwicky (1964, 1965) “Type V” SN classification. Its spectrum, dominated by narrow emission lines of H, He I, and Fe II, suggested a maximum expansion velocity of only $`2000`$ km s-1. Along with only SN 1987A in the Large Magellanic Cloud (LMC; e.g., Gilmozzi et al. 1987; Sonneborn, Altner, & Kirshner 1987), SN 1993J in M81 (Aldering, Humphreys, & Richmond 1994; Cohen, Darling, & Porter 1995), and SN 1997bs in M66 (Van Dyk et al. 2000), SN 1961V has a possible progenitor star identified on pre-SN images.

The progenitor of SN 1961V may have been one of the most extraordinary stars known. Identified as a $`m_{\mathrm{pg}}18`$ mag star in many photographs of NGC 1058 prior to the explosion (Bertola 1964; Zwicky 1964; Klemola 1986), it was extremely luminous, at $`M_{\mathrm{pg}}^012`$ mag for a distance of 9.3 Mpc (determined using Hubble Space Telescope \[HST\] observations of Cepheids; Silbermann et al. 1996), leading to an initial (almost certainly erroneous) theoretical mass estimate of $`2000`$ $`M_{\mathrm{}}`$ (Utrobin 1984). Nonetheless, with somewhat more likely mass estimates of $``$100–200 $`M_{\mathrm{}}`$ (Goodrich et al. 1989), it was among the most luminous and massive known individual stars in any galaxy.

Based on an optical spectroscopic detection of SN 1961V in 1986 and reanalysis of older data, Goodrich et al. (1989) postulate that SN 1961V was not a genuine supernova (defined as the explosion of a star at the end of its life), but rather was the giant eruption of a luminous blue variable (LBV) star, much like an exaggerated version of $`\eta `$ Carinae (see Davidson & Humphreys 1997 for a review of $`\eta `$ Car). Goodrich et al. suggest that SN 1961V actually survived the eruption, and that the progenitor faded because of the formation of optically thick dust in the ejecta. They further predict that at the site of SN 1961V should be a quiescent S Dor-type, hot Of/WN star enshrouded by dust and possibly a dense wind; if its circumstellar extinction $`A_V5`$ mag and its bolometric correction BC $`4`$ mag, the star should now have $`V27`$ mag. Alternatively, if $`A_V4`$ mag, as in the case of $`\eta `$ Car, and BC $`3`$ mag, then $`V25`$–26 mag. Finally, $`V`$ could be as bright as 22–23 mag if the star has $`A_V=4`$–5 mag and an active, optically thick wind with BC $`0`$ mag.

Motivated by this prediction, Filippenko et al. (1995; hereafter F95) attempted to recover the surviving star via multi-band imaging using HST, utilizing the radio position obtained by Cowan, Henry, & Branch (1988). Their pre-refurbishment WFPC $`VRI`$ (flight system filters F555W, F702W, and F785LP) images, accompanied by precise astrometry obtained from deep ground-based optical data, led them to associate SN 1961V with one star in a cluster of three stars, all within $`2\mathrm{}`$ ($`90`$ pc) of each other, as the possible eruption survivor. Their most likely candidate (“Object 6”), shown in Figure 1 of the present paper222This is Figure 2 in F95; it is given again here because it came out so poorly in F95, both in the original version and in the Erratum., has $`V25.6`$, $`VR1.0`$, and $`RI0.9`$ mag, roughly the colors of a mid-K-type supergiant; if dereddened by reasonable amounts of interstellar and circumstellar reddening, the colors are consistent with those of O-type stars. F95 point out that unaberrated images, particularly in bluer bands, are required to better isolate the possible survivor.

Based on its fading radio emission recently detected in new observations using the Very Large Array (VLA), Stockdale et al. (2001) argue that the current radio properties of SN 1961V are consistent with those of some “peculiar” Type II radio SNe. In particular, they compare SN 1961V with SN 1986J in NGC 891 (Weiler, Panagia, & Sramek 1990), which was originally classified as Type V (Rupen et al. 1987), but is generally considered to be a prototypical Type IIn SN (see Schlegel 1990 and Filippenko 1997 for discussion of this SN subtype). Stockdale et al. contend that it is less likely, based on its radio properties, that SN 1961V was similar to $`\eta `$ Car or other LBVs. (However, there is a lack of radio observations of bright LBVs with ages similar to that of SN 1961V, preventing a secure determination of the true nature of SN 1961V.) In other words, they argue there should be no survivor, only a very old SN or young SN remnant emerging forty years after explosion.

In this paper we further investigate whether SN 1961V survived a LBV superoutburst by exploiting archival HST WFPC2 images that contain the SN site. We recover “Object 6” of F95, and we show that an additional faint and very red object is consistent with the radio position of SN 1961V and may represent the possible surviving star. We consider it less likely to be a fading, very old SN, but additional observations are necessary to further clarify this issue.

## 2. Analysis of Archival HST WFPC2 Images

Two HST programs, GO-5446 and GO-9042, have imaged NGC 1058, the host galaxy of SN 1961V; the datasets are publicly available in the HST archive. The former program obtained on 1994 September 8 a pair of 80-s exposures with the F606W filter, while the latter program obtained on 2001 July 3 two pairs of 230-s exposures, one pair with the F450W filter and the other with the F814W filter. These bandpasses roughly correspond to $`V`$, $`B`$, and $`I`$. Using the radio position of SN 1961V from Stockdale et al. (2001; it is coincident with the position measured by Cowan et al. 1988, to within the errors), and correcting for geometric distortions on the WF3 chip, we isolate the site of the SN for both sets of observations.

Figures 2 and 3 show the SN site on the F814W and F450W images, respectively. An error circle is plotted at the radio position, based on the astrometric information in the headers of the F450W and F814W images, with a (conservative) radius of $`1\stackrel{}{\mathrm{.}}5`$ (the observations by GO-9042 employed the fine-lock mode, which results in about this level of astrometric uncertainty; see the discussions of HST astrometry in F95 and Van Dyk et al. 1999b). Figure 4 shows the site in the F606W image; this exposure was obtained in “gyro” mode and likely has less reliable astrometry than the F814W and F450W images, so we impose the error circle of Figures 2 and 3 onto Figure 4. The source numbering is from Figure 1, as in F95. Object 8 is situated off the F450W and F814W images, but is detected on the F606W image.

To measure the brightnesses of the various sources through each of the three bands, we employed version 1.1 of the package HSTphot (Dolphin 2000a,b). We followed the “recipe” in the HSTphot online333http://www.noao.edu/staff/dolphin/hstphot/ manual, initially adopting a $`4\sigma `$ detection threshold and using the following tasks in sequential order: mask, crmask, coadd, getsky, hotpixels, and hstphot. As Dolphin (2000b) has shown, HSTphot produces results that are quite consistent with those obtained from both the DAOPHOT and DoPHOT packages, while accounting for WFPC2 point-spread function variations and charge-transfer effects across the chips, zeropoints, aperture corrections, etc., automatically within one package. (We have also conducted limited tests of HSTphot vs. DAOPHOT and find very good agreement in the results.) Table 1 gives the results of our photometry in the flight system bandpasses. Not all sources were detected in all three bands, given the relatively low signal-to-noise ratio (S/N) of these images. Detection limits ($`3\sigma `$) are $`m_{\mathrm{F450W}}25.3`$, $`m_{\mathrm{F606W}}25.4`$, and $`m_{\mathrm{F814W}}25.3`$ mag.

In order to further analyze the archival images, we have transformed our flight system magnitudes into Johnson-Cousins $`BVI`$ magnitudes, in an analogous manner to that of F95, i.e., via synthetic photometry of normal stars of a wide span of spectral types and luminosity classes, obtained by applying the STSDAS package SYNPHOT to the Bruzual Spectral Synthetic Atlas. In Table 1 we list the $`BVI`$ magnitudes, wherever possible, for the objects in the SN 1961V environment. In two cases, we used the detection limit at F606W to set a lower limit on the $`I`$ magnitude; for Object 8, we could only approximate the transformation.

Despite the passband differences between the two studies (F606W and F814W here vs. F555W and F785LP in F95), the $`V`$ and $`I`$ magnitudes of Objects 1–10 agree fairly well overall. Our $`V`$-band magnitudes are brighter by $`0.08`$ mag, but with a rather large dispersion, 0.34 mag. Our $`I`$-band magnitudes are fainter by 0.04 mag, with a dispersion of 0.20 mag. The average errors in our $`V`$ and $`I`$ magnitudes are 0.18 and 0.17 mag, respectively; for the measurements in F95, they are 0.14 and 0.20 mag, respectively. Thus, the $`I`$-band magnitudes agree to within the errors, but the dispersion in our $`V`$ magnitudes is about a factor of two larger than the uncertainties in either study. We cannot account for this discrepancy entirely through differences in bandpass or S/N; it may arise from underestimates of the uncertainty in the aberrated WFPC photometry. For example, the largest disagreements are for Object 5 (which is clearly extended), Object 7 (which is in a crowded environment with a high background), and Object 10 (which also may well be extended).

In Figures 5 and 6 we show the color-magnitude diagrams for the objects in the SN 1961V environment. As pointed out by F95, SN 1961V appears to have occurred in a region containing several massive supergiant stars. Objects 4, 7, and 9 are consistent with being early-type supergiants. Objects 1 and 10 are also of early type and are possible supergiants, but we cannot be certain because substantial disagreement exists between their respective $`BV`$ and $`VI`$ colors. Object 5 is consistent with either an intermediate-type supergiant or an early-type supergiant reddened by $`A_V0.8`$ mag. We have no color information for Object 8, but F95 identified it as an intermediate-type supergiant. Both Objects 2 and 6 have quite similar colors and brightnesses, consistent with late-type supergiants. An additional faint, red object is detected only in the F814W image, just to the southeast of Object 6, and was not identified by F95. We list this in Table 1 as “Object 11” and represent it with limits in Figure 6 (its $`I`$-band brightness and $`VI`$ color limits are based on the 3$`\sigma `$ detection limit of the F606W image).

Figure 3 also hints at an underlying population of fainter blue stars, extending from Object 5 northeast to Object 8, and to the west around Object 10. Such a background of blue stars may have also been detected by Utrobin (1987) and Goodrich et al. (1989). Figure 4 hints at extended nebular emission, primarily H$`\alpha `$ falling within the F606W bandpass, that appears to trace these blue stars. The S/N in both the F450W and F606W images is simply not sufficient to better detect this population and nebulosity.

## 3. Discussion

As can be seen in Figure 2, Objects 2, 5, 6, and 11 are within the positional error circle for SN 1961V and, therefore, can be considered viable candidates for its survivor or remnant. Object 5 is quite diffuse and fairly blue; barely visible in Figure 2, it is not detected as a source in the F814W image. F95 identify it as a possible O6 I–III star, although it appears here as a possible supergiant intermediate in type, or perhaps of early type and somewhat reddened. Even if SN 1961V were a fading Type II SN, its color should be fairly red, even without local interstellar or circumstellar reddening, since its optical spectral energy distribution should be dominated by H$`\alpha `$ emission (e.g., Long, Blair, & Krzeminski 1989; Fesen 1993; Chevalier & Fransson 1994). We therefore do not further consider Object 5 as a possible candidate.

Both Objects 2 and 6 are quite red, so they are more consistent with our expectations for either an old SN or reddened LBV star. They are roughly equidistant from the radio position (Cowan et al. 1988; Stockdale et al. 2001) on the WFPC2 images. However, central to the argument by F95 that SN 1961V survived the superoutburst as Object 6 is its near-coincidence with the radio position projected onto the WFPC images. This projection is based on the application of an independent, accurate astrometric grid, established from deep ground-based imaging, to the WFPC chip on which the SN environment lies, with the knowledge that the relative astrometry within a given chip is excellent. This projection falls within our error circle and may, in fact, better represent the actual SN position on the HST images. Object 2 is more distant from this position; hence, we no longer consider it to be a candidate. Furthermore, to within the errors, both Objects 2 and 6 have not changed in brightness or color between the WFPC and WFPC2 exposures. It is possible that they are both LBVs, reddened by similar amounts and remaining at relatively the same brightness; it is more likely, however, that both Objects 2 and 6 are similar late-type supergiants. We therefore also now consider it less likely that Object 6 is the SN 1961V survivor.

Another candidate for the survivor or very young remnant of SN 1961V may be Object 11. It is faint ($`I24.3`$ mag), but also quite red ($`BI1.0`$ mag and $`VI1.1`$ mag). Unlike the possibly early stellar-type Object 10, which is actually fainter in $`I`$, Object 11 is not detected in the F606W and F450W images. From the independent astrometric grid, the radio position actually lies just to the southeast, within a few pixels of Object 6 (Fig. 1), and is certainly within the WFPC2 error circle. However, it is virtually coincident with Object 11. A hint of this source can be seen in the deconvolved F785LP image (Fig. 1). It is seen in each of the raw F814W images before cosmic-ray rejection and coaddition and is detected by HSTphot in the coadded image at the 5$`\sigma `$ level.

Could the brightness and color of Object 11 be consistent with a decades-old SN? If SN 1961V is currently an old radio SN, as Stockdale et al. (2001) suggest, its optical luminosity, particularly at H$`\alpha `$, might still be relatively high, since late-time SN radio and optical emission appear to be correlated (e.g., Chevalier & Fransson 1994). SN 1961V has not dropped off precipitously in radio brightness at very late times, in contrast with SN 1957D (Cowan, Roberts, & Branch 1994; Long, Winkler, & Blair 1992) or SN 1980K (Montes et al. 1998), so we also might not expect the optical emission to have dropped off completely. Therefore, an old SN 1961V could be bright enough to detect with HST, and may be fairly red.

We can explore this further by comparing SN 1961V with two old Type II radio SNe which are seen in HST images. SN 1986J is in publicly available HST archival images, obtained by GO-9042 on 2001 July 3 with the F450W and F814W filters (see Van Dyk et al. 1999b, their Figure 4, for a finder chart). Again, applying HSTphot to these images containing SN 1986J, we obtain $`m_{\mathrm{F450W}}=23.13\pm 0.06`$ ($`B23.4`$ mag) and $`m_{\mathrm{F814W}}=20.39\pm 0.01`$ ($`I20.4`$ mag), corresponding to a color $`BI3.0`$ mag. We can also compare with the Type II-L SN 1979C in M100, which Van Dyk et al. (1999a) detect in HST images from 1996. We transform the brightness of SN 1979C in the F439W, F555W, and F814W bands to $`B23.2`$ mag, $`V22.1`$ mag, and $`I21.0`$ mag; hence, $`BI2.2`$ and $`VI1.1`$ mag. Both SNe are quite red, but the brightness of SN 1986J is evolving much more slowly than that of SN 1979C (compare $`I20.4`$ mag for SN 1986J with Gunn $`i20.2`$ mag in 1986; Rupen et al. 1987).

The host galaxies for SNe 1961V and 1986J are in the same group, at a distance of 9.3 Mpc (Silbermann et al. 1996). If we assume $`A_V=1.5`$ mag for SN 1986J (Leibundgut et al. 1991), we find that $`M_I^010.3`$ mag. For SN 1979C, at 16.1 Mpc and $`A_V0.5`$ mag (Van Dyk et al. 1999a), $`M_I^010.5`$ mag. (If either SN is experiencing higher intrinsic circumstellar extinction, then these absolute magnitudes will be even brighter.) In contrast, for Object 11, assuming circumstellar $`A_V4`$–5 mag, we derive $`M_I^07.9`$ to $`8.5`$ mag; assuming less extinction, of course, means that Object 11 would be intrinsically fainter.

Making a direct comparison of SN 1961V in the optical with these old radio SNe, of course, is difficult, since both SN 1986J and SN 1979C are about 20 years younger than SN 1961V 444Although not discovered until 1986, the explosion date of SN 1986J was possibly sometime between 1978 and 1983 (Rupen et al. 1987).; they could still significantly fade over the next two decades. Also, SNe II-L and IIn behave differently in both the optical and the radio. Unfortunately, no optical data on old radio SNe as old as SN 1961V are available, other than the ground-based emission-line data for SN 1957D from 1987 through 1988 (Long et al. 1989). Nonetheless, although SNe 1986J and 1979C are indeed quite red, they both appear to be intrinsically brighter than Object 11, making it less likely that Object 11 is an old SN.

Another distinct possibility is that no positional coincidence actually exists — SN 1961V is not Object 11 or any of the other objects detected in the WFPC2 images. This implies that SN 1961V may have faded below detectability in all of these bands (specifically, it must have $`I25.3`$ mag). If SN 1961V is an old SN, then it would have intrinsically faded or circumstellar extinction would have significantly increased since the 1986 spectroscopic observations by Goodrich et al. (1989). (In fact, if this is true, then the SN must also have faded significantly between 1986 and the 1991 WFPC observations by F95.) Even if SN 1961V is an $`\eta `$ Car-like LBV that has survived a superoutburst, it also may no longer be detectable: A similar, younger case of a possible $`\eta `$ Car-like LBV, SN 1997bs in M66 (Van Dyk et al. 2000), has faded considerably and become bluer in recent HST Snapshot images (Li et al. 2002), rather than redder and brighter, as predicted for increasing amounts of dust.

## 4. Conclusions

Ultimately, we conclude that SN 1961V is most likely Object 11, and that it has survived its eruption. The brightness and color limits for Object 11, shown in Figure 6, are consistent with an early-type supergiant star with $`A_V1`$ mag — specifically, the $`V`$ magnitude ($`27`$–25 mag) predicted by Goodrich et al. (1989) for circumstellar $`A_V4`$–5 mag and BC $`3`$ to $`4`$ mag of a quiescent LBV surrounded by optically thick dust and a stellar wind. Its current absolute brightness is more consistent with a low-luminosity LBV, such as R71 in the LMC (Humphreys & Davidson 1994), possibly in quiescence, than with a decades-old “peculiar” Type II SN.

Clearly, full resolution of the nature of SN 1961V still eludes us. As Stockdale et al. (2001) point out, further multi-wavelength observations are needed. Ideally, one might obtain spectra of the candidates: Gruendl et al. (2002) attempted to detect the SN with the echelle spectrograph on the Kitt Peak 4-m telescope, but failed to do so. They conclude that they simply may have misplaced the spectrograph slit. However, if SN 1961V is indeed Object 11, an alternative explanation is that its relative faintness and redness may have made detection very difficult.

An upcoming HST Cycle 11 program, GO-9371 (PI: Y.-H. Chu), will obtain spectra of the SN 1961V environment within $`1^{\prime \prime }`$ of the SN position using STIS, to find signatures of LBV ejecta nebulae, old SNe, and SN remnants. They will also obtain WFPC2 $`VRI`$ images. If the slit is centered on Object 11, and the exposures are of sufficient depth, these spectra might provide the definitive answer. If the broad-band imaging is also sufficiently deep and at the highest possible resolution, the brightnesses and colors of the candidate and other stars in the environment could be more accurately measured. SN 1961V is also a GTO target using the newly installed ACS at the F475W, F625W, and F775W bands, complementing the GO program.

In addition, deep, high-resolution, near-infrared imaging may be necessary, since we might expect old, cooling SNe and dust-enshrouded hot LBVs to have different infrared colors. For example, $`\eta `$ Car has $`JH0.9`$ mag and $`HK1.5`$ mag (Whitelock et al. 1983), while at late times ($``$100 d) the Type IIn SN 1998S (also a radio SN) has $`JH0.4`$ mag and $`HK0.8`$ mag due to dust in the circumstellar medium (Fassia et al. 2000). Unfortunately, published data are not available for any SN IIn (and, with the exception of SN 1987A, for any other SN of any type) in the infrared at very late times (i.e., several years to decades after explosion), but we might expect that these SNe get redder through additional dust formation. Although this is not a rigorous test, the infrared detection of SN 1961V would also help determine its nature. The red color of Object 11 implies that it should be a relatively bright infrared source for detection with a large-aperture telescope under excellent seeing conditions, or with HST/NICMOS and SIRTF.

Finally, as an important aside, if SN 1961V, SN 1986J, the old “Type V” SN classification, and Type IIn SNe bear any connection, and if the progenitor of SN 1961V did survive, then the “explosion” mechanism is not core collapse in all SNe IIn. That is, the SN IIn subclass, which is already known to span a very broad range of properties (e.g., Filippenko 1997), may include a number of SN “impostors.” The occurrence of SN 1961V and its possible cousins \[SN 1954J (Smith, Humphreys, & Gehrz 2001), SN 1997bs (Van Dyk et al. 2000), SN 1999bw (Filippenko, Li, & Modjaz 1999), and SN 2000ch (Filippenko 2000)\], and their apparent resemblance to $`\eta `$ Car and LBV superoutbursts, add an intriguing twist to the evolution of very massive stars.

We thank A. J. Barth for producing Figure 1, which was originally published by Filippenko et al. (1995). The work of A.V.F.’s group at UC Berkeley is supported by NSF grant AST-9987438, as well as by NASA grant AR-08754 from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. A.V.F. is grateful for a Guggenheim Foundation Fellowship.