text:
# Application of interactive parallel visualization for commodity-based clusters using visualization APIs

## 1 Introduction

We are interested in large scale high-performance scientific parallel visualization with real-time interaction. More precisely, we want to support the interactive extraction of insight from large sets of data. Also, in order not to miss important details in the data, we would like to develop high resolution visualization. Our interest is motivated by today’s technological advances, which stimulate and facilitate the development of increasingly complicated mathematical models. Visualization with real-time interaction is essential for better analysis and understanding of the masses of data produced by such models, or by devices such as Magnetic Resonance Imaging (MRI) and laser-microscope scanners. Some of the applications that we are interested in and work on are given on Figure 1. Others come from high-energy physics, climate modeling, etc. Currently the size of the data sets that we are dealing with is around $`1`$ GB.

We develop an inexpensive and, as the results show, efficient application. The approach is as follows:

* Do remote parallel rendering to a large tiled display (see for example Figure 1 where we show displays composed of $`4`$ tiles).
* Use commodity-based clusters connected with high speed network.
* Extend and combine already existing APIs such as Chromium, Open Inventor, VTK, etc.

The APIs considered are open source. Open Inventor (see ) is a library of C++ objects and methods for building interactive 3D graphics applications, VTK (see ) is another widely used API for visualization and image processing, and Chromium (see ) is an OpenGL  interface for cluster visualization to a tiled display. Chromium is based on WireGL  and provides a scalable display technology. The latest Chromium paper  reports on the implementation of components (stream processing units, or SPUs) that can be used in sort-last parallel graphics applications. More on the sort-last parallel rendering approach can be found in  and .

There are feasible alternatives to our approach. For example, the tiled display can be replaced with a high resolution flat panel driven by a single workstation, the visualization clusters can be replaced with fast sequential machines, and the process of extending and combining already existing APIs can be replaced with building the entire visualization system from scratch. The IBM’s T221 display with its maximum resolution of $`3840\times 2400`$ brings to the scientist’s desktop a resolution that is high enough for many applications. For such cases this is a preferred alternative and we have it as an option in our implementation. Similarly, for “small” data sets we prefer sequential processing. Extending and combining visualization APIs is a popular and in many cases preferred development approach since one can easily leverage already existing and powerful APIs. Examples are the extension of VTK to ParaView (see ) and parallel VTK (see ), VisIt (see ), etc. Building an entire visualization system from scratch may be efficient for very specific requirements, but in general is expensive and time consuming (see for example NASA’s long term project ParVox ). Our goals, apart from developing the application outlined above, also include studying and identifying the singularities of the Chromium API. We concentrate on Chromium’s sort-first rendering, a model inherited, and further developed from the WireGL API.

The parallel model that we consider is MIMD (multiple instructions – multiple data). MIMD parallel models are common in the design of parallel visualization software for clusters of workstations. A fundamental framework is when cluster nodes process separate parts of a global scene and their output is composited and rendered to a tiled display. Providing MIMD model visualization systems with efficient user interaction has become a task of great interest. See for example . The continuous interest also prompted the development of a new cluster rendering utility toolkit (CRUT) which was reported in . CRUT is a glu-like toolkit that will facilitate the development of user interaction within Chromium, and hence the development of applications and visualization API extensions like the one that we developed. We pursue a Master/Slave paradigm: we declare one of the cluster nodes as GUI Master, use the Master to intercept the sequential user input, and broadcast that input in a user defined protocol to the other nodes (here called Slave nodes). This yields a system where the GUI Master sends to the Slave nodes instructions of “how and when” to redraw their part of the scene. The user interface in this case is part of the visualization software. Another implementation is to have the GUI Master separate from the visualization software. The interaction in this case would be through a GUI window that will simulate a “parallel interaction device” that takes a sequential user input and broadcasts it to the cluster nodes (through sockets). The details are given in the following sections.

The article is organized as follows. In Section 2 we describe our model framework and its implementation in extending Open Inventor and VTK. Next (in Section 3) we discuss issues that are important for the development of high performance visualization on clusters of workstations. Also, we provide some of our performance results. The last section (Section 4) summarizes the results of this work.

## 2 Parallel GUI model and implementation

Parallel GUI for a MIMD parallel visualization model requires the presence of a Master node that will synchronize with the other nodes the redrawing of the consecutive scenes according to a sequential user input. We consider two variations of parallel GUI. Both have similar visualization pipelines (see Figure 2). The parallel GUI (ParaMouse) gets the sequential user input (through mouse, keyboard, etc.) and broadcasts it to the OpenGL applications through sockets. Every application is visualizing a separate part of a global scene. The applications are bound together by MPI communications. The OpenGL calls that the applications make are intercepted by Chromium and sent through the network to the visualization servers, which composite the input and render it to a tiled display.

Chromium’s parallelization model (denoted in the article by CPM) is represented with the pseudo-code:

```

  glXMakeCurrent(getDisplay(), getNormalWindow(),
                 getNormalContext());
  if (clearFlag)
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
  glBarrierExecCR(MASTER_BARRIER);

  do sequential OpenGL rendering of the local scene

  glBarrierExecCR(MASTER_BARRIER);
  if (swapFlag)
    glXSwapBuffers(getDisplay(), getNormalWindow());
  else
    glXSwapBuffers(getDisplay(), CR_SUPPRESS_SWAP_BIT);

```

where the display obtained is composited if clearFlag is $`1`$ for all processors, swapFlag is $`1`$ for processor with rank $`0`$, and $`0`$ for the rest, and tiled if clearFlag and swapFlag are $`1`$ only for processor with rank $`0`$.

To extend Open Inventor and VTK (or any API) to support interactive cluster visualization within the above framework we did the following:

* Apply the CPM to the APIs rendering method(s).
* Implement the ParaMouse parallel GUI or extend the API’s GUI using the Master-Slave concept.

To implement the first step for Open Inventor we extended the SoXtRenderArea::redraw() method. For VTK this step is implemented by David Thompson, Sandia National Laboratory . We implemented in Open Inventor the Master-Slave concept by extending the GUI that ivview provides. Function main in ivview was extended by implementing the pseudo-code:

```

  declare processor with rank 0 as Master;
  if (Master)
    run GUI as implemented in ivview;
  else
    listen for instructions from the Master;

```

We extended the callback functions triggered by the devices that Open Inventor’s GUI handles (mouse, keyboard, etc). The callback functions responding to user input have MPI\_Isend of the invoked user input in the Master node to the Slave nodes. The Slave nodes are in “listening” mode (MPI\_Recv or MPI\_Irecv while animating), which is given with the pseudo-code below, and upon receiving data, representing instructions in an internally defined protocol, they call the action for the invoked input and redraw their part of the global scene (if necessary).

```

  // Initialize Inventor
  SoDB::init();
  SoNodeKit::init();
  SoInteraction::init();

  // Build the Inventor’s objects and scene graphs
  SoSeparator *root = new SoSeparator;
  root->ref();
  readScene(root, files);

  // Create a Slave node ExaminerViewer
  SoExaminerViewer *viewer = new SoExaminerViewer(crrank);
  viewer->setSceneGraph(root);

  // Chromium initialization
  crctx = crCreateContextCR(0x0, visual);
  crMakeCurrentCR(crwindow, crctx);
  glBarrierCreateCR( MASTER_BARRIER, crsize);

  glEnable(GL_DEPTH_TEST);
  viewer->mainLoop();

```

The SoExaminerViewer class is based on the Open Inventor’s SoXtExaminerViewer class. The difference is that SoExaminerViewer does not have Xt window interface function calls and the rendering is with the CPM.

The user interface is through the ivview window (see Figure 3, left), which is blank and used only for the user input. The scene is drawn in separate windows/tiles. The example from Figure 3 (right) shows a display with $`4`$ tiles. There are $`4`$ applications running, each of which visualize a sphere.

In VTK’s Master-Slave model we extend the vtkXRenderWindowInteractor class by:

* Implementing “pure” X Windows GUI (no GL calls).
* Adding MPI send/receive calls for mouse and keyboard events.

For the ParaMouse interface we added in VTK a new interactor style,
called vtkInteractorStylePMouse.

## 3 Observations and performance results

We did our implementation and testing on a Beowulf cluster with $`4`$ nodes, each node with $`2`$ Pentium III processors, running at 1 GHz. Every node has Quadro2 Pro graphics card. The nodes are connected into a local area network with communications running through $`100`$ Mbit/sec fast Ethernet or $`1`$ Gbit/sec fiber optic network. More about the general performance of this particular machine can be found in .

The experience that we had in developing interactive parallel visualization is summarized as follows:

1. For the parallel model considered performance is problem and interaction specific.
For example, depending on the data and the user interaction, the entire global scene may be mapped to a single tile of the display. Also, the Chromium bucketing strategy  will not work for scenes composited of non-localized consecutive polygons.
2. The GUI communications time is negligible compared to the visualization time, which leads us to the scalability results reported in . This statement is supported by our next observation.
3. We can send approximately $`5,000`$ “small” (less then 100 bytes) messages/sec (see ).
4. Chromium automatically minimizes data flow, except for geometry flow (see below).
5. It is advisable to keep small static scenes in display lists (see below).

The network data flow is usually the bottleneck in the visualization pipeline that we consider. Chromium provides several techniques to minimize it. They are: simplified network protocol, bucketing, and state tracking (see ). None of these techniques however are intended for the automatic minimization of the geometry flow, which usually is the most expensive component. For example in animation the same objects are drawn without (or with minimal) change from frame to frame. Nevertheless the scene is transmitted over the network for every frame, creating an enormous bottleneck, as shown in the next example. Data flow minimization techniques that exploit spacial or temporal coherence of consecutive frames (see for example ) are not supported within the current framework.

The following example illustrates items $`4`$ and $`5`$ from the above observations. We use a small static scene composed of spheres with spikes, as the ones shown on Figure 3, right. Each sphere is composed of $`11,540`$ triangles. Sequentially, one sphere is drawn by the Open Inventor at a rate of $`195`$ frames/sec, i.e. $`2,251,569`$ triangles/sec. For $`2`$ spheres the rate is $`117`$ frames/sec or $`2,709,130`$ triangles/sec, etc. The speed of visualizing $`2`$ spheres on $`2`$ processors ($`2`$ tiles) is approximately $`6`$ frames/sec. The composited visualization ($`2`$ processors, $`1`$ tile) is approximately $`3`$ frames/sec. Runs with various problem sizes and parallel visualization configurations give similar results in favor of the sequential execution.

The enormous difference is due to the fact that in the first case the scene resides in the graphics card memory, while in the second the same scene is transmitted over the network for every frame. Table $`1`$ gives more results related to the network’s bandwidth bottleneck and the overall performance of the system. The results are for different applications using the $`100`$ Mbit and the $`1`$ Gbit network. The first one, mouse brain, comes from medical science and is a 3D surface of a mouse brain. The size of the data is $`16`$ MB. The $`100`$ Mbit network gets saturated and the frame rate of $`0.46`$ frames/second is expected, since the data traffic, although dependent on user interaction and locality of the scene polygons, is proportional to the data size. The $`0.46`$ frames/second translates to $`2.17`$ seconds per frame. The network transfer takes $`1.75`$ seconds. The next application, fluid flow, is an isosurface extracted from a fluid flow simulation data. The beetle head is also an isosurface. It was extracted from the X-ray computed microtomography data of an Alaskan spruce bark beetle. With performance mainly depending on the data size we observe that doubling the data size reduces the frame rate two times. Another application type that we tested is ray tracing volume visualization. We applied it to X-ray computed microtomography data of a rock sample of size $`1`$ GB. This is an example where a substantial part of the visualization is done in the CPU and only the result, in terms of OpenGL primitives, is sent through the network. On the $`100`$ Mbit network we get $`12.5`$ seconds per frame. $`10.1`$ seconds are spent in the ray tracing algorithm (run on $`4`$ cluster nodes with dual processor on each node) and $`2.4`$ seconds in transfer and rendering of the OpenGL commands issued in the ray tracing algorithm.

Switching to the $`1`$ Gbit network approximately doubles the performance results. Similar improvement was observed in  for various numerical analysis applications. The ray casting relies mostly on the CPU’s performance and switching to the $`1`$ Bbit network did not show any speed up.

We used SGI’s pmchart to monitor the network traffic for the different applications described above and in Table $`1`$. The network gets saturated for the applications discussed. The network traffic in and out of every cluster node looks like the one given on Figure 4.

The non-automatic mechanism that Chromium provides for acceleration (minimization) of the geometry flow is display lists. They are supported by sending the lists to each rendering server, and thus guaranteeing their presence on the server when they are called. Display lists in Open Inventor are created for the parts of the scene that have as root SoSeparator nodes with field renderCaching turned ON or AUTO (for more information see , pages $`224`$-$`227`$). Using display lists we get a speed of $`90`$ frames/sec for visualizing $`2`$ spheres by $`2`$ processors on composited display, and $`80`$ frames/sec on tiled display. Different parallel and display configurations illustrate the same order of improvement. These results are comparable to the sequential ones. The trade-offs are that the process is not automatic and that the display lists are broadcast to all rendering servers and thus replicating the scene in each node. Work on the automatic optimization of the geometry data flow can be found in .

## 4 Conclusions

We developed an application for interactive parallel visualization on tiled displays for commodity-based clusters. We used Chromium’s parallel rendering and extended popular visualization APIs, such as Open Inventor and VTK. The most original part of this work is the combined use of Open Inventor and Chromium. We gave implementation details and described our experience in the combined use of the Chromium’s rendering technology with Open Inventor and VTK. A general conclusion is that fast sequential visualization often relies on graphics acceleration hardware, data sampling methods, static scenes, and data size within the hardware limitations, while the visualization developed does not require expensive graphics hardware, provides high resolution, facilitates well the visualization of time-varying scenes, and can handle large data sets. We demonstrated a low cost (time, money, effort, etc.) of development. Our benchmarks quantitatively demonstrated the bottleneck that the network’s bandwidth presents in the sort-first rendering architecture. This bottleneck makes the sort-first architecture more appealing for time-varying data sets, commodity-based clusters with slow (or no) graphics acceleration, and visualization algorithms that rely on the CPU’s speed, such as ray casting.

Acknowledgments. We would like to thank Beverly Tomov from Cold Spring Harbor Laboratory, NY, for her attentive editing and remarks. We thank the anonymous referees whose constructive remarks improved the presentation. We also thank Keith Jones from the BNL’s Environmental Research and Technology Division for providing the data produced by the synchrotron-based computed microtomography at the Brookhaven National Synchrotron Light Source.


format: text

text:
# Embedding spherical spacelike slices in a Schwarzschild solution

## I Spherical spacelike slices in a Schwarzschild spacetime

It is clear that many spherical spacelike slices can be embedded in a given extended Schwarzschild spacetime. In the negative mass Schwarzschild (or in flat space-time) all the slices will have topology $`R^3`$. In the positive mass Schwarzschild one has a range of topologies. There are slices which start at one spacelike infinity, run through the middle and out the other end. These have topology $`S^2\times R^1`$. There are slices which start at one or other of the $`R=0`$ singularities and go out to one or other of the infinities. These have topology $`R^3`$. Finally, there are slices which start and end on one or other of the singularities. These have topology $`S^3`$. There are variants of the infinite ones, where instead of going to spacelike infinity, the slice remains within the horizon and runs more or less along one of the $`R`$ = constant lines, with $`R<2m`$. These slices are asymptotically cylindrical. Inside the horizons, the lines of constant $`R`$ are spacelike, while outside they are timelike. Therefore inside the horizons a spacelike spherical slice can have an oscillatory areal radius with as many maxima and minima as one wishes, while outside the horizon the area must monotonically increase. These features are, of course, reflected in the nature of the spacetime in which any given slice can be embedded.

Any spacelike slice embedded in a solution of the Einstein equations will have an intrinsic geometry, given by a three-metric $`g_{ab}`$, and an extrinsic curvature, given by a symmetric tensor $`K^{ab}`$. These are not independent, they must satisfy the constraints

$$^{(3)}K_{ab}K^{ab}+\left(g_{ab}K^{ab}\right)^2=0,$$
(1)
$$_a\left(K^{ab}g^{ab}g_{cd}K^{cd}\right)=0,$$
(2)
where $`^{(3)}`$ is the scalar curvature of $`g_{ab}`$. Eq.(1) is called the Hamiltonian constraint, Eq.(2) is called the momentum constraint.

Given a spherical slice in a spherical space-time then both the three-geometry and the extrinsic curvature are spherically symmetric. There are several ‘natural’ choices of coordinates that are used to write down the spherical three-metric. The one we favor is the ‘proper distance gauge’ where one writes the three-metric as

$$ds^2=dl^2+R^2(l)d\mathrm{\Omega }^2,$$
(3)
where $`d\mathrm{\Omega }^2=d\theta ^2+\mathrm{sin}^2\theta d\varphi ^2`$ is the round two-metric on a sphere, $`R(l)`$ is the areal (Schwarzschild) radius of the isometry two-spheres (expressed as a function of $`l`$) and $`l`$ is the proper distance between surfaces of constant $`R`$. Any spherical three-geometry can be expressed in this form. If the topology of the slice is $`S^2\times R`$, $`l(\mathrm{},\mathrm{})`$, if the topology is $`R^3`$, $`l[0,\mathrm{})`$, and if the topology is $`S^3`$, $`l[a,+b]`$ where $`a+b`$ is the proper distance between the points where $`R0`$.

A spherically-symmetric symmetric two-tensor has only two independent components. We can write the extrinsic curvature in a form consistent with spherical symmetry as gom

$$K^{ab}=n^an^bK_l+(g^{ab}n^an^b)K_R,$$
(4)
where $`K_l`$ and $`K_R`$ are two scalars and $`n^a`$ is the outward-pointing unit normal to the two-surfaces of constant $`R`$, $`n^a=(1,0,0)`$ in our preferred gauge.

The Hamiltonian and momentum constraints can now be written in the special case of spherical symmetry as gom

$$K_R\left[K_R+2K_l\right]\frac{1}{R^2}\left[2RR^{\prime \prime }+R^21\right]=0,$$
(5)
and

$$K_R^{}+\frac{R^{}}{R}\left[K_RK_l\right]=0,$$
(6)
where  represents the derivative with respect to $`l`$. There is a first integral of the constraints gom

$$m=\frac{R^3K_R^2}{2}+\frac{R}{2}\left[1R^2\right],$$
(7)
where $`m`$ is the Schwarzschild mass. This is the Misner-Sharp, Hawking et al mass formula. An immediate consequence of Eq.(7) is that if we define

$$M=\mathrm{max}\frac{R}{2}\left[1R^2\right],$$
(8)
where the maximum is taken over the whole three-geometry then

$$mM.$$
(9)
This is our key equation, it is an algebraic relation between a global spacetime quantity, $`m`$, and a quantity, $`M`$, which depends only on the three-geometry.

In this article we will show that essentially the converse of condition (9) holds. More precisely, given a three-geometry in the proper distance gauge, one can always compute the quantity $`M`$ as defined by Eq.(8). Let us assume that it is finite. This only requires that $`Rl`$ as $`l\pm \mathrm{}`$. We will show that this spacelike spherical geometry can be embedded in any Schwarzschild solution whose Schwarzschild mass satisfies

$$m>M.$$
(10)
The quantity $`M`$ is really a coordinate independent object. $`R`$ is the area of the two-dimensional isometry spheres and $`2R^{}/R`$ is the mean extrinsic curvature of the isometry two-spheres as embedded surfaces in the three-geometry.

## II Embedding slices in a Schwarzschild solution

Let us assume that we are given a spherical three metric, say in the form of Eq.(3). It could be a complete manifold or just a patch. We wish to consider this as the metric of a hypersurface in a spherical vacuum solution of the Einstein equations. Therefore we need to find a spherical extrinsic curvature $`(K^{ab})`$ expressed, say, in the form of Eq.(4), such that the combination $`(g_{ab},K^{ab})`$ satisfy both the Hamiltonian and momentum constraints, Eqs. (1) and (2) (or, rather, Eqs. (5) and (6)). A priori, at least, this seems to be a reasonable task. We have two free functions, $`K_R`$ and $`K_l`$, and we need to pick them so as to satisfy two scalar equations. However, we do not have a completely free choice. Given the metric, we can compute $`M`$ via Eq.(8) and if we succeed in embedding this in a Schwarzschild solution of mass $`m`$, we must have that $`mM`$. In addition to being a necessary condition, it is almost a sufficient condition. More precisely, we can prove:

Proposition: Given a spherical Riemannian three-geometry with finite $`M`$ as defined by Eq.(8) and any (extended) vacuum Schwarzschild solution whose mass $`m`$ satisfies $`m>M`$, then a spacelike slice can be found in this spacetime which is isometric to the given three-geometry.

Proof: We start off with the mass expression, Eq.(7), and write it as

$$K_R=\sqrt{\frac{2m}{R^3}\frac{1}{R^2}\left[1R^2\right]}.$$
(11)
The quantity under the square root is positive definite since $`m>M`$ and so the equation makes sense. We can choose either the positive or the negative root. With either choice we get a well-defined function $`K_R(l)`$ which does not change sign. It may well be that if $`R0`$ then $`|K_R|\mathrm{}`$. All this does is reflect the fact that the slice goes into the Schwarzschild singularity.

Given that $`K_R`$ is bounded away from zero we now can solve the Hamiltonian constraint, Eq.(5), for $`K_l`$. More precisely, we rewrite Eq.(5) as

$$2K_RK_l+\frac{2m}{R^3}\frac{2R^{\prime \prime }}{R}=0$$
(12)
and manipulate this to get

$$K_l=\frac{1}{K_R}\left[\frac{R^{\prime \prime }}{R}\frac{m}{R^3}\right].$$
(13)
This is clearly well-defined and finite except possibly, again, as $`R0`$. It is very straightforward to show that $`K_R`$ from Eq.(11) and $`K_l`$ from Eq.(13) satisfy the momentum constraint, Eq.(6).

We have only shown that we can construct spherically symmetric initial data that satisfies the vacuum Einstein constraints. However, the Einstein evolution equations allow us to propagate this data so as to construct at least a patch of spacetime (which will be spherically symmetric) and which satisfies the vacuum Einstein equations. In turn, Birkhoff’s theorem guarantees that this must be part of the Schwarzschild solution with given mass $`m`$.

## III A concrete example

We wish to embed the ‘static’ (moment-of-time-symmetry) slice from one Schwarzschild solution (of mass $`m_1`$) in another Schwarzschild solution (of mass $`m_2`$). We assume $`m_2>m_1`$. The three-metric we are given, written in Schwarzschild coordinates is

$$ds^2=\frac{dR^2}{1\frac{2m_1}{R}}+R^2d\mathrm{\Omega }^2.$$
(14)
To convert to proper distance gauge we would need to integrate

$$l(r)=_{2m_1}^r\frac{dR}{\sqrt{1\frac{2m_1}{R}}}$$
(15)
but there is no real need to do so, all we really use is

$$R^{}=\frac{dR}{dl}=\sqrt{1\frac{2m_1}{R}}\frac{R}{2}\left[1R^2\right]=m_1.$$
(16)
When this is substituted into Eq.(11) we get

$$K_R=\sqrt{\frac{2(m_2m_1)}{R^3}}.$$
(17)
We choose the negative root to get the slice in the upper half plane.

We can use the Hamiltonian constraint, Eq.(5), remembering that the scalar curvature of the metric given by Eq.(14) vanishes, to get

$$K_l=\frac{K_R}{2}=\sqrt{\frac{m_2m_1}{2R^3}}.$$
(18)
It is easy to confirm that this choice of extrinsic curvature satisfies the momentum constraint, Eq(6).

We know that maximal slices which run from one end to the other of a Schwarzschild solution cannot approach the singularity closer than $`R=3m/2`$ bom . If we relax the maximal condition, and replace by the requirement that the three-scalar-curvature be nonnegative we get no such restriction, we can approach the singularity as closely as we wish.

If we let $`m_10`$ we get the well-known flat slice of the Schwarzschild solution. Now $`R0`$ so the flat slice runs into the singularity and the extrinsic curvature is given by

$$K_R=2K_l=\sqrt{\frac{2m_2}{R^3}}.$$
(19)
We can also embed the ‘static’ slices of the negative mass Schwarzschild solution in flat spacetime. These slices have intrinsic metric

$$ds^2=\frac{dR^2}{1+\frac{2m_1}{R}}+R^2d\mathrm{\Omega }^2,$$
(20)
and extrinsic curvature

$$K_R=2K_l=\sqrt{\frac{2m_1}{R^3}}.$$
(21)
This slice is regular everywhere except at the origin.

## IV The significance of $`m`$ $`>M`$

We seek the maximum of $`\frac{R}{2}[1R^2]`$. If we wish to embed this slice in a negative mass Schwarzschild solution (or in flat spacetime) we need $`\frac{R}{2}[1R^2]<0`$ over the entire slice. This implies $`R^{}1`$. The areal radius starts off at $`R=0`$ and monotonically increases, there cannot be a throat. Not only that but it must increase rapidly. It must look like a trumpet. This is completely in agreement with the requirement that the topology of slices embedded in the negative mass Schwarzschild must be $`R^3`$.

In general, however, the maximum of $`\frac{R}{2}[1R^2]`$ will be positive and therefore this slice can only be embedded in a positive mass Schwarzschild solution. For example, if the area is not monotonic then the function $`R(l)`$ has either a local maximum or minimum (or both). Thus the slice has at least one point where $`R^{}=0`$, and if the areal radius equals $`R_0`$ at that point then $`\frac{R}{2}[1R^2]=R_0/2`$ there. Thus we have that $`MR_0/2`$ and a necessary condition that this slice be embeddable in a Schwarzschild solution of mass $`m`$ is that $`2mR_0`$. This means that if we have an extremum of the area we can only embed this slice in a positive mass Schwarzschild solution and also that all maxima and minima of the areal radius must occur inside the horizon.

If we have a local maximum of $`R(l)`$, i.e., a point where $`R^{}=0,R^{\prime \prime }0`$, this point is also a local maximum of $`\frac{R}{2}[1R^2]`$. This is easy to show, all one needs to note is that both $`R/2`$ and $`[1R^2]`$ have a maximum at that point. No equivalent statement can be made about a point where $`R(l)`$ is a local minimum. One could have a maximum or minimum of $`\frac{R}{2}[1R^2]`$ at that point or nothing at all. It is also worth noting that the maximum of $`\frac{R}{2}[1R^2]`$ need not occur at an extremum of $`R(l)`$. It may occur at a ‘large’ value of $`R`$ where $`R^{}`$ may be small without ever going to zero. One can then choose an $`m>M`$ such that this point is outside the horizon.

## V The marginal case: what happens when $`m=M`$

If we have a spherical slice in a Schwarzschild solution we know that $`mM=\mathrm{max}\frac{R}{2}[1R^2]`$. Conversely, we have shown that if we have a spherically symmetric Riemannian three metric we can embed it in a Schwarzschild solution with mass $`m`$ if $`m>M`$. In this section we would like to discuss the issue of when and if one can embed a given spherical slice in a Schwarzschild solution satisfying $`m=M`$. Very useful tools for analysing spherical spacetimes are the optical scalars, $`\mathrm{\Theta }_\pm `$ (see e.g. gom ; mom ). These are defined by

$$\mathrm{\Theta }_\pm =\frac{2}{R}\left(R^{}\pm RK_R\right).$$
(22)
These are the divergences of the future pointing and past pointing outward radial light rays from the isometry spheres. Both are defined so as to be positive on flat space and near infinity. Note that the definition of extrinsic curvature used agrees with Wald w  and not with Misner, Thorne, and Wheeler mtw .

The mass expression, Eq.(7), can be rewritten as

$$\frac{2m}{R}=1\frac{\mathrm{\Theta }_+\mathrm{\Theta }_{}R^2}{4}.$$
(23)
If $`m0`$ then $`\mathrm{\Theta }_+`$ and $`\mathrm{\Theta }_{}`$ must always remain positive. If $`m>0`$ and if $`2m/R>1`$, i.e., if we are inside the horizon, then one or other of the optical scalars must be negative. In the upper half of the extended Schwarzschild solution it is $`\mathrm{\Theta }_+`$ which is negative and in the lower half it is $`\mathrm{\Theta }_{}`$ which is negative. Eq.(11) is double valued because we can choose either sign when we take the square root. This means that we have two slices with the same intrinsic geometry. Given a solution to Eqs.(5, 6) and multiply $`K_R`$ and $`K_l`$ by a minus sign we again have a solution with the same three-geometry. These two solutions are just reflections of each other about the $`t=0`$ plane. The slice with $`K_R>0`$ crosses the horizon below $`t=0`$ and either crashes into the past singularity or continues through while the slice with $`K_R<0`$ tends to live in the upper half plane.

Say we are given a spherical three-geometry with finite $`M`$ and try to embed this in a Schwarzschild spacetime with $`m=M`$. This means that we no longer, from Eq.(11), get that $`K_R>0`$. Rather we get that $`K_R=0`$ at the point(s) where $`M`$ achieves its maximum. At such a point, since $`\left(\frac{R}{2}[1R^2]\right)`$ is at its maximum, we know that

$$0=\frac{}{l}\left(\frac{R}{2}[1R^2]\right)=\frac{R^2R^{}}{4}\times \frac{2}{R^2}\left[1R^22RR^{\prime \prime }\right].$$
(24)
Therefore at the point(s) where $`K_R=0`$ we are guaranteed that either $`R^{}=0`$ or $`\frac{2}{R^2}\left[1R^22RR^{\prime \prime }\right]=0`$. From Eq.(5) it is clear that we need the second expression (which is nothing else but the three scalar curvature, $`^{(3)}`$, of the three-geometry) to vanish whenever $`K_R=0`$. Hence we need to distinguish between the case where the maximum of $`\frac{R}{2}[1R^2]`$ occurs at a point where $`R^{}=0`$ (which is bad, except if simultaneously $`^{(3)}=0`$ there, which cannot be guaranteed) and the case where the maximum occurs at a point where $`R^{}0`$ (which is good because $`^{(3)}=0`$ there). This question is trivial in the case where $`M0`$ because there cannot be any point(s) with $`R^{}=0`$ in the geometry and so we automatically have that $`^{(3)}=0`$ at the point where $`M`$ achieves its maximum.

It is also clear that we are unable to use Eq.(13) directly to evaluate $`K_l`$ at the point(s) where $`K_R=0`$. However, we can return to the momentum constraint Eq.(6) and rewrite it as

$$K_l=\frac{RK_R^{}}{R^{}}+K_R.$$
(25)
This shows us that $`K_l`$ is finite and well defined at the point(s) where $`K_R=0`$ so long as $`R^{}0`$ at those points. This is consistent with Eq.(13) because the term in square brackets in this equation is essentially $`^{(3)}`$ so Eq.(13) becomes $`K_l=0/0`$ and use of de l’Hospital’s rule gives us a finite value. Therefore the only case we need to worry about is when the maximum of $`M`$ coincides with $`R^{}=0`$. We cannot, in general, hope to solve the constraints with $`m=M`$ in this situation. Nevertheless, there are special cases.

It is clear that if there exists a point where $`K_R=0`$ and $`R^{}=0`$ then both optical scalars, $`\mathrm{\Theta }_\pm `$, vanish simultaneously. In the Schwarzschild solution the only place where this can happen is at the bifurcation ‘point’, where both horizons cross. Note again that $`R^{}=0`$ is only possible with $`m>0`$. If the given three-geometry has the property that that the maximum of $`\frac{R}{2}[1R^2]`$ occurs at a point where $`R^{}=0`$ and we try $`m=M`$, then clearly $`K_R=0`$ and $`R^{}=0`$ at this point; if this slice is to be embedded in this Schwarzschild solution it must pass through the bifurcation point. This places both local and global restrictions on the geometry of the slice. This point where $`R^{}=0`$ must be a local minimum of the area, i.e., $`R^{\prime \prime }0`$. More than that, we require that the scalar curvature vanishes at that point, hence $`2R^{\prime \prime }=1/R`$. Further, it must be the global minimum. The area cannot oscillate. If we had a local maximum of $`R`$ then the value of $`\frac{R}{2}[1R^2]`$ at this point would be larger than the value at the ‘throat’, where it is supposed to be the maximum. Therefore we can only have one minimum.

The behaviour of $`K_R`$ in the neighbourhood of the point where $`K_R=0`$ is worth noting. Let us define $`f(l)=\frac{R}{2}[1R^2]`$ and assume $`f=f_0=m`$ at $`l=l_0`$, $`f_0^{}=0`$, and $`f_0^{\prime \prime }<0`$. We have $`K_R^2=\frac{2m2f}{R^3}`$. Hence we have $`(K_R^2)_0=0`$, $`(K_R^2)_0^{}=0`$, and $`(K_R^2)_0^{\prime \prime }=(\frac{2f^{\prime \prime }}{R^3})_0>0`$. Now make a Taylor expansion of $`K_R^2`$ around $`l=l_0`$ in terms of $`x=ll_0`$. In general we will get $`K_R^2(x)=A+Bx+Cx^2+\mathrm{}`$ with $`A=B=0`$ and $`C>0`$. Thus we get $`K_R(x)=\sqrt{C}x+\mathrm{}`$. Therefore $`K_R`$ will pass through zero at $`l=l_0`$ with nonzero slope and so must change sign.

There is a further special case. The requirement that $`f(l)`$ be a maximum at $`l=l_0`$ only requires $`f_0^{\prime \prime }0`$. If $`f_0^{\prime \prime }=0`$ we would have $`(K_R^2)_0=0`$, $`(K_R^2)_0^{}=0`$, and $`(K_R^2)_0^{\prime \prime }=(\frac{2f^{\prime \prime }}{R^3})_0=0`$. Now the requirement that we have a maximum of $`f`$ forces $`f_0^{\prime \prime \prime }=0`$ and the first nontrivial derivative must be at fourth order. In this case the Taylor expansion of $`K_R^2`$ can only start at $`x^4`$, i.e. $`K_R^2=Dx^4+\mathrm{}`$ and $`K_R=\sqrt{D}x^2+\mathrm{}`$. In this case $`K_R=0`$ and $`K_R^{}=0`$ at $`l=l_0`$. Following from this not only does $`K_R`$ vanish at $`l=l_0`$ but also, from Eq.(25), we have that $`K_l=0`$ at $`l_0`$.

To summarize: If the maximum of $`\frac{R}{2}[1R^2]`$ occurs at a point where $`R^{}0`$ we can choose $`m=M`$. This condition of $`R^{}0`$ is trivially satisfied if $`M<0`$. If the maximum occurs at a point where $`R^{}=0`$ we can only choose $`m=M`$ if this point becomes the bifurcation point. This places further requirements (both local and global) on the three-geometry.

## VI Spherical Superspace and the thick sandwich theorem

Wheeler identified the configuration space of canonical general relativity as being the space of all spacelike three-geometries, which he called ‘superspace’. A trajectory of a solution of the Einstein equations in this configuration space corresponds to the sequence of Riemannian three-geometries generated by a foliation. This is not a unique curve, rather we have ‘a spray of geodesics’ in the language of DeWitt. Every slicing of a spacetime generates a different sequence of three-geometries and by changing the slicing one changes the curve through superspace. It is interesting to ask what fraction of superspace do all the curves corresponding to one given spacetime pass through. This set of three-geometries which is the union of all the solution curves is exactly the same as the set of all spacelike three-geometries that can be embedded in a given spacetime.

Let us define spherical superspace $`(SS)`$ as the space of all spherically symmetric spacelike three-geometries. Now consider a Schwarzschild solution $`(S(m))`$, which is defined by its mass $`m`$. There will be many spacelike spherical three-geometries that can be embedded in $`S(m)`$. We define $`B_S(m)`$ as the collection of all spherically symmetric spacelike three-geometries that can be embedded in a given $`S(m)`$. Obviously $`B_S(m)SS`$. The analysis so far tells us a great deal about the relationship between $`B_S(m)`$ and $`SS`$.

It is clear that any single $`B_S(m)`$ with given mass $`m`$ covers a large fraction of spherical superspace. We have a measure on spherical spacelike three-geometries given by $`M=\mathrm{max}\frac{R}{2}\left[1R^2\right]`$ and we have shown that all three-geometries which satisfy $`M<m`$ belong to $`B_S(m)`$. We also know that any three-geometry which satisfies $`M>m`$ cannot belong to $`B_S(m)`$. Thus $`B_S(m)`$ is defined by a single algebraic condition. The only uncertainty is about those three-geometries on the boundary of $`B_S(m)`$ i.e., those metrics which satisfy $`M=m`$. If $`m0`$ we know that $`B_S(m)`$ is closed; every geometry which satisfies $`M=m`$ can be embedded in the appropriate Schwarzschild solution. This is because $`M0R^{}0`$ so, as discussed in the previous section, we have no difficulty in solving for $`K_l`$ at the point where $`K_R=0`$. If $`m>0`$ we have that $`B_S(m)`$ is neither open nor closed. We have that most of the geometries on the boundary can be embedded. However, we have a relatively small class of geometries, those for which the maximum of $`M`$ occurs at a local maximum of the area, i.e., a point where $`R^{}=0`$, which cannot be embedded in the Schwarzschild solution with mass $`m=M`$.

Further, we have a nested structure on spherical superspace. Given two Schwarzschild solutions with masses $`m_1`$ and $`m_2`$ with $`m_1<m_2`$ then every spherical slice that can be embedded in the first solution can also be embedded in the second one. This means that $`B_S(m_1)B_S(m_2)`$. The foliation freedom in a given solution to the Einstein equations is represented by the fact that one can choose an arbitrary lapse (the shift freedom can be ignored because we are considering the geometries rather than three-metrics). Starting with spherical data and maintaining the spherical symmetry means that we must have a spherical lapse. Therefore the foliation choice is represented by one single spherical function. However, as we have seen in Eq.(3), the freedom in spherical three-geometries is also represented by one single spherical function. Thus it should not come as too much of a surprise that any one Schwarzschild solution should cover so much of spherical superspace.

Since this kind of counting argument works so well in the spherical case, it is interesting to try it in more general situations. Let us define axial superspace as the set of all axially symmetric spacelike geometries. The foliation freedom of axially symmetric slicings of an axially symmetric spacetime is represented by a single axially symmetric function. However, to give a general axially symmetric three-geometry (say written in the Brill conformal gauge b ) we have to give two axially symmetric functions. Therefore we expect that a given axially symmetric spacetime will visit only a subset of axially symmetric superspace, corresponding more-or-less to the square root of the whole. In the general, unsymmetric, case the foliation freedom is represented by a single arbitrary function while the general three-geometry is represented by three arbitrary functions.

Spherical geometries are also useful in that they offer a simple counterexample to a longstanding idea in canonical gravity, the thick-sandwich approach. Wheeler introduced the concept of the ‘thick sandwich theorem’. The plan was to choose an initial and final point in configuration space (in this case a pair of three-geometries) and that the equations of motion (the dynamical part of Einstein equations) would give a natural path joining these two points. In other words to try and find a unique four manifold, satisfying the Einstein equations, which fills in between the two given three-geometries. The analysis in this article shows how ill posed this ‘thick sandwich’ concept is.

Choose any two spherical three-geometries. Evaluate $`M`$ for each geometry. Let them be $`M_1`$ and $`M_2`$ respectively. Pick any Schwarzschild solution with mass $`m`$ which satisfies $`m>\mathrm{max}(M_1,M_2)`$. Both of the given three-geometries can simultaneously be embedded in the given Schwarzschild solution. This generates a ‘filling’ (which satisfies the Einstein equations) between the given slices which is highly nonunique. Further, it is impossible to lift this ambiguity by demanding that the given three-geometries are, with respect to any measure, close to one another. If we find any Schwarzschild solution that interpolates between them, all Schwarzschild solutions with larger mass will also do the job.

###### Acknowledgements.

K.R. acknowledges the support of the TEMPUS project JEP-12249-97. We wish to thank Edward Malec for much help.


text:
# Nuclear Magnetic Resonance in a Ferromagnet-Semiconductor Heterostructure

## Abstract

We report the observation of nuclear magnetic resonance (NMR) in a ferromagnet-semiconductor heterostructure in the presence of a spin-polarized current. Spin-polarized electrons injected from a metallic ferromagnet generate a large nuclear spin population in a GaAs quantum well by dynamic polarization. The characteristic time for the polarization process is approximately 20 sec, and the nuclear polarization can persist for several minutes after the current is turned off. Resonant depolarization is observed in the presence of an AC magnetic field or when the injection current is modulated at the NMR frequency.

A detailed understanding of the hyperfine interaction is a prerequisite for the coherent control of electron and nuclear spin systems in semiconductor heterostructures. Electrical spin injectionMolenkamp\_Ohno\_Ploog\_Jonker\_IMEC  provides a simple means of introducing spin-polarized carriers into a semiconductor, but it has proven more difficult to achieve the degree of tunability over the combined electron and nuclear spin systems demonstrated in optical pumping experiments.Kikkawa:Awschalom:allOptNMR ; Salis:OptNMR-PRL ; Epstein:Awschalom:FPP\_in\_n-GaAs ; Stephens:Awschalom:Imaging  For example, Poggio et al.Poggio:Awschalom:gated\_NMR  have recently shown that it is possible to electrically control the spatial profile of optically-pumped nuclear polarization within a single quantum well. One of the advantages of optical pumping is the efficiency of dynamic nuclear polarization (DNP), in which a non-equilibrium nuclear polarization is generated by the spin relaxation of electrons.Overhauser ; Lampel ; Paget:GaAsNuclearSpinCoupling ; NucMag:Abragam  In principle, the DNP process can also be driven by a spin-polarized current.Johnson:ElectricalDNP-proposal  This was inferred in recent experiments on Fe/AlxGa1-xAs heterosructures,Strand:Crowell:electrical-DNP  although an explicit dynamical signature of DNP was not demonstrated.

In this letter we report the direct observation of nuclear magnetic resonance (NMR) in a GaAs quantum well (QW) in the presence of a spin-polarized current. We demonstrate that DNP can be partially supressed by resonant depolarization of the nuclear spins. The depolarization can be achieved either by the traditional approach of applying a resonant AC magnetic field or by modulating the spin-polarized current at the NMR frequency. The time dependence of the DNP process is measured, indicating a characteristic build-up time of $`20`$ sec and persistence of nuclear polarization up to eight minutes after the spin-polarized current has been turned off.

The heterostructure is grown by molecular beam epitaxy on a p+ GaAs (100) substrate and consists of p-Al0.1Ga0.9As/QW/n-Al0.1Ga0.9As/Fe/Al. The QW comprises 100 Å GaAs between intrinsic Al0.1Ga0.9As barriers, and the Fe film is 50 Å thick. The metal layers are grown in situ at $`0^{}`$ C, and a $`3\times 10^{13}`$ atoms/cm2 $`\delta `$-doped layer of Si is grown 25 Å from the Fe/AlxGa1-xAs interface to create a thin Schottky tunnel barrier. The devices are fabricated into bars with an active region 80 $`\mu `$m wide (see schematic drawing in Fig. 1), and the experiments are carried out in a magneto-optical cryostat at temperatures ranging from 2 K to 40 K.

The measurements are performed in the Voigt geometry with the applied field $`𝐇_0`$ in the plane of the QW. Electroluminescence (EL) is collected along the growth direction, indicated by $`𝐤`$ in Fig. 1, which coincides with the direction of the AC magnetic field $`𝐇_1`$ generated by a four-turn coil placed on top of the sample. The EL is due to the recombination of electrons tunneling from the Fe layer with heavy holes flowing from the substrate. The circular polarization of the EL, referred to as the ELP signal, is equal to the component of electron spin polarization along $`𝐤`$ at the time of recombination.

Typical ELP data as a function of applied field are shown in Fig. 1. All of the measurements under discussion here were conducted with $`𝐇_0`$ nearly parallel to the $`[01\overline{1}]`$ direction, which corresponds to a magnetic hard axis. As $`H_0`$ decreases below 500 Oe, the magnetization $`𝐌`$ rotates away from $`[01\overline{1}]`$ towards the $`[011]`$ (easy) direction. The non-zero angle between the injected spin $`𝐒_0𝐌`$ and $`𝐇_0`$ leads to precession of the spin in the QW, resulting in a steady-state component of $`𝐒`$ along $`𝐤`$ and hence the observed ELP signal. The ELP hysteresis loops shown in Fig. 1 can be modeled extremely well provided that the total field $`𝐁=𝐇_0+𝐁_N`$, where $`𝐁_N`$ is assumed to be the hyperfine field due to nuclei polarized by DNP. As discussed in Ref. Strand:Crowell:electrical-DNP, , the magnitude of $`𝐁_N`$ depends on the injection current density and the angle $`\theta `$ between $`𝐒_0`$ and $`𝐇_0`$ and reaches several kilogauss for current densities $`10`$ A/cm2 and $`\theta <\pi /2`$.

The large nuclear magnetic field is due to the long spin-relaxation time for nuclei and the large difference between the electron and nuclear magnetic moments.NucMag:Abragam  Since the nuclear spin-lattice relaxation time can be many minutes,OO ; Kikkawa:Awschalom:allOptNMR  we expect the field $`𝐁_N`$ to build up on laboratory time scales when the current is turned on and then decay slowly when the current is turned off. Figure 2 shows the ELP signal as the device current is switched on and off for variable amounts of time. In this case, the current is off for ten minutes prior to the start of data collection at $`t=0`$. The LED is then turned on and the ELP signal is measured for three minutes. The increase in the ELP signal is due to the increase in total internal magnetic field and the consequent increase in the angle of precession of the injected spin polarized electrons. To probe the lifetime of the nuclear polarization in the absence of continuous pumping, the current is turned off for a time $`t_{\mathrm{off}}`$ and then on again for three minutes. Figure 2 shows a series of these measurements with $`t_{\mathrm{off}}`$ varying from 1 min to 9 min. The first data point taken after each “off” period indicates the degree of nuclear polarization that has been retained. By repeating this procedure with increasing $`t_{\mathrm{off}}`$, we observe a decay time of approximately 5 min. Although both the build-up and decay times for the ELP signal are long, it is difficult to make a direct quantitative mapping between the magnitude of the ELP signal and the actual nuclear polarization. In particular, the relationship between the ELP signal and $`B_N`$ is very non-linear. The ELP signal increases rapidly at small nuclear fields and then saturates,Strand:Crowell:electrical-DNP  and this is responsible in part for the large difference between the build-up and decay times for the ELP signal.

Although the nuclear polarization can relax slowly through coupling to the lattice, the most effective way of destroying DNP is to apply an oscillating magnetic field $`𝐇_1`$ at the NMR frequency along a direction perpendicular to $`𝐇_0`$. The transverse field couples directly to the nuclear spin system, leading to depolarization of the nuclei that are in resonance. $`H_1`$ is generated by a 1 cm diameter, four turn coil driven by the sinusoidal output of a function generator. By sweeping the frequency of $`H_1`$ at fixed $`H_0`$, depolarization signatures for each of the three isotopes in the GaAs QW \[75As ($`\gamma =0.731`$ kHz/Oe), 69Ga ($`\gamma =1.025`$ kHz/Oe), and 71Ga ($`\gamma =1.302`$ kHz/Oe)\]CRCgamma  are observed, as seen in Fig. 3(a) for $`H_0=330`$ Oe. In each case, the linewidth of the depolarization signature is 10 - 15 kHz, with no systematic dependence on magnetic field. The linewidth at these fields ($`<500`$ Oe) is therefore not limited by inhomogeneous dephasing.

Frequency sweeps were conducted at several different magnetic fields. The peaks of the depolarization features for all three isotopes are plotted versus $`H_0`$ in Fig. 3(b). Linear fits of $`\nu =\gamma H_0`$ were obtained for each isotope, yielding gyromagnetic ratios $`\gamma _{{}_{}{}^{75}\mathrm{As}}=0.71\pm 0.1`$ kHz/Oe, $`\gamma _{{}_{}{}^{69}\mathrm{Ga}}=1.00\pm 0.04`$ kHz/Oe, and $`\gamma _{{}_{}{}^{71}\mathrm{Ga}}=1.27\pm 0.04`$ kHz/Oe, in agreement with the accepted values. The deviations from the fits are smaller than the limit set by trapped flux in the magnet ($`10`$ Oe).

The sign of the depolarization features warrants comment, since the dips in Fig. 3(a) correspond to an increase in the magnitude of the ELP signal, which is negative in this field range. We have confirmed that the features in Fig. 3(a) correspond to a decrease in the nuclear polarization by comparing with the observed dependence of the ELP signal on current density.Strand:Crowell:electrical-DNP  The reason for this non-intuitive result is that for sufficiently large nuclear fields, the average angle through which the electrons precess during their lifetime becomes larger than $`\pi /2`$. An increasing nuclear field therefore results in a decrease in the magnitude of the ELP signal because the electrons are actually precessing further away from the measurement axis. Although a complete quenching of $`B_N`$ would result in nearly total suppression of the ELP signal, only a single isotope is depolarized at each resonance. We have confirmed that for low $`B_N`$, corresponding to applied fields less than 60 Oe, resonant depolarization results in a decrease in the ELP magnitude as expected.

Since the nuclear polarization in these devices is generated by the injected electrons, an alternative approach to depolarization is resonant modulation of the device current. In this case, the coil is not used ($`H_1=0`$), and the DC bias voltage is set above the threshold for light emission. An AC voltage is then added so that the average bias voltage corresponds approximately to the value used in the measurements discussed above. In this case, the bias is being used only to modulate the injection current. Given the square well design, we cannot control the relative position of the nuclear and electronic spins as achieved recently in parabolic quantum wells.Poggio:Awschalom:gated\_NMR  Fig. 4 shows the ELP signal as a function of bias voltage frequency for $`H_0=330`$ Oe, $`V_{DC}=2.1`$ V, and $`V_{AC}=0.5`$ VPP. The principal resonances are clearly visible, as can be seen by comparing the data taken under alternating device current (lower curve) with the data taken under AC magnetic field (upper curve). In the case of alternating device current, however, there are many other peaks in the data, and frequency sweeps at different fields indicate that these also correspond to nuclear resonances. Some of the additional peaks can be mapped to sub-harmonics of the primary resonances, such as $`\nu =\gamma H_0/3`$. Other peaks occur at either harmonics of the primary resonances, such as $`\nu =2\gamma H_0`$, or their sub-harmonics, such as $`\nu =2\gamma H_0/3`$. Similar harmonic and subharmonic resonances have also been observed in quantum well systems with quadrupolar interactions in which a strong electric field gradient is modulated with a gate voltage.Poggio:Awschalom:gated\_NMR

Sub-harmonics of the primary resonances are never observed with an AC magnetic field, and we attribute them to the non-sinusoidal time dependence of the injection current as the bias voltage is modulated. This is a consequence of the diode structure of the device. The origins of harmonics of the fundamental resonances are less clear. These peaks are also observed in AC magnetic field experiments at high amplitudes ($`H_1>0.4`$ Oe). The observation of the $`\mathrm{\Delta }m=\pm 2`$ peaks in the AC magnetic field experiments at high amplitudes indicates that in this experiment they originate from some intrinsic non-linearity. Dipole-dipole coupling can lead to both $`\mathrm{\Delta }m=\pm 2`$ and $`\mathrm{\Delta }m=\pm 3`$ transitions,Kalevich:Fleisher:2spin\_resonance\_GaAs  and may play a more significant role because of the large nuclear polarization in these experiments.

The spectroscopic measurements presented here unambiguously demonstrate dynamic nuclear polarization by electrical spin injection. The ability to modulate the nuclear polarization with a spin-polarized current provides a unique means of controlling the coupled electron-nuclear spin system. This capability, in combination with the large nuclear polarization achieved at low fields, provides a new approach to NMR spectroscopy of heterostructures.

We thank S. A. Crooker and M. S. Conradi for helpful discussions. This work was supported by ONR, the DARPA/ONR SPINS program, and the University of Minnesota MRSEC (NSF DMR-0212032). A.F.I. and B.D.S. thank 3M for support.


text:
# Fingerprinting Hysteresis

## I Introduction

The conventional methods He et al. (1992); Proksh and Moskowitz (1994); Hedja et al. (1994) to characterize magnetic interactions in hysteretic systems, such as the $`\delta M`$ method Che and Bertram (1992); El-Hilo et al. (1992), utilize isothermal remanent magnetization (IRM) and dc demagnetization remanence (DCD) curves based on the Wohlfarth relation Wohlfarth (1958). Recently, FORC diagrams Pike et al. (1999); Katzgraber et al. (2002) have been introduced to study hysteretic systems. Their extreme sensitivity has helped to “fingerprint” several experimental systems as well as theoretical models ranging from geological samples and recording media to paradigmatic models of random magnets and spin glasses Katzgraber et al. (2002).

In this work we perform numerical simulations of random magnets (and spin glasses) in order to test the predictive power of FORC diagrams by comparing to a histogram of up- and down-switching fields of the underlying switching units along the major hysteresis loop.

The aforementioned re-parametrization of the major hysteresis loop (switching-field histogram) displays the information carried by the major loop in a more comprehensive way and provides a good comparison to the FORC diagram. We find, that the major-loop behavior predicts the minor-loop behavior captured by FORC diagrams well. We present a comparison of both distributions and discuss some differences between them. We argue that switching-field histograms are useful to study hysteretic systems in more detail than with conventional methods due to their simplicity and ease to compute.

## II Model & Algorithm

The Hamiltonian of the random-field Ising model (RFIM) is given by Ji and Robbins (1992)

$$=\underset{i,j}{}J_{ij}S_iS_j\underset{i}{}h_iS_iH\underset{i}{}S_i.$$
(1)
Here $`S_i=\pm 1`$ are Ising spins on a square lattice of size $`N=L^3`$ in three dimensions with periodic boundary conditions. The interactions between the spins are uniform ($`J_{ij}=1`$) and nearest-neighbor, and $`H`$ represents the externally applied field. Disorder is introduced into the model by coupling the spins to site-dependent random fields $`h_i`$ drawn from a Gaussian distribution with zero mean and standard deviation $`\sigma _\mathrm{R}`$.

We simulate the zero-temperature dynamics of the RFIM by changing the external field $`H`$ in small steps starting from positive saturation. After each field step we compute the local field $`f_i`$ of each spin:

$$f_i=\underset{j}{}J_{ij}S_jHh_i.$$
(2)
A spin is unstable if it points opposite to its local field, i.e., $`f_iS_i<0`$. We then flip a randomly chosen unstable spin and update the local fields at neighboring sites. This procedure is repeated until all spins are stable.

For the rest of this work we set $`L=50`$ ($`N=125000`$ spins) and $`\sigma _\mathrm{R}=5.0`$, unless otherwise specified. The different figures are calculated by averaging over 5000 disorder realizations in order to reduce finite-size effects.

## III FORC Diagrams

In oder to calculate an FORC diagram, a family of First Order Reversal Curves (FORCs) with different reversal fields $`H_\mathrm{R}`$ is measured, with $`M(H,H_\mathrm{R})`$ denoting the resulting magnetization as a function of the applied and reversal fields. Computing the mixed second order derivative Della Torre (1999); Pike et al. (1999)

$$\rho (H,H_\mathrm{R})=\frac{1}{2}[^2M/HH_\mathrm{R}]$$
(3)
and changing variables to $`H_\mathrm{c}=(HH_\mathrm{R})/2`$ and $`H_\mathrm{b}=(H+H_\mathrm{R})/2`$, the local coercivity and bias respectively, yields the “FORC distribution” $`\rho (H_\mathrm{b},H_\mathrm{c})`$. FORC diagrams resemble the commonly known Preisach diagrams Preisach (1935); Mayergoyz (1986), yet they are model-independent and therefore more general.

Figure 1 shows an FORC diagram of the RFIM at high disorder strength ($`\sigma _\mathrm{R}=5.0>\sigma _{\mathrm{crit}}2.16`$) Perkovic et al. (1999). Note the vertical ridge at $`H_\mathrm{c}1`$ reminiscent of domain-wall nucleation Drossel and Dahmen (1998). A vertical cross-section of the ridge ($`H_\mathrm{c}=1`$) mirrors the distribution of the applied random fields, because these can be viewed as a distribution of random biases acting on the spins when $`\sigma _\mathrm{R}\sigma _{\mathrm{crit}}`$. We have tested this in detail by selecting the random fields from a box distribution. The resulting FORC diagram is qualitatively similar to the Gaussian case, yet a vertical cross-section of the ridge is box-shaped. This is not evident by studying the major hysteresis loop for different disorder-distribution shapes and illustrates the advantages of the FORC method over conventional approaches for studying hysteretic systems.

## IV Switching-Field Histograms

In order to test the predictive power of FORC diagrams, we simulate the RFIM with the zero-temperature dynamics described in Sec. II and store the up- and down-switching fields of the spins along the major hysteresis loop. We then create a histogram of the number of flipped spins for a given pair of up- ($`H_{}`$) and down-switching fields ($`H_{}`$). By changing the variables to the coercivity \[$`H_\mathrm{c}=(H_{}H_{})/2`$\] and bias \[$`H_\mathrm{b}=(H_{}+H_{})/2`$\] of each spin, we obtain a distribution of the coercivities and biases of the spins in the system along the major hysteresis loop.

Figure 2 shows the switching-field histogram (SFH) for the RFIM. One can see a close resemblance with the corresponding FORC diagram presented in Fig. 1. In order to better compare FORC diagram and SFH, in Fig. 3 we present the absolute difference between both diagrams.

Even though the SFH and the FORC diagram of the RFIM differ slightly (see Fig. 3), the main characteristics representing the underlying physical properties of the model are the same (vertical ridge representing multi-domain nucleation). It is interesting that a zeroth-order reversal curve (the major hysteresis loop) contains possibly all the necessary information to reconstruct the first order reversal curves of the system (the FORC diagram).

The differences found between the FORC diagram and the SFH could be due to numerical error in the derivatives of the FORCs because noise is amplified in numerical derivatives considerably. In addition, the differences could be attributed to either hysteron correlations or the failure of the (simple) Preisach picture of hysterons. The latter would imply that a generalization of “classical” hysterons is required.

Figure 2 also illustrates how the re-parametrization of the major hysteresis loop in terms of an SFH shows more details about the microscopic structure of the system. The gained information is similar to the information provided by an FORC diagram, yet the computation of an SFH is considerably faster than calculating an FORC diagrams (generally $`10^2`$ times faster) and involves no numerical derivatives of the data, thus reducing numerical error.

## V Frustrated Systems

As the random-field Ising model is a random magnet with no frustration, we also calculate the FORC diagram and SFH for the 3D Edwards-Anderson Ising spin glass Binder and Young (1986) (EASG). Due to frustration, a spin can flip more than twice along the full hysteresis loop. With the current definition of the SFH this is not taken into account and differences to an FORC diagram are expected.

The Hamiltonian of the Edwards-Anderson Ising spin glass is given by Eq. (1) where the $`J_{ij}`$ are nearest-neighbor interactions chosen according to a Gaussian distribution with zero mean and standard deviation unity, and $`h_i=0`$ $`i`$. $`H`$ represents the externally applied magnetic field and periodic boundary conditions are applied. For the simulations we use the zero-temperature algorithm presented in Sec. II. Frustration is introduced by the random signs of the interactions $`J_{ij}`$.

Figure 4 shows the FORC diagram of the EASG. One can see a pronounced ridge along the $`H_\mathrm{c}`$-axis together with an asymmetric feature at small coercivities. The underlying details of the EASG FORC diagram have been discussed elsewhere Katzgraber et al. (2002).

In Fig. 5 the SFH of the EASG is shown. Note that the asymmetry present in the FORC diagram in Fig. 4 is lost. The weight of the asymmetric part of the FORC diagram shifts to the ridge at $`H_b=0`$.

Although some of the features in the FORC diagram (Fig. 4) of the EASG are missing in the corresponding SFH (Fig. 5), the horizontal ridge reminiscent of the underlying reversal-symmetry of the Hamiltonian Katzgraber et al. (2002) is conserved. In particular, by comparing the FORC diagram and SFH one can study the effects of frustration on the hysteretic behavior of a spin glass.

## VI Conclusions

By re-parameterizing the major hysteresis loop with a switching-field histogram we show that for systems with no frustration (random-field Ising model) the SFH closely resembles the FORC diagram. Small differences can be attributed to numerical error in the calculation of an FORC diagram, hysteron correlations, or the breakdown of the hysteron picture.

SFHs show more details about the system than the major hysteresis loop and are considerably faster to computer than FORC diagrams. Therefore they are an efficient alternative in order to study the microscopic distributions of coercivity and bias of the switching units.

Because the switching fields of the underlying microscopic switching units have to be recorded for the computation of an SFH, the method is in general limited to numerical studies of hysteretic systems. Experimental applicability might be possible with synthetic particulate samples Liu et al. (2002) where the individual switching units can be traced during the magnetic field sweep.

We also present results on the (frustrated) Edwards-Anderson Ising spin glass. We show that there are clear differences between the FORC diagram and the SFH because SFHs do not take into account multiple switching events of the spins, a hallmark of spin glasses. We suggest these differences can be used to quantify the effects of frustration in FORC diagrams.

The FORC and SFH methods promise to be powerful tools to “fingerprint” hysteretic systems. Still, the breadth of information they provides remains to be understood fully.

## Acknowledgments

We would like to thank F. Pázmándi and R. T. Scalettar for discussions.


text:
# Network dependence of strong reciprocity

## 1 Homo oeconomicus versus homo reciprocans

The assumption of self-interest as a motivation for social and economic behavior is widely used as a guiding principle for social modeling. In a game theory context the idea of maximization of self-interest leads to the notion of (noncooperative) Nash equilibrium. A strategy is a Nash equilibrium if no player can improve his payoff by changing his strategy, when the strategies of the other players are fixed.

Given any environment situation, in a Nash equilibrium solution, each player tries to maximize his gains regardless of what happens to the other players. It is the rational expectations attitude of what has been called the Homo oeconomicus, a notion which is at the basis of many theoretical economics constructions. Whether this is a realistic notion when applied to human societies is an important issue. Experiments have been carried out and, in many cases, when played by human players, games have outcomes very different from the Nash equilibrium points. An interesting case is the ultimatum game. A simplified version of this game is the following:

One of the players (the proposer $`P`$) receives 100 coins which he is told to divide into two non-zero parts, one for himself and the other for the other player (the responder $`R`$). If the responder accepts the split $`\left(R_0\right)`$, it is implemented. If the responder refuses $`\left(R_1\right)`$, nothing is given to the players. Consider, for example, a simple payoff matrix corresponding to two different proposer offers $`\left(P_0\text{ and }P_1\right)`$

|  | $`R_0`$ | $`R_1`$ |
| --- | --- | --- |
| $`P_0`$ | $`a,c`$ | $`0,0`$ |
| $`P_1`$ | $`b,b`$ | $`0,0`$ |

(1)
with $`ac`$, $`a+c=2b`$ (for example $`a=99,c=1,b=50`$).

The unique Nash equilibrium is $`(P_0,R_0)`$, corresponding to the payoffs $`(a,c)`$. However, when the game is played with human players, such greedy proposals are most often refused, even in one-shot games where the responder has no material or strategic advantage in refusing the offer. Based on this and similar results in other situations (public goods games, etc), Bowles and Gintis  developed the notion of strong reciprocity (Homo reciprocans) as a better model for human behavior. Homo reciprocans would come to social situations with a propensity to cooperate and share but would respond to selfish behavior on the part of others by retaliating, even at a cost to himself and even when he could not expect any future personal gains from such actions. This should be distinguished from cooperation in a repeated game or reciprocal altruism or other forms of mutually beneficial cooperation that can be accounted for in terms of self-interest.

The same authors, in collaboration with a group of anthropologists, conducted a very interesting “ultimatum game experiment” in many small-scale societies around the world. Homo oeconomicus is rejected in all cases and consistently different results are obtained in different societies, the players’ behavior being strongly correlated with existing social norms and the market structure in their societies. This and other experiments   strongly suggest that human decision problems involve a mixture of self-interest and a background of (internalized) social norms  .

Strong reciprocity is a form of altruism  in that it benefits others at the expense of the individual that exhibits this trait. Monitoring and punishing selfish agents or norm violators is a costly (and dangerous) activity without immediate direct benefit to the agent that performs it. It would be much better to let others do it and to reap the social benefits without the costs.

Strong reciprocator agents contribute more to the group than selfish ones and they sustain the cost of monitoring and punishing free-riders. For this reason it was thought that the strong reciprocity trait could not invade a population of self-interested agents, nor could it be maintained in a stable population equilibrium. To counter this belief, Bowles and Gintis  developed a simple (mean-field type) model that might apply to the structure of the small hunter-gatherer bands of the late Pleistocene. Taking the view that the strong reciprocity trait has a genetic basis, this would be a period long enough to account for a significant development in the modern human gene distribution. The model would give an evolutionary explanation of the phenomenon. Of course, if instead of gene-based, strong reciprocity is culturally inherited, emergence and (or) modification of this trait could be much faster.

Because I intend to explore the influence of the social (network) structure on the evolution of strong reciprocity, I will start by discussing a simplified version of the Bowles-Gintis model. The main simplification is that migration in and out of the evolving group to an outside pool of agents is not considered. The consideration of these migrations may be of interest for a realistic picture of the hunter-gatherer bands of the Pleistocene, but not for the general picture of strong reciprocity in a wider society. By simplifying and somewhat enlarging the punishment scenario (beyond ostracism) of the Bowles-Gintis model and framing it as a replicator one-dimensional map, a clear view is obtained of its dynamical aspects.

## 2 Emergence of strong reciprocity. The Bowles-Gintis model

One considers a population of size $`N`$ with two species of agents, one denoted reciprocators (R-agents) and the other self-interested (S-agents). In a public goods activity each agent can produce a maximum amount of goods $`q`$ at cost $`b`$ (with goods and costs in fitness units). The benefit that an S-agent takes from shirking public goods work is the cost of effort $`b\left(\sigma \right)`$, $`\sigma `$ being the fraction of time the agent shirks. The following conditions hold

$$b\left(0\right)=b,b\left(1\right)=0,b^{^{}}\left(\sigma \right)<0,b^{^{\prime \prime }}\left(\sigma \right)>0$$
(2)
Furthermore $`q\left(1\sigma \right)>b\left(\sigma \right)`$ so that, at every level of effort, working helps the group more than it hurts the worker.

For $`b\left(\sigma \right)`$ one chooses

$$b\left(\sigma \right)=\frac{2}{2\sigma 1+\sqrt{1+4/b}}\frac{2}{1+\sqrt{1+4/b}}$$
(3)
which satisfies the constraints (2).

R-agents never shirk and punish each free-rider at cost $`c\sigma `$, the cost being shared by all R-agents. For an S-agent the estimated cost of being punished is $`s\sigma `$, punishment being ostracism or some other fitness decreasing measure. $`s`$ is the weight given by an S-agent to the punishment probability. It may or may not be the same as the actual fitness cost of punishment. Each S-agent chooses $`\sigma `$ (the shirking time fraction) to minimize the function

$$B\left(\sigma \right)=b\left(\sigma \right)+sf\sigma q\left(1\sigma \right)\frac{1}{N}$$
(4)
$`f`$ being the fraction of R-agents in the population, $`f\sigma `$ is the probability of being monitored and punished. The last term is the agent’s share of his own production. The value $`\sigma _S`$ that minimizes $`B\left(\sigma \right)`$ is

$$\sigma _S=\mathrm{max}(\mathrm{min}(\frac{1}{2}\sqrt{\frac{1}{4}+\frac{1}{b}}+\frac{1}{\sqrt{sf+\frac{q}{N}}},1),0)$$
(5)
The contribution of each species to the population in the next time period is proportional to its fitness given by

$$\begin{array}{ccc}\pi _S^{^{}}\left(f\right)\hfill & =\hfill & q\left(1\left(1f\right)\sigma _S\right)b\left(\sigma _S\right)\gamma f\sigma _S\hfill \\ \pi _R^{^{}}\left(f\right)\hfill & =\hfill & q\left(1\left(1f\right)\sigma _S\right)bc\left(1f\right)\frac{N\sigma _S}{Nf}\hfill \end{array}$$
(6)
for S- and R-agents. The baseline fitness is zero, that is, $`\pi _{S,R}=\mathrm{max}(\pi _{S,R}^{^{}},0)`$.

The first term in both $`\pi _S^{^{}}`$ and $`\pi _R^{^{}}`$ is the benefit arising from the produced public goods and the second term the work effort. The last terms represent the fitness cost of punishment for S-agents and the cost incurred by R-agents.

$`\gamma =1`$ corresponds to ostracism from the group, other values to general coercive measures affecting the fitness of S-agents. The last term in $`\pi _R^{^{}}`$ emphasizes the collective nature of the punishment. Notice however the improbable heavy punishing burden put on reciprocators when in small number.

Finally one obtains111Here replicator dynamics is used for the population evolution. Notice that Bowles and Gintis use a different (incremental) dynamics. a one-dimensional map for the evolution of the fraction of R-agents

$$f_{\text{new}}=f\frac{\mathrm{\Pi }_R\left(f\right)}{\left(1f\right)\mathrm{\Pi }_S+f\mathrm{\Pi }_R\left(f\right)}$$
(7)
Figs. 1 and 2 display this map, as well as $`\sigma _S\left(f\right)`$, $`\mathrm{\Pi }_R\left(f\right)`$ and $`\mathrm{\Pi }_S\left(f\right)\mathrm{\Pi }_R\left(f\right)`$ for two different values of $`\gamma `$, the other parameters being the same. They show the general behavior of the map in Eq.(7). If $`\gamma `$ (the fitness impact of punishment) is large enough, the map has an unstable fixed point $`A`$ at $`f_A`$ and a left-stable one $`B`$ at $`f_B`$. Between $`f_B`$ and $`1`$ there is a continuum of marginally stable fixed points. For smaller $`\gamma `$ the region between $`f_A`$ and $`f_B`$ (where $`\mathrm{\Pi }_S\mathrm{\Pi }_R`$ is negative) disappears and only the marginally stable fixed points remain. In both cases the asymptotic behavior corresponds either to $`f=0`$ (and $`\sigma _S=1`$) or to $`f`$ between $`0`$ and $`1`$ but $`\sigma _S=0`$. That is, in this second case, both reciprocators and shirkers remain in the population but shirkers choose not to shirk because the minimum of $`B\left(\sigma \right)`$ is at $`\sigma _S=0`$.

For an initial $`f`$ smaller than $`f_A`$ the fraction of reciprocators falls very rapidly to zero. This reflects the (maybe unrealistic) fact that in this case a very small number of reciprocators has to carry the burden of punishing very many shirkers.

Hence, from the point of view of intragroup dynamics, either reciprocators are completely eliminated from the population or they remain in equilibrium with a probably large number of shirkers, which do not shirk for fear of being punished. Therefore intragroup dynamics, by itself, cannot explain how the reciprocator attitude might have become a dominant trait. However when very many groups are considered, for example assembled at random from a pool containing both reciprocators and shirkers  , then only the groups that contain at the start a fraction $`f`$ greater than $`f_A`$ will have in the end a nonzero fitness. In all others, S-agents invade the population and suffer a “tragedy of the commons” situation with final zero fitness. Therefore from an intergroup dynamics perspective the groups with reciprocators tend to dominate and impose an above average predominance of the reciprocator trait.

Although the model, together with intergroup dynamics, explains why strong reciprocity is an adaptive trait , the marginally stable nature of the (above $`f_B`$) fixed points also suggest that the shirker trait is never eliminated and will remain in the population.

Small independent groups assembling and disassembling is a likely scenario for the development of the reciprocator trait. In this sense the hunter-gatherer bands of the Pleistocene might have indeed provided the appropriate environment for the evolution of the trait, whether gene-coded or culturally-inherited.

It is well known that group size affects monitoring in public goods provision . Therefore, a natural question is what happens when, later on, the Pleistocene reciprocators and their fellow shirkers become imbedded into a larger society. Monitoring and punishment of shirkers by reciprocators necessarily looses its global collective nature. Once monitoring looses its global nature, it becomes the business of the neighbors of the shirker. In addition to the individual cost of monitoring and (or) punishing free-riders, such punishing requires an amount of force that, in particular, insures the effectiveness of the punishment and on the other hand puts the punisher safe from direct retaliation from the violator. This is one of the reasons for the creation of central authorities for this purpose. However if central authorities have enough force to implement punishment without retaliation, they are at times quite ineffective at monitoring. Also laws and central authorities, on the role of reciprocators play, a role in the control of serious offenses, but not so much on the day to day monitoring of public goods work. Therefore in a large society the nature of the control performed by the neighbors is certainly going to play a role on the evolution of the reciprocator trait.

If the trait is genetically encoded, maybe the wider societies developed by modern man had no time to make significant changes on its structure. However if it is (at least in part) culturally inherited then a much shorter time scale may be involved. What about the big city tales of a guy being mugged in full daylight while a crowd of passersby moves along quite indifferent to the event? Is it the ($`1f`$) remnants of non-reciprocators in the population or are we watching the emergence of Homo Oeconomicus in his full glory? Or is it something else?

## 3 Network dependence of strong reciprocity

To explore the possible effect of the social network structure on the evolution of strong reciprocity I will consider a agent-based model, which later on will be interpreted in a mean field sense similar to the model in Section 2.

As before one considers R-agents and S-agents and the monitoring function performed by R-agents is kept at the neighbors level. However punishment is only implemented if at least two neighbors are willing to do so. It is the same as to say that punishing a norm-violator cannot be an individual action, but requires a minimal social power and consensus. The need to be close to monitor and the need for agreement of at least two neighboring reciprocators to implement punishment, immediately suggests that the structure of the network is going to play a role on the evolution of the group. The following is the mathematical coding of this idea:

As before one has two agent species (S-agents and R-agents), the fraction of R-agents being $`f`$. The agents are placed in a network where, on average, each agent is connected to $`k`$ other agents. $`k`$ is called the degree of the network. To the whole population of dimension $`N`$ one associates 3 $`N`$dimensional vectors, $`Wk`$, $`Pu`$, $`Cpu`$. $`Wk`$ is called the work vector, $`Pu`$ the punishment vector and $`Cpu`$ the cost of punishment vector.

The link structure of the network is chosen as in the $`\beta `$model of Watts and Strogatz . Starting from a regular ring structure where each node is symmetrically connected to its $`k`$ closest neighbors, each link is examined in turn and, with probability $`\beta `$, replaced by a random link to some other node in the network.

At time zero, $`fN`$ R-agents and $`\left(1f\right)N`$ S-agents are placed at random in the network. The local neighborhood of agent $`i`$, that is the set of other agents connected to $`i`$, is denoted $`\mathrm{\Gamma }_i`$. The entries of the vectors $`Wk`$, $`Pu`$, $`Cpu`$ are then computed as follows:

\# For the $`Wk`$ vector

R-agents; $`Wk\left(i\right)=1`$

S-agents; $`Wk\left(i\right)=\frac{n_R\left(i\right)}{k}`$, where $`n_R\left(i\right)`$ is the number of R-agents connected to this S-agent, $`n_R\left(i\right)=\mathrm{\#}\{j:jR,j\mathrm{\Gamma }_i\}`$

\# For the $`Pu`$ vector

R-agents; $`Pu\left(i\right)=0`$

S-agents; $`Pu\left(i\right)=n_P\left(i\right)\left(1Wk\left(i\right)\right)`$, where $`n_P\left(i\right)`$ is the number of pairs of R-agents in $`\mathrm{\Gamma }_i`$ which are also neighbors among themselves, $`n_P\left(i\right)=\mathrm{\#}\{(j,k):(j,k)R,(j,k)\mathrm{\Gamma }_i,j\mathrm{\Gamma }_k\}`$ and $`\left(1Wk\left(i\right)\right)`$ is the shirking fraction.

\# For the $`Cpu`$ vector

R-agents; $`Cpu\left(i\right)=_{kS}n_C(i,k)\left(1Wk\left(k\right)\right)`$ where $`n_C(i,k)`$ is the number of times that the agent $`i`$ is in a R-pair punishing an S-agent $`k`$, $`n_C(i,k)=\mathrm{\#}\{(i,j):kS,(i,j)R,(i,j)\mathrm{\Gamma }_k,j\mathrm{\Gamma }_i\}`$

S-agents; $`Cpu\left(i\right)=0`$

Summarizing:

Each reciprocator, on detecting an S-agent $`k`$, looks for another reciprocator in his own neighborhood also connected to $`k`$. If he finds one he punishes $`k`$ by an amount proportional to the fraction of shirking. An S-agent may be punished several times by all different pairs of reciprocators in his neighborhood.

The amount of work that an S-agent does is inversely proportional to the number of reciprocators in his neighborhood. However lack of communication between neighboring reciprocators may make the probability of punishment much smaller.

The (average) fitness of R-agents and S-agents is

$$\pi _R^{^{}}=\frac{q}{N}\underset{all}{}Wk\left(i\right)\frac{b}{fN}\underset{iR}{}Wk\left(i\right)\frac{c}{fN}\underset{iR}{}Cpu\left(i\right)$$
(8)
$$\pi _S^{^{}}=\frac{q}{N}\underset{all}{}Wk\left(i\right)\frac{b}{\left(1f\right)N}\underset{iS}{}Wk\left(i\right)\frac{\gamma }{\left(1f\right)N}\underset{iS}{}Pu\left(i\right)$$
(9)
The baseline fitness is zero, that is

$$\pi _{R,S}=\mathrm{max}(\pi _{R,S}^{^{}},0)$$
(10)
Once the fitness is computed the replicator equation

$$f_{\text{new}}=f\frac{\pi _R}{f\pi _R+\left(1f\right)\pi _S}$$
(11)
is applied and a new cycle starts with a new random distribution, on the network, of $`Nf_{\text{new}}`$ R-agents and $`N\left(1f_{\text{new}}\right)`$ S-agents.

Running this agent model for several values of $`\beta `$ and, in each case, for random initial $`f_0`$’s one finds two separate regions in the $`(f_0,\beta )`$ plane (Fig.3). In region 1 the evolution drives $`f`$ towards zero as well as the overall fitness $`\pi `$ (Example in Fig.4a)

$$\pi =f\pi _R+\left(1f\right)\pi _S$$
(12)
In region 2 there is an asymptotic nonzero value for $`f`$ and for the fitness (Example in Fig.4b).

As $`\beta `$ increases it becomes less likely to have a stable nonzero $`f`$. the origin of this effect is clear. Although $`\beta `$rewiring maintains the average degree of the network, the probability of two neighbors of an agent to be themselves neighbors decreases. Therefore it becomes increasingly difficult for reciprocators to find local consensus for the punishment of S-agents.

The average probability of two R-neighbors of a network node in S to be themselves neighbors, is called the (relative) clustering coefficient,

$$C_R=\overline{\frac{n_P\left(i\right)}{\left(\begin{array}{c}\mathrm{\#}\left\{\mathrm{\Gamma }_iR\right\}\\ 2\end{array}\right)}}$$
(13)
$`\left(\begin{array}{c}\mathrm{\#}\left\{\mathrm{\Gamma }_iR\right\}\\ 2\end{array}\right)`$ being the maximum possible number of links between the R-neighbors of S-agent $`i`$. The network clustering coefficient is related to the notion of transitivity used in the sociological literature.

For the $`\beta `$rewiring model, the clustering may be estimated from the number $`\mathrm{\Phi }`$ of shortcuts which in this case is proportional to $`\beta `$ .

$$C_\beta (\mathrm{\Phi },k)=\frac{\frac{3}{4}\left(1\mathrm{\Phi }\right)^2\left(k\frac{2}{3}\right)\left(1\mathrm{\Phi }\right)}{k1}$$
(14)
Therefore a mean field version of the agent model may be written as follows

$$\mathrm{\Pi }_R^{^{}}=q\left(1\left(1f\right)\sigma _S\left(f\right)\right)b\left(\sigma _S\left(f\right)\right)\gamma fC_\beta (\mathrm{\Phi },fk)\sigma _S\left(f\right)$$
(15)
$$\mathrm{\Pi }_S^{^{}}=q\left(1\left(1f\right)\sigma _S\left(f\right)\right)bc\left(1f\right)\frac{fk}{2}C_\beta (\mathrm{\Phi },fk)\sigma _S\left(f\right)$$
(16)
Notice the term $`fk`$ in $`C_\beta (\mathrm{\Phi },fk)`$ and in the cost of punishment term in $`\mathrm{\Pi }_S^{^{}}`$. It reflects the fact that neighborhood relations for reciprocators are to computed on their subnetwork of size $`fN`$.

$`b\left(\sigma \right)`$ is as in Eq.(3) with $`\sigma _S`$ being computed to minimize

$$B\left(\sigma \right)=b\left(\sigma \right)+sfC_\beta (\mathrm{\Phi },fk)\sigma q\left(1\sigma \right)\frac{1}{N}$$
(17)
This mean field version gives results identical to the agent-based model. Clustering appears therefore as the determining network parameter driving the evolution of the reciprocator trait.

## 4 Conclusions

1 - With a structure of small groups with collective monitoring of the agents’ activities, the fitness difference between groups with a sizable amount of reciprocators and those where they have disappeared, makes the emergence of the strong reciprocity trait a likely event.

However rather than being completely invaded by reciprocators, maintenance of a certain amount of self-interested types is also likely, which only cooperate for fear of being punished. If, at a later stage, the social structure changes, they may be a source of instability and invade the population.

2 - In a large population, monitoring of the public goods behavior of the agents cannot be a fully collective activity, rather being the chore of those in close contact with the free-riders. Punishment of free-riders also requires a certain amount of local consensus among reciprocators. Therefore the clustering nature of the society may play an important role in the maintenance and evolution of the reciprocator trait.

Maybe the indifferent passersby that let the poor guy be mugged are not yet homo oeconomicus. Maybe they are just reciprocators in the middle of strangers with whom they do not communicate nor trust. A clustering problem.

3 - Culturally-inherited traits may have a much faster dynamics than gene-based ones. Modern societies are “small worlds” in the sense of short path lengths but not necessarily in the sense of also maintaining a high degree of clustering. Therefore if the reciprocator trait has a high cultural component, it may very well happen that, eventually, we will see homo oeconomicus leaving the benches of economy classes for a life on the streets.