text:
# Large scale correlations in normal non-Hermitian matrix ensembles

## 1 Introduction

A matrix $`M`$ is called normal if it commutes with its Hermitian conjugated: $`[M,M^{}]=0`$. The both matrices can be simultaneously diagonalized, the eigenvalues being complex numbers. The partition function of normal matrices has the general form

$$Z_N=_{\mathrm{normal}}𝑑\mu (M)e^{\frac{1}{\mathrm{}}\mathrm{tr}W(M,M^{})}.$$
(1)
Here $`\mathrm{}`$ is a parameter, and the measure of integration over normal $`N\times N`$ matrices is induced from the flat metric on the space of all complex matrices. Like in other random matrix models, the large $`N`$ limit of interest implies $`N\mathrm{}`$, $`\mathrm{}0`$, while $`N\mathrm{}`$ stays finite.

In a particular case, when Laplacian of the potential $`W(z,\overline{z})`$ is a constant in a big domain of a complex plane, i.e.,

$$W=MM^{}+V(M)+\overline{V}(M^{}),$$
(2)
where $`V(z)`$ and $`\overline{V}(z)`$ are holomorphic functions, the normal matrix ensemble is equivalent to the ensemble of general complex matrices. It generalizes the Ginibre-Girko Gaussian ensemble . In this, perhaps the most interesting case for applications, the model bears some formal similarities with the model of two Hermitian random matrices  and the matrix quantum mechanics in a singlet sector . Unlike models of few Hermitian matrices, the normal matrix model is integrable for a general class of potentials, not only of the form (2).

Applications and studies of matrix ensembles with complex eigenvalues are numerous. A large list of references can be found in recent papers . New applications to diffusion limited growth models (Laplacian growth) , complex analysis - and Quantum Hall effect  were found recently.

Despite of a comprehensive literature on the model of one and two Hermitian matrices, interest to the normal matrix ensemble, first introduced in  and further studied in , as well as to a model of general complex matrices with potential (2), just starts to grow. In this paper we revisit the ensemble of normal random matrices to calculate the large $`N`$ limit of correlation functions under condition that separation between arguments is much more than an average distance between eigenvalues (macroscopic, or smoothed correlations). Short scale (microscopic) correlations in the ensemble with the potential (2) are well studied (see e.g.  for a review).

At large $`N`$ the eigenvalues of the random matrices are distributed within a domain (with sharp edges) of a complex plane with a density proportional to $`\mathrm{}^1\mathrm{\Delta }W`$, where $`\mathrm{\Delta }`$ is the 2D Laplace operator. We assume that the support of eigenvalues is a single connected bounded domain $`D`$, and that the boundary is a Jordan curve. Analytical properties of this curve determines the correlation functions.

We will show that correlation functions are expressed through the objects of the Dirichlet boundary problem of the domain complimentary to the support of eigenvalues. Namely, the two-point correlations are expressed through the Dirichlet Green function, while higher multi-point correlations are expressed through the Neumann jump on the boundary, through the Bergman kernel and through the curvature of the boundary. The objects similar to the correlation functions previously appeared in studies of thermal fluctuations in classical confined Coulomb plasma  (see also  for a review) and in recent studies of integrable structure of the Dirichlet boundary problem in . If the potential $`W`$ is such that the support of eigenvalues collapses to a cut (or cuts), the normal matrix model reproduces known large $`N`$-limit features of the Hermitian matrix model.

We compute two-, three- and four-point connected functions of density of eigenvalues in the leading order in $`N`$, in a manner that the algorithm of computing higher correlation functions becomes clear. Connected density correlators are localized at the edge of the support of eigenvalues and show some universal features. They depend on the potential only through a shape of the support of the eigenvalues and boundary values of a finite number of derivatives of the potential. In the literature on Hermitian random matrices the universal character of correlations has been emphasized in Refs. . The two-point functions are distinguished since they depends on the support of the eigenvalues only. The two-point function has been previously computed in Refs.  in the course of studies of confined Coulomb plasma.

Other aspects, including integrable structure of the model (discussed already in ), semiclassical properties of biorthogonal polynomials, critical points (a boundary with cusps), disconnected supports of eigenvalues, correlations at short distances, etc., will be addressed elsewhere.

## 2 Preliminaries

### The measure of normal matrices.

The measure in (1) is induced from the flat metric $`\delta M^2=\text{tr}(\delta M\delta M^{})`$ in the space of all complex matrices. Formally, one can write $`d\mu (M)=\delta ([M,M^{}])dMdM^{}`$, where the delta-function selects the subspace of normal matrices. To introduce coordinates in this subspace, one uses of the decomposition $`M=UZU^{}`$ , where $`U`$ is a unitary matrix and $`Z=\text{diag}(z_1,\mathrm{},z_N)`$ is a diagonal matrix of eigenvalues . The measure is then given by

$$d\mu (M)=\frac{d\mu _0(U)}{N!\text{Vol}(𝒰(N))}|\mathrm{\Delta }_N(z)|^2\underset{i=1}{\overset{N}{}}d^2z_i,$$
(3)
here $`d^2zdxdy`$ for $`z=x+iy`$, $`d\mu _0`$ is the Haar measure on the unitary group $`𝒰(N)`$, and $`\mathrm{\Delta }_N(z)=det(z_j^{i1})_{1i,jN}=_{i>j}^N(z_iz_j)`$ is the Vandermonde determinant.

If $`A(M)`$ is any invariant function (i.e., a symmetric function of eigenvalues) of a matrix $`M`$ (and $`M^{}`$), then the mean value

$$A=\frac{𝑑\mu (M)A(M)e^{\frac{1}{\mathrm{}}W}}{𝑑\mu (M)e^{\frac{1}{\mathrm{}}W}}$$
(4)
is expressed through the integral over eigenvalues

$$A=\frac{1}{N!Z_N}A(z)|\mathrm{\Delta }_N(z)|^2\underset{j=1}{\overset{N}{}}\left(e^{\frac{1}{\mathrm{}}W(z_j,\overline{z}_j)}d^2z_j\right),$$
(5)
where the partition function is

$$Z_N=\frac{1}{N!}|\mathrm{\Delta }_N(z)|^2\underset{j=1}{\overset{N}{}}\left(e^{\frac{1}{\mathrm{}}W(z_j,\overline{z}_j)}d^2z_j\right).$$
(6)
### Potential.

For notational simplicity, we shall write simply $`W(z)`$ instead of $`W(z,\overline{z})`$. We assume that $`W`$ is a real-valued function with (at least local) minimum at the origin and set $`W(0)=0`$ for convenience. We also assume that the integral (6) converges and that $`W`$ is a regular function in both variables at the origin. We shall also set

$$\sigma (z)=\frac{1}{\pi }_z_{\overline{z}}W(z)=\frac{1}{4\pi }\mathrm{\Delta }W,$$
(7)
and assume that $`\sigma (z)>0`$.

A special interesting case  occurs if the potential is harmonic in some big domain around the origin. Then in this domain $`\sigma (z)=\text{constant}`$, say, set to be $`\pi ^1`$, and

$$W=|z|^2+V(z)+\overline{V}(\overline{z}),$$
(8)
where $`V(z)`$ is a holomorphic function.

### The measure of complex matrices.

If the potential is chosen in the form (8), then the same measure, up to a numerical factor, appears in complex random matrices. In this case the relevant decomposition reads: $`M=U(Z+R)U^{}`$, where $`U`$, $`Z`$ are again unitary and diagonal matrix respectively, and $`R`$ is an upper triangular matrix. The measure (3) acquires a multiplicative factor $`_{ij}dR_{ij}e^{R_{ij}^2}`$ (we use the fact that $`\text{tr}MM^{}=\text{tr}ZZ^{}+\text{tr}RR^{}`$). As a result the representation (5) holds .

### Coherent states of particles in magnetic field.

Most of the known matrix ensembles represent coherent states of $`N`$ fermions. Fermonic representations are well known in the case of Hermitian matrices (see, e.g., ), where fermions live on a line and are confined by a potential.

Complex matrix ensembles also enjoy a fermionic representation . In this case fermions are situated on a plane and occupy the lowest energy level of a strong magnetic field $`B(z)`$, the magnetic field is not necessarily uniform. It is related to the potential by $`B(z)=2\pi \sigma (z)`$. The coherent state of the fully occupied lowest level is

$$\mathrm{\Psi }(z_1,\mathrm{},z_N)=\frac{1}{\sqrt{N!}}\mathrm{\Delta }_N(z)e^{_{j=1}^N\frac{1}{2\mathrm{}}W(z_j)},Z_N=|\mathrm{\Psi }(z_1,\mathrm{},z_N)|^2d^2z_i.$$
(9)
So $`Z_N`$, in this picture, is the normalization factor of the wave function.

### Example: Gaussian model.

In the case

$$W(z)=\pi \sigma |z|^2+2e(t_1z+t_2z^2),\sigma >0,$$

(the Ginibre-Girko ensemble) the partition function can be found explicitly even for finite $`N`$ (see e.g. ):

$$Z_N(\sigma ,t_1,t_2)=Z_N^{(0)}(\pi \sigma )^{\frac{1}{2}(N^2N)}(\pi ^2\sigma ^24t_2\overline{t}_2)^{\frac{1}{2}N^2}\mathrm{exp}\left(\frac{N}{\mathrm{}}\frac{t_1^2\overline{t}_2+\overline{t}_1^2t_2+\pi \sigma |t_1|^2}{\pi ^2\sigma ^24t_2\overline{t}_2}\right),$$
(10)
where $`Z_N^{(0)}=\mathrm{}^{\frac{1}{2}N(N+1)}\pi ^N_{j=1}^{N1}j!`$ is the partition function for $`W=|z|^2`$. Correlation functions in this case are expressed through Hermite polynomials and are known explicitly for any $`N`$.

### Correlation functions of traces.

In this paper, we are interested in correlation functions of products of $`n`$ traces, $`_{i=1}^n\text{tr}f_i(M)`$. Clearly, they are expressed through $`n`$-point correlation functions of the density of eigenvalues,

$$\rho (z)=\mathrm{}\underset{i=1}{\overset{N}{}}\delta ^{(2)}(zz_i).$$

We obviously have:

$$\mathrm{}^n\underset{i=1}{\overset{n}{}}\text{tr}f_i(M)=\rho (z_1)\mathrm{}\rho (z_n)f_1(z_1)\mathrm{}f_n(z_n)\underset{j=1}{\overset{n}{}}d^2z_j.$$
(11)
So, the density correlation functions carry the necessary information.

It is customary to deal with connected part of a correlation function. In case of the 2-point function it is

$$\rho (z_1)\rho (z_2)_{\mathrm{conn}}=\rho (z_1)\rho (z_2)\rho (z_1)\rho (z_2).$$

As $`N\mathrm{}`$, the $`n`$-point correlation function of densities is $`O(1)`$ while the connected part of the $`n`$-point function is $`O(N^{22n})`$.

### The method of functional derivatives.

Correlation functions are the linear responses to a small variation of the potential. Variation of the partition function (6) over a general potential $`W(z)`$ inserts $`\frac{1}{\mathrm{}}_i\delta ^{(2)}(zz_i)`$ into the integral. Then

$$\rho (z)=\mathrm{}^2\frac{\delta \mathrm{log}Z_N}{\delta W(z)},\rho (z_1)\rho (z_2)_{\mathrm{conn}}=\mathrm{}^2\frac{\delta \rho (z_1)}{\delta W(z_2)},$$
(12)
and, more generally,

$$\rho (z_1)\mathrm{}\rho (z_n)_{\mathrm{conn}}=\mathrm{}^{2n}\frac{\delta ^n\mathrm{log}Z_N}{\delta W(z_1)\mathrm{}\delta W(z_n)}.$$

We shall use this method in the following version. Set $`\delta W(z)=\epsilon g(z)`$, where $`g`$ is an arbitrary smooth function on the plane and $`\epsilon 0`$. Then, in the first order in $`\epsilon `$,

$$\mathrm{}^2\delta \text{tr}f(M)=\epsilon \text{tr}f(M)\text{tr}g(M)_{\mathrm{conn}}.$$
(13)
It is often convenient to consider correlations of the potential

$$\phi (z)=\mathrm{log}|zz^{}|^2\rho (z^{})d^2z^{}=\mathrm{}\text{ tr}\left(\mathrm{log}(zM)(\overline{z}M^{})\right),$$
(14)
and of the current field

$$J(z)\phi (z)=\mathrm{}\text{tr}\left(\frac{1}{zM}\right),$$
(15)
rather than correlations of density. The potential is the Bose field or the loop field of a (collective) field theory of the matrix model (more accurately, $`\phi `$ is a negative mode part of a Bose field). The field theory is proved to be successful approach to Hermitian matrix ensembles . Elsewhere we will develop this approach for the non-Hermitian ensembles. The potential is harmonic outside the support of eigenvalues except at infinity where it behaves as $`2N\mathrm{}\mathrm{log}|z|`$. The current is holomorphic outside the support of eigenvalues.

In order to obtain the correlations of the potentials, one has to vary the partition function by $`W(z)W(z)+\epsilon \mathrm{log}|z\zeta |^2`$ where $`\zeta `$ is a parameter. We denote this particular deformation of the potential by $`\delta _\zeta `$:

$$\delta _\zeta W(z)=\epsilon \mathrm{log}|z\zeta |^2.$$
(16)
Under this variation the correlation function changes by insertion of the field $`\phi (\zeta )`$:

$$\mathrm{}^2\delta _\zeta A=\epsilon A\phi (\zeta )_{\mathrm{conn}}.$$
(17)
This is the linear responce relation used in the Coulomb gas theory .

While varying the potential it is important to distinguish a harmonic variation of the potential $`W(z,\overline{z})W(z,\overline{z})+V(z)+\overline{V}(\overline{z})`$, where $`V(z)`$ is a holomorphic function. This variation does not change $`\mathrm{\Delta }W`$. To implement a harmonic variation one may extend the potential by adding a harmonic function $`WW+2e_{k1}t_kz^k`$ and apply the operator $`D(z)=_{k1}\frac{z^k}{k}_{t_k}`$ used in Refrs. -. Then correlators of holomorphic parts of the potential $`\varphi (z)=\mathrm{}\text{tr}\mathrm{log}(zM)`$ are

$$\begin{array}{c}\varphi (z)=\mathrm{}N\mathrm{log}z\mathrm{}^2D(z)\mathrm{log}Z_N,\hfill \\ \\ \underset{i=1}{\overset{n}{}}\varphi (z_i)_{\mathrm{conn}}=(1)^n\mathrm{}^{2n}\underset{i=1}{\overset{n}{}}D(z_i)\mathrm{log}Z_N,n2.\hfill \end{array}$$
(18)
### Dirichlet boundary problem.

We list some elements of the external Dirichlet boundary problem, which are extensively used below. More details can be found in .

Let $`D`$ be a closed connected domain of the complex plane bounded by a smooth curve. Given a real analytic function $`f(z)`$ in a vicinity of the boundary, we may restrict it to the boundary of $`D`$. The external Dirichlet problem is to find a harmonic function in the exterior of $`D`$, whose value on the boundary is $`f(z)`$. We call this harmonic function the harmonic extension of $`f(z)`$ to the exterior and denote it by $`f^H(z)`$. It is given by the

* Dirichlet formula:

$$f^H(z)=\frac{1}{2\pi }_n^{}G(z,z^{})f(z^{})|dz^{}|,$$

where $`_n^{}`$ is the normal derivative at the boundary with respect to the second variable, with the normal vector being directed to the exterior of the domain $`D`$.

The Green function $`G(z,z^{})`$ for the exterior problem is a harmonic function everywhere outside $`D`$ including infinity except the point $`z=z^{}`$, where it has a logarithmic singularity: $`G(z,z^{})\mathrm{log}|zz^{}|`$ as $`zz^{}`$. If $`z^{}\mathrm{}`$, then $`G(z,\mathrm{})\mathrm{log}|z|`$. The Green function is symmetric in $`z,z^{}`$ and vanishes on the boundary. In particular, the harmonic extension of $`\mathrm{log}|z\zeta |`$ in $`z`$ is

$$(\mathrm{log}|z\zeta |)^H=\mathrm{log}|z\zeta |G(z,\zeta )+G(z,\mathrm{}).$$
(19)
If the point $`\zeta `$ happens to be inside, $`G(z,\zeta )`$ in this formula is understood to be null.
* The Green function can be expressed through a conformal map $`w(z)`$ from $`𝐂D`$ (the exterior of $`D`$) onto the exterior of the unit disk:

$$G(z,z^{})=\mathrm{log}\left|\frac{w(z)w(z^{})}{1\overline{w(z)}w(z^{})}\right|.$$
(20)
This formula does not depend on the normalization of the map. It is convenient to fix $`w(z)`$ by the condition that it sends infinity to infinity and the coefficient in front of the leading term as $`z\mathrm{}`$ is real positive.
* Neumann external jump is the difference between normal derivatives of a smooth function at the boundary and its harmonic extension. The Neumann external jump operator $`𝖱`$ acts as follows:

$$f(z)(𝖱f)(z)=_n^{}(f(z)f^H(z)),zD.$$
(21)
The upperscript indicates that the derivative is taken in the exterior of the boundary. As it follows from the Dirichlet formula, the Neumann jump is an integral operator on the boundary curve with the kernel given by normal derivative of the Green function in both arguments: $`(𝖱f)(z)=_nf(z)+\frac{1}{2\pi }_n_n^{}G(z,z^{})f(z^{})|dz^{}|`$. In fact this integral is not yet well-defined since the kernel has a second order pole on the contour. The operations of taking normal derivative and contour integration do not commute. The above formula has to be understood as $`_n_n^{}G(z,z^{})f(z^{})|dz^{}|`$ Alternatively the formula is understood as the principal value integral:

$$(𝖱f)(z)=_nf(z)+\frac{1}{2\pi }\text{P.V.}_n_n^{}G(z,z^{})(f(z)f(z^{}))|dz^{}|.$$
(22)
* Hadamard formula describes deformation of the Dirichlet Green function under deformation of the domain.

A change of the boundary can be characterized by its normal displacement $`\delta n(z)`$, such that $`\delta n(z)`$ is a continuous function on the boundary (see Fig. 1). The Hadamard formula  expresses the deformation of the Green function through the Green function itself:

$$\delta G(z_1,z_2)=\frac{1}{2\pi }_D_nG(z_1,\xi )_nG(z_2,\xi )\delta n(\xi )|d\xi |.$$
(23)
A corollary of the Hadamard formula is the variation of the boundary value of a harmonic function under variation of the boundary. It reads

$$\delta f^H(z)=(𝖱f)(z)\delta n(z),zD.$$
(24)

## 3 Large $`N`$ limit

The large $`N`$ limit is understood as $`N\mathrm{}`$, $`\mathrm{}0`$ while $`\mathrm{}N`$ is kept finite. The expansion in $`N^1`$ is then equivalent to the expansion in $`\mathrm{}`$. We call it a semiclassical limit.

### Semiclassical density and the support of the eigenvalues.

To elaborate the semiclassical limit, we have to find the maximum of the integrand $`|\mathrm{\Psi }_N(z_1,\mathrm{},z_N)|^2`$ in (9). At finite $`N`$, it is given by the conditions $`\mathrm{log}|\mathrm{\Psi }_N|/z_i=0`$ for every $`i`$:

$$\mathrm{}\underset{j=1,i}{\overset{N}{}}\frac{1}{z_iz_j}+_{z_i}W(z_i)=0,$$

or by the equation

$$_z\left(\phi _0(z)+W(z)\right)=_{\overline{z}}\left(\phi _0(z)+W(z)\right)=0,zD$$
(25)
in the limit. Here $`\phi _0(z)=\phi (z)=\mathrm{log}|zz^{}|^2\rho (z)(z^{})d^2z^{}`$ and $`\rho (z)`$ are expectation values of the potential and the density in the leading semiclassical approximation. The last equation holds inside the support of eigenvalues, $`D`$. It does not hold outside it. In the semiclassical approximation, $`D`$ is a domain with a well-defined sharp boundary determined by the potential $`W`$. An assumption that the support is a compact connected domain bounded by a Jordan curve implies restrictions on the potential. The restrictions do not reduce the number of parameters in the potential but rather ranges of their variation. We do not discuss them here. The assumption is valid, for example, for small perturbations of the Gaussian potential.

The solution for density is obtained by applying $`_{\overline{z}}`$ to eq. (25):

$$\rho (z)=\{\begin{array}{cc}\sigma (z)\hfill & \text{if}zD\hfill \\ & \\ 0\hfill & \text{ if}z𝐂D.\hfill \end{array}$$
(26)
The function $`\sigma `$ is introduced in (7). Consequently, $`\phi _0(z)={\displaystyle _D}\mathrm{log}|zz^{}|^2\sigma (z^{})d^2z^{}`$. This function is harmonic in the exterior of the domain except at infinity where it has a logarithmic singularity. Eq. (25) means that inside $`D`$, including the boundary, it is equal to $`W`$ plus a constant. Since $`W(0)=0`$, the constant is $`\phi _0(0)={\displaystyle _D}\mathrm{log}|\xi |^2\sigma (\xi )d^2\xi `$, and so $`\phi _0(z)+W(z)=\phi _0(0)`$ for $`zD`$. Moreover, according to (25),

$$_n\left(\phi _0(z)+W(z)\right)=0,zD,$$
(27)
so both tangential and normal derivatives of $`\phi _0+W`$ at the boundary vanish.

Since $`\phi _0(z)N\mathrm{}\mathrm{log}|z|^2`$ is harmonic outside, the harmonic extension of $`W`$ is

$$W^H(z)=\phi _0(0)\phi _0(z)+N\mathrm{}\left(\mathrm{log}|z|^2\left(\mathrm{log}|z|^2\right)^H\right),$$
(28)
whence it follows that

$$𝖱\left(W(z)+N\mathrm{}\mathrm{log}|z|^2\right)=0.$$
(29)
This condition means, in other words, that the domain $`D`$ is such that the function $`_zW`$ on its boundary is the boundary value of an analytic function in $`𝐂D`$. The latter, together with the normalization condition $`_D\mathrm{\Delta }Wd^2z=4\pi \mathrm{}N`$, determines the shape of the support of eigenvalues.

For the potential of the form (8) with $`V(z)=_{k1}t_kz^k`$ this condition is somewhat more explicit. In this case, the shape of the domain is determined, by the relations

$$\frac{1}{\pi k}_Dz^k𝑑z=t_k,\frac{1}{\pi }_Dd^2z=\mathrm{}N.$$

such that$`\pi \mathrm{}N`$ is the area of $`D`$ and $`\pi kt_k`$ are harmonic moments of the domain complementary to $`D`$. The problem is thus equivalent to the inverse problem of 2D potential theory.

Using the relation $`W(z)=\phi _0(0)\phi _0(z)`$, it is easy to find the value of $`\mathrm{}^2\mathrm{log}|\mathrm{\Psi }_N|^2`$ at the saddle point, which we denote by $`F_0`$:

$$F_0=_D_D\mathrm{log}\left|\frac{1}{z}\frac{1}{z^{}}\right|\sigma (z)\sigma (z^{})d^2zd^2z^{}.$$
(30)
The leading asymptote of the partition function as $`\mathrm{}0`$ is therefore $`Z_Ne^{F_0/\mathrm{}^2}`$. By inspection one can check that the average density (26) can be also obtained by variation of (30): $`\rho (z)=\delta F_0/\delta W(z)`$.

The function $`F_0`$ plays an important role in the complex analysis. It generates conformal maps from the exterior of $`D`$ onto the unit disk and gives a formal solution to the Dirichlet boundary problem. See Sec. 4 of  for details.

### Variation of the support of eigenvalues.

Let us examine the change of the support of eigenvalues $`D`$ under a small change of the potential at fixed $`N`$. In general, we can write: $`_{\delta D}f(z)d^2z=_D\delta n(z)f(z)|dz|`$ for any function $`f`$, where $`\delta D`$ stands for a strip between the domains $`D(W+\delta W)`$ and $`D(W)`$ and $`\delta n`$ has the same meaning as in the previous section (Fig. 1).

The variation of the saddle point condition is most conveniently found from (27):

$$\delta (_n(W+\phi _0))=\delta n_n^2(W+\phi _0)+_n\delta (W+\phi _0)=0.$$

where $`z`$ is on the boundary. Since both tangential and normal derivatives of $`W+\phi _0`$ vanish on the boundary, one writes $`_n^2(W+\phi _0)=\mathrm{\Delta }(W+\phi _0)=\mathrm{\Delta }W=4\pi \sigma `$. This gives the variation of the boundary:

$$\delta n(z)=\frac{1}{4\pi \sigma (z)}(𝖱\delta W).$$
(31)
Here $`𝖱\delta W=_n(\delta W(\delta W)^H)`$ is the external Neumann jump operator defined by (21). (We used the fact that $`\delta \phi _0`$ is, up to a constant, the harmonic extension of $`\delta W`$ to the exterior domain, as is seen from (28).) Note that if the variation of the potential is harmonic outside $`D`$ including infinity the domain does not change.

Let us check that our result meets the requirement that $`N`$ stays constant under a variation of the potential. Since

$$\delta n(z)\sigma (z)|dz|=\frac{1}{4\pi }_n\delta W(z)|dz|=\frac{1}{4\pi }_D\mathrm{\Delta }\delta W(z)d^2z$$

it is easy to see that

$$\delta \rho (z)d^2z=\delta n(z)\sigma (z)|dz|+_D\delta \sigma (z)d^2z=0,$$

so that $`N`$ is kept constant under this variation.

Consider special variations of the potential of the form $`\delta _\zeta W=\epsilon \mathrm{log}|z\zeta |^2`$. Then one is able to express $`𝖱\delta _\zeta W`$ through the Dirichlet Green function (cf. (19)). It is convenient to introduce the modified Green function:

$$𝒢(z,\zeta )=G(z,\zeta )G(z,\mathrm{})G(\mathrm{},\zeta ).$$
(32)
It is easy to see that

$$\delta _\zeta n(z)=\frac{1}{2\pi \sigma (z)}_n𝒢(z,\zeta ).$$

(The last term in $`𝒢`$ vanishes under the normal derivative. It is included for the symmetry and future convenience.)

### Variation of the boundary under variation of the size of the matrix.

If one varies the size of the matrix, $`NN+\delta N`$, keeping the potential fixed, the support of eigenvalues also changes its shape. To find how it grows, we note that $`\delta W=0`$ means that $`\delta \phi _0(z)=\delta \phi _0(0)`$ for all $`zD`$. Plugging here the integral representation of $`\phi _0`$, we conclude that $`\mathrm{log}|z^1\xi ^1|\delta n(\xi )\sigma (\xi )|d\xi |=0`$ for all $`z`$ in $`D`$, with the variation of the normalization condition being $`\mathrm{}\delta N=\delta n(z)\sigma (z)|dz|`$. These conditions are met with

$$\delta n(z)=\frac{\mathrm{}\delta N}{2\pi \sigma (z)}_nG(z,\mathrm{}),$$
(33)
where $`G(z,\mathrm{})=\mathrm{log}|w(z)|`$. This describes the interface dynamics in Laplacian growth models, known as Darcy’s law (cf. ).

## 4 Connected correlation functions in the first non-vanishing order in $`\mathrm{}`$.

### Connected two-point function.

In order to obtain the two-point function one has to vary the one point function $`\rho (z)=\mathrm{\Theta }(z;D)\sigma (z)`$, where $`\mathrm{\Theta }(z;D)`$ is the characteristic function of the domain $`D`$: $`\mathrm{\Theta }(z;D)=1`$ for $`zD`$ and $`\mathrm{\Theta }(z;D)=0`$ for $`zD`$. The variation reads

$$\delta \rho =\delta \sigma \mathrm{\Theta }(z;D)+\sigma \delta \mathrm{\Theta }(z;D).$$
(34)
The second term is localized on the contour. Let $`\delta (z;D)`$ be a $`\delta `$-function located on the boundary of $`D`$, defined by the condition that $`f(z)\delta (z;D)d^2z=_Df(z)|dz|`$ for any smooth function $`f`$. It is clear that $`\delta \mathrm{\Theta }(z;D)=\delta n(z)\delta (z;D)`$. Using the relation (31) between a variation of the potential and deformation of the domain, we write:

$$4\pi \delta \rho =\mathrm{\Delta }\delta W\mathrm{\Theta }(z;D)+(𝖱\delta W)\delta (z;D).$$

To keep track of the singular boundary terms, it is helpful to integrate the variation of density with some reasonable function $`f`$ on the plane. This is equivalent to calculating $`\delta \text{tr}f(M)`$ instead of $`\delta \rho `$. Setting $`\delta W=\epsilon g`$ we have:

$$\frac{4\pi }{\epsilon }\delta \rho fd^2z=_D\mathrm{\Delta }gfd^2z+_Df(𝖱g)|dz|.$$

The result is symmetric with respect to $`fg`$. With a help of the Green formula, it can be expressed through the Bergmann kernel (22):

$$\begin{array}{c}\mathrm{}^2\text{tr}f(M)\text{tr}g(M)_{\mathrm{conn}}=\frac{1}{\pi }e_D_zf(z)_{\overline{z}}g(z)d^2z+\\ \\ +\frac{1}{8\pi ^2}f(z)_n_n^{}G(z,z^{})g(z^{})|dz||dz^{}|.\end{array}$$
(35)
Choosing $`f(z)=\mathrm{log}|z_1z|^2`$, $`g(z)=\mathrm{log}|z_2z|^2`$, we find the pair correlation function of the Bose field $`\phi (z)`$:

$$\frac{1}{2\mathrm{}^2}\phi (z_1)\phi (z_2)_{\mathrm{conn}}=𝒢(z_1,z_2)\mathrm{log}\frac{|z_1z_2|}{r},$$
(36)
where $`𝒢`$ is introduced in (32) and $`\mathrm{log}r=lim_z\mathrm{}(\mathrm{log}|z|+G(z,\mathrm{}))`$ is Robin’s constant ($`r`$ is the exterior conformal radius of the domain $`D`$, see e.g. ). The function in the r.h.s. is harmonic outside the domain. If one of the points, say $`z_1`$, is located inside, one sets the corresponding Green functions in (32) to be zero. In particular, if both points are inside, the correlation function is just

$$\frac{1}{2\mathrm{}^2}\phi (z_1)\phi (z_2)_{\mathrm{conn}}=\mathrm{log}\frac{|z_1z_2|}{r}.$$
(37)
This result is valid for well separated points ($`|z_1z_2|^2>>N\mathrm{}/\sigma `$).

Taking holomorphic or antiholomorphic derivatives of (36), we find pair correlations of currents:

$$\mathrm{}^2J(z_1)J(z_2)_{\mathrm{conn}}=\frac{1}{(z_1z_2)^2}+2_{z_1}_{z_2}G(z_1,z_2),$$
(38)
$$\mathrm{}^2J(z_1)\overline{J}(z_2)_{\mathrm{conn}}=\pi \delta ^{(2)}(z_1z_2)+2_{z_1}_{\overline{z}_2}G(z_1,z_2)$$
(39)
(here it is implied that both points are outside). These results resemble the two-point functions of the Hermitian 2-matrix model found in ; they were also obtained in  in the study of thermal fluctuations of a confined 2D Coulomb gas.

Outside the domain formulas (36), (38) describe correlations also at merging points away from the boundary. In particular, the mean square fluctuation of the current is

$$\mathrm{}^2J^2(z)_{\mathrm{conn}}=\frac{1}{6}\{w;z\},z𝐂D.$$
(40)
Here

$$\{w;z\}=\frac{w^{\prime \prime \prime }(z)}{w^{}(z)}\frac{3}{2}\left(\frac{w^{\prime \prime }(z)}{w^{}(z)}\right)^2=6\underset{z^{}z}{lim}\left(2_z_z^{}G(z,z^{})\frac{1}{(zz^{})^2}\right)$$
(41)
is the Schwarzian derivative of the conformal map $`w(z)`$.

These formulas show that there are local correlations in the bulk as well as strong long range correlations at the edge. (See  for a similar result in the context of classical Coulomb systems). At the same time further variation of the pair density correlation function suggests that, starting from $`n=3`$, the connected $`n`$-point density correlations vanish in the bulk in all orders of $`\mathrm{}`$ (in fact they are exponential in $`1/\mathrm{}`$). All the leading contribution (of order $`\mathrm{}^{2n2}`$) comes from the boundary.

Note that the result (35) is universal in the sense that it depends on the potential only through the form of the domain $`D`$. The universality holds for any connected correlation function of two traces.

The kernel in the boundary term in (35) is the absolute value of the Bergman kernel  of the domain $`𝐂D`$ at the boundary. Presumably, this result can be generalized to the more complicated case of non-connected supports of eigenvalues, with the boundary term being expressed through the Bergman kernel of the Schottky double of the Riemann surface $`𝐂D`$. For the Hermitian one-matrix model, where the support of eigenvalues shrinks to a number of cuts on the real axis, a similar result was recently obtained in .

### Connected three-point function.

The 3-point function can be obtained by further varying eq. (35). Let us first transform the r.h.s., to bring it to the form convenient for the varying. Using the Green formula, one rewrites the two point function as an integral over the entire complex plane plus the integral over the exterior of the domain

$$4\pi \mathrm{}^2\text{tr}f(M)\text{tr}g(M)_{\mathrm{conn}}=_𝐂f\mathrm{\Delta }gd^2z+_{𝐂D}(gg^H)\mathrm{\Delta }fd^2z.$$

The variation of the first term is zero since it does not depend on the contour. The variation of the second term consists of two parts: the one coming from variation of the boundary and another one from variation of the integrand. The first part actually vanishes because the integrand equals zero on the boundary, where $`g=g^H`$. The variation of the integrand is $`{\displaystyle _{𝐂D}}(\delta g^H)\mathrm{\Delta }fd^2z={\displaystyle (\delta g^H)(𝖱f)|dz|}`$ (where we substituted $`fff^H`$ under the Laplace operator and again used the Green formula). Next, we use the Hadamard formula (23) to compute the $`\delta g^H`$, the response of the harmonic extension of the function $`g`$ to a small change of the domain, taken on the boundary of the initial domain (24). The result is

$$4\pi \mathrm{}^2\delta \text{tr}f(M)\text{tr}g(M)_{\mathrm{conn}}=(𝖱f)(𝖱g)\delta n|dz|.$$

Now, let us redefine $`f=f_1`$, $`g=f_2`$ and plug (31) for the $`\delta n(z)`$ with a function $`f_3`$: $`\delta n={\displaystyle \frac{\epsilon }{4\pi \sigma }}𝖱f_3`$. Using (17), we obtain the correlation function of three traces:

$$\underset{j=1}{\overset{3}{}}\text{tr}f_j(M)_{\mathrm{conn}}=\frac{\mathrm{}^4}{16\pi ^2}\frac{|dz|}{\sigma (z)}\underset{j=1}{\overset{3}{}}𝖱f_i(z).$$
(42)
The answer is non-universal, i.e., it depends explicitly on the boundary value of the Laplacian of the potential. Note that if at least one of the functions $`f_j`$ is harmonic in $`𝐂D`$, then the correlation function vanishes (in this leading order in $`\mathrm{}`$ of course).

Alternatively, one can apply the Hadamard formula directly to the 2-point correlation function of the potentials $`\phi `$ (36). We have

$$\phi (z_1)\phi (z_2)\phi (z_3)_{\mathrm{conn}}=\frac{\mathrm{}^2}{\epsilon }\delta _{z_3}\phi (z_1)\phi (z_2)_{\mathrm{conn}}=2\frac{\mathrm{}^4}{\epsilon }\delta _{z_3}\left[𝒢(z_1,z_2)+\mathrm{log}r\right]=$$

$$=\frac{\mathrm{}^4}{\epsilon \pi }_n𝒢(z_1,\xi )_n𝒢(z_2,\xi )\delta _{z_3}n(\xi )|d\xi |=$$

$$=\frac{\mathrm{}^4}{2\pi ^2}\frac{|d\xi |}{\sigma (\xi )}_n𝒢(z_1,\xi )_n𝒢(z_2,\xi )_n𝒢(z_3,\xi ).$$

The result agrees with the formula for the third order derivatives of the tau-function obtained in  within a different approach (also based on the Hadamard variational formula).

For example, in the case $`W=|z|^2`$, the support is a disk of the radius $`R=\sqrt{\mathrm{}N}`$. The conformal map is simply $`w(z)=z/R`$ and the above formula gives

$$\underset{j=1}{\overset{3}{}}\phi (\lambda _j)_{\mathrm{conn}}=2\mathrm{}R^2e(\frac{1}{(\lambda _1\overline{\lambda }_2R^2)(\lambda _1\overline{\lambda }_3R^2)}+[12],+[13]).$$

Finite $`N`$ formulas for the general $`n`$-point correlations reviewed in  should be able to be used to reclaim the same result.

### Variation of contour integrals.

One may proceed in the same way to find the connected $`n`$-point correlation functions. Starting from $`n=4`$, however, one encounters technical difficulties.

When passing from 2-point to 3-point functions we transformed integral over the boundary into a bulk integral since the latter is easier to vary. Passing to 4-point functions, one needs to vary the contour integral in (42) which does not seem to be naturally representable in a bulk form. We have to learn how to vary contour integrals. Here are general rules.

Consider the contour integral of the general form $`F(f(z),_nf(z))𝑑s`$ where $`ds=|dz|`$ is the line element along the boundary curve and $`F`$ is any fixed function. Calculating the linear response to the deformation of the contour, one should vary all items in the integral independently and add the results. There are four elements to be varied: the support of the integral $``$, the $`_n`$, the line element $`ds`$ and the function $`f`$. By variation of the $``$ we mean integration of the old function over the new contour. This gives $`\delta n_nFds.`$ The change of the slope of the normal vector results in $`\delta \frac{}{n}=_s(\delta n)\frac{}{s}.`$ The rescaling of the line element gives $`\delta ds=\kappa \delta nds`$, where $`\kappa (z)`$ is the local curvature of the boundary curve. The curvature is $`\kappa =d\theta /ds`$ where $`\theta `$ is the angle between the outward pointing normal vector to the curve and the $`x`$-axis. The formula $`\kappa (z)=_n\mathrm{log}\left|w(z)/w^{}(z)\right|`$ expresses the curvature through the conformal map. This formula is useful in calculating the 5-point function. We do not attempt to do this here. Another element is the Laplace operator on the boundary in terms of normal ($`_n`$) and tangential ($`_s`$) derivatives: $`\mathrm{\Delta }=_n^2+_s^2+\kappa _n`$.

At last, we have to vary the function $`f`$ if it explicitly depends on the contour. In particular, if this function is the harmonic extension of a contour-independent function on the plane, its variation on the boundary is given by (24).

### Connected four-point function.

Let us elaborate the case $`n=4`$. In accordance with (17), we find $`_{j=1}^4\text{tr}f_j(M)_{\mathrm{conn}}`$ from the response of the r.h.s. of eq. (42) to the variation of the potential $`\delta W=\epsilon f_4`$, so that

$$\delta \sigma =\frac{\epsilon }{4\pi }\mathrm{\Delta }f_4,\delta n=\frac{\epsilon }{4\pi \sigma }𝖱f_4.$$
(43)
To get the result, we apply the above rules to the contour integral (42). Note that in this particular case a change of the slope of the normal vector gives no contribution since the functions under the normal derivative do not change along the boundary. Along this way we obtain the following result:

$$\begin{array}{c}\frac{64\pi ^3}{\mathrm{}^6}\underset{j=1}{\overset{4}{}}\text{tr}f_j(M)_{\mathrm{conn}}=\frac{|dz|}{\sigma ^2}\underset{i=1}{\overset{4}{}}\mathrm{\Delta }f_i\underset{k=1,i}{\overset{4}{}}𝖱f_k\\ \\ \frac{|dz|}{\sigma ^2}(_n\mathrm{log}\sigma +2\kappa )\underset{k=1}{\overset{4}{}}𝖱f_k\\ \\ |dz|[\frac{𝖱f_1𝖱f_2}{\sigma }_n\left(\frac{𝖱f_3𝖱f_4}{\sigma }\right)^H+[13]+[23]].\end{array}$$
(44)
The first two terms are explicitly symmetric with respect to all permutations of $`(1,2,3,4)`$. The third one is seemingly not, but in fact it is, as is clear from the Green formula.

In the 4-point correlation function the first term in the r.h.s. of (44) vanishes. For $`\sigma (z)=1/\pi `$, we get

$$\begin{array}{c}\underset{j=1}{\overset{4}{}}\phi (\lambda _j)_{\mathrm{conn}}=\frac{\mathrm{}^6}{8\pi ^2}|dz|_n|dz^{}|_n^{}𝒢(z,z^{})\times \\ \\ \times (_n𝒢(\lambda _1,z)_n𝒢(\lambda _2,z)_n^{}𝒢(\lambda _3,z^{})_n^{}𝒢(\lambda _4,z^{})+[13]+[23])\\ \\ \frac{\mathrm{}^6}{2\pi }|dz|\kappa (z)\underset{j=1}{\overset{4}{}}_n𝒢(\lambda _j,z),\end{array}$$
(45)
where $`𝒢`$ is defined by (32) and $`\lambda _j`$’s are assumed to be outside.

## 5 Discussion

We have shown that large scale properties of the normal matrix ensemble are obtained from the analytical properties of a curve on the complex plane. The curve bounds a semiclassical support of eigenvalues and is determined by the relation $`W(z)=\mathrm{}\text{ tr}\frac{1}{zM}`$.

Large scale correlation functions of the ensemble are objects of the Dirichlet boundary problem for the non-compact exterior domain complementary to the compact domain $`D`$. The two point function is essentially the Dirichlet Green function, the higher order functions are related to the deformation of the Green function under deformations of the curve. They are determined by successive applications of the Hadamard formula and are therefore expressed through the Neumann jump on the curve and through the Bergman kernel.

We expect that other objects of matrix ensembles, such as the genus ($`1/N`$) expansion of the partition function are recorded in the analytic properties of the curve. In particular, we expect that $`F_1`$ \- a genus 1 correction to the partition function $`\mathrm{log}Z_N`$ (a correction of order $`\mathrm{}^0`$) is related to the determinant of the Laplace operator in the exterior domain $`𝐂D`$.

## Acknowledgments

We acknowledge discussions with O.Agam, J.Ambjorn, E.Bettelheim, O.Bohigas, A.Boyarsky, A.Caceres, L.Chekhov, A. Gorsky, V.Kazakov, I.Kostov, Y.Makeenko, A.Marshakov, M.Mineev-Weinstein, O.Ruchayskiy, R.Theodorescu and P.Zinn-Justin. We are indebted to P.J.Forrester who drew our attention to Refs. . A.Z. thanks B.Jancovici for useful remarks. P.W. was supported by grants NSF DMR 9971332 and MRSEC NSF DMR 9808595. P.W. thanks S.Ouvry for the hospitality in LPTMS at Universite de Paris Sud at Orsay, where the paper was completed. The work of A.Z. was supported in part by RFBR grant 01-01-00539, by grant INTAS-99-0590 and by grant 00-15-96557 for support of scientific schools.


format: text

text:
# Ownership and Trade from Evolutionary Games

## 1 Introduction

In the last three decades behavioral ecology has made strides in understanding the origin and utility of behavioral traits in animal societies (Maynard Smith 1988; Wilson 1975). Much of this success has come within the analytic framework of evolutionary game theory. Evolutionary games, which mimic the dynamics of Darwinian natural selection, differ fundamentally from classical games (Maynard Smith 1982). While the latter focus on strategic interactions between rationally calculating agents, evolutionary games model repeated interactions between adaptive but otherwise non-thinking agents.

As one tends to believe humans think, it is not surprising that classical game theory has established a dominant position in law and economics (Baird, Gertner, and Picker 1994) while evolutionary game theory has attracted much less attention.111This is not to suggest the absence of a strong evolutionary tradition in law. Scholars who have written about the evolution of law include Boyd and Richerson (1985); Cooter and Kornhauser (1980); Elliot (1985); Ellickson (1991); Epstein (1980); Hirshleifer (1982); Huang and Wu (1994); Johnston (1996); Priest (1977); Rubin (1977); and Yee (1998). A notable exception to this trend is Hirshleifer, who pointed out parallels between evolved animal behaviors (and their game theory models) and economically efficient human practices (Hirshleifer 1977, 1980, 1982). Hirshleifer proposed three metallic norms in particular: the “Golden Rule of communal sharing,” the “Silver Rule of private rights,” and the “Iron Rule of dominance.” Each, he asserted, has evolved because they have sufficient socioeconomic advantages.

This article develops an evolutionary game model of property ownership and trade. To start, it extends and interprets a well-known evolutionary model of animal territoriality as a model of human property ownership. That one could do this is perhaps not surprising. The more surprising result is that trade – the bilateral transfer of property for renumeration – then emerges as a strategy which is evolutionarily preferred over permanent ownership without trade. In other words, traders are evolutionarily superior to owners who do not trade. The identification that trade is an evolutionary stable strategy in evolutionary game theory is the main contribution of this article.

The swallowtail butterfly, Papilio zelicaon, provides a provacative example of what may underlie animal territorial behavior222 Other species, such as the hamadryas baboon, have been documented to exhibit similar behavior. (Maynard Smith 1988, p. 214). Because the swallowtail lives in low density populations, one might expect matchmaking—the finding of a sexual partner—to be a problem. This problem is solved in swallowtail society by “hilltopping.” Males establish territories at or near the tops of hills and wait for virgin females, who instinctively seek out hilltops to mate. Since there are typically more males than hilltops, most males are relegated to disadvantageous positions lower down the slopes, where they attempt to waylay females on their way uphill. Although the lower altitude males sometimes succeed, hilltop males mate most. Curiously, despite their enviable estates, hilltoppers are seldom challenged by intruders. On occasion when an intruder does confront a hilltopper, the visitor tends to retreat after a brief, mild contest.

How is hilltop occupancy negotiated? In experiments, Larry Gilbert333Because a caterpillar virus wiped out his butterfly population before Gilbert completed his studies, Gilbert’s original study, the one described in Maynard Smith (1988), was never published (L. E. Gilbert, private communication to Yee). In any case, R. Lederhouse describes similar observations of a closely related butterfly species (see Scriber, Tsubaki, and Lederhouse 1995). tested two alternative explanations: (A) Potential intruders are intimidated from invading because hilltoppers are physically strongest, and because this fact is perceived by the intruders; or (B) Swallowtails defer to prior possessors, that is, whoever first stakes out a hilltop is granted a socially-established privilege to keep it (not unlike in the Anglo-American doctrine of adverse possession).

In a series of experiments with pairs of randomly selected male butterflies, Gilbert ruled out A, and found support for the prior possessor theory, B, as follows. He convinced each butterfly of a pair that it was the sole occupier of the same hilltop by letting it be the hilltopper on alternate days. On its days off, Gilbert kept the sidelined butterfly unaware of its counterpart’s existence by confining it in dark room. After a couple of weeks, each male clearly acted as the hilltop’s rightful proprietor, chasing away all comers, who invariably retreated without much protest. When Gilbert finally released both males to the same hilltop on the same day, an abnormally prolonged contest between the two “proprietors” ensued, lasting many minutes and causing serious injury to each contestant. As a result, Gilbert concluded that deference to prior hilltop possessors is an instinctive trait of swallowtails.

Maynard Smith (1982) and others hypothesized that such instincts evolved by natural selection of the fittest and constructed evolutionary game theories modeling the evolutionary processes. In these models, deference to possession—a dispute resolution strategy based on pre-existing status—is evolutionarily preferred over an always-fight strategy, which costs too much, and a never-fight strategy, which yields too little.

Ownership conventions in human societies range from the simplest unspoken norms, such as not cutting ahead of somebody else in a grocery store queue to much more involved rights bundles expressed in the Common Law. A broad range of viewpoints (Rose 1985), ranging from Locke’s labor-mixing theory of property444 In a nutshell, Locke’s view is that one owns one’s body and, by extension, the fruits of his body’s labor. Hence, first useful possession establishes a property right (Epstein 1979). to law and economics555The normative law and economics view may be summarized as follows. Property rights may be thought of as a bundle of strands of primitive rights. Particular strands—tailored to the situation—are granted to encourage social-welfare-optimizing use and investment. A benefit-cost balance determines whether a particular strand should be granted enforcement. Costs include enforecement costs (Posner 1992). may be called upon to justify them. Legal property rights can be enforced by various combinations of liability and injunctive remedies (e.g. Calabresi and Melamed 1972).

Whatever the justifications, most ownership rights bundles consist of two primitive strands: (a) Possession, the right to occupy or possess what one owns, and (b) Trade, the right to buy and sell ownership. By constructing evolutionary game theory models, I will illustrate how evolutionary forces can serve to establish these two strands as stable strategies. The trade model shows that those who trade are evolutionarily preferred over Possessors who don’t. What is new here is that evolutionary forces are enough by themselves to establish Possession and, given Possession, the practice of trade.

Section 2 reviews Maynard Smith’s (1982) construction of evolutionary stable strategy (ESS) and argues that ESSes have interpretation as social norms. Sections 3 and 4 present two models corresponding, respectively, to the two strands of ownership. In the Section 3 model, possession is an ESS; in Section 4, possession with the right to trade is an ESS.

## 2 Evolutionary Stable Strategies

Social benefits come only at a price: for every benefit accrued, there must be set of behavioral constraints or obligations to be fulfilled. Social norms—prevalent responses to recurring social situations—are the reciprocal constraints enabling the social benefits.

Evolutionary game theory provides a quantitative dynamical theory of how such social constraints on behavior can emerge from anarchy.666Axelrod (1986) argues that social norms cannot survive unless there is a secondary norm enforcing the first norm. This implies the need for a tertiary norm to enforce the secondary norm, and so forth ad infinitem. (See also Martinez Coll and Hirshleifer 1991; Lomborg 1996.) In contrast, the elegance of evolutionary stable strategies (as exemplified by the models herein) is that, within the context of the game, they do not require such an infinite hierarchy of supporting norms. The basic ingredients of evolutionary game models for our purposes are:

* a self-contained community of disputants interacting in repeated, random pairwise encounters;
* who, at each encounter, select from a predetermined777In more ambitious formulations the strategy space is permitted to evolve via mutations in analogy to genes in biology. menu of strategies, say, $`\{\alpha ,\beta ,\gamma ,\mathrm{}\}`$ which may be either pure or mixed.
* If player $`\mathrm{\#}1`$ chooses strategy $`\alpha `$ and $`\mathrm{\#}2`$ chooses $`\beta `$, the payoff to player $`\mathrm{\#}1`$ is denoted $`w_{\alpha \beta }`$, which is determined entirely by $`\alpha `$ and $`\beta `$.
* A round consists of many random encounters for each player. After each round, the community undergoes “natural selection,” in which strategies replicate in proportion to how far above average their scores in the just-completed round were; those with below average scores die off in proportion to how far below average they were.

Consequently, an adaptive population evolves according to a set of coupled first-order differential equations—one corresponding to each strategy—in the same spirit as the Malthusian predator-prey relations. Solving the equations yields a phase diagram in strategy space, which typically has several fixed points (Friedman 1991). Some of these fixed points are attractive: populations always evolve into them so that they are evolutionarily stable against subversive mutations. Each attractive fixed point in a phase diagram represents an “evolutionary stable strategy (ESS).”

In other words, an ESS is a possible mode of behavior that might lock-in. Thus, an ESS can be interpreted as a social norm.

For purposes of this article, it will be sufficient to state a criteria defining ESS(es) given the game’s defining payoff matrix.888From Maynard Smith (1982). Other authors (e.g. Zeeman 1981) have pointed out that these criteria expressed in terms of payoff matrix elements do not always recover all the aforementioned attractive fixed points. For our models, the technical distinction is immaterial. To this end, for a given game strategy $`\alpha _{}`$ is an ESS if, starting from a status quo where $`\alpha _{}`$ is the norm, it is not possible for insurgents to achieve higher payoffs with a renegade strategy, say, $`\gamma `$. Algebraically, this means $`\alpha _{}`$ is an ESS if either

$$w_{\alpha _{}\alpha _{}}>w_{\gamma \alpha _{}}\gamma $$
(1)
or

$$w_{\alpha _{}\alpha _{}}=w_{\gamma \alpha _{}}\mathrm{and}w_{\alpha _{}\gamma }>w_{\gamma \gamma }\gamma .$$
(2)
Condition (2) is sufficient to establish an ESS because, even if $`w_{\gamma \alpha _{}}=w_{\alpha _{}\alpha _{}}`$, invaders behaving according to $`\gamma `$ cannot successfully gain a foothold if they perform so poorly against each other that they prevent themselves from becoming a sizable fraction of the population. Note that whether a strategy is an ESS or not depends on the strategy space, the set of competing strategies. What is an ESS in one strategy space may cease to be an ESS if the strategy space is expanded to include other strategies. Also, games can (and usually do) have more than one ESS.

Of direct relevance to us is the Hawk-Dove game (Maynard Smith 1973), whose payoffs are depicted in Figure 1. In this game, two equally matched parties vie for the same asset, worth $`\$V`$ to each. (Imagine two randomly paired strangers fighting over a parking space at a crowded shopping mall.) Suppose only two strategies, Hawk (H) and Dove (D), are available. Hawks always fight and, since both parties are equal, in a fight each Hawk has only a one-half chance of winning the asset. Fighting has an expected total cost to each participant of $`\$h`$. In addition to the expected injury, $`h`$ contains all other costs including the expected energy expenditure and any risk-bearing costs. Doves retreat when confronted by a Hawk. If two Doves meet, a random one of the two Doves retreats and leaves the other to the spoils.

In neoclassical game theory (e.g. Baird, Gertner, and Picker 1994), $`\frac{1}{2}V>h`$ corresponds to the Prisoner’s Dilemma while $`\frac{1}{2}V<h`$ corresponds to the Chicken game. In the Prisoner’s Dilemma, $`H`$ is called “defection” and defection is the unique Nash equilibrium strategy for both players. In the Chicken game, the unique Nash equilibrium is for one player to be the chicken (D) and the other to be the hawk (H).

When $`\frac{1}{2}V>h`$ (the Prisoner’s Dilemma case), H is an ESS and D is not by virtue of Criterion (1). In other words, aggressiveness is evolutionarily preferred when the rewards outweigh the costs of fighting. Note that this does not mean H optimizes social welfare. In fact, an all-Dove population maximizes social welfare.

When the possibility of serious injury is sufficiently large, $`\frac{1}{2}V<h`$ (the Chicken game case). In this case, neither pure strategies $`H`$ or $`D`$ are ESSes. It turns out that the only ESS (within the H-D strategy subspace) is a mixed strategy $`\alpha _{}=pH+(1p)D`$, where $`p\frac{1}{2}V/h`$. This mixed strategy can be achieved in two ways. Either the population is homogeneous and at an encounter each individual exercises H with probability $`p`$, or the population is comprised of fraction $`p`$ Hawks and $`(1p)`$ Doves. In either case, the more potential costs exceed potential gains the less evolutionarily attractive hawkishness is.

The material in this section was established by Maynard Smith (1982) and others. While the Hawk-Dove game has provided insights into animal behavior, it is too simple to allow for more sophisticated human strategies. Humans act like neither Hawks, Doves, nor mixtures thereof. Rather, we have more sophisticated options. The next section turns to one of them.

## 3 Possession as an ESS

“Finder’s keepers” and “first come, first serve” are not only basic thumb rules in playground citizenship, they are powerful norms that have been recognized by the courts and applied widely in such varied settings as adverse possession, abandoned property, fisheries, wildlife, seabed minerals, groundwater rights, intellectual property, debt collection, oil and gas, pollution permits, the radio frequency spectrum, satellite orbits, and ownership of wartime spoils.

To model this norm as an ESS, introduce the “Possessor (P)” strategy:

$$P\{\begin{array}{cc}\mathrm{H}\hfill & \text{if current owner;}\hfill \\ \mathrm{D}\hfill & \text{if current intruder.}\hfill \end{array}$$

The Possessor strategy models the practice of “possession”. Unlike Hawks and Doves, Possessors observe convention based on their status; their behavior depends on whether they are the owner or intruder.

Figure 2 depicts the Hawk-Dove-Possessor game payoffs. If $`f=\frac{1}{2}`$, the Hawk-Dove-Possessor game reduces to the Hawk-Dove-Bourgeois game (Maynard Smith 1988; Hirshleifer 1982, pp. 22-23). What distinguishes the Hawk-Dove-Possessor game from the Bourgeois game is the parameter $`f`$. $`f`$ is the expected fraction of confrontations in which a disputant anticipates she will be in the role of Owner. We will show that P is the unique ESS for all $`f`$ in the interval $`(0,1)`$. This implies that the Possessorship strategy is robust to the wealth (or dearth) of ownership opportunities available – all one requires is a nonzero chance ($`f0`$) to be an owner.

The payoffs of the Hawk-Dove-Possessor game are motivated as follows. First, the parameter $`f`$ has the following interpretation. Imagine a neighborhood in which everyone possesses a separate (but otherwise identical) plot of land. At the same time, everyone is wandering the neighborhood interested in obtaining possession of additional plots. Suppose that each agent wanders randomly throughout the neighborhood so that, in the repeated disputes which arise, each disputant is in the role of owner fraction $`f`$ the time and intruder the other $`1f`$. For instance, if there were $`N`$ disputants each owning one plot of land, in a round-robin tournament every disputant would intrude $`N1`$ times and be reciprocally intruded upon $`N1`$ times. In this case, $`f=\frac{1}{2}`$ because every disputant is an owner half the time and an intruder the other half.

Also of interest is when land is scarce and not everyone can own a plot. Then the probability $`f`$ of being an owner in a random encounter is no longer $`\frac{1}{2}`$, but will be some positive number less than that. To obtain an expression for $`f`$, let $`n`$ denote the number of available plots, and $`N`$ the total number of disputants. Assume $`Nn`$ and that $`N`$ is very large. Assume everyone has an equal chance to be an owner and nobody owns more than one plot of land at any time. Then the chance to be an owner is $`g=\frac{n}{N}`$. The total number of encounters a plot-owner endures in a round-robin tournament where every disputant (owners and non-owners alike) attempts to intrude once on every alien plot is $`(N1)+(n1)`$ because the owner defends his land once against the other $`N1`$ disputants and intrudes once on the other $`n1`$ plot owners. Hence, ignoring the higher order effect of what happens when property transfers resulting from encounters lead to a disputant temporarily owning more than a single plot of land, the fraction of encounters where an owner plays the role of owner is

$$\frac{N1}{N1+n1}.$$

The overall fraction $`f`$ of time a random disputant plays the role of owner in an encounter is

$`f`$ $`=`$ $`\left(\mathrm{chance}\mathrm{to}\mathrm{be}\mathrm{an}\mathrm{owner}\right)\times \left({\displaystyle \genfrac{}{}{0pt}{}{\mathrm{fraction}\mathrm{of}\mathrm{encounters}}{\mathrm{in}\mathrm{owner}^{}\mathrm{s}\mathrm{role}\mathrm{if}\mathrm{owner}}}\right)`$
$`=`$ $`g\times \left({\displaystyle \frac{N1}{N1+n1}}\right)`$
$`=`$ $`g\times \left({\displaystyle \frac{1\frac{1}{N}}{1+\frac{n}{N}\frac{2}{N}}}\right).`$
In the limit where the number $`N`$ of disputants is large and the ratio $`g=n/N`$ is finite ($`n`$ may or may not be large depending on the value of $`g`$, but $`n`$ is strictly less than $`N`$),

$$\underset{\genfrac{}{}{0pt}{}{N\mathrm{}}{g=n/N}}{lim}f\frac{g}{1+g}.$$

In this limit, $`0f\frac{1}{2}`$ because $`0<g<1`$. In particular, $`f=0`$ when the chance of land ownership is zero ($`g=0`$). When land is scarce ($`n<<N`$), the chance of land ownership is small ($`g<<1`$) and so is $`fg`$. On the other hand, when the chance of land ownership is almost certain999When land is plentiful relative to the number of disputants ($`n>>N`$), the chance $`g`$ for a disputant to be an owner is given by a complex combinatoric function if disputants are allowed to own more than one plot of land; the expression $`g=n/N`$ is valid only when ownership is capped at one plot per person. ($`g1`$), then $`f=\frac{1}{2}`$.

In the Hawk-Dove-Possessor game, at each encounter a disputant can either behave like a Hawk, a Dove, or a Possessor. Figure 2 depicts the payoffs of the Hawk-Dove-Possessor game. In deriving Figure 2, payoffs have been averaged over many encounters. Since disputants are owners $`f`$ of the time, and intruders $`1f`$ of the time, the payoff to disputant $`\mathrm{\#}1`$ is the weighted average of two conditional payoffs:

$$w=fw_{|\mathrm{owner}}+(1f)w_{|\mathrm{intruder}}.$$

For example,

$$w_{PH}=f(\frac{1}{2}Vh)+(1f)\times 0=(\frac{1}{2}Vh)f,$$

$$w_{HP}=fV+(1f)(\frac{1}{2}Vh)=\frac{1}{2}Vh+(\frac{1}{2}V+h)f,$$

and

$$w_{PP}=fV+(1f)\times 0=fV.$$

When $`\frac{1}{2}V<h`$ (the Chicken game case), P is the only ESS of the Hawk-Dove-Possessor game. (The reader is invited to confirm this using Conditions (1) and (2).) In fact, one can go further and prove there is no mixed ESS so that P is the unique ESS if the strategy space is H-D-P.) As this is true for any value of $`f(0,1)`$, it takes only a tiny (any nonzero) chance $`f`$ of ownership to establish P as an ESS. The fact that P is not an ESS exactly at $`f=0`$ is moot because $`f=0`$ means there is no property to be owned ($`n=0`$).

In defining the Hawk-Dove-Possessor payoff matrix, I assumed that owners and intruders value the disputed plots equally, that it costs owners and intruders the same to engage, and that owners and intruders each have equal chances of winning a fight. In real life, due to informational advantages of being an owner, owners probably value and can defend their properties more and better than intruders. So it is likely that possessorship is even more preferred than in this stylized model. What the Hawk-Dove-Possessor model shows is that possessorship is evolutionarily stable despite ignoring all the likely advantages an owner has over an alien intruder.

While P is an ESS for the Chicken game, it does not resolve the Prisoner’s Dilemma. When $`\frac{1}{2}V>h`$ (the Prisoner’s Dilemma case), $`H`$ remains the only pure-strategy ESS. When the expected gain from fighting exceeds expected losses, aggressiveness is evolutionarily preferred over dovishness and possessorship.

Possessor is an ESS for the Chicken game because it provides a predictable and costless mechanism for resolving disputes. As in the game, possession in the real world serves to establish an asymmetric dispute resolution mechanism in an otherwise symmetric situation. This role is consistent with the prolific public notice requirements usually associated with possession, including title registration requirements and the adverse possession doctrine, under which owners can be divested of property for failing to protest against dispossession in a timely manner.

The utilitarian value of possession was recognized early on in the common law. In Pierson v. Post101010 Pierson v. Post. 3 Cai. R. 175 (N.Y. Sup. Ct. 1805) a famous wild-fox case from the nineteenth century, a hunter, Post, had a fox in his gunsight but before he could fire an interloper killed the fox and ran off with the prize. The indignant Post sued the interloper on the theory that his pursuit of the fox established his right to have it. The court disagreed. It held that possession requires a clear act putting the world on notice that the “pursuer has an unequivocal intention of appropriating the animal to his individual use.” (Id. at 178.) Gaining property rights over a wild animal requires either establishing physical control over it or mortally wounding it.

### 3.1 Why not Anti-Possessor?

I asserted at the beginning of Section 3 that the virtue of Possessorship is dispute resolution via an objective symmetry-breaking criterion. Yet, Possessorship is not the only symmetry-breaking mechanism one can imagine. In particular, consider its mirror image, the “Anti-Possessor (AP) strategy:”

$$AP\{\begin{array}{cc}\mathrm{D}\hfill & \text{if current owner;}\hfill \\ \mathrm{H}\hfill & \text{if current intruder.}\hfill \end{array}$$

One might reasonably guess that Anti-Possessorship would serve just as well as Possessorship as a symmetry-breaking device. AP indeed is an ESS in a game consisting of Hawks, Doves, and Anti-possessors.

In animal societies, AP is rare, but not totally unheard of. The social spider Oecibus civitas lives together in groups, but each constructs its individual web. If a spider is driven from its web, it may dart into the web of a neighbor. If the neighbor is in residence, it does not expel the intruder, but instead darts out into somebody else’s web. As a result, dislocating a single spider may trigger a chain reaction—a game of musical webs (Maynard Smith 1982, p. 96).

It is apparent why P is more common than AP in both nature and culture. AP is disfavored due to nomadicy costs. In an AP culture, citizens are forced to alternate as owners and intruders in an unending tag-team match. Due to relocation costs, AP is not as viable as P.

## 4 Trade as an ESS

While possession retention is the first strand of property rights, the trading is the second strand. Trade is efficient because if an intruder values a property at $`\$V`$ while the current owner values it at $`\$v<\$V`$, then both benefit if the owner sells it to the intruder for $`\$x`$ where

$$v<x<V.$$

Beyond restricting $`x`$ to the interval $`(v,V)`$, the model here does not have anything to say about whether the transaction price $`\$x`$ is paid in money or another asset of value $`\$x`$. The model also does not specify how the value of $`x`$ is determined – whether by case-by-case negotiation or by an exogenously given recipe. In this section, I will simply show that trading is an ESS for any $`x(v,V)`$. The point is that trade is evolutionarily preferred over stoic possessorship (P).

Suppose in every encounter, the two disputants do not value the disputed property equally because, for example, they always have slightly different reasons for wanting the property. Let $`\$V`$ and $`\$v`$, with $`v<V`$, denote the two mismatched valuations. Define a “Trader (T)” as a Possessor who is willing to sell or buy for $`\$x`$ when dealing with a fellow Trader. In particular, when both owner and intruder of a particular encounter are Traders, and the intruder values the property more (e.g., at $`\$V`$) than the owner, he will want to purchase the property from the owner for $`\$x`$, where $`v<x<V`$. Symbiotically, a Trader-owner who values the property less (e.g., at $`\$v`$) than a Trader-intruder will readily agree to sell the property to the latter for $`\$x`$. In other words,

$$T\{\begin{array}{cc}\mathrm{P}\hfill & \text{if counterpart is not T;}\hfill \\ \{\begin{array}{cc}\mathrm{sell}\mathrm{for}x\hfill & \text{if }v\text{-valuing owner;}\hfill \\ \mathrm{buy}\mathrm{for}x\hfill & \text{if }V\text{-valuing intruder;}\hfill \\ \mathrm{P}\hfill & \text{otherwise;}\hfill \end{array}\hfill & \text{if counterpart is T.}\hfill \end{array}$$

Figure 3 depicts the payoffs of the Hawk-Dove-Possessor-Trader game (“Trader game” for short). In deriving the payoff matrix, I have assumed that disputants are owners fraction $`f`$ of the time, and intruders $`1f`$ of the time. Additionally, each encounter is between a disputant valuing the property at $`\$V`$ and one valuing it at $`\$v`$. Thus, the payoff to disputant $`\mathrm{\#}1`$ is the weighted average of four conditional payoffs:

$$w=\frac{f}{2}(w_{|\genfrac{}{}{0pt}{}{V\mathrm{valuing}}{\mathrm{owner}}}+w_{|\genfrac{}{}{0pt}{}{v\mathrm{valuing}}{\mathrm{owner}}})+\frac{1f}{2}(w_{|\genfrac{}{}{0pt}{}{V\mathrm{valuing}}{\mathrm{intruder}}}+w_{|\genfrac{}{}{0pt}{}{v\mathrm{valuing}}{\mathrm{intruder}}}).$$

For instance,

$$w_{TT}=\frac{f}{2}(V+x)+\frac{1f}{2}(Vx+0)=\frac{1}{2}V(\frac{1}{2}f)x.$$

The ESSes of the Trader game depends the relative values of the cost of fighting $`h`$, the average spoils of sharing without fighting $`\overline{V}\frac{V+v}{4}`$, and the incremental expected per party gain from trading $`\frac{Vv}{2}`$. As depicted in Figure 4, when the cost of fighting is small (Region I), H is an ESS and T is not. On the other hand, when the cost of fighting is comparably large (Region III), T is an ESS and H is not.

When fighting cost is intermediate (Region II), both H and T are ESSes. If $`V=v`$, then $`\mathrm{\Delta }=0`$ and Region II is subsumed by Region I. Hence, Region II owes its existence to a valuation discrepancy between owner and intruder. In Region II, hawkishness is viable because the cost of fighting is less than the expected spoils $`\overline{V}`$. At the same time, the cost of fighting exceeds the expected harm from not Trading, so Trading is also viable. As it is, if everyone is trading in Region II, a disputant is better off to trade. On the other hand, if everyone is a Hawk in Region II, a disputant is advised to also be a Hawk.

These results follow from verifying the necessary and sufficient ESS Criteria (1) and (2). The intuition for why Region III ($`h>\overline{V}`$) is a trade-only region is as follows. In this case, within the H-D-P strategy subspace P is the ESS (for the same reasons why P is the ESS in the Hawk-Dove-Possessor game). T fares no worst than P against Hawks and Doves since Traders and Possessors act the same way against Hawks and Doves. The remaining issue is why T fares better than P within the P-T subspace. The reason is because, in trades the seller-owner always nets a surplus of $`\$(xv)`$ while the buyer-intruder nets a surplus of $`\$(Vx)`$. Because both parties gain and nobody loses in a trade, Traders fare better than Possessors and T is the only pure-strategy ESS.

Clearly, Region III corresponds to bartering in ancient societies and commerce in modern ones. When the cost of fighting is high, trade evolves.

Likewise, Region I is reminiscent of nations fighting over land, oil fields, and governance structures. When the expected value of the asset far exceeds the cost of war, hawkishness emerges.

Region II, when fighting cost is intermediate, is the most interesting. In this case, either H or T can emerge and, once one of them does, it locks in as the ESS. However, either H or T has an equal ex ante chance of emerging. For instance, Region II may help explain some litigation practices in the United States, but whether it does or not is difficult to assess because determining the parties’ costs and benefits in these situations is difficult. Suppose costs are such that Region II applies in some legal disputes, such as bankruptcy or tort cases. Then bankruptcy cases may be more likely settled out of court (a T strategy) while tort litigation may go to trial more frequently (an H strategy) simply because historically a T strategy has emerged as normal practice in bankruptcy while an H strategy is the norm in tort litigation. Once an ESS is established, the strategy (T or H) becomes the modus operandi when that kind of dispute arises even though the other strategy could fare equally well if everyone were to adopt it instead.

Finally, Figure 4 highlights the cost $`h`$, which is exogenously given, as the critical determinant of whether T or H or both are evolutionarily favored. As law has the power to change $`h`$ by exacting extra penalties and fees on disputants, an implication of the Trading game is that law can inspire or hinder the evolution of trade by its assessment of penalties on hawkishness or trading.

## 5 Concluding Remarks

The Hawk-Dove-Possessor game and the Hawk-Dove-Possessor-Trade game provide two messages. The first message is that deference by intruders to owners is evolutionarily preferred over non-status-based behavior. This is because prior possession provides a ready asymmetrizing criterion enabling resolution of possession conflicts. Individuals who avoid transactions costs by resolving conflicts based on a cultural asymmetrizing criterion are better off than those who do not.

Can one do better than mere possession? The second message answers this question affirmatively. A Trade strategy trumps Possessorship with no trade. Trade is the ability to buy and sell according to what optimizes personal gain; trading does not occur unless both parties gain. Accordingly, traders always benefit from trade and, so, are evolutionarily preferred. Analysis of the Hawk-Dove-Possessor-Trade game shows that those who trade are evolutionarily preferred over Possessors who don’t. That voluntary exchange tends to improve social welfare is not new; it is already well known from traditional Walrasian analysis. What is new here is that, without the help of rationality or utility maximization, evolution is enough by itself to generate trade.

This raises the question, If trade is evolutionarily preferred over possession without trade, why don’t Gilbert’s butterflies in Section 1 trade?111111To be more precise, there is no evidence to suggest that butterflies trade. However, we cannot rule out the possibility that they do trade. Other species exhibit a rudimentary form of trade biologists call “mutualism.” Cleaning fish, which eat debris off the bodies of bigger fish, trade their cleaning services for food (Chapter 13, Maynard Smith 1982). Why doesn’t an animal with extra food barter away some in exchange for a future meal? To be sure, this does occur in a rudimentary way when animal family groups hunt together in cooperative packs and share prey. But it is probably fair to say that animals do not trade with strangers unless the exchange is immediately mutually gratifying (as in the cleaning fish example).

The main barrier to animal trade is probably logistical and traces back to the exchange value $`x`$ in the Hawk-Dove-Possessor-Trade game. The existence of $`x`$ tacitly presumes a form of book keeping, either mental accounting or the use of money. Money has not evolved in the animal kingdom and, absent money, anonymous trading is difficult because unrelated animals cannot keep track of who owes who how much. Hence, the failure of animals to evolve a monetary system probably has obstructed the emergence of animal trade from animal territoriality instincts. But why has money not evolved in animal kingdoms? Perhaps, like law or other sociopolitical structures, money is higher up in the evolutionary path, and it emerges only along with human (or prehistoric human) traits such as language and agriculture.

Beyond property rights, can evolutionary games explain the origin of other human social norms? An affirmative answer requires two ingredients: a theory of the origin of social norms, and an identification of the specific evolutionary games corresponding to each norm. This article has illustrated how to realize the second of these two ingredients by constructing the games corresponding to the practice of Possession and Trade. Let’s now turn to the first.

For an evolutionary theory of social norms, it is unlikely that human social norms are biologically evolved like butterfly territoriality is. Rather, the allusion to biological evolution must be metaphorical. In biological evolution, genes undergo natural selection in an ecology of water, geology, food, and weather. To survive, genes and their organic manifestations must adapt, mutate, and establish parasitic, predatory, or mutualistic niches in the food chain. Natural selection of the fittest weeds out inappropriate genes. In the long run only genes which manage to stake out biological niches in the food chain reproduce and propagate.

Norms are the “genes” undergoing natural selection. Through their embodiment in the behaviors of their human carriers, norms compete against each other for social popularity against a backdrop of cultural traditions and legal and political institutions – themselves manifestations of successful ideas. Not unlike between genes in biological evolution, the competition between norms for social acceptance and influence is a life and death struggle. Those that don’t successfully adapt across time and sociopolitical barriers die.

This notion that norms are like genes is not new. In The Selfish Gene, Dawkins (1976) proposes that social ideas, what he calls “memes,” are a nonorganic form of life. His examples of memes include tunes, catch-phrases, taboos, and architectural fashions. In Dawkins’ view, the fundamental characteristics of life are replication and evolution. In biological life, genes serve as the fundamental replicators. In human culture, memes are the fundamental replicators. Both genes and memes evolve by mutation-coated replication and natural selection of the fittest. Analogous to how genes encode the essence of biological life, Dawkins regards memes as genetic carriers of a memotic life form—what we know as human culture. This theme has found popularity in several fields (Waldrop 1992; Wilson 1975; Yee 1997).

Epstein’s theory (1980) of evolutionary norms goes beyond Dawkins’ by incorporating a sociobiological dimension. The thrust of Epstein’s hypothesis is that human beings who abide by certain rules of conduct are more likely to survive, reproduce, and pass on both biological genes and ideological memes to their children. Over generations, natural selection operating at a combined biological and cultural level leads to ingrained behavioral traits. Epstein listed four categories of law which he thinks has evolutionary roots: (i) prohibition against violence except in self-defense; (ii) first possession as the root of title; (iii) obligations of parents to their offspring; (iv) promissory obligations.

A cross between Dawkins’ memes and Epstein’s theory of sociobiological enhancement is compelling. In this picture, less successful individuals and groups within a population must imitate the behavior of their more successful peers in order to successfully compete for resources. Accordingly, the more above average an individual is, the more others copy his behavior. As a result of peer mimicry, the population establishes and self-enforces over time standards of normal behavior. Normal behavior may either be time-independent or it may cycle through a range of behaviors. This picture is naturally reconciled with evolutionary games because the evolutionary process is essentially a scenario of replication dynamics based on survival of the fittest. Any process which favors iterated, merit-based growth of some subgroups at the expense of peers—such as Darwinian evolution or proportional group learning—can be described by evolutionary games.

The appeal of evolutionary games is that participants do not have to be endowed with superhuman characteristics like unfailing Bayesian rationality. Even butterflies and baboons are qualified to play. All that is asked is that the parties learn by trial and error, incorporate what they learn in future behavior, and die if they don’t.

This article also raises another question, What is the connection between norms and formal law. The Possessor and Trade ESSes suggest only that these ownership conventions may emerge as social norms – they do not draw any connection between norms and formal law. Indeed, Ellickson (1991) suggests that not only are fundamental norms like neighborly cooperation pervasive but they exist independently of and oblivious to legal standards. Based on a case study of ranchers in Shasta County, California, Ellickson conjectures that the omnipresent threat of reputation-damaging gossip is an effective means of enforcing social norms in the Shasta County ranching community without the need of laws. So why laws?

If law is not necessary to sustain Trade in Regions II or III of Figure 4, is the (only) role of law to guide society into (or out of) these two Regions by twiddling cost $`h`$ with the exaction of penalties?


text:
# A Primordial Origin of the Laplace Relation Among the Galilean Satellites

> Understanding the origin of the orbital resonances of the Galilean satellites of Jupiter will constrain the longevity of the extensive volcanism on Io, may explain a liquid ocean on Europa, and may guide studies of the dissipative properties of stars and Jupiter-like planets. The differential migration of the newly formed Galilean satellites due to interactions with a circumjovian disk can lead to the primordial formation of the Laplace relation $`n_\mathrm{𝟏}\mathbf{}\mathrm{𝟑}n_\mathrm{𝟐}\mathbf{+}\mathrm{𝟐}n_\mathrm{𝟑}\mathbf{=}\mathrm{𝟎}`$, where the $`n_i`$ are the mean orbital angular velocities of Io, Europa, and Ganymede, respectively. This contrasts with the formation of the resonances by differential expansion of the orbits from tidal torques from Jupiter.

The orbital resonances among the Galilean satellites of Jupiter have led to sustained dissipation of tidal energy to produce astounding volcanos on Io and, probably, to maintain a liquid ocean on Europa. Understanding the origin of these resonances will constrain the history of the satellites and their formation scenarios. The two models proposed are an assembly of the resonances through differential tidal expansion of the orbits from tides raised on Jupiter (?, ?) and a primordial origin of unspecified assembly and subsequent evolution away from exact resonance (?, ?). Here we discuss a means of assembling the resonances during the formation process in support of a primordial origin.

The currently observed orbital resonances at the 2:1 mean motion commensurabilities involving Io-Europa and Europa-Ganymede are such that the resonance variables $`\theta _1=\lambda _12\lambda _2+\varpi _1`$ and $`\theta _3=\lambda _22\lambda _3+\varpi _2`$ librate about $`0^{}`$ and $`\theta _2=\lambda _12\lambda _2+\varpi _2`$ librates about $`180^{}`$, all with small amplitude. Here the $`\lambda _i`$ are mean orbital longitudes and the $`\varpi _i`$ are the longitudes of periapse numbered consecutively for Io, Europa, and Ganymede, respectively. The simultaneous libration of $`\theta _2`$ and $`\theta _3`$ means the Laplace angle $`\theta _5=\theta _2\theta _3=\lambda _13\lambda _2+2\lambda _3180^{}`$ with small libration amplitude, leading to $`n_13n_2+2n_3=0`$, where $`n_i=\dot{\lambda }_i`$ are the mean angular velocities with the dot indicating time differentiation. This last libration condition is called the Laplace relation because Laplace first understood its stability. The geometry is such that conjunctions of Io and Europa occur when Io is near the periapse of its orbit and Europa is near its apoapse, and conjunctions of Europa and Ganymede occur when Europa is near its periapse. Ganymede can be anywhere in its orbit at the conjunctions with Europa, because the resonance variable $`\theta _4=\lambda _22\lambda _3+\varpi _3`$ circulates through all angles. The existence of this group of resonances with small amplitudes of libration implies an assembly of originally random orbits by dissipative processes.

A currently popular means of resonance assembly is by differential orbital expansion due to gravitational tides raised on Jupiter, principally by Io, and damping of the libration amplitudes by tidal dissipation within the satellites, mostly in Io (?, ?). Io’s orbit expands more rapidly than Europa’s, and capture of Io and Europa into the 2:1 orbital resonances is certain for a reasonably wide range of initial orbital eccentricities. Continued orbital expansion tends to drive the pair deeper into the resonance (i.e., driving $`n_12n_2`$ to smaller values) and thereby to increase the orbital eccentricities. But dissipation of tidal energy within the satellites tends to decrease the eccentricities, so an equilibrium can be reached where the eccentricities and the ratio of the mean motions are almost constant as the pair of satellites moves out together locked in the resonances (?, ?). Europa eventually encounters the 2:1 mean motion commensurability with Ganymede, where capture into the Laplace relation as the resonance variable $`\theta _3`$ is trapped into libration about $`0^{}`$ has a probability of about 0.9 (?). $`\theta _4`$ is circulating because the free eccentricity (?) of Ganymede’s orbit exceeds the relatively small forced eccentricity in the resonance (?). The librations are damped by the tidal dissipation in Io. Variations of this scenario have been proposed by Malhotra and colleagues (?, ?, ?). The creation of the Laplace relation by differential tidal expansion of the orbits depends on the dissipation function $`Q_J`$ of Jupiter being sufficiently small ($`10^6`$) for non-negligible expansion of the orbits to occur, and the maintenance of the current configuration in equilibrium would require an even smaller $`Q_J`$ ($`2\times 10^4`$).

Greenberg, motivated by the difficulty of accounting for a sufficiently small $`Q_J`$ necessary for the tidal origin or equilibrium described above and the observed positive $`\dot{n}_1`$, proposed a primordial origin of the Laplace relation (?, ?). He described evolution from much deeper in the resonances (i.e., with $`n_12n_2`$ approaching zero and thereby forcing much larger orbital eccentricities). The path from deep resonance as dissipation decreased the eccentricities is characterized by stable librations of the resonance angles about centers much different from the $`0^{}`$ or $`180^{}`$ currently observed. The high dissipation rates at the high eccentricities imply rapid decay of the resonances ($`<10^7`$ year time scale), and the only way to preserve them for $`4.6\times 10^9`$ years was to have a sufficiently low $`Q_J10^6`$ to allow at least an episodic heating and cooling of Io, albeit not as low as would be required by the equilibrium model. Now we would be observing the system on its way out of the resonance, which lowers the eccentricities and, hence, the heating rates. Eventually Io cools and the dissipation function $`Q_I`$ of Io increases to the point at which the Jupiter torque “wins,” pushing the system back to higher eccentricities to repeat the cycle (?, ?). It is not certain that Io could cool rapidly enough in this scenario, but $`Q_J10^6`$ could delay the evolution from deep resonance sufficiently to prevent its destruction (?). Major problems with this scenario were that (i) there was no mechanism for placing the system in the resonances at the time of satellite formation and (ii) evolution through the high-eccentricity phases could not be followed with the analytic theory. Here, we describe a mechanism that addresses both of these problems and adds support to Greenberg’s hypothesis.

Constraints on the formation of the Galilean satellites of Jupiter have been established by theoretical analysis and the Galileo spacecraft observations. The nearly coplanar and nearly circular orbits of the four large satellites imply formation within a dissipative, equatorial disk of gas and solid particles formed, perhaps, by one of several processes during the formation of Jupiter (?). Ganymede and Callisto are about 50% ice and 50% rock by mass (?). Europa has a layer of water and ice 80 to 170 km thick (?). The large amount of relatively volatile material in Ganymede and Callisto implies that the disk temperature remained sufficiently cool during the accretion process that icy particles could persist. Callisto may not be fully differentiated (?, ?), a state that requires an accretion time scale $`10^6`$ years to keep interior temperatures below 273 K (?, ?).

Accretion models for the Galilean satellites have typically assumed an initial disk with sufficient solid material to form the satellites and with gas added to bring the mixture up to solar composition — the so-called minimum mass sub nebula (MMSN) (?). But a careful evaluation of the processes in such a massive disk (?) indicate that temperatures would generally be too high and accretion rates too fast to allow formation of the satellites with their observed properties. In particular, it would be difficult if not impossible to keep Callisto only partially differentiated. Therefore, it is likely that the satellites formed near the end of Jupiter’s formation when the flux of material into Jupiter through an accretion disk was greatly reduced (?). Theoretical calculations (?) show that Jupiter would limit this flux by opening a gap in the disk of material surrounding the Sun when the planet’s mass approached its current value. Hydrodynamical simulations (?, ?, ?, ?) show that gas and particles continue to trickle through the gap in two streams entering at the inner and outer Lagrange points to form an accretion disk around Jupiter. The surface mass density of the accretion disk resulting from this scenario is sufficiently less than the MMSN to allow the slow accumulation of the satellites in a cool environment.

By modeling processes in this disk, Canup and Ward (?) have found a scenario of satellite accretion that appears to satisfy all of the observational and theoretical constraints imposed by the satellites’ current properties. The important characteristics of the model, for our purposes, are that the formation time scales of the individual Galilean satellites are similar and that the satellites migrate toward Jupiter while the disk persists. The migration is the result of the gravitational interaction of the satellites with spiral density waves generated in the disk material at mean motion commensurabilities with the satellites (Lindblad and corotation resonances). The satellites are not sufficiently massive to open gaps in the disk and, thus, they execute the so-called Type I drift (?), whose time scale is inversely proportional to both the disk surface mass density and the satellite mass.

The key ingredient of the Canup-Ward model that allows an alternative history of the Galilean satellite resonant system is that Ganymede’s orbit converges on those of Europa and Io because of its faster Type I migration induced by its greater mass. Callisto is left behind because of its smaller mass. Various choices of the viscosity parameter $`\alpha `$ and opacity of the disk lead to many different migration rates for the satellites, but all reasonable choices lead to converging orbits of the inner three satellites \[figs. 7c, 8c, and 9c of (?)\] because of Ganymede’s more rapid migration. Any other model of formation will also suffice for our purposes here, provided that the formation time scales of all the satellites are comparable and that the disk persists long enough for differential migration to lead to converging orbits (?).

We demonstrate (Fig. 1) that the Laplace relation can be established primordially during the final stages of satellite accretion by allowing Ganymede’s orbit to migrate inward on the $`10^5`$ year time scale of the Canup-Ward model (?), whereas Europa and Io initially migrate more slowly on a $`2\times 10^5`$ year time scale. (The very rapid migration rates implied for disk densities corresponding to the MMSN would most likely preclude the capture of the satellites in any resonances during differential migration.) The satellites start in an initial configuration where the orbits are separated by distances greater than those appropriate to the 2:1 mean motion commensurability. The numerical calculations are done both by the Wisdom-Holman (?) symplectic integrator contained in the SWIFT package (?) and in a center-of-mass cartesian coordinate system with a Bulirsch-Stoer algorithm. The algorithms are modified as described in (?) to accommodate the migration and eccentricity damping, and both algorithms yield similar results. Europa is first captured into the 2:1 mean motion resonances with Ganymede and Io is subsequently picked up in the 2:1 resonances with Europa, simultaneously establishing the Laplace relation (?).

The orbital eccentricities of the satellites are damped by interactions with the gas disk \[e.g., (?)\] and eventually by tidal dissipation in the satellites \[e.g., (?)\]. In the case of disk interactions, the ratio of the inverse time scales is given by (?)

$$\frac{1/\tau _e}{1/\tau _a}=\frac{|\dot{e}/e|}{|\dot{a}/a|}=9\frac{C_e}{C_a}\left(\frac{H}{a}\right)^230,$$
(1)
where $`H/a`$ is the ratio of the scale height of the disk $`H`$ to the semimajor axis of the orbit $`a`$. $`H/a0.1`$ is assumed constant in the disk, where the numerical value follows from the conditions in the Canup-Ward model (?). $`C_a3`$ is the coefficient of the expression for $`\dot{a}/a`$ (?, ?), and $`C_e0.1`$ is the coefficient in the expression for $`\dot{e}/e`$ (?). For the integrations, we assumed that the time scales for migration in $`a`$ and in eccentricity damping are constant and, hence, the ratio of the inverse time scales is fixed, with $`|\dot{e}/e|/|\dot{a}/a|=30`$ from Eq. 1 assumed for all three satellites. The eccentricity damping causes the eccentricities to reach constant asymptotic values (Fig. 1), which are reduced if the ratio of the inverse time scales in Eq. 1 is increased. The evolution is reasonably independent of the rate of migration up to a certain point. An example where the satellite migration rates were 30 times larger led to a very similar evolution.

Although all of the resonance variables, including the Laplace angle $`\theta _5`$, are librating in the primordial configuration of the satellites (Fig. 1), the geometry of the system differs from the current configuration, where all resonance variables, except $`\theta _4`$, librate about either $`0^{}`$ or $`180^{}`$, as described above. The centers of libration of $`\theta _5`$ and $`\theta _2`$ are near $`80^{}`$ instead of $`180^{}`$ and that of $`\theta _1`$ is near $`30^{}`$ instead of $`0^{}`$ (Fig. 2). $`\theta _4`$ would normally librate about $`180^{}`$ for smaller eccentricities, but here it is near $`10^{}`$. Only $`\theta _3`$ has the same libration center near $`0^{}`$ that it has now. Conjunctions would no longer occur near the lines of apsides of the orbits. This is consistent with the locus of stable libration centers for $`\theta _5`$ branching symmetrically from $`180^{}`$ as the system is driven deeper into the resonances than the current configuration (?). The occurrence of libration centers of both $`\theta _2`$ and $`\theta _5`$ near $`80^{}`$ in some of our calculations is verification of the symmetric branches.

As the satellites move closer to Jupiter and the disk dissipates, tidal dissipation in the satellites eventually dominates the eccentricity damping. At the same time, the dissipation of tidal energy in Jupiter from tides raised by the satellites results in the torques that tend to push the satellites outward thereby tending to increase the eccentricities \[e.g., (?)\]. Because the current orbital eccentricities are smaller than the primordial asymptotic values (Fig. 1), the tendency for the eccentricities to decrease due to dissipation in the satellites must dominate the effect of the tides raised on Jupiter that tend to increase the eccentricities—at least at the end of migration. The tidal time scale for eccentricity damping is $`|e/\dot{e}|6.5\times 10^4Q_I/f`$ years for Io, which is derived from expressions given in (?). Here $`Q_I30\mathrm{to}\mathrm{\hspace{0.17em}100}`$ for rocky material and $`f`$ is a factor to account for an inhomogeneous structure of Io, which can be as high as 13 for a two layer model consisting of a solid mantle over a molten interior (?). The tidal time scale for eccentricity damping is probably more than 6 times as long as the migration time scale of the Canup-Ward model, and possibly 100 times as long as the disk eccentricity damping time scale if Eq. 1 is correct. This means that tides will not take over the eccentricity damping until after the migration has essentially stopped.

After migration has stopped, we assumed that the dissipation in Io is reducing its eccentricity on the same time scale as disk interaction did during the migration ($`6.7\times 10^3`$ years) (Figs. 3 and 4). We maintained this time scale for the damping calculation to keep the computer time within reasonable bounds. The algorithm is such that forces are applied to Io that would cause the assumed eccentricity damping if Io were isolated but, because it must drag down the eccentricities of all the satellites locked in the resonances, the actual time scale is much longer. The system naturally relaxes to the current configuration through a series of states with offset libration centers, much like Greenberg described (?). The eccentricities at the end of the calculation have almost reached the values observed currently, and the libration centers have migrated to the current geometry from their positions at the end of migration (Fig. 2)—even to the detail that $`\theta _4`$ has ended up circulating just as we observe now.

The eccentricity of Ganymede, $`e_3`$, evolves in an interesting way in conjunction with the behavior of $`\theta _4`$ (Figs. 3 and 4). Initially, the center of libration of $`\theta _4`$ shifts from near $`10^{}`$ to near $`180^{}`$ in about $`10^4`$ years while $`e_3`$ decreases. At $`1.5\times 10^4`$ years, a sudden jump in $`e_3`$, perhaps from an encounter with a secondary resonance (?), results in $`\theta _4`$ changing to circulation (angle increasing or decreasing without bound) where it remains until $`t=5\times 10^4`$ years. Here $`e_3`$ suddenly becomes smaller, perhaps from another change in the effect of the secondary resonance, and $`\theta _4`$ again librates about $`180^{}`$. After $`10^4`$ years, $`\theta _4`$ switches from libration to permanent circulation due to the slow reduction in the forced eccentricity while the free eccentricity remains near $`0.0016`$. In another calculation with a different migration rate, the resonance variable $`\theta _4`$ shows similar behavior but without the sudden changes in $`e_3`$ and with different timing of events. In both cases, $`\theta _4`$ ends up circulating with $`e_30.0016`$. This is about the current value of Ganymede’s free eccentricity, but there may be circumstances where less excitation of the free eccentricity leaves $`\theta _4`$ librating.

A primordial origin of the Laplace relation relieves Jupiter of providing sufficient tidal torque over history to expand the orbits enough for the completely tidal origin. However, it has been argued that even if the origin of the resonances is primordial, the average $`Q_J`$ cannot be far from the bounds $`6\times 10^4<\text{ }Q_J<\text{ }2\times 10^6`$ established by (?, ?, ?, ?, ?). The current dissipation in Jupiter must still be sufficient to prevent the destruction of the resonances from the tidal dissipation in Io and might be sufficiently high to preserve the system in equilibrium (?, ?). In contrast, could it be that we are witnessing the further decay of the Laplace resonance due to the dissipation in Io (?, ?)? Because the energy dissipated must come from the orbits, a continuing decay of the Laplace configuration would result in an increasing $`n_1`$. Measures of $`\dot{n}_1=\ddot{\lambda }_1`$, based on 300 years of precise ephemerides of the Galilean satellites coupled with current precise observations, imply either an acceleration of Io’s mean motion, indicating further decay of the resonant configuration (?, ?), or insufficient deceleration (?) from Jupiter torques on Io to accommodate the high heat flux from Io in an equilibrium eccentricity configuration. The heat flux is now estimated to be about $`3\mathrm{W}/\mathrm{m}^2`$ (?). Although the measurements of $`\dot{n}_1`$ remain controversial, we should keep in mind that a heat flux from Io of $`3\mathrm{W}/\mathrm{m}^2`$ would require a current $`Q_J2\times 10^4`$ to maintain the system in equilibrium, which is a factor of 3 lower than the minimum historical average $`Q_J=6\times 10^4`$ from the proximity of Io to Jupiter (?). But there is no reason to believe that the current $`Q_J`$ should be constrained much by a bound on its historical average (?).

A primordial origin of the Laplace relation among the Galilean satellites, created by differential migration of the satellites in the gaseous disk, with the help of tidal dissipation, reproduces every detail of the current configuration. Changing assumptions will change details of this scenario such as the asymptotic values of the eccentricities during the migration; a different damping rate may change the circulation of $`\theta _4`$. But as long as the satellites exist simultaneously in the disk and the more rapid migration of Ganymede forces the inner-three orbits to converge and capture into the Laplace resonance results (?), the qualitative features of this alternative early history will persist.

* We thank R. Canup and W. Ward for sending us an early preprint of their paper on Galilean satellite origin. This research was supported in part by NASA PGG grants NAG5 3646 and NAG5 11666 and SSO grants NAG5 7177 and NAG5 12087.