text:
# Partonic Effects in Heavy Ion Collisions at RHIC

## 1 Introduction

In heavy ion collisions at RHIC, a quark-gluon plasma is expected to be created in the initial stage. Many observables have been suggested as possible signals for the quark-gluon plasma. These include enhanced production of dileptons  and strange hadrons , suppression of $`J/\psi `$ production , Hanbury-Brown-Twiss interferometry between pions , elliptic flow of various particles , quench of jet production , and fluctuations in the net charge , net baryon  as well as total baryon numbers . However, to establish these observables as unambiguous signals for the quark-gluon plasma requires that they are not affected appreciably by the final hadronic dynamics in the collisions. To study quantitatively these signals, we have developed a multiphase transport model  that includes both the dynamics of initial partonic and final hadronic matters. With this model, we have studied not only the rapidity and transverse momentum distributions of various charged particles  but also their elliptic flows . Furthermore, we have investigated multistrange baryon production , $`J/\psi `$ suppression  and pion interferometry  in these collisions. We find that the model describes very well the experimental data on the charged particle rapidity distributions and transverse momentum spectra. However, these results are more sensitive to the hadronic interactions than the partonic interactions. On the other hand, both the charged particle elliptic flow and the two-pion correlation function are affected by partonic interactions, and to explain the observed large elliptic flow and the measured pion correlation function in the experiments requires a large parton-parton scattering cross section in the initial partonic matter.

## 2 Multiphase Transport Model

In the AMPT model, the initial conditions are obtained from the HIJING model  by using a Woods-Saxon radial shape for the colliding nuclei and including the nuclear shadowing effect on parton production via the gluon recombination mechanism of Mueller-Qiu . After the colliding nuclei pass through each other, the Gyulassy-Wang model  is then used to generate the initial space-time information of partons. In the default HIJING model, these minijet partons are allowed to lose energy via the gluon splitting mechanism and transfer their energies to the nearby strings associated with initial soft interactions. Such a jet quenching is replaced in the AMPT model by explicitly taking into account parton-parton collisions via Zhang’s Parton Cascade (ZPC) . At present, only gluon elastic scatterings with a default cross section of 3 mb are included, so the partons do not suffer any inelastic energy loss as they traverse the dense matter. After partons stop interacting, they combine with their parent strings and, after an average proper formation time of 0.7 fm/c, are converted to hadrons using the Lund string fragmentation model . Dynamics of the resulting hadronic matter is then described by a relativistic transport model (ART) , which has been improved to include baryon-antibaryon production from meson-meson interactions and their annihilation .

The parameters in the AMPT model are determined by fitting the experimental data from Pb+Pb collisions at center-of-mass energy $`\sqrt{s}=17A`$ GeV at SPS . Specifically, to describe the measured net baryon rapidity distribution, we have included in the Lund string fragmentation model the popcorn mechanism for baryon-antibaryon production with equal probabilities for baryon-meson-antibaryon and baryon-antibaryon configurations. Also, to account for both the pion yield and the enhanced kaon yield in this reaction, we have adjusted the two parameters in the momentum splitting function used in the Lund string fragmentation model, following the expectation that the string tension is modified in the dense matter formed in the initial stage of heavy ion collisions.

Since the probability for minijet production is very small in heavy ion collisions at SPS energies, the partonic stage does not play any role in these collisions. We find that final-state hadronic scatterings reduce the proton and antiproton yields, but increase kaon and antikaon production by about 20%. In contrast, the kaon yield in the default HIJING model is smaller than our final result by about 40%.

## 3 Rapidity and Transverse Momentum Distributions

Results on the charged particle rapidity distributions in central Au+Au collisions at RHIC are shown in the left four panels of Fig. 1 for center-of-mass energies of $`56A`$ GeV (dashed curves) and $`130A`$ GeV (solid curves) together with data from the PHOBOS collaboration . The measured total charged particle multiplicities at mid-pseudorapidity at both energies are well reproduced by our model. Also, the $`\overline{p}/p`$ ratio predicted by the AMPT model is consistent with that measured in experiments . The transverse momentum spectra of pions, kaons, and protons in Au+Au collisions at $`\sqrt{s}=130A`$ GeV are shown in the right panel of Fig. 1. The results from the the AMPT model are seen to also agree with the experimental data from the PHENIX collaboration .

Furthermore, the AMPT model describes very well the centrality dependence of the charged particle rapidity distribution measured by the BRAHMS collaboration  as shown in Fig. 2 by dashed curves. Compared to the default HIJING model given by solid curves, hadronic scatterings are responsible for the broader rapidity distributions seen in the experiments.

Without hadronic interactions, there is a significant increase in the numbers of total charged particles, pions, protons, and antiprotons at midrapidity. The kaon number is, on the other hand, reduced slightly. As a result, the ratios of $`\overline{p}/p`$ and $`K^+/\pi ^+`$ in the absence of final-state hadronic interactions are $`0.80`$ and $`0.13`$, respectively, instead of $`0.66`$ and $`0.18`$ from the default AMPT model, which are more consistent with the data. Although the default HIJING  gives a total charged particle multiplicity at midrapidity that is consistent with the PHOBOS data, including hadronic scatterings would reduce its prediction appreciably.

Effects of partonic dynamics on final hadronic observables can be studied in the AMPT model by turning off the partonic cascade. We find that this leads to less than $`5\%`$ change in the final charged particle yields in Au+Au collisions at $`\sqrt{s}=130A`$ GeV. The multiplicity distribution of hadrons is thus not very sensitive to parton elastic scatterings.

## 4 Elliptic Flow

Elliptic flow in heavy ion collisions measures the anisotropy of particle momentum distributions in the plane perpendicular to the beam direction. It results from the initial spatial anisotropy in non-central collisions and is thus sensitive to the properties of the hot dense matter formed during the initial stage of heavy ion collisions . In transport models based only on the parton cascade, the elliptic flow has been shown to be sensitive to the parton scattering cross section, and a large value can be obtained with a large cross section .

Results on the charged particle elliptic flow from the AMPT model for Au+Au collisions at center-of-mass energy of $`130A`$ GeV are shown in Fig. 3 as a function of both the charged particle multiplicity (left panel) and transverse momentum (right panel). Theoretical values from the default AMPT model (no string melting) are much smaller than that measured in the experiments by the STAR collaboration . This is mainly due to the small fraction of energy that is carried by minijet partons and the lack of transverse collective motion of the strings. As a result, the elliptic flow is significantly reduced after minijet partons combine with strings and fragment to hadrons.

Since the initial energy density of produced matter in central Au+Au collisions at RHIC is more than an order of magnitude higher than the critical energy density ($`1`$ GeV/fm3) expected for the formation of the quark-gluon plasma, keeping strings in the high energy density region thus underestimates the partonic effects in these collisions. To model the melting of strings into partons in regions of high energy density, we have extended the AMPT model to allow the initially excited strings to fragment into hadrons using the LUND fragmentation model and then convert these hadrons to their constituent quarks but with bare masses, i.e., a meson is converted to a quark and an anti-quark, while a baryon is converted to three quarks. We further assume that quarks are produced isotropically in the rest frame of a hadron and start to interact only after a proper formation time given by the inverse of the hadron transverse mass. Scatterings among these quarks are then treated using the parton cascade ZPC. After the quarks stop interacting, we model the hadronization by combining the nearest quark and antiquark into a meson and three quarks into a baryon with the same flavor. The resulting hadrons are given an additional formation time of 0.7 fm$`/c`$ in their rest frame and then imported to the ART hadronic transport model to take into account their rescatterings.

With strings converted to partons, the initial energy originally stored in strings also contributes to the parton dynamics. This leads to a larger elliptic flow and also a stronger dependence on the parton scattering cross section, making it possible to determine the strength of partonic interactions from the final elliptic flow. As shown in Fig. 3, a large parton scattering cross section of 6-10 mb is needed to account for the observed dependence of charged particle elliptic flow on the total charged particle multiplicity and their transverse momentum.

## 5 Pion Interferometry

Particle interferometry based on the Hanbury-Brown Twiss (HBT) effect can provide information not only on the spatial extent of an emission source but also on its expansion velocity and emission duration . In particular, the long emission time as a result of the quark-gluon plasma to hadronic matter phase transition in the initial stage of heavy ion collisions is expected to lead to a much larger radius parameter along the direction of the total momentum of detected two particles than that perpendicular to both this direction and the beam direction.

From the emission function, which is given by the space-time and momentum distribution of particles at freezeout in the AMPT model, we can evaluate the correlation function $`C(𝐐,𝐊)`$ of two identical pions including their Coulomb interaction. Here, $`𝐐`$ and $`𝐊`$ are the relative and total momenta of the two pions. The six-dimension correlation function is usually shown as a function of the invariant momentum ($`Q_{\mathrm{inv}}=\sqrt{Q^2}`$) or as a function of the projection of the relative momentum $`𝐐`$ in the “out-side-long” ($`\mathrm{o}\mathrm{s}\mathrm{l}`$) system , defined by the beam direction ($`Q_{\mathrm{long}}`$), the direction along the total transverse momentum of the two pions ($`Q_{\mathrm{out}}`$), and the direction orthogonal to the above two directions ($`Q_{\mathrm{side}}`$).

In the left panel of Fig. 4, we show the one-dimensional projections of the correlation function for midrapidity pions with transverse momentum $`125<p_{}<225`$ MeV$`/c`$ in Au+Au collisions at $`\sqrt{s}=130A`$ GeV. Also shown is the measured $`\pi ^{}`$ correlation function by the STAR collaboration without correcting the effect due to Coulomb interaction . In evaluating the one-dimensional projection of the correlation function onto one of the $`Q_{\mathrm{out}}`$, $`Q_{\mathrm{side}}`$, $`Q_{\mathrm{long}}`$ axes, we have integrated the other two $`𝐐`$ components over the range $`035`$ MeV/$`c`$. Dash-dotted curves in the figure are from the default AMPT model with a parton scattering cross section of $`\sigma _\mathrm{p}=3`$ mb, while other curves are from the extended AMPT model with string melting but different values of parton scattering cross section. It is seen that with string melting both the width of the $`Q_{\mathrm{inv}}`$ correlation function and its height decrease with increasing $`\sigma _\mathrm{p}`$. To reproduce the measured one-dimensional correlation functions from the STAR collaboration, we need a parton scattering cross section of about 10 mb in the extended AMPT model with string melting.

In the right panel of Fig. 4, we show the distribution of the emission source from the AMPT model for midrapidity pions with $`125<p_{}<225`$ MeV$`/c`$ in the $`x_{\mathrm{out}}x_{\mathrm{side}}`$ and $`x_{\mathrm{out}}t`$ spaces. It is seen that the emission source shows a large halo around a central core. The halo consists of not only pions from decays of long-lived resonances such as the $`\omega `$ but also thermal pions. The latter contribution to the halo becomes increasingly important when the collective expansion velocity of the source is large as a result of increasing parton cross section. The emission source also shows a strong correlation between $`x_{\mathrm{out}}`$ and $`t`$, with the width of its $`x_{\mathrm{out}}`$ distribution increasing with the emission time $`t`$.

The size of emission source can be estimated from the curvature of the correlation function at $`𝐐=0`$, i.e, the second moments of the emission function from the AMPT model. In the left panel of Fig. 5, we show the dependence of the source radii of midrapidity pions ($`0.5<y<0.5`$) on the average transverse mass $`m_\mathrm{T}`$ of the pion pair from the AMPT model with string melting and a parton cross section of 10 mb. At low transverse mass, all three radii $`R_{\mathrm{out}}^0`$, $`R_{\mathrm{side}}^0`$, and $`R_{\mathrm{long}}^0`$ are large and have values between 10 and 25 fm/c. The emission duration $`\mathrm{\Delta }\tau `$ also lasts a long time of more than 25 fm/c. As the average pion transverse mass increases, both the emission duration and the source radii, particularly $`R_{\mathrm{long}}^0`$, are found to decrease. The ratio $`R_{\mathrm{out}}^0/R_{\mathrm{side}}^0`$ has a value of about 1.2 for most values of $`m_\mathrm{T}`$ and is smaller than the prediction based on the hydrodynamical model with freezeout treated via the hadronic transport model .

Empirically, radii of the emission source are obtained from fitting the measured correlation function by a Gaussian function in the relative momentum of the pion pair. Results from such a fit to the correlation function obtained from the AMPT model with string melting and $`\sigma _\mathrm{p}=10`$ mb are shown in the right panel of Fig. 5 as functions of the average pion transverse mass. The fitted radius parameters $`R_{\mathrm{out}}`$, $`R_{\mathrm{side}}`$, and $`R_{\mathrm{long}}`$ are about a factor of two to three smaller than corresponding source radii evaluated directly from the emission function, indicating that the emission source from the AMPT model deviates appreciably from a Gaussian one. Furthermore, the ratio $`R_{\mathrm{out}}/R_{\mathrm{side}}`$ of the fitted Gaussian radius parameters is much closer to one than that from the source radii obtained from the emission function.

## 6 Conclusions

Using a multiphase transport model, which includes both initial partonic and final hadronic interactions, we have studied the rapidity distributions of charged particles such as protons, antiprotons, pions, and kaons in heavy ion collisions at RHIC. With the model parameters constrained by data from central Pb+Pb collisions at $`\sqrt{s}=17A`$ GeV at SPS, theoretical results on the total charged particle multiplicity at midrapidity in central Au+Au collisions at $`\sqrt{s}=56A`$ and $`130A`$ GeV agree quite well with the data from the PHOBOS collaboration. Also, the predicted $`\overline{p}/p`$ ratio is consistent with that from the STAR and BRAHMS collaborations. Furthermore, the centrality dependence of the charged particle rapidity distribution measured by the BRAHMS collaboration is well reproduced. These hadronic observables are, however, less sensitive to the initial partonic interactions than to the final hadronic interactions. On the other hand, partonic effects are important for both the elliptic flow of charged particles and the two-pion correlation function. To reproduce the measured elliptic flow and two-pion correlation function by the STAR collaboration, we need not only to convert the initial strings to partons but also to use a large parton-parton scattering cross section. Our results thus demonstrate the possibility of studying the dynamics of partonic matter in heavy ion collisions at RHIC. To further test these partonic effects, it is necessary to study their dependence on the initial parton distribution, the inelastic scatterings in the partonic stage, and the mechanism for melting the initial strings as well as for hadronization. Such an improved study based on the AMPT model will be very useful in finding the signals for the quark-gluon plasma formed in heavy ion collisions at RHIC.

## Acknowledgments

This talk is based on work supported by the National Science Foundation under Grant Nos. PHY-9870038 and PHY-0098805, the Welch Foundation under Grant No. A-1358, and the Texas Advanced Research Program under Grant No. FY99-010366-0081.


format: text

text:
# Threshold temperature for pairwise and many-particle thermal entanglement in the isotropic Heisenberg model

\[

## Abstract

We study the threshold temperature for pairwise thermal entanglement in the spin-1/2 isotropic Heisenberg model up to 11 spins and find that the threshold temperature for odd and even number of qubits approaches the thermal dynamical limit from below and above, respectively. The threshold temperature in the thermodynamical limit is estimated. We investigate the many-particle entanglement in both ground states and thermal states of the system, and find that the thermal state in the four-qubit model is four-particle entangled before a threshold temperature.

PACS numbers: 03.65.Ud, 75.10.Jm

\]

Entanglement is a nonlocal quantum correlation and entangled states constitute indeed a valuable resource in quantum information processing . Entanglement in systems of interacting spins  as well as in systems of indistinguishable particles  has been investigated. In particular entanglement in both the ground state  and thermal state  of a spin-1/2 Heisenberg spin chain have been analyzed in the literature. The intriguing issue of the relation between entanglement and quantum phase transition  have been addressed in a few quite recent papers .

The state we will consider is the thermal equilibrium state in a canonical ensemble. The system state at finite temperature is given by the Gibb’s density operator $`\rho _T=\mathrm{exp}\left(H/kT\right)/Z,`$ where $`Z=`$tr$`\left[\mathrm{exp}\left(H/kT\right)\right]`$ is the partition function, $`H`$ the system Hamiltonian, $`k`$ is Boltzmann’s constant which we henceforth will take equal to 1, and $`T`$ the temperature. As $`\rho _T`$ represents a thermal state, the entanglement in the state is called thermal entanglement. It is important to stress that, although the central object of statistical physics, the partition function, is determined by the eigenvalues of $`H`$ only, thermal entanglement properties generally require in addition the knowledge of the energy eigenstates. In some models such as the isotropic Heisenberg qubit rings the pairwise thermal entanglement can be completely determined by the partition function .

Let us briefly review the main results of Ref. . We consider a physical model of a ring of $`N`$ qubits interacting via the isotropic Heisenberg Hamiltonian

$$H=J\underset{i=1}{\overset{N1}{}}S_{i,i+1}+JS_{N,1},$$
(1)
where $`S_{j,j+1}=`$$`\frac{1}{2}\left(1+\stackrel{}{\sigma }_i\stackrel{}{\sigma }_{i+1}\right)`$ is the swap operator between qubit $`i`$ and $`j`$, $`\stackrel{}{\sigma }_i=(\sigma _{ix},\sigma _{iy},\sigma _{iz})`$ is the vector of Pauli matrices, and $`J`$ is the exchange constant. In Ref.  we proved that there is no thermal entanglement between two nearest-neighbour qubits in the ferromagnetic ($`J<0`$) isotropic Heisenberg rings at any temperature. So we only consider the antiferromagnetic case and the exchange constant $`J`$ is then assumed to be 1. A direct relation is established between the concurrence quantifying two-qubit entanglement and a macroscopic thermodynamical function, the internal energy $`U`$:

$$C(N)=\mathrm{max}(0,\frac{U}{N})=\mathrm{max}(0,\frac{1}{NZ}\frac{Z}{\beta }),$$
(2)
where $`\beta =1/T`$ and $`C(N)`$ refers to the concurrence for two nearest-neighbor qubits of the $`N`$-qubit ring. The entanglement is uniquely determined by the partition function of the system. Note that there is a slight difference between the expression of the concurrence and the one given in Ref.  due to a constant term added in Eq.(1) for convenience of numerical calculations. The reason why the entanglement is determined only by the eigenvalues of the system is that there are many symmetries in the isotropic Heisenberg model such as translational invariance and rotational symmetry along $`x,y,z`$ directions, et al. If we apply a nonzero magnetic field along $`z`$ direction we only have rotational symmetry along $`z`$ direction. In this case we can not determine the entanglement by eigenvalues alone.

From Eq.(2) the concurrence for the ground state (GS) with energy $`E_{GS}`$ is then followed by

$$C_{GS}(N)=\mathrm{max}(0,\frac{E_{GS}}{N}),$$
(3)
which is obtained in Ref.  for even-number qubits. But now we know that it is valid for any number of qubits. From the available results  of the ground-state energy we can determine the concurrence directly. For the case of even-number qubits the concurrence is given in Ref. . Here in Table I we give the concurrence of ground states for $`N`$ qubits up to $`N=11`$. From the table we find that the entanglement increases with the increase of odd number $`N`$. On the contrary the entanglement decreases with the increase of even $`N`$. We also see that the ground state is always entangled except the case of $`N=3`$.

Since $`U/T0`$ and $`C_{GS}(3)=0`$ there is no entanglement at any temperature for the case of $`N=3`$. For other cases, $`C_{GS}(N)>0`$ $`(N3),`$ and the concurrence (2) is a monotonic decreasing function of the temperature. It is then obvious that there exists a $`N`$-dependent threshold temperature $`T_{th}(N)`$ after which the thermal entanglement disappears. The threshold temperature is determined by the equation $`U(T_{th})=0`$. Therefore the problem left is to examine the threshold temperature and we address it in this brief report. We will determine the threshold temperature both analytically and numerically.

The number of qubits we will consider is up to $`N=11`$. Analytical results can be obtained for very small rings. For $`N=2`$, the eigenvalues of $`H`$ is easily found to be -2(1) and 2 (3), where the number in the parenthesis is the degeneracy. For $`N=3`$ the eigenvalues of $`H`$ are 0(4) and 3(4). For $`N=4,`$ by direct diagonalization, the eigenvalues are given by -2(1), 0 (3), 2(7) and 4(5) (see later discussions). From the eigenvalues the partition functions $`Z(N)`$ $`(N=2,3,4)`$ simply follow as

$`Z(2)`$ $`=`$ $`e^{2\beta }+3e^{2\beta },`$ (4)
$`Z(3)`$ $`=`$ $`4+4e^{3\beta },`$ (5)
$`Z(4)`$ $`=`$ $`3+e^{2\beta }+7e^{2\beta }+5e^{4\beta }.`$ (6)
The substitution of Eq.(6) into (2) leads to

$`C(2)`$ $`=`$ $`{\displaystyle \frac{1}{Z(2)}}\mathrm{max}(0,e^{2\beta }3e^{2\beta }),`$ (7)
$`C(3)`$ $`=`$ $`{\displaystyle \frac{1}{Z(3)}}\mathrm{max}(0,4e^{3\beta })=0,`$ (8)
$`C(4)`$ $`=`$ $`{\displaystyle \frac{1}{2Z(4)}}\mathrm{max}(0,e^{2\beta }7e^{2\beta }10e^{4\beta }),`$ (9)
which are explicit analytical expressions for the concurrence. The threshold temperatures are directly obtained from the above equation: $`T_{th}(2)=4/\mathrm{ln}33.641,`$ $`T_{th}(3)=0,`$ and $`T_{th}(4)=1.726728`$.

For $`N5`$ we use numerical method to determine the threshold temperature. We mention that some analytical results of the eigenvalue problem for $`N=5,6`$ and $`7`$ can be found in Ref. . Bethe’s Ansatz can be used for the calculations of eigenvalues. Here we use the exact diagonalization method to calculate the eigenvalues. The obvious symmetry $`[H,S_z]=0`$ allows to decompose the Hilbert space of the system into a set of orthogonal subspaces. Here $`S_z=1/2_{i=1}^N\sigma _{iz}`$ is the collective operator for $`N`$ qubits. The dimension of the subspace is much smaller than the original Hilbert space, and then we can diagonalize the Hamiltonian in the subspaces.

The numerical results for the threshold temperature are given in Table II and III, from which we find that threshold temperature increases (decreases) as $`N`$ increases for the case of odd (even) qubits. We also observe that the threshold temperature converges very quickly. Then we can estimate the threshold temperature in the thermodynamic limit $`T_{th}(\mathrm{})\mathrm{1.\hspace{0.17em}5903}`$ from the numerical results of these finite systems. So we have the following relation

$$T_{th}(N)|_{N\text{ odd}}<T_{th}(\mathrm{})<T_{th}(N)|_{N\text{ even}},$$
(10)
which implies that the threshold temperatures for $`N`$odd and $`N`$ even approach the thermal dynamical limit from below and above, respectively.

We have considered only pairwise thermal entanglement of a multiqubit system. Another type of entanglement is the $`N`$-particle entanglement which involves all $`N`$ particles. As far as we know there are no discussions on many-particle thermal entanglement. In the following we will consider the many-particle entanglement in the ground states of the Heisenberg models up to 11 qubits and then study the many-particle thermal entanglement within the four-qubit Heisenberg model.

For a three-qubit pure state, there exists a good measure of three-particle entanglement. Here we use the state preparation fidelity $`F`$ of a $`N`$-qubit state $`\rho `$ to investigate the many-particle entanglement. The state $`\rho `$ can be either pure or mixed. The fidelity is defined as

$$F(\rho )=\mathrm{\Psi }_{\text{GHZ}}|\rho |\mathrm{\Psi }_{\text{GHZ}},$$
(11)
where $`|\mathrm{\Psi }_{\text{GHZ}}=1/\sqrt{2}(|00\mathrm{}0+|11\mathrm{}1)`$ is the $`N`$-particle Greenberger-Horne-Zeilinger (GHZ) state.

The sufficient condition for $`N`$-particle entanglement is given by

$$F(\rho )>1/2.$$
(12)
We have a freedom to choose other GHZ states such as $`|\mathrm{\Psi }_{\text{GHZ}}=1/\sqrt{2}(|00\mathrm{}01\pm |11\mathrm{}10)`$, etc. By local unitary operations we can transfer these states to the original GHZ state and the operations do not change the entanglement. Detail discussions on this sufficient condition can be found in Ref..

Let us study the many-particle entanglement in the ground states of the Heisenberg models for the number of qubits up to 11. It is known that  the ground state which has $`S_z=0`$, is non-degenerate if $`N`$ is even, and the ground state which has $`S_z=\pm 1/2`$, is four-fold degenerate if $`N`$ is odd. Then the ground state is a pure (mixed) state for even (odd) $`N`$. We numerically find the ground state of the system, from which we calculate the fidelity. The form of the GHZ state is choose in order that the fidelity is as large as possible. The numerical results is given in Table IV. First we observe that the fidelity of four-qubit ground state is larger than one half, hence it is a many-particle entangled state. The fidelity is less than one half for $`N`$ from 2 to 11 except $`N=4`$. If we consider even or odd number $`N`$, the fidelity decreases as $`N`$ increases. From the analytical expressions of the ground states for $`N=3`$ and $`N=4`$ (see below), we can analytically obtain the fidelity $`F=1/6`$ for $`N=3`$ and $`F=2/3`$ for $`N=4`$, respectively.

Having found that the ground state of the four-qubit Heisenberg model is a 4-particle entangled state, we expect that the thermal state will be also a 4-particle entangled state for small temperature. We also expect a threshold temperature before which the ground state is 4-particle entangled. To study the many-particle thermal entanglement we need to know all the eigenvalues and eigenstates. For convenience following notations are used

$`|n_0`$ $`=`$ $`𝒯^{n1}|1_0,|1_0=|1000,`$ (13)
$`|n_1`$ $`=`$ $`𝒯^{n1}|1_1(n=1,2,3,4),|1_1=|1100,`$ (14)
$`|n_2`$ $`=`$ $`𝒯^{n1}|1_2(n=1,2),|1_2=|1010.`$ (15)
Here $`𝒯`$ is a unitary cyclic right shift operator defined by its action on the basis $`𝒯|m_1m_2m_3m_4=|m_4m_1m_2m_3`$ and $`|m_1m_2m_3m_4=|m_1|m_2|m_3|m_4`$.

A direct diagonalization of the four-qubit Hamiltonian gives all eigenvectors and eigenvalues,

$`E_0`$ $`=`$ $`4,|\mathrm{\Psi }_0=|0000,`$ (16)
$`E_{15}`$ $`=`$ $`4,|\mathrm{\Psi }_{15}=|1111,`$ (17)
$`E_k`$ $`=`$ $`(2+t_k+t_k^1),|\mathrm{\Psi }_k={\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{4}{}}}t_k^n|n_0,`$ (18)
$`E_{k+4}`$ $`=`$ $`(2+t_k+t_k^1),|\mathrm{\Psi }_{k+4}={\displaystyle \frac{1}{2}}{\displaystyle \underset{n=1}{\overset{4}{}}}t_k^n\mathrm{\Lambda }_x|n_0,`$ (19)
$`E_9`$ $`=`$ $`2,|\mathrm{\Psi }_9={\displaystyle \frac{1}{2\sqrt{3}}}\left({\displaystyle \underset{n=1}{\overset{4}{}}}𝒯^{n1}|1_12{\displaystyle \underset{n=1}{\overset{2}{}}}𝒯^{n1}|1_2\right),`$ (20)
$`E_{10}`$ $`=`$ $`0,|\mathrm{\Psi }_{10}={\displaystyle \frac{1}{\sqrt{2}}}(|1010|0101),`$ (21)
$`E_{11}`$ $`=`$ $`2,|\mathrm{\Psi }_{11}={\displaystyle \frac{1}{\sqrt{2}}}(|1100|0011),`$ (22)
$`E_{12}`$ $`=`$ $`2,|\mathrm{\Psi }_{12}={\displaystyle \frac{1}{\sqrt{2}}}(|1001|0110),`$ (23)
$`E_{13}`$ $`=`$ $`2,|\mathrm{\Psi }_{13}={\displaystyle \frac{1}{2}}(|1100+|0011|1001|0110),`$ (24)
$`E_{14}`$ $`=`$ $`4,|\mathrm{\Psi }_{14}={\displaystyle \frac{1}{\sqrt{6}}}\left({\displaystyle \underset{n=1}{\overset{4}{}}}𝒯^{n1}|1_1+{\displaystyle \underset{n=1}{\overset{2}{}}}𝒯^{n1}|1_2\right),`$ (25)
where the parameters $`t_k=\mathrm{exp}(i\pi k/2)(k=1,2,3,4)`$, and the operator $`\mathrm{\Lambda }_x=\sigma _x^4=\sigma _x\sigma _x\sigma _x\sigma _x`$.

These eigenvalues and eigenstates completely determine the thermal state $`\rho _T`$. We let our GHZ state be $`|\mathrm{\Psi }_{\text{GHZ}}=1/\sqrt{2}(|1010+|0101)`$. Then from Eq.(11) and all the eigenvalues and eigenstates we get

$$F(\rho _T)=\frac{e^{4\beta }/3+2e^{2\beta }/3}{3+e^{2\beta }+7e^{2\beta }+5e^{4\beta }}.$$
(26)
It is easy to verify that $`F=2/3`$ in the limit of $`T0`$ as we expected.

In Fig.1 we give the numerical result for the state preparation fidelity against temperature. For comparison we also plot the concurrence. First we observe that there exists a threshold temperature $`T_{th}`$ for the state preparation fidelity before which the thermal state $`\rho _T`$ is four-particle entangled state. The threshold temperature $`T_{th}0.83`$. Since the condition $`F>1/2`$ is a sufficient condition we can not say that the thermal state has no four-particle entanglement after $`T_{th}`$. From the figure we can see that the thermal state has both pairwise entanglement and four-particle entanglement before the threshold temperature. But both entanglement are not maximal.

In conclusion we investigated the threshold temperature for the pairwise thermal entanglement in the isotropic Heisenberg model with finite number of qubits up to $`N=11`$. We find that the threshold temperatures for odd and even number of qubits approach the thermal dynamical limit from below and above, respectively. Although we only make calculations of the threshold temperatures for small spin rings we can well estimate the threshold temperature in the thermodynamical limit as $`T_{th}(\mathrm{})\mathrm{1.\hspace{0.17em}5903}`$ since $`T_{th}(N)`$ converges very quickly as $`N`$ increases. By using one sufficient condition for many-particle entanglement, the state preparation fidelity is larger than one half, we investigated the many-particle entanglement of both ground states and thermal states in the Heisenberg models. We find that the ground state of the four-qubit model is a four-particle entangled state, and also find that there exists a threshold temperature $`T_{th}0.83`$ before which the thermal state is a four-particle entangled state. It is interesting to see that pairwise entanglement and four-particle entanglement (not maximal) coexists in this system.

Acknowledgments The authors thanks for the helpful discussions with Paolo Zanardi, Irene D’Amico, Barry C Sanders, and Prof. V.E.Korepin, H. Fu, A. I. Solomon, and Guenter Mahler. This work has been supported by the European Community through grant IST-1999-10596 (Q-ACTA).


text:
# Multi-object searching algorithm using subgrouped oracles

## Abstract

We present in this work, if a set of well organized suboracles is available, a quantum algorithm for multiobject search with certainty in an unsorted database of $`N`$ items. Depending on the number of the objects, the technique of phase tunning is included in the algorithm. If one single object is to be searched, this algorithm performs a factor of two improvement over the best algorithm for a classical sorted database. While if the number of the objects is larger than one, the algorithm requires slightly less than $`\mathrm{log}_4N`$ queries, but no classical counterpart exists since the resulting state is a superposition of the marked states.

Quantum mechanical algorithms are much more powerful than their classical counterparts becsause they are quantum physical processes possessing the unique feature of quantum parallelism due to superpositions and entanglements of quantum states. The quantum parallelism potentially can bring about an exponential speedup in running time over the classical computation processing. A typical example is the factorizing algorithm discovered by Shorshor  . The running time required in Shor’s algorithm is proportional to $`O(n^2)`$ for the factorization of a number registered in $`n`$ bits, while the running time for the classical algorithm would take $`O(\mathrm{exp}(n^{1/3}))`$. This is the first quantum algorithm that reduces a classical NP problem to a P problem, showing a quantum computer, if it can ever be built, can transcend any classical one. One the other hand, Grovergrover1 grover2  proposed the quantum searching algorithm, which performs a quadratic speedup over its classical counterpart. If there is an unsorted database of $`N`$ items, and out of which only one marked item satisfies a given condition, then using Grover’s algorithm one will find the object in $`O(\sqrt{N})`$ quantum mechanical steps, instead of $`O(N)`$ classical steps. The Grover algorithm can be extended to find $`M`$ objects out of a database of size $`N`$. It would then required $`O(\sqrt{N/M})`$ iterations to search the superposition of these $`M`$ marked states. The superposition of the marked states is an unique feature that only quantum algorithm can search and therefore no classical algorithms can ever do. It has been shown that Grovers original algorithm is optimal optimal1 optimal2 optimal3 , in the sense it requires the minimal queries to undertake the search task. Although the quantum searching algorithm did not provide an exponential speedup, it indeed greatly improved in reducing the exhausting effort needed in the classical counterpart. Up to date the searching problem remains an NP problem.

Grover’s algorithm utilizes a designed unitary operator, the Grover operator, successively on an uniform superposition of all the possible states to enhance the probability of the marked state. The Grover operator is a product of two operators with specific functions, namely, the selective inversion operator $`G_\tau `$ and the inverse-about-average operator $`G_{\overline{0}}`$. Expressed mathematically, the Grover operator is given by $`G=G_{\overline{0}}`$ $`G_\tau `$. Iterating $`G`$ about $`(\pi /4)\sqrt{N}`$ times on the initial state and then measuring, one will obtain the marked state. Alternatively, the searching probelm can be rephrased in terms of the action of an oracle. To carry out the search, one requires an oracle $`f(x)`$ with output $`0`$ or $`1`$ and, during the processing, is assigned to find a specific marked item $`\tau `$ corresponding to the query $`f(x)=1`$, for $`x=\tau `$ and $`f(x)=0`$, for $`x\tau `$. In the Grover algorithm each action of $`G_\tau `$ needs one oracle call, so the quantum searching task consume at least $`(\pi /4)\sqrt{N}`$ queries.

Instead of using one single oracle in Grover’s algorithm for the search in an unsorted database, Patelpatel  recently suggested that if a set of factorized oracles is available, one can locate a desired item in an unsorted database using $`O(\mathrm{log}_4N)`$ queries, which is a factor of two improvement over the best search algorithm for a classical sorted database. Patel’s scenario, in terms of the action of oracle calls, is to search a marked state

$$|\tau =|\tau ^{(n/2)}\mathrm{}|\tau ^{(2)}|\tau ^{(1)}$$
(1)
satisfying the factorizable function

$$f(\tau )=f_{n/2}(\tau ^{(n/2)})\mathrm{}f_2(\tau ^{(2)})f_1(\tau ^{(1)})=1,$$
(2)
provided that the state $`|\tau `$ is registered in $`n`$($`=\mathrm{log}_2N`$) qubits while each of the factorized state $`|\tau ^{(j)}`$, $`j=1,\mathrm{},n/2`$, is in two qubits. The outline of this factorized search algorithm is arranged to successively find the factorized state $`|\tau ^{(j)}`$ in the corresponding four-dimensional subspace in a single iteration according to the corresponding factorized oracle $`f_j(x)`$, $`x\{0,1\}^2`$, so the total queries are $`n/2`$. The factorized search algorithm is successful in searching one single marked state since the single state is natually factorizable. It turns out to be, however, that this algorithm will not be applicable to a multi-solution searching problem. As mentioned, if $`M`$ marked states are to be searched, a quantum algorithm is designed to enhance the probability amplitued of the superposition of these objects. The superposition is an entangled state, so the factorized search algorithm developed by Patel fails to solve the milti-solution searching problem. In this work, we are thus motivated to show how to search the superposition of the $`M`$ marked states in $`O(\mathrm{log}_4N)`$ iterations provided that a specific set of subgrouped oracles is available.

Suppose we have an unsorted database of $`N`$ items at hand and according to its size we need $`n`$($`=\mathrm{log}_2N`$) qubits to register the state in the database. The $`n`$ qubits are divided into $`\eta +1`$ subgroups in which all are of two qubits except one $`n_0`$-qubit subgroup, so $`n=2\eta +n_0`$, where $`n_0`$ is chosen depending on the number of the desired items $`M`$ and will be defined in what follows. The marked states then can be expressed as

$`|\tau _j`$ $`=`$ $`|\tau _j^{(\eta +1)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)}`$ (3)
$`=`$ $`|\tau _j^{(\eta +1)}\mathrm{}|\tau _j^{(2)}|\tau _j^{(1)},\text{ }j=1,\mathrm{},M,`$
where

$$\tau _j^{(1)}\{1,0\}^{n_0},\text{ }n_0=\mathrm{log}_2(4M)2,$$
(4)
$$\tau _j^{(k)}\{1,0\}^2,\text{ }k=2,\text{ }\mathrm{},\text{ }\eta +1.$$
(5)
The ”floor” symbol $``$ $``$ in (4) denotes the largest integer smaller than the value in it. Although the ordering of the $`n_0`$-qubit factorized states can be arbitrary specified, we have designated them as $`\tau _j^{(1)}`$ in the first subspace for convenience. As mentioned, our purpose is to outlet the superposition of the marked states, namely,

$$|\tau =\frac{1}{\sqrt{M}}\underset{j=1}{\overset{M}{}}|\tau _j=\frac{1}{\sqrt{M}}\underset{j=1}{\overset{M}{}}|\tau _j^{(\eta +1)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)}.$$
(6)
The initial state is prepared in the uniform superposition of all the possible states, given by

$$|s=W_n|\overline{0}=\frac{1}{\sqrt{N}}\underset{y\{0,1\}^n}{}|y$$
(7)
where $`W_n`$ is the $`n`$-qubit Walsh-Hadamard transformation and $`|\overline{0}`$ is the standard state with $`0`$ in every qubit.

By expression (3) we realize that the global, $`2^n`$-dimensional Hilbert space $``$ for the search is the tensor product of $`\eta +1`$ subspaces, expressed as $`=_2^\eta _{n_0}`$, where the subspace $`_{n_0}`$ is $`2^{n_0}`$-dimensional and the subspace $`_2`$ is four-dimensional. The scenario of the algorithm contains $`\eta +1`$ sequential stages. It begins with finding the $`M`$ factorized states $`|\tau _j^{(1)}`$ in the subspace $`_{n_0}`$ in just one searching step. Knowing that Grover’s algorithm promises to search $`M`$ objects with certainty in a database containing $`4M`$ items by applying the Grover operator once, we therefore choose $`n_0=\mathrm{log}_2(4M)`$ for the first subgrouped qubits. Even if $`\mathrm{log}_2(4M)`$ is not an integer, we can also search the objects with certainty in just one searching step by choosing $`n_0=\mathrm{log}_2(4M)`$ and using the technique of phase tunninghsieh and li . In the $`2^{n_0}`$-dimensional subspace $`_{n_0}`$, we apply the subgrouped Grover operator $`G_1`$, given by

$$G_1=(I+(e^{i\varphi }1)|s^{(1)}s^{(1)}|)(I+(e^{i\varphi }1)|\tau ^{(1)}\tau ^{(1)}|),$$
(8)
where

$$|\tau ^{(1)}=\frac{1}{\sqrt{M}}\underset{j=1}{\overset{M}{}}|\tau _j^{(1)},\text{ }|s^{(1)}=\frac{1}{\sqrt{2^{n_0}}}\underset{y\{0,1\}^{n_0}}{}|y,$$
(9)
and the phase angle, according to Hsieh and Lihsieh and li , is given by

$$\varphi =\varphi (M)=2\mathrm{sin}^1(\sqrt{\frac{2^{n_0}}{4M}}).$$
(10)
Notice that if $`\mathrm{log}_2(4M)`$ is an integer, $`\varphi =\pi `$ and $`G_1`$ reduces to the original Grover operator defined in $`_{n_0}`$, while if $`\mathrm{log}_2(4M)`$ is not an integer, the choice $`n_0=\mathrm{log}_2(4M)`$ leads to the smallest deviation of the phase $`\varphi `$ form $`\pi `$. The oracle used in the first subgrouped Grover operator $`G_1`$ should be

$`f_1(y)`$ $`=`$ $`1,\text{ }y=\tau _j^{(1)},\text{ }j=1,\mathrm{},M,`$ (11)
$`f_1(y)`$ $`=`$ $`0,\text{ }y\tau _j^{(1)},\text{ }j=1,\mathrm{},M.`$
In the first stage the searching operator that we apply on the initial state $`|s`$ is

$$S_1=I_2^\eta G_1,$$
(12)
where $`I_2`$ denotes the identity operator in the subspace $`_2`$. The outlet at the end of this stage then becomes

$$|t_1=S_1|s=(\frac{1}{\sqrt{2^{nn_0}}}\underset{y\{0,1\}^{nn_0}}{}|y)|\tau ^{(1)}.$$
(13)
The outlet state $`|t_1`$ can be rewritten, for convenience in the depiction of the next stage,

$$|t_1=(\frac{1}{\sqrt{2^{nn_{02}}}}\underset{y\{0,1\}^{nn_02}}{}|y)|s^{(2)},$$
(14)
where

$$|s^{(2)}=\frac{1}{\sqrt{4M}}\underset{j=1}{\overset{M}{}}(\underset{x\{0,1\}^2}{}|x)|\tau _j^{(1)}$$
(15)
is located in the subspace $`_2_{n_0}`$. Now, observing that in the subspace $`_2_{n_0}`$ the state $`|s^{(2)}`$ is a superposition of $`4M`$ nonvanishing orthonormal states $`|x\tau _j^{(1)}`$, $`x\{0,1\}^2`$, $`j=1,\mathrm{},M`$, and in them the $`M`$ factorized marked states $`|\tau _j^{(2)}\tau _j^{(1)}`$ are embedded. We then once again can apply a second subgrouped Grover operator $`G_2`$ to search these factorized marked states with certainty in just one searching step. The second subgrouped Grover operator $`G_2`$ is defined in the subspace $`_2_{n_0}`$, given by

$$G_2=(I2|s^{(2)}s^{(2)}|)(I2|\tau ^{(2)}\tau ^{(2)}|),$$
(16)
where

$$|\tau ^{(2)}=\frac{1}{\sqrt{M}}\underset{j=1}{\overset{M}{}}|\tau _j^{(2)}\tau _j^{(1)}.$$
(17)
The oracle implemented in $`G_2`$ is of course defined by

$`f_2(y)`$ $`=`$ $`1,\text{ }y=\tau _j^{(2)}\tau _j^{(1)},\text{ }j=1,\mathrm{},M,`$ (18)
$`f_2(y)`$ $`=`$ $`0,\text{ }y\tau _j^{(2)}\tau _j^{(1)},\text{ }j=1,\mathrm{},M.`$
The searching operator that we apply in the second stage then is

$$S_2=I_2^{(\eta 1)}G_2,$$
(19)
and the outlet at the end of the second stage is

$`|t_2`$ $`=`$ $`S_2|t_1=S_2S_1|s`$
$`=`$ $`({\displaystyle \frac{1}{\sqrt{2^{nn_{02}}}}}{\displaystyle \underset{y\{0,1\}^{nn_02}}{}}|y)|\tau ^{(2)}`$
$`=`$ $`({\displaystyle \frac{1}{\sqrt{2^{nn_{04}}}}}{\displaystyle \underset{y\{0,1\}^{nn_04}}{}}|y)|s^{(3)},`$
where the state $`|s^{(3)}`$ is written

$$|s^{(3)}=\frac{1}{\sqrt{4M}}\underset{j=1}{\overset{M}{}}(\underset{y\{0,1\}^2}{}|x)|\tau _j^{(2)}\tau _j^{(1)},$$
(21)
which locates in the third subspace $`_2^2_{n_0}`$ and is also a superposition of the $`4M`$ nonvanishing orthonormal states $`|x\tau _j^{(2)}\tau _j^{(1)}`$, $`x\{0,1\}^2`$, $`j=1,\mathrm{},M`$, with the $`M`$ factorized marked states $`|\tau _j^{(3)}\tau _j^{(2)}\tau _j^{(1)}`$ embedded in them. So the scenario obvisouly can proceed staightfoward to the final statge following the same way as depicted from the first to the second stage. That is, in a next stage we add two more qubits to the former subgrouped qubits, so a next subspace is the one enlarged by multiplying the dimensions of the former subspace by four, and apply a searching operator like those given by (12) and (19) and a corresponding subgrouped oracle like those expressed in (11) and (18).

As a result, we can thus derive the general formula for the present searching algorithm. Excluding the first stage, where as we can apply expressions (8)-(13), in the $`k`$-th stage, $`2k\eta +1`$, the searching operator is

$$S_k=I_2^{(\eta k+1)}G_k,$$
(22)
and the $`k`$-th subgrouped Grover operator $`G_k`$ is defined in the subspace $`_2^{(k1)}_{n_0}`$, given by

$$G_k=(I2|s^{(k)}s^{(k)}|)(I2|\tau ^{(k)}\tau ^{(k)}|),$$
(23)
where

$$|\tau ^{(k)}=\frac{1}{\sqrt{M}}\underset{j=1}{\overset{M}{}}|\tau _j^{(k)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)},$$
(24)
and

$$|s^{(k)}=\frac{1}{\sqrt{4M}}\underset{j=1}{\overset{M}{}}(\underset{y\{0,1\}^2}{}|x)|\tau _j^{(k1)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)}.$$
(25)
The corresponding subgrouped oracle required in the $`k`$-th stage is defined by

$`f_k(y)`$ $`=`$ $`1,\text{ }y=\tau _j^{(k)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)},\text{ }j=1,\mathrm{},M,`$ (26)
$`f_k(y)`$ $`=`$ $`0,\text{ }y\tau _j^{(k)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)},\text{ }j=1,\mathrm{},M.`$
The outlet state at the end of the $`k`$-th stage then is

$`|t_k`$ $`=`$ $`S_k|t_{k1}=S_k\mathrm{}S_2S_1|s`$
$`=`$ $`({\displaystyle \frac{1}{\sqrt{2^{nn_02(k1)}}}}{\displaystyle \underset{y\{0,1\}^{nn_02(k1)}}{}}|y)|\tau ^{(k)}.`$
As $`k=\eta +1`$, i.e., as we arrive at the final stage, we will reach the final state at the end,

$`|t_{\eta +1}`$ $`=`$ $`S_{\eta +1}|t_\eta =S_{\eta +1}\mathrm{}S_2S_1|s`$
$`=`$ $`|\tau ^{(\eta +1)}={\displaystyle \frac{1}{\sqrt{M}}}{\displaystyle \underset{j=1}{\overset{M}{}}}|\tau _j^{(\eta +1)}\mathrm{}\tau _j^{(2)}\tau _j^{(1)}`$
$`=`$ $`|\tau ,`$
which is exactly the object that we search for. It should be mentioned that a projected / measurement operator, although not shown, has been combined within the searching operator used in each of the stages to make sure the eliminated states do not take part in the subsequent stages of the algorithm. To summarize, totally $`\eta +1`$ queries are required in our subgrouped searching algorithm since in each stage of the scenario only one query is utilized. The number $`\eta +1`$, being equal to $`(nn_0+2)/2`$, is in fact smaller than $`n/2`$ for $`M>1`$. As $`M=1`$, our algorithm will reduce to the factorized algorithm presented by Patel. A circuit representation of the present algorithm is shown in Fig. 1 to schematically depict the whole scenario.

In this work, we have shown a subgrouped algorithm using only $`O(\mathrm{log}_4N)`$ queries in searching $`M`$ objects out of an unsorted database of $`N`$ items, if in advance a set of subgrouped oracles is available. The subgrouped oracles required in the algorithm are expressed in (11) and (26), while the searching operator are defined by (12) and (22), with the aides of (8)-(10) and (23)-(25), respectively. The first searching operator $`S_1`$ is more particular than the others because in it the technique of phase tunning is necessary for the action of the subgrouped Grover operator $`G_1`$ especially when $`\mathrm{log}_2(4M)`$ is not an integer. Using this subgrouped searching algorithm, one will require $`O(\mathrm{log}_4N)`$ queries in the search of multi-object. The resulting state at the end of the algorithm is the superposition of the $`M`$ marked states, which is an unique feature of the quantum algorithm and can not be performed by any classical algorithm. Difficulties, however, may arise from how to physically implement the black boxes associated with the subgrouped oracles, but these are not of concern in the complexity analysis of this work and should be studied in a furture work.


text:
# Connectedness applied to closure spaces and state property systems11footnote 1Published as Aerts, D., Deses, D. and Van der Voorde, A., 2001, “Connectedness applied to closure spaces and state property systems”, Journal of Electrical Engineering, 52, 18 - 21.

## 1 Introduction

In  a physical entity is represented by a mathematical model called a state property system. This model contains a complete lattice of properties of the physical entity. In  it is shown that the lattice can viewed as the lattice of closed sets of a closure space. We introduce the concept of classical property of the entity, and show that these correspond exactly to the clopen (open and closed) subsets of the associated closure space. Using the concept of connectedness for closure space we decompose the state property system into smaller ones which will be ‘completely quantum mechanical’ (no classical properties) and another one which will be ‘totally classical’. Finally we introduce a way to extract the classical properties of the entity, and interpret this in term of the closure space. Let us first introduce the basic definitions and concepts.

###### Definition 1.

A triple $`(\mathrm{\Sigma },,\xi )`$ is called a state property system if $`\mathrm{\Sigma }`$ is a set, $``$ is a complete lattice and $`\xi :\mathrm{\Sigma }𝒫()`$ is a function such that for $`p\mathrm{\Sigma }`$, $`0`$ the minimal element of $``$ and $`(a_i)_i`$, we have:

$$0\xi (p)\text{ and }ia_i\xi (p)_ia_i\xi (p)$$

and for $`a,b`$ we have: $`a<br\mathrm{\Sigma }:a\xi (r)\text{ then }b\xi (r)`$.
If $`(\mathrm{\Sigma },,\xi )`$ is a state property system then its Cartan map is the mapping $`\kappa :𝒫(\mathrm{\Sigma })`$ defined by :

$$\kappa :𝒫(\mathrm{\Sigma }):a\kappa (a)=\{p\mathrm{\Sigma }|a\xi (p)\}$$

The physical interpretation of this mathematical structure (introduced in ) is the following. Considering an entity $`S`$, the set $`\mathrm{\Sigma }`$ consists of states of $`S`$ while the set $``$ consists of properties of $`S`$. These two sets are linked by means of a function $`\xi :\mathrm{\Sigma }𝒫()`$ which maps a state $`p`$ to the set $`\xi (p)`$ of all properties that are actual in state $`p`$. This means that the statement ‘a state $`p`$ makes the property $`a`$ actual’ is mathematically expressed by the formula: $`a\xi (p)`$.

###### Definition 2.

A closure space $`(X,)`$ consists of a set $`X`$ and a family of subsets $`𝒫(X)`$ satisfying the following two conditions:

$$\mathrm{}\text{ and }(F_i)_i_iF_i$$

The closure operator corresponding to the closure space $`(X,)`$ is defined as

$$cl:𝒫(X)𝒫(X):A\{F|AF\}$$

The following theorem shows how we can associate with each state property system a closure space and vice versa.

###### Proposition 1.

If $`(\mathrm{\Sigma },,\xi )`$ is a state property system then $`F(\mathrm{\Sigma },,\xi )=(\mathrm{\Sigma },\kappa ())`$ is a closure space. Conversely, if $`(\mathrm{\Sigma },)`$ is a closure space then $`G(\mathrm{\Sigma },)=(\mathrm{\Sigma },,\xi ^{})`$, where $`\xi ^{}(x)=\{F|xF\}`$, is a state property system.

This is a consequence of the fact that there is an categorical equivalence between state property systems and closure spaces, as proven in  .

## 2 Super selection rules and classical properties

Superposition states in quantum mechanics are those states that do not exist in classical physics and hence their appearance is one of the important quantum characteristics. This concept can be traced back within this general setting, by introducing the idea of ‘super selection rule’. Two properties are separated by a super selection rule iff there do not exist ‘superposition states’ related to these two properties.

###### Definition 3.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$. For $`a,b`$ we say that $`a`$ and $`b`$ are separated by a super selection rule, and denote $`a`$ ssr $`b`$, iff for $`p\mathrm{\Sigma }`$ we have: $`ab\xi (p)a\xi (p)\mathrm{or}b\xi (p)`$

###### Lemma 1.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$ and its corresponding closure system $`=\kappa ()`$. For $`a,b`$ we have:

$$assrb\kappa (ab)=\kappa (a)\kappa (b)\kappa (a)\kappa (b)$$

###### Proof.

This is an easy verification. ∎

We are ready now to introduce the concept of a ‘classical property’.

###### Definition 4.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$. We say that a property $`a`$ is a ‘classical property’, if there exists a property $`a^c`$ such that $`aa^c=I`$, $`aa^c=0`$ and $`a`$ ssr $`a^c`$.

Remark that for every state property system $`(\mathrm{\Sigma },,\xi )`$ the properties $`0`$ and $`I`$ are classical properties. Note also that if $`a`$ is a classical property, we have for $`p\mathrm{\Sigma }`$ that $`a\xi (p)a^c\xi (p)`$ and $`a\xi (p)a^c\xi (p)`$. This follows immediately from the definition of a classical property. From the previous definition and lemma 1 one can prove the following lemma.

###### Lemma 2.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$. If $`a`$ is a classical property, then $`a^c`$ is unique and is a classical property. We call it the complement of $`a`$. Further the following is satisfied:

$$(a^c)^c=a,a<bb^c<a^c,\kappa (a^c)=\kappa (a)^C$$

###### Definition 5.

A closure space $`(X,)`$ is called connected if the only clopen (i.e. closed and open) sets are $`\mathrm{}`$ and $`X`$.

We can show now that the subsets that make closure systems disconnected are exactly the subsets corresponding to classical properties.

###### Proposition 2.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$ and its corresponding closure space $`(\mathrm{\Sigma },\kappa ())`$. For $`a`$ we have: $`a\text{ is classical }\kappa (a)`$ is clopen.

###### Proof.

From the previous lemmas it follows that if $`a`$ is classical, then $`\kappa (a)`$ is clopen. So now consider a clopen subset $`\kappa (a)`$ of $`\mathrm{\Sigma }`$. This means that $`\kappa (a)^C`$ is closed, and hence that there exists a property $`b`$ such that $`\kappa (b)=\kappa (a)^C`$. We clearly have $`ab=0`$ since there exists no state $`p\mathrm{\Sigma }`$ such that $`p\kappa (a)`$ and $`p\kappa (b)`$. Since $`\mathrm{\Sigma }=\kappa (a)\kappa (b)`$ we have $`ab=I`$. Further we have that for an arbitrary state $`p\mathrm{\Sigma }`$ we have $`a\xi (p)`$ or $`b\xi (p)`$ which shows that $`a`$ ssr $`b`$. This proves that $`b=a^c`$ and that $`a`$ is classical. ∎

This means that the classical properties correspond exactly to the clopen subsets of the closure system.

###### Corollary 1.

Let $`(\mathrm{\Sigma },,\xi )`$ be a state propery system. T.F.A.E.
(1) The properties $`0`$ and $`I`$ are the only classical ones.
(2) $`F(\mathrm{\Sigma },,\xi )=(\mathrm{\Sigma },\kappa ())`$ is a connected closure space.

###### Definition 6.

A state property system $`(\mathrm{\Sigma },,\xi )`$ is called a ‘pure nonclassical state property system’ if the properties $`0`$ and $`I`$ are the only classical properties.

###### Proposition 3.

Let $`(\mathrm{\Sigma },)`$ be a closure space. T.F.A.E.
(1) $`(\mathrm{\Sigma },)`$ is a connected closure space.
(2) $`G(\mathrm{\Sigma },)=(\mathrm{\Sigma },,\xi )`$ is a pure nonclassical state property system.

###### Proof.

Let $`(\mathrm{\Sigma },)`$ be a connected state property system. Then $`\mathrm{}`$ and $`\mathrm{\Sigma }`$ are the only clopen sets in $`(\mathrm{\Sigma },)`$. Since the Cartan map associated to $`\xi `$ is given by $`\kappa :𝒫(\mathrm{\Sigma }):FF`$, we have $`\kappa (\mathrm{})=\mathrm{}`$ and $`\kappa (\mathrm{\Sigma })=\mathrm{\Sigma }`$. Applying proposition 2, we find that $`\mathrm{}`$ and $`\mathrm{\Sigma }`$ are the only classical properties of $``$. Conversely, let $`G(\mathrm{\Sigma },)=(\mathrm{\Sigma },,\xi )`$ be a pure nonclassical state property system. Then by corollary 1, $`(\mathrm{\Sigma },)=FG(\mathrm{\Sigma },)`$ is a connected closure space. ∎

## 3 Decomposition theorem

As for topological spaces, every closure space can be decomposed uniquely into connected components. In the following we say that, for a closure space $`(X,)`$, a subset $`AX`$ is connected if the induced subspace is connected. It can be shown that the union of any family of connected subsets having at least one point in common is also connected. So the component of an element $`xX`$ defined by $`K_{\mathrm{𝐂𝐥𝐬}}(x)=\{AX|xA,A\text{ connected }\}`$ is connected and therefore called the connection component of $`x`$. Moreover, it is a maximal connected set in $`X`$ in the sense that there is no connected subset of $`X`$ which properly contains $`K_{\mathrm{𝐂𝐥𝐬}}(x)`$. From this it follows that for closure spaces $`(X,)`$ the set of all distinct connection components in $`X`$ form a partition of $`X`$. In the following we will decompose state property systems similarly into different components.

###### Proposition 4.

Let $`(\mathrm{\Sigma },,\xi )`$ be a state property system and let $`(\mathrm{\Sigma },\kappa ())`$ be the corresponding closure space. Consider the equivalence relation on $`\mathrm{\Sigma }`$ given by: $`pqK_{\mathrm{𝐂𝐥𝐬}}(p)=K_{\mathrm{𝐂𝐥𝐬}}(q)`$ with equivalence classes $`\mathrm{\Omega }=\{\omega (p)|p\mathrm{\Sigma }\}`$. If $`\omega \mathrm{\Omega }`$ we define:

$`\mathrm{\Sigma }_\omega `$ $`=`$ $`\omega =\{p\mathrm{\Sigma }|\omega (p)=\omega \}`$
$`s(\omega )`$ $`=`$ $`s(\omega (p))=a,\text{ such that }\kappa (a)=\omega (p)`$
$`_\omega `$ $`=`$ $`[0,s(\omega )]=\{a|0as(\omega )\}`$
$`\xi _\omega `$ $`:`$ $`\mathrm{\Sigma }_\omega 𝒫(_\omega ):p\xi (p)_\omega `$
then $`(\mathrm{\Sigma }_\omega ,_\omega ,\xi _\omega )`$ is a state property system.

###### Proof.

Since $`_\omega `$ is a sublattice (segment) of $``$, it is a complete lattice with maximal element $`I_\omega =s(\omega )`$ and minimal element $`0_\omega =0`$. Let $`p\mathrm{\Sigma }_\omega `$. Then $`0\xi (p)`$. So $`0\xi (p)_\omega =\xi _\omega (p)`$. If $`a_i\xi _\omega (p),i`$, then $`a_i_\omega `$ and $`a_i\xi (p),i`$. Hence $`a_i_\omega \xi (p)=\xi _\omega (p)`$. Finally, let $`a,b_\omega `$ with $`a<_\omega b`$ and let $`r\mathrm{\Sigma }_\omega `$. If $`a\xi _\omega (r)`$, then $`a_\omega `$ and $`a\xi (r)`$, thus $`b_\omega `$ and $`b\xi (r)`$. So $`b\xi _\omega (r)`$. Conversely, if $`a,b_\omega `$ and $`r\mathrm{\Sigma }_\omega :a\xi _\omega (r)b\xi _\omega (r)`$ then we consider a $`q`$ such that $`a\xi (q)`$ ($`q`$ must be in $`\mathrm{\Sigma }_\omega `$ by definition of $`_\omega `$). Then $`a\xi _\omega (q)`$ implies that $`b\xi _\omega (q)`$. So $`b\xi (q)`$ and $`a<b`$. Thus $`a<_\omega b`$. ∎

Moreover we can show that the above introduced state property systems $`(\mathrm{\Sigma }_\omega ,_\omega ,\xi _\omega )`$ have no proper classical properties, and hence are pure nonclassical state property systems.

###### Proposition 5.

Let $`(\mathrm{\Sigma },,\xi )`$ be a state property system. If $`\omega \mathrm{\Omega }`$, then $`(\mathrm{\Sigma }_\omega ,_\omega ,\xi _\omega )`$ is a pure nonclassical state property system.

###### Proof.

If $`a`$ is classical element of $`_\omega `$, then $`\kappa (a)`$ must be a clopen set of the associated closure space $`(\mathrm{\Sigma }_\omega ,\kappa (_\omega ))`$ which is a connected subspace of $`(\mathrm{\Sigma },\kappa ())`$. Hence there are no proper classical elements of $`_\omega `$. ∎

###### Proposition 6.

Let $`(\mathrm{\Sigma },,\xi )`$ be a state property system. If we introduce the following :

$`\mathrm{\Omega }`$ $`=`$ $`\{\omega (p)|p\mathrm{\Sigma }\}`$
$`𝒞`$ $`=`$ $`\{s(\omega _i)|\omega _i\mathrm{\Omega }\}`$
$`\eta `$ $`:`$ $`\mathrm{\Omega }𝒫(𝒞):\omega =\omega (p)\xi (p)𝒞`$
then $`(\mathrm{\Omega },𝒞,\xi )`$ is a state property system.

###### Proof.

First we remark that $`\eta `$ is well defined because if $`\omega (p)=\omega (q)`$, then $`\xi (p)𝒞=\xi (q)𝒞`$. Indeed, if $`s(\omega _i)\xi (p)`$ then $`p\kappa (s(\omega _i))=cl(\omega _i)`$ in the corresponding closure space $`(\mathrm{\Sigma },\kappa ())`$. Since $`cl(\omega _i)`$ is not connected we have that $`K_{\mathrm{𝐂𝐥𝐬}}(p)=\omega (p)=\omega (q)cl(\omega _i)`$ so $`qcl(\omega _i)=\kappa (s(\omega _i))`$ and $`s(\omega _i)\xi (q)`$. Now, since $`𝒞`$ is a sublattice of $``$ it is a complete lattice with $`1_𝒞=1`$ and $`0_𝒞=0`$. By definition $`𝒞`$ is generated by its atoms $`\{s(\omega )|\omega \mathrm{\Omega }\}`$. Clearly $`0\eta (\omega (p))`$ because $`0\xi (p)`$. If $`a_i\eta (\omega (p))=\xi (p)𝒞,i`$, then $`a_i\xi (p)𝒞=\eta (\omega (p))`$. Finally, let $`a,b𝒞`$ with $`a<_𝒞b`$. Let $`\omega (p)\mathrm{\Omega }`$ with $`a\eta (\omega (p))`$. Thus $`a\xi (p)`$. $`a<_𝒞b`$ implies $`a<b`$. So we have $`b\xi (p)𝒞=\eta (\omega (p))`$. Conversely, let $`a,b𝒞`$ and assume that $`p\mathrm{\Sigma }:a\eta (\omega (p))b\eta (\omega (p))`$. Then we have $`p\mathrm{\Sigma }:a\xi (p)b\xi (p)`$. Thus $`a<b`$ and $`a<_𝒞b`$. ∎

###### Proposition 7.

$`(\mathrm{\Omega },𝒞,\eta )`$ is a totally classical state property system, in the sense that the only quantum segments (i.e. segments with no proper classical elements) are trivial, i.e. $`\{0,s(\omega )\}`$.

###### Proof.

Suppose $`[0,a]`$ is a quantumsegment of $`𝒞`$, then in the corresponding closure space $`(\mathrm{\Sigma },\kappa ())`$ the subset $`\kappa (a)`$ is connected hence $`\kappa (a)\omega `$ for some $`\omega \mathrm{\Omega }`$, hence $`a<s(\omega )`$. Since $`s(\omega )`$ is an atom, $`a=s(\omega )`$. Thus $`[0,a]=\{0,s(\omega )\}`$. ∎

###### Corollary 2.

The closure space associated with $`(\mathrm{\Omega },𝒞,\eta )`$ is a totally disconnected closure space.

Summarizing the previous results we get:

###### Proposition 8.

Any state property system $`(\mathrm{\Sigma },,\xi )`$ can be decomposed into: a number of pure nonclassical state property systems $`(\mathrm{\Sigma }_\omega ,_\omega ,\xi _\omega ),\omega \mathrm{\Omega }`$ and a totally classical state property system $`(\mathrm{\Omega },𝒞,\eta )`$

## 4 The classical part of a state property system

In this section we want to show how it is possible to extract the classical part of a state property system. First of all we have to define the classical property lattice related to the entity $`S`$ that is described by the state property system $`(\mathrm{\Sigma },,\xi )`$.

###### Definition 7 (Classical property lattice).

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$. We call $`𝒞^{}=\{_ia_i|a_i\mathrm{is}\mathrm{a}\mathrm{classical}\mathrm{property}\}`$ the classical property lattice corresponding to the state property system $`(\mathrm{\Sigma },,\xi )`$.

###### Proposition 9.

$`𝒞^{}`$ is a complete lattice with the partial order relation and infimum inherited from $``$ and the supremum defined as follows: for $`a_i𝒞^{}`$, $`_ia_i=_{b𝒞^{},a_ibi}b`$.

Remark that the supremum in the lattice $`𝒞^{}`$ is not the one inherited from $``$.

###### Proposition 10.

Consider a state property system $`(\mathrm{\Sigma },,\xi )`$. Let $`\xi ^{}(q)=\xi (q)𝒞^{}`$ for $`q\mathrm{\Sigma }`$, then $`(\mathrm{\Sigma },𝒞^{},\xi ^{})`$ is a state property system which we shall refer to as the classical part of $`(\mathrm{\Sigma },,\xi )`$.

###### Proof.

Clearly $`0\xi ^{}(p)`$ for $`p\mathrm{\Sigma }`$. Consider $`a_i\xi ^{}(p)i`$. Then $`a_i\xi (p)𝒞^{}i`$, from which follows that $`_ia_i\xi (p)𝒞^{}`$ and hence $`_ia_i\xi ^{}(p)`$. Consider $`a,b𝒞^{}`$. Let us suppose that $`ab`$ and consider $`r\mathrm{\Sigma }`$ such that $`a\xi ^{}(r)`$. This means that $`a\xi (r)𝒞^{}`$. From this follows that $`b\xi (r)𝒞^{}`$ and hence $`b\xi ^{}(r)`$. On the other hand let us suppose that $`r\mathrm{\Sigma }:a\xi ^{}(r)\mathrm{then}b\xi ^{}(r)`$. Since $`a,b𝒞^{}`$, this also means that $`r\mathrm{\Sigma }:a\xi (r)\mathrm{then}b\xi (r)`$. From this follows that $`ab`$. ∎

Since $`(\mathrm{\Sigma },𝒞^{},\xi ^{})`$ is a state property system, it has a corresponding closure space $`(\mathrm{\Sigma },\kappa (𝒞^{}))`$. In order to check some property of this space we introduce the following concepts.

###### Definition 8.

Let $`(X,)`$ be a closure space and $``$. $``$ is called a base of $`(X,)`$ iff $`F:B_i:F=B_i`$. $`(X,)`$ is called weakly zero-dimensional iff there is a base consisting of clopen sets.

###### Proposition 11.

The closure space $`(\mathrm{\Sigma },\kappa (𝒞^{}))`$ corresponding to the state property system $`(\mathrm{\Sigma },𝒞^{},\xi ^{})`$ is weakly zero-dimensional.

###### Proof.

To see this recall that $`a`$ is classical iff $`\kappa (a)`$ is clopen in $`(\mathrm{\Sigma },\kappa ())`$, hence $`\kappa (𝒞^{})`$ is a family of closed sets on $`\mathrm{\Sigma }`$ which consists of all intersections of the clopen sets of $`(\mathrm{\Sigma },\kappa ())`$. ∎

In general the lattice of $`(\mathrm{\Sigma },𝒞^{},\xi ^{})`$ does not need to be atomistic, hence it is different from the totally classical state property system $`(\mathrm{\Omega },𝒞,\eta )`$ associated with $`(\mathrm{\Sigma },,\xi )`$.

Didier Deses has presented this subject at SCAM 2001. He is a research assistant of the fund for scientific research flanders and Phd. student at the free university of Brussels, his supervisor is Prof. Eva Colebunders. His main interest lies in General and Categorical Topology.


text:
# On the Role of External Constraints in a Spatially Extended Evolutionary Prisoner’s Dilemma Game

## Introduction

The successful applications of game theory in the area of economics and political decisions initiated its increasing development after the second world war vNM . Originally, the game theory is devoted to find the optimal strategy for a given game between two intelligent players. The straightforward developments involve the generalization toward the iterated games of $`n`$ players with assuming local interactions among the spatially distributed players. The spatial evolutionary prisoner’s dilemma games (SEPDG) has attracted a particular attention because of its applicability in the human and behavior sciences as well as in biology HS ; sigmund ; msmith ; axelrod ; weibull . Nowadays the prisoner’s dilemma game is considered to be the metaphor for studying the emergence of cooperation among selfish individuals. The emerging cooperation appears to be crucial at many transitions in evolution MSSz . The first numerical investigations have shown that the cooperation can be maintained by very simple strategies in the iterated games axelrod . Very recently it is demonstrated that the players can be as simple as bacteriophages (viruses that infect bacteria)TC ; NS99 .

In these systems the players wish to maximize their individual income coming from games with other players. The prisoner’s dilemma game is a simple version of the two-player matrix games where the players’ income depend on their simultaneous choice between two options. Following the widely accepted expressions each player can choose defection or cooperation with the other player. The defector reaches the highest payoff $`t`$ (called temptation to defect) against the cooperator, which receives then the lowest reward $`s`$ (called sucker’s payoff). For mutual cooperation \[defection\] each player receives the same payoff $`r`$ (reward for mutual cooperation) \[$`p`$ (punishment)\]. The game is symmetric in the sense that player’s income is independent of the player itself, it depends only on their choice. The mentioned payoff values satisfy the inequalities $`t>r>p>s`$ and $`2r>t+s`$. These assumptions provide the largest total payoff for the mutual cooperators. Comparing to this situation the defector reaches extra income against the cooperator whose loss exceeds the defector’s benefit. Consequently, the choice of defection can be interpreted as an exploiting behavior. These are the main features for which the prisoner’s dilemma games are used to study the emergence of mutual cooperation, altruism and ethic norms among selfish individuals axelrod ; NMS .

The rational players should defect as this choice provides the larger income, independently of the partner’s decision. However, this situation creates a dilemma for intelligent players as mutual cooperation would result in higher income for each of them than mutual defection does.

In the iterated round-robin prisoner’s dilemma games we can introduce some simple evolutionary processes without assuming intelligent players (who are capable to find the best strategy if it exists). These games are started from an initial set of strategies, which defines the player’s decision in the knowledge of their previous choices. The evolutionary process is devoted to model the Darwinian selection principle among $`n`$ (selfish) players whose total income comes from $`n1`$ games within a given round. In the simplest evolutionary models the worst player adopts the winner’s strategy round by round.

The numerical simulations have demonstrated the crucial role of the so-called “tit for tat” strategy in the emergence of mutual cooperation. Despite of its simplicity the “tit for tat” strategy won the computer tournaments conducted by Axelrod axelrod . The “tit for tat” strategy cooperates in the first step and then always repeats his co-player’s previous decision. This strategy cooperates forever with all the other so-called nice strategies which never defect first. Furthermore, its defection and cooperation can be interpreted as a punishment and forgiveness when reacting to the previous decision for other strategies. The most remarkable feature of this strategy is that it is capable to sustain the mutual cooperation among themselves in the presence of defectors.

Early numerical investigations have also indicated the importance of local interactions because it favors the formation of cooperating colonies. In the simplest models the players are distributed on a lattice and the interaction (the games between two player as well as the strategy adoption) is limited to a given neighborhood. Evidently, the short range interactions enhance the role of fluctuations at the same time. These models can be well investigated by sophisticated methods of non-equilibrium statistical physics.

For the numerical investigation of the spatial effects Nowak and May NM  have introduced an SEPDG model, which is equivalent to a two-state cellular automaton. Each lattice site can be in one of the two states $`C`$ and $`D`$, representing the two simple strategy “always cooperate” and “always defect” respectively. The income for a given player (site) comes from games with its neighbors (and also with itself in some version of the model). According to the cellular automaton rule the players modify their strategy simultaneously in discrete time steps. Namely, each player adopts the best strategy found in its neighborhood. The step by step visualization of the strategy distribution on a two-dimensional lattice exhibits different spatio-temporal patterns (homogeneous and coexisting strategies, transitions between these states, competing interfacial invasions, etc.) depending on the payoff matrix, which is characterized by a single parameter. In these models the randomness is restricted to the initial states. In a subsequent work Nowak et al. NBM  have extended the former models by allowing irrational strategy adoptions with some probability. The simulations indicated that the randomness favors the spreading of $`D`$ strategies. These results have initiated systematic numerical investigations of many stochastic cellular automata KD96 ; KD98 ; KDK ; CO .

The study of spatio-temporal patterns observed in nature, however, requires continuous time description HG ; NM . Moreover, it is difficult to analyze the above mentioned stochastic cellular automata in the framework of generalized mean-field approximation, which is often used in non-equilibrium physics. To reduce the technical difficulties Szabó and Tőke have suggested a simplified dynamics epdg2s . The systematic investigations of this model have justified that when tuning the model parameters the stationary state undergoes two consecutive phase transitions which belong to the directed percolation (DP) universality classepdg2s ; CO . Very recently this SEPDG model has been extended by allowing three strategies for the players epdg3s . In the present work this three-strategy model will be compared with its simplified version. During the model descriptions and discussion, our attention will be focused on the motivations, the elementary processes and their consequences as well as on the universal features relating the SEPDGs to the area of complex systems.

## Spatial evolutionary model with three strategies

In the present spatial evolutionary prisoner’s dilemma game the players are located on the sites $`𝐱=(i,j)`$ of a square lattice, where $`i,j=1,\mathrm{},L`$. To avoid the undesired boundary effects we assume periodic boundary conditions. Each player follows one of the three strategies: $`D`$ defects always; $`C`$ cooperates unconditionally; $`T`$ accommodating to the partner’s strategy chooses defection against $`D`$ and cooperation with $`C`$ and $`T`$. In fact the name $`T`$ refers to the strategy “tit for tat” which first cooperates and later repeats the partner’s previous decision. Consequently, after the first step the decisions of these two strategies are equivalent against $`C`$, $`D`$ and themselves. The consequences of the different first decisions become irrelevant if the strategy changes (defined below) are rare comparing to the frequency of games. At the site $`𝐱`$ the player’s strategy is denoted by a three-component unit vector whose possible values are

$$𝐬(𝐱)=\left(\begin{array}{c}1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right)$$
(1)
corresponding to the $`D`$, $`C`$, and $`T`$ strategies respectively. At a given time the state of the whole system is described by the variables $`𝐬(𝐱)`$.

For each player the total payoff comes from the games with its four nearest neighbors. Using the above formalism the total payoff $`m(𝐱)`$ for the player at site $`𝐱`$ is given as

$$m(𝐱)=\underset{\delta 𝐱}{}𝐬^{}(𝐱)\mathrm{𝐌𝐬}(𝐱+\delta 𝐱)$$
(2)
where $`𝐬^{}(𝐱)`$ is the transpose matrix of $`𝐬(𝐱)`$ and the summation runs over the four nearest neighbors ($`\delta 𝐱`$). Accepting the simplified payoff matrix suggested by Nowak and May NM  $`𝐌`$ is given by the following expression:

$$𝐌=\left(\begin{array}{ccc}0& b& 0\\ 0& 1& 1\\ 0& 1& 1\end{array}\right)$$
(3)
where the only free parameter $`b`$ ($`1<b<2`$) measures the temptation to defect. In the above mentioned notation the present payoff matrix corresponds to the choices: $`r=1`$, $`p=0`$, $`t=b`$, and $`s=\epsilon `$ in the limit $`\epsilon 0`$.

To model the Darwinian selection rule the players are allowed to modify their strategy. In the simplest case the system evolution is governed by random sequential updates. It means that a randomly chosen player (e.g. at site $`𝐱`$) adopts one of its neighboring strategy, $`𝐬(𝐱+\delta 𝐱)`$, if $`m(𝐱+\delta 𝐱)>m(𝐱)`$ and this elementary process is iterated many times.

Here it is worth mentioning that a state consisting only of $`C`$ and $`T`$ strategies leads to a uniform payoff distribution \[$`m(𝐱)=4`$\] and the above dynamics leaves this state unchanged. An example of a more complicated situation is given in Figure 1. The payoffs associated with the three different strategies, $`D`$, $`C`$ and $`T`$, are explicitly given.

The reader can easily check that inside a $`D`$ region the defectors receive zero payoffs. The same is true for a solitary $`T`$ surrounded by defectors. In the absence of $`C`$ strategies, however, two (or more) neighboring $`T`$ strategies will invade the $`D`$ territories because their mutual cooperation gives them some incomes, while the defectors’ payoff remain zero.

In the presence of $`C`$ strategies, however, the above situation becomes quite different as the exploitation provides large incomes for the defectors. As a result, the defectors can invade the neighboring $`C`$ or $`T`$ sites for some configurations. This process dominates the time evolution for small $`T`$ and large $`C`$ concentrations as illustrated in a ternary diagram (see Figure 2). Note that the trajectories are two dimensional projections of a many dimensional space. Accordingly, there can be crossing of trajectories. As the average defector’s payoff decreases with the $`C`$ concentration, sooner or later the $`TD`$ invasion processes will govern the system evolution and, finally, all the $`D`$ strategies extinct. Figure 2 shows clearly that the ratio of $`C`$ and $`T`$ strategies in the final (frozen) state depends on the initial state.

It is emphasized that in the absence of $`T`$ strategies the defectors will dominate the present system in the final state. It is not evident as in Figure 1 one can find many $`D`$-$`C`$ pairs where $`C`$ beats $`D`$. In general, these pairs are located along the horizontal and vertical straight fronts separating the $`D`$ and $`C`$ domains. The random sequential invasions, however, makes the smooth fronts irregular and this situation generally prefers the $`DC`$ invasion to the opposite one. As a result, the “sharp” $`D`$ fronts cut the $`C`$’s domains into small pieces and finally all the $`C`$s will be eliminated.

The reader can easily recognize that in most of the $`C`$-$`D`$ (or $`T`$-$`D`$) confrontations the direction of dominance is not affected by the value of $`b`$ within the prescribed region ($`1<b<2`$). The systematic analysis shows that there is only one situation when the value of $`b`$ becomes important. Namely, if a defector has a payoff of $`2b`$ while its $`C`$ (or $`T`$) counterpart has 3. In this case, $`D`$ wins if $`b>3/2`$, otherwise $`D`$ will be invaded. These types of elementary processes, however, do not modify the system behavior drastically epdg3s , therefore the subsequent investigations will be focused on the case $`b>3/2`$.

The above dynamical rules introduce some noises (irrational choices) in the system evolution. Now an additional (superimposing) noisy term is introduced by allowing the appearance of mutants with probability $`P`$. In fact, the effect of two different external constraints (mutation mechanisms) will be studied in models A and B as a function of $`P`$.

### Model A

In the first model, the above evolutionary rule is modified as follows. Each randomly chosen player adopts with probability $`P`$ a randomly chosen strategy among the two other strategies. With probability $`1P`$ it follows the old rule.

This model can describe the behavior of those biological and economical problems where the appearance of mutants cannot be neglected msmith ; HS . The main feature of this model is that this mutation mechanism does not allow the extinction of any strategy.

### Model B

In the second model the mutation mechanism is restricted to the adoption of $`C`$ strategies epdg3s . In other words, the randomly chosen player adopts the $`C`$ strategy with probability $`P`$, otherwise it adopts one of its neighboring strategy if this neighbor has higher income. Note that in this case the extinction of the $`D`$ and/or $`T`$ strategies is permitted.

Model B is devoted to describe the effect of an external constraint which enforces the cooperative behavior naively by supporting an unconditional cooperation. Such a phenomenon can be observed in human societies in which any kind of social pressure enforces the $`D`$ and $`T`$ players to choose the $`C`$ strategy. Furthermore, a $`T`$ player surrounded by only cooperating strategies ($`C`$ or $`T`$) is motivated to adopt the $`C`$ strategy also because of its convenience. In fact, playing $`C`$ is simpler than playing $`T`$, which requires the knowledge of the previous decision of your neighbors.

## Mean-field approximation

In the classical mean-field approximation the system is described by the strategy concentrations which satisfy the normalization condition $`c_D(t)+c_C(t)+c_T(t)=1`$. In this approach the average payoffs are given as:

$`m_D`$ $`=`$ $`bc_C,`$
$`m_C`$ $`=`$ $`c_C+c_T,`$ (4)
$`m_T`$ $`=`$ $`c_C+c_T.`$
For model A, the time dependent concentrations satisfy the following equations of motion:

$`\dot{c}_D`$ $`=`$ $`{\displaystyle \frac{P}{2}}(c_C+c_T2c_D)(1P)c_D(c_C+c_T),`$
$`\dot{c}_C`$ $`=`$ $`{\displaystyle \frac{P}{2}}(c_T+c_D2c_C)\pm (1P)c_Dc_C,`$ (5)
$`\dot{c}_T`$ $`=`$ $`{\displaystyle \frac{P}{2}}(c_D+c_C2c_T)\pm (1P)c_Dc_T,`$
where the upper (lower) signs are valid if $`m_D<m_C=m_T`$ ($`m_D>m_C=m_T`$). In these expressions the first terms describe the effect of external constraint, the second terms come from the Darwinian selection mechanism.

After some algebraic manipulations one can easily get the following stationary solution (for $`P<1`$):

$`c_D`$ $`=`$ $`{\displaystyle \frac{1+P/2\sqrt{1P+9P^2/4}}{2(1P)}},`$
$`c_C`$ $`=`$ $`c_T={\displaystyle \frac{1c_D}{2}}.`$ (6)
Here all the three strategies are present for arbitrary values of $`P`$. Notice that the concentrations of $`C`$ and $`T`$ strategies are the same due to the symmetries of Eqs. (5). In the limit $`P0`$, however, the concentration of $`D`$ strategy vanishes. Evidently, the concentration of the three strategies becomes equal when the evolution is governed exclusively by the mutation ($`P=1`$).

For model B the corresponding equations of motion are similar to those given by Eqs. (5), the differences appear in the first terms proportional to $`P`$. Namely,

$`\dot{c}_D`$ $`=`$ $`Pc_D(1P)c_D(c_C+c_T),`$
$`\dot{c}_C`$ $`=`$ $`+P(c_T+c_D)\pm (1P)c_Dc_C,`$ (7)
$`\dot{c}_T`$ $`=`$ $`P2c_T\pm (1P)c_Dc_T,`$
where the average payoff values are given by Eqs. (4) and the conditions of validity of the upper and lower signs are defined as above. The analytical solution of these equations predicts strikingly different behavior in the stationary state epdg3s , that is, for $`0<P<1/2`$

$`c_D`$ $`=`$ $`{\displaystyle \frac{12P}{1P}},`$
$`c_C`$ $`=`$ $`{\displaystyle \frac{P}{1P}},`$ (8)
$`c_T`$ $`=`$ $`0,`$
while the system goes to the absorbing state ($`c_C=1`$ and $`c_D=c_T=0`$) for $`P>1/2`$. The most surprising result is the extinction of $`T`$ strategy if $`P>0`$.

We have to emphasize the non-analytical behavior in the limit $`P0`$. As illustrated in the upper plot of Figure 3, without the mutation ($`P=0`$) the system evolves toward either a homogeneous $`D`$ state ($`c_D=1`$) or a mixed state composed of $`C`$ and $`T`$ strategies with a ratio depending on the initial conditions. However, the homogeneous $`D`$ state is unstable against $`T`$ invasions, therefore in the close vicinity of this state some small perturbations can drive the system toward the state of $`C`$+$`T`$. Conversely, this mixed state becomes unstable at a given concentration (where $`m_D=m_C=m_T`$) against small perturbations increasing $`c_C`$ and $`c_D`$ simultaneously. In other words, the system evolves toward the $`D`$ dominance when the state is positioned on the right hand side of dashed line (see the upper plot in Figure 3 as a results of fluctuations. This feature explains why the system is so sensitive to applied external constraints.

Figure 3 illustrates that for model A the mutation drives the (concentration) trajectories away from the boundaries. On the right hand side of the dashed line ($`m_D>m_C=m_T`$), $`c_C`$ and $`m_D`$ decrease while $`c_D`$ and $`c_T`$ increase until one crosses the dashed line. On the left hand side all the initial states tend toward the only fixed point given by Eq. (6). For model B, however, there is no fixed point on the left hand side. In this region the external insertion of $`C`$ strategies increases the value of $`c_C`$ until $`m_D`$ becomes larger than $`m_C=m_T`$ and then the $`D`$ invasion drives the system toward the fixed point defined by Eq. (8). During the $`D`$ invasion the external constraint can compensate only the loss of $`C`$ strategies. Consequently, the $`T`$ strategies die out exponentially fast.

Notice that the variation of $`b`$ leaves the fixed points unchanged, but modifies only the slope of the dashed line separating the two regions mentioned above.

Within the framework of mean-field theory, the extinction of $`T`$ strategies is a consequence of the fact that here $`m_T=m_C`$ \[see Eq. (4)\] in contrary to the spatially extended case, as illustrated in Figure 1.

## Monte Carlo simulations

Systematic Monte Carlo simulations have been performed on a square lattice consisting of $`L\times L`$ sites with periodic boundary conditions, $`L`$ varying from 200 to 1500. The larger sizes were used in the vicinity of the critical points. Each run started from a random initial state. During the simulations we have monitored the number of players playing a given strategy ($`N_\alpha `$; $`\alpha =D`$, $`C`$ or $`T`$) and the payoffs related to a given strategy. After some relaxation time we have determined the average concentrations

$$c_\alpha =N_\alpha /L^2$$
(9)
and the fluctuations

$$\chi _\alpha =L^2(N_\alpha /L^2c_\alpha )^2$$
(10)
by averaging over a sampling time varying between $`10^4`$ and $`10^6`$ Monte Carlo steps (MCS) per sites. The results obtained respectively for model A and B are the following.

### Results for Model A

Figure 4 shows a typical strategy distribution for the stationary state at a small value of $`P`$. In contrary to the mean-field prediction \[see Eqs. (6)\] the system is dominated by the $`T`$ strategies. The randomly inserted $`D`$ and $`C`$ strategies form small islands. Occasionally the larger $`C`$ islands are occupied by $`D`$s, however, a consecutive $`T`$ invasion will eliminate the larger $`D`$ territories and maintains the $`T`$ dominance. At the same time this process prevents the formation of large $`C`$ islands inside a $`T`$ domain.

One can observe in Figure 5 that when increasing the value of $`P`$, the concentration of $`D`$ and $`C`$ strategies increases monotonously. In the limit $`P1`$, the strategy distribution on the lattice tends toward a random (uncorrelated) one $`c_D=c_C=c_T=1/3`$ in agreement with the classical mean-field theory \[see Eq. (6)\]. In this case, instead of the neighbor invasions, the system evolution is ruled by the stochastic mutation mechanism.

As shown in Figure 5, the Monte Carlo data agree remarkably well with the results of the pair approximation. This pair approximation is considered as a generalized mean-field theory taking the nearest-neighbor correlations explicitly into account. The details of this calculation are available in many previous works epdg2s ; epdg3s ; ST . The good agreement refers to the absence of long-range correlations which is observable in the “homogeneous” strategy distribution (see Figure 4). It is worth mentioning that the pair approximation is capable to describe the dominance of $`T`$ strategies in the limit $`P0`$.

### Results for Model B

In order to visualize the relevant differences between the two models at small $`P`$ values the strategy distribution for model B is displayed in Figure 6. When comparing the corresponding snapshots (Figures 4 and 6) the reader can easily recognize the most striking differences. Namely, the appearance of a strongly correlated spatial structure for model B. In this case the formation of large $`C`$ domains inside the sea of $`T`$ strategies is not prevented by the random appearance of $`D`$ mutants as happened in the previous case. The large $`C`$ domains (white areas in Figure 6), however, are unprotected against the $`D`$ invasion. Figure 6 shows some $`D`$ domains (black areas) invading the $`C`$’s territories. These $`D`$ domains are “strip-like” because their territories are invaded simultaneously by the $`T`$ strategies. This invasion process is similar also for larger values of $`P`$ (but $`0<P<P_{c1}`$, see later), and only the average invasion velocity changes. On the other hand, the randomly inserted $`C`$ strategies survive and accumulate in the $`T`$ domains. Consequently, far behind the $`T`$-$`D`$ invasion front the $`T`$’s territory will be occupied by the externally inserted $`C`$s and then this area becomes unprotected against the $`D`$ invasion. Sooner or later this area will be invaded by $`D`$s and the above process repeats itself. This means that the cyclic invasion maintains a self-organizing domain structure. Here we have to emphasize that this cyclic (rock-scissors-paper game like) dominance is provided by this external constraint.

Similar processes are observed in the forest-fire models GK ; DS  introduced by Bak et al. BCT  to model the phenomenon of self-organized criticality. In these models each cell can be in one of the following three states: non-burning tree, burning tree and ash. The dynamics are governed by cyclic dominance, similarly to our model B. Note that the consequence of cyclic invasion with three (or more) states are studied in Lotka-Volterra models Lotka ; Volterra ; AD ; ST  and in cyclically dominated voter models TI ; Tain93 ; T94 ; FKB ; SSM .

In model B the transition from the $`T`$ to $`C`$ state introduces a characteristic length and time unit, both proportional to $`1/P`$. In other words, this length unit is characteristic to the typical (linear) size of the $`T`$+$`C`$ domain, and the time unit corresponds to periodic time of cyclic invasion processes at a given site.

When increasing the value of $`P`$, the typical size of $`T`$+$`C`$ domains decreases and the concentration of $`D`$’s increases. It is found that the $`T`$ strategies die out if $`P>P_{c1}=0.1329(1)`$. Figure 7 shows a typical snapshot in the vicinity of this critical value. In this case the external support is sufficiently strong to maintain small $`C`$ clusters inside the $`D`$ domains. The most remarkable feature of this snapshot is that the $`T`$’s form non-uniformly distributed small (isolated) colonies. The observation of time evolution of configuration shows that these $`T`$ colonies walk randomly, they can extinct spontaneously, a single colony can split into two, or two colonies can merge. This phenomenon is analogous to the branching annihilating random walks (BARW) exhibiting a critical transition when varying the control parameters CT96 ; CT98 . The corresponding critical transitions, both for our model B and for BARW, belongs to the so-called directed percolation universality class Janssen ; Grassb .

For $`P>P_{c1}`$, the concentration of $`D`$ decreases monotonously if $`P`$ is increased and vanishes at $`P=P_{c2}=0.3678(1)`$. This extinction process is similar to the previous one, i.e. it also belongs to the DP universality class. The similarity in the correlations is recognizable in the spatial distribution of the extincting strategies when comparing the snapshots displayed in Figures 7 and 8.

For $`P>P_{c2}`$, any initial state evolves toward the absorbing state where all the players follow the $`C`$ strategy.

The results of our systematic investigations are summarized in Figure 9. Systematic numerical investigations in the close vicinity of the critical points show that the vanishing concentrations follow the same power law behavior. Namely,

$`c_T`$ $`=`$ $`(P_{c1}P)^\beta ,`$
$`c_D`$ $`=`$ $`(P_{c2}P)^\beta ,`$ (11)
in the limits $`P_{c1}P0`$ and $`P_{c2}P0`$ respectively and $`\beta =0.57(3)`$ in both cases epdg3s . Within the statistical error this value of the exponent $`\beta `$ agrees with the one of the 2+1 dimensional directed percolation BFM ; JFD .

As expected, these critical transitions are accompanied with the divergence of concentration fluctuations, i.e.

$`\chi _T`$ $`=`$ $`(P_{c1}P)^\gamma ,`$
$`\chi _D`$ $`=`$ $`(P_{c2}P)^\gamma ,`$ (12)
in the vicinity of the corresponding critical points. The numerical fitting yields $`\gamma =0.37(9)`$ in agreement with the DP values BFM ; JFD ; Janssen ; Grassb .

Despite the same universal behavior there is a remarkable difference between the two extinction processes. The second extinction process (at $`P=P_{c2}`$) results in a frozen (time independent) absorbing state. Conversely, the transition at $`P=P_{c1}`$ is an example where the extinction of $`T`$ strategies happens on a fluctuating background. In other words, the properties of the absorbing state (frozen or fluctuating) do not affect the critical behavior of our model.

As demonstrated in Figure 9 the results of Monte Carlo simulations are reproduced qualitatively well by the pair approximation epdg3s . The striking differences are related to the long-range correlations accompanying the critical transitions at $`P=P_{c1}`$ and $`P_{c2}`$. Due to the strongly correlated domain structure, illustrated in Figure 6, the largest deviation can be observed for small $`P`$ values. We note that the concentration fluctuations, defined by Eq. (10), also diverge in the limit $`P0`$. Unfortunately, in this particular case, we could not deduce a reliable value for the exponent $`\gamma `$ because of the significant size effects. Further systematic analyzes are required to clarify what happens in this limit.

## Conclusions

We have studied quantitatively the effect of external constraints on the emergence of cooperation in an evolutionary prisoner’s dilemma game with three possible strategies (cooperation, defection and tit for tat). In the present spatial model the players are distributed on a square lattice and their interactions are restricted to nearest neighbors. The Darwinian selection rule is modeled by the adoption of the neighboring successful strategies. This evolutionary process is superimposed by two types of mutation mechanisms (external constraints) whose strength is characterized by a control parameter $`P`$.

The choice of these three possible strategies yields non-analytical behavior in the limit $`P0`$ for both the mean-field approximation and Monte Carlo simulation. The time-dependent predictions of mean-field theory are sensitive to the small perturbations.

According to the Monte Carlo simulations, in the absence of external constraint the system tends toward a frozen state composed from $`C`$ and $`T`$ strategies whose ratio depends on the initial concentrations. For both types of external constraints (models A and B) the system evolves toward a stationary state independently of the initial condition, and the defector concentration vanishes linearly as $`P0`$. In the limit $`P0`$, however, model A and B will exhibit different ratio of $`C`$ and $`T`$ strategies. This difference is related to the appearance of self-organizing patterns for model B. The present investigation indicates that such a society of strategies (or species) are very sensitive to the type of external supports (or the ration of mutation rates).

The measure of mutual cooperation can be well characterized by the average payoff whose maximum (4) can be reached only in the absence of defectors. Figure 10 compares the Monte Carlo results for the models A and B. Surprisingly, for weak external support (small $`P`$) the average payoff is larger for model A than for model B. In contrary to the naive expectation, the weak support of defenseless cooperators results in opposite consequence. Namely, this mechanism feeds the defectors and simultaneously prevents their elimination by the retaliatory ($`T`$) strategies.

Examples from the political and economical world justifies the above conclusions. In general, the exploiters are preferred by the governmental support for the defenseless layer of a society. The most dangerous effect is the reduction in the $`T`$ type population which can maintain the mutual cooperation against the exploiters. From the view point of cooperation, it is better to help those individuals who are able to prevent themselves against the exploitation.

Evidently, for sufficiently large $`P`$ values the random insertion of $`C`$ strategies can provide their dominance. In this case the A type external support is preferred to the B one if we wish to improve cooperation. Above a threshold value this type of external constraints yields a homogeneous $`C`$ state which is defenseless against any defector appearing occasionally in a real system. Further systematic research is required to clarify what happens in those models where the mutation mechanism is characterized by three independent control parameters.

The present study confirms that the $`T`$ strategy is able to prevent the spreading of defection in the spatial models. We have to emphasize, however, that according to the simplest mean-field theory, $`T`$ dies out if the external support is of B type. Consequently, the defectors will dominate those systems where the mean-field theory is exact (e.g. infinite range of interaction, or randomly chosen partnership). In these mean-field like systems, the games between the “parent“ and “its offspring” is not emphasized (they are not neighbors), which is an advantage for the defectors comparing to spatially extended models. In the light of this feature our investigations imply many interesting questions related to the transition from the “short range” spatially extended systems to the “long range” of mean-field like ones.

#### Acknowledgments.

This work was supported by the Hungarian National Research Fund under Grant No. T-23552 and by the Swiss National Foundation.