text:
# A high-efficiency quantum non-demolition single photon number resolving detector

(October 01, 2003)

## Abstract

We discuss a novel approach to the problem of creating a photon number resolving detector using the giant Kerr nonlinearities available in electromagnetically induced transparency. Our scheme can implement a photon number quantum non-demolition measurement with high efficiency ($``$99%) using less than 1600 atoms embedded in a dielectric waveguide.

In recent years we have seen signs of a new technological revolution, caused by a paradigm shift to information processing using the laws of quantum physics. One natural architecture for realising quantum information processing (QIP) technology would be to use states of light as the information processing medium. There have been significance developments in all optical QIP following the recent discovery by Knill, Laflamme and Milburn that passive linear optics, photo-detectors, and single photon sources can be used to create massive reversible nonlinearitiesKLM . Such nonlinearities are an essential requirement for optical QIP and many communication applications. These nonlinearities allow efficient gate operations to be performed. In principle, fundamental operations such as the nonlinear sign shift and CNOT gates have been demonstrated experimentallyPittman03 ; OBrien03 ; Gasparoni04 . However, such operations are relatively inefficient (they have a probability of success significantly less than 50%) and so are not scalable, due primarily to the current state of the art in single photon sources and detectors. Good progress is being made on the development of single photon sourcesyamamoto1 ; yamamoto2 . Absorptive single photon resolution detection is possiblemiller ; rehacek ; waks03 , with efficiencies up to $``$90% (visible spectrum) and $``$30% (infrared, microwaves). However, true universal optical QIP will require significant further improvements in detector efficiencies, which will likely require a drastic change of approach to detection technologywaks03 ; james02 ; iammoglu02 .

In this letter, we propose an implementation of the quantum non-demolition (QND) single-photon detection scheme originally described by Imoto, Haus and Yamamotoimot85 , with the required optical nonlinearity provided by the giant Kerr effect achievable with AC Stark-shifted electromagnetically induced transparency (EIT)Imamoglu96 . We show below that the scheme uses about 1600 EIT atoms and a weak pulse in the probe mode to achieve an error probability less than 1%. The effect of the QND measurement in turn means that signal photons are not destroyed and can be reused if requiredgrang98 . Furthermore, for a signal mode in a superposition state (such as a weak coherent state) the number-resolving QND measurement projects the signal mode into a definite number stateMilburn 1984 , so the detector can be used as a heralded source of number eigenstates. We focus on EIT as an example here since it has already been used successfully to demonstrate cross-Kerr nonlinearities at low light levels in rubidiumharr99 ; harris04 . From the perspective of realising our detector application, we concentrate on potential condensed-matter mechanisms of EIT. However, clearly any system capable of producing a comparable form and strength of Kerr interaction can be used, including optical fibersli04 , silica whispering-gallery micro-resonatorskipp04  and cavity QED systemsgrang98 ; kimble95 .

Before we begin our detailed discussion of the EIT detection scheme, we first consider the photon number QND measurement using a cross-Kerr nonlinearityimot85 ; howe00 . The Kerr Hamiltonian has the canonical form

$$H_{\text{QND}}=\mathrm{}\chi a^{}ac^{}c,$$
(1)
where the signal (probe) mode has the creation and destruction operators given by $`a^{},a`$ ($`c^{},c`$) respectively, and $`\chi `$ is the strength of the nonlinearity. If the signal field contains $`n_a`$ photons and the probe field is in an initial coherent state with amplitude $`\alpha _c`$, the cross-Kerr optical nonlinearity causes the combined system to evolve as

$$|\mathrm{\Psi }(t)_{out}=e^{i\chi ta^{}ac^{}c}|n_a|\alpha _c=|n_a|\alpha _ce^{in_a\chi t}.$$
(2)
We observe immediately that the Fock state $`|n_a`$ is unaffected by the interaction, but the coherent state $`|\alpha _c`$ picks up a phase shift directly proportional to the number of photons $`n_a`$ in the $`|n_a`$ state. If we measure this phase shift using a homodyne measurement (depicted

schematically in Fig. 1, we can infer the number of photons in the signal mode $`a`$. The homodyne apparatus allows measurement of the quadrature operator $`\widehat{x}(\varphi )ce^{i\varphi }+c^{}e^{i\varphi }`$, with an expected result $`\widehat{x}(\varphi )=2\text{Re}\left[\alpha _c\right]\mathrm{cos}\delta +i2\text{Im}\left[\alpha _c\right]\mathrm{sin}\delta `$, where $`\delta =\varphi +n_a\chi t`$. For a real initial $`\alpha _c`$, a highly efficient homodyne measurement of the momentum quadrature $`Y\widehat{x}(\pi /2)`$ would yield the signal $`Y=2\alpha _c\mathrm{sin}\left(n_a\chi t\right)`$ with a variance of one, thus giving a signal-to-noise ratio of $`\text{SNR}_Y=2\alpha _c\mathrm{sin}\left(n_a\chi t\right)`$. If the input in mode $`a`$ is either the Fock state $`|0`$ or $`|1`$, the respective output states of the probe mode $`c`$ are the coherent states $`|\alpha _c`$ or $`|\alpha _ce^{i\chi t}`$. Using the momentum quadrature measurement, the probability of misidentifying one of these states for one another is then $`P_{\text{error}}=\frac{1}{2}\text{erfc}(\text{SNR}_Y/2\sqrt{2})`$. A signal-to-noise ratio of $`\text{SNR}_Y=4.6`$ would thus give $`P_{\text{error}}10^2`$. To achieve the necessary phase shift we require $`\alpha _c\mathrm{sin}\left(\chi t\right)2.3`$, which can be achieved in a number of ways dependent upon the range of values available for $`\alpha _c`$ and $`\chi t`$. For example, we could choose $`\alpha _c2.3`$ with $`\chi t`$ small and satisfy the above inequality; alternatively we could choose $`\chi t=\pi /2`$ with $`\alpha _c=2.3`$. The particular regime chosen depends on the strength of the Kerr nonlinearity achievable in the physical system.

Figure 2 generalizes the detector shown schematically in Fig. 1 to the case where the polarization of the input state is resolved into different paths by a polarizing beamsplitter. In general we may wish to apply different phase shifts to the two distinct polarizations, but in the case shown in Fig. 2 an identical phase shift is applied to each path, so that the detector is insensitive to the polarization of the input state. This is a particularly useful approach when the efficiency of the EIT system and/or the optical propagation path (e.g., as provided by a photonic crystal waveguide optimized for either TE or TM modes) is polarization-dependent.

We now address the generation of the cross-Kerr nonlinearity required to perform the QND measurement. We consider a model (depicted in Fig. 3) of the nonlinear electric dipole interaction

between three quantum electromagnetic radiation fields with angular frequencies $`\omega _a`$, $`\omega _b`$, $`\omega _c`$ and a corresponding four-level $`𝒩`$ atomic systembeau03 . Mode $`b`$ should be thought of as a pump or coupling field: by choosing the correct conditions, we can factor both the coupling field and the atom out of the evolution of the atom-field system, creating an effective cross-Kerr nonlinear interaction between modes $`a`$ and $`c`$. The effective vacuum Rabi frequency for each mode is defined as $`\left|\mathrm{\Omega }_k\right|^2=(\sigma _k/\eta _k𝒜)A_k\mathrm{\Delta }\omega _k/8\pi `$, where $`\sigma _k3\lambda _k^2/2\pi `$ is the resonant atomic absorption cross section at wavelength $`\lambda _k2\pi c/\omega _k`$cohe92 , $`\eta _k`$ is the refractive index of the waveguide material, $`𝒜`$ is the effective laser mode cross-sectional area, $`A_kf_ke^2\omega _k^2/2\pi ϵ_0m_ec^2`$ for a transition with oscillator strength $`f_k`$, and $`\mathrm{\Delta }\omega _k`$ is the bandwidth of the profile function describing the adiabatic interaction of a pulsed laser field with a stationary atomblow90 ; chan02 ; domo02 .

It is difficult to achieve a substantial vacuum Rabi frequency using free-space fieldsvane00 , but encapsulating one or more atoms in a waveguide (such as a line defect in a photonic crystal structure) allows field transversality to be maintained at mode cross-sectional areas on the order of $`𝒜(\lambda /3\eta )^2`$. Consider, then, a two-dimensional photonic crystal waveguide constructed from diamond thin film ($`\eta =2.4`$), with nitrogen-vacancy color centers fabricated in the center of the waveguide channelhemm01 ; shah02 . The optical transition at 637 nm in NV-diamond has an oscillator strength of approximately 0.12, within a factor of three of rubidium, which has been used successfully to demonstrate cross-Kerr nonlinearities at low light levelsharr99 ; harris04 . An EIT transmission window of about 8 MHz has been observed experimentallyhemm01 , so a pulse with $`\mathrm{\Delta }\omega _k/2\pi 5`$ MHz should propagate through this window with negligible loss. The corresponding vacuum Rabi frequency is therefore $`\mathrm{\Omega }3.6`$ MHz.

We consider a number $`N`$ of $`𝒩`$ atoms, fixed and stationary within a cylinder that is narrow but long compared to the optical wavelengths, with the three frequency modes of the system driven by Fock states containing $`n_a`$, $`n_b`$, and $`n_c`$ photons, respectively. If the durations of the three pulse envelope functions are long compared to the lifetime of atomic level $`|2`$, the evolution of the amplitude where all of the atoms are in the ground state $`|1`$ is simply given by $`|1,n_a,n_b,n_ce^{iWt}|1,n_a,n_b,n_c`$. In general, $`W`$ is complex (see Eq. (139) in Ref. beau03 ); other states contribute to the full evolution and there is photon absorption and loss in the system. However, for our detector we require the probability of even single photon loss from mode $`a`$ to be very small, with a real $`W`$ preserving the norm of $`|1,n_a,n_b,n_c`$. For this, we assume that the laser frequencies $`\omega _a`$ and $`\omega _b`$ are both precisely tuned to the corresponding atomic transition frequencies, so that the EIT Raman resonance condition is satisfied. Then, in the weak-signal regime, we assume that the decoherence rate $`\gamma _k`$ is dominated by spontaneous emission from atomic level $`|k`$, and that $`|\mathrm{\Omega }_a|\gamma _2`$, so that $`W`$ is given by

$$W=\frac{N\left|\mathrm{\Omega }_a\right|^2\left|\mathrm{\Omega }_c\right|^2n_an_c}{\nu _c\left|\mathrm{\Omega }_b\right|^2n_b+i\left(\gamma _4\left|\mathrm{\Omega }_b\right|^2n_b+\gamma _2\left|\mathrm{\Omega }_c\right|^2n_c\right)},$$
(3)
where $`\nu _c\omega _c\omega _{43}`$. In principle, in this regime we must have $`\nu _c|\mathrm{\Omega }_b|^2n_b\gamma _4|\mathrm{\Omega }_b|^2n_b+\gamma _2|\mathrm{\Omega }_c|^2n_c`$ to obtain a nearly real $`W`$ and a low residual absorption. As we shall show below, in practical cases where $`|Wt|1`$, this constraint can be substantially relaxed. For NV-diamond, $`\gamma _2^1=2\times 25`$ ns beve01 ; beve02 , and the spin decoherence lifetime is 0.1 ms shah02 ; hemm01 , so for $`N10000`$, dephasing can be neglected and $`|\mathrm{\Omega }_a|/\gamma _21`$.

Under these conditions, the state $`|1,n_a,n_b,n_c`$ simply acquires a phase-shift which is the basis for the emergence of the approximate cross-Kerr non-linearitybeau03 . In general, when the pump and probe fields are intense coherent states (parameterized by $`\alpha _b`$ and $`\alpha _c`$, respectively), the evolution of the state $`|1,n_a,\alpha _b,\alpha _c`$ can be approximated as $`|1,n_a,\alpha _b,\alpha _ce^{in_a(\theta i\kappa )}`$beau03 , where the angle $`\theta `$ and the residual absorption $`\kappa `$ is defined for the case $`\left|\mathrm{\Omega }_b\right|^2=\left|\mathrm{\Omega }_c\right|^2`$ by

$$\theta i\kappa =\frac{N\left|\mathrm{\Omega }_a\right|^2t}{\nu _c\left|\alpha _b\right|^2+i\left(\gamma _4\left|\alpha _b\right|^2+\gamma _2\left|\alpha _c\right|^2\right)}.$$
(4)
This is equivalent to a damped evolution generated by the cross-Kerr Hamiltonian of (1), with $`\theta =\chi t`$.

What values of $`\theta `$—and, therefore, $`\text{SNR}_Y`$—are achievable? To establish an estimate we need to make several assumptions about the physical system and its geometry. We assume that the interaction region (where the light and $`𝒩`$ atoms interact) is encapsulated within the photonic crystal waveguide described above, and that the pulses have weakly super-Gaussian profiles so that the bandwidth-interaction time product is $`\mathrm{\Delta }\omega _kt3\pi `$, giving $`\left|\mathrm{\Omega }_a\right|^2t81\eta \gamma _2/8\pi `$. Suppose now that the largest phase shift that can be applied in practice (without significantly distorting the signal pulse) is $`\theta _{\text{max}}`$, and that the corresponding value of $`\alpha _c`$ needed to obtain a given signal-to-noise ratio is therefore $`\alpha _c=\text{SNR}_Y/2\theta _{\text{max}}`$ (if $`\theta _{\text{max}}1`$). Using Eq. (4), we can calculate the minimum number of atoms and the corresponding minimum detuning $`\nu _c`$ needed to generate this phase shift, by taking the real part and solving for $`\nu _c`$ explicitly. As $`\nu _c`$ is required to be real we find

$`N_{\text{min}}`$ $`=`$ $`{\displaystyle \frac{2\theta _{\text{max}}}{\left|\mathrm{\Omega }_a\right|^2t}}\left(\gamma _4\left|\alpha _b\right|^2+\gamma _2\left|\alpha _c\right|^2\right)`$ (5)
$`=`$ $`{\displaystyle \frac{\gamma _2\text{SNR}_Y^2}{2\theta _{\text{max}}\left|\mathrm{\Omega }_a\right|^2t}}\left({\displaystyle \frac{\gamma _4\left|\alpha _b\right|^2}{\gamma _2\left|\alpha _c\right|^2}}+1\right)`$
with

$$\nu _{c\text{min}}=\frac{N_{\text{min}}\left|\mathrm{\Omega }_a\right|^2t}{2\left|\alpha _b\right|^2\theta _{\text{max}}}=\frac{\gamma _4\left|\alpha _b\right|^2+\gamma _2\left|\alpha _c\right|^2}{\left|\alpha _b\right|^2}$$
(6)
When we choose $`N=N_{\text{min}}`$ and $`\nu _c=\nu _{c\text{min}}`$, we find that $`\kappa =\theta `$, so that the state $`|1,n_a,\alpha _b,\alpha _c`$ evolves according to $`|1,n_a,\alpha _b,\alpha _ce^{(1+i)n_a\theta }`$. Therefore, when $`\theta 1`$, the residual absorption of a signal photon in mode $`a`$ can also be made intrinsically small, even though the detuning is not large compared to the absorption linewidthbeau04b .

Figure 4 shows the minimum number of NV-diamond color centers as a function of the maximum single-photon phase angle for three different values of the error probability $`P_{\text{error}}`$. We note that $`\gamma _4=\gamma _2`$ for the optical transitions in NV-diamond, and for convenience we have chosen $`n_b=10n_c`$, requiring a minimum detuning $`\nu _{c\text{min}}/\gamma _2=1.1`$ in all three cases. Note that the minimum number of atoms needed to obtain a given phase shift decreases as the phase shift increases, because the constraint that the signal-to-noise ratio remain constant allows the values of $`|\alpha _b|`$ and $`|\alpha _c|`$ to decrease as $`\theta _{\text{max}}`$ increases. As an example, we choose $`P_{\text{error}}=0.01`$ and $`\theta _{\text{max}}=0.01`$ radians, requiring $`n_c=5.6\times 10^4`$ to maintain $`\text{SNR}_Y=4.6`$. Therefore, $`N1600`$ is sufficient to achieve the desired phase shift, resulting in a residual absorption of less than 1%.

There is considerable flexibility in the engineering design parameter space for this implementation of the QND detector. For example, if we choose $`\theta _{\text{max}}=0.1`$, following the design procedure outlined above for $`P_{\text{error}}=0.01`$ leads to a residual absorption of almost 10% for $`N160`$ and $`n_c560`$. However, if we increase the number of atoms to 800 and the detuning to $`\nu _c=11\gamma _2`$, then the absorption is reduced to 1%. Note also that the detector can perform a QND measurement on a Fock state with $`n_a>1`$ with single-photon resolution. For example, as $`n_a`$ increases from 1 to 2, the phase shift $`n_a\theta `$ doubles, and the SNR also increases from 4.6 to 9.2 for constant $`\alpha _c`$. The detector sensitivity improves until the phase shift becomes so large that one of two fundamental limits is reached: either the SNR decreases below the 1% error threshold, or the strong nonlinear interaction begins to significantly distort the pulse profile of the signal Fock state.

In summary, we have presented a scheme for a highly efficient photon number quantum non-demolition detector (with single-photon resolution) based on the cross-Kerr nonlinearity produced by an EIT condensed matter system with approximately 1600 color centers. We have explored several different operating regimes, and we have examined in detail the performance of the detector for an NV-diamond photonic crystal waveguide system. In particular, we have shown that efficient detection is possible with small phase shifts, which will likely be necessary to ensure that the EIT optical nonlinearity doesn’t distort the pulse envelope of the signal state. Future modelling will address detector performance for pulsed Fock state profile functions, and much experimental work remains to be done to implement such a detector. For example, fabricating EIT atomic or molecular systems into a dielectric waveguide is challenging but feasible. A method for orienting the color center spins uniformly in such a condensed-matter system must be found, and spatial hole-burning techniques will be needed to overcome the effects of inhomogeneous broadening on the transparency of a condensed matter mediumturu02 . The overall efficiency of the detector is likely to be limited by the efficiency of the homodyne measurement of the phase shift, which will depend on the degree to which the homodyne detector can be spatiotemporally mode-matched to a single-photon signal. Nevertheless, EIT provides us with the best known candidate mechanism for the implementation of the original QND proposal by Imoto, Haus, and Yamamotoimot85 , and even a weak nonlinearity could allow efficient reuse of resources in linear optics quantum computation schemes.

Acknowledgments: This work was supported by the European Project RAMBOQ. KN acknowledges support in part from MPHPT, JSPS, Asahi-Glass and Japanese Research Foundation for Opto-Science and Technology research grants. WJM acknowledges support in the form of a JSPS fellowship.


format: text

text:
# Analysis of Norms Game in networked societies

## Abstract

Norms, defined as generally accepted behaviour in societies without central authority (and thus distinguished from laws), are very powerful mechanism leading to coherent behaviour of the society members. This paper examines, within a simple numerical simulation, the various effects that may lead to norm formation and stability. The approach has been first used by Axelrod, who proposed two step model of norm and meta-norm enforcement. We present here an extension and detailed analysis of the original work, as well as several new ideas that may bear on the norm establishment mechanisms in societies. It turns out that a relatively simple model for simulated norm enforcement predicts persistent norm breaking even when it is associated with high punishment levels. The key factors appear to be the combination of the level of penalty for breaking the norm and proximity of norm enforcers. We also study a totally different mechanism of norm establishment, without meta-norms but using instead the direct bonus mechanism to norm-enforcers.

## 1 Definitions of the Norms Game variants

The original work of Axelrod , later reprinted and expanded, as Chapter 3 of  has described a possible scenario for modelling norm formation, called the Norms Game (NG), and its higher level variant, the Meta Norms Game (MNG). The games are based on extended Prisoner’s Dilemma game, including additional action by the rest of community: *punishment* of a *defector* (NG), and additional possibility of punishment of *for failing to punish the defector* (MNG). The idea behind the model is to see if such enforcement may lead to formation of stable cooperative society, and whether simple enforcement has to be augmented by meta-enforcement, to achieve stability of the system.

### 1.1 The classic formulation of Norms Game.

The model presented by Axelrod is based on an $`N`$-person PD game. Within the game, each player $`i`$ has several "action opportunities" each turn, during which he can either play safe, in which case he gets no benefit, or defect in which case he gets a treasure of $`T>0`$, while players $`ji`$ get extra harm $`H<0`$. The willingness of an agent $`i`$ to commit a defection is called boldness $`b_i`$. For obvious reasons, the only stable strategy in such game is that of total defection — each "safe player" is at distinct disadvantage with the defectors.

The next step is the introduction of norm enforcement activities (the norm being understood as the neutral, safe behaviour).

If the player $`i`$ does not defect, other players do not take action. If, on the other hand, $`i`$ does break the norm, other players $`j`$ may, with probability $`v_j`$ (vengefulness), punish him for misbehaviour. In the original treatment, the punishment value $`P`$ was assumed greater than the treasure $`T`$. At the same time, each punishing agent pays punishment costs $`E`$ (enforcement). The punishment for defection constitutes the proper Norms Game: society uses active measures to enforce the norms.

After each round of interactions between agents, Axelrod used simple population adjustment, giving higher population ratios to those players with higher overall benefits. He has also used a small "mutation" rate to introduce small random shifts in strategies. Results of simulations were presented after 100 generations. In Axelrod’s simulations two out of five runs ended up in boldness $`b1`$ and $`v0`$, while other simulations ended with various levels of vengefulness, with very small boldness. The pure norm enforcing behaviour, depending on the simulation history or initial conditions led to either almost fully defecting society or to one pacified by vengefulness. Axelrod’s argument was that pure NG is too weak to ensure norm respecting behaviour.

Introducing Meta Norms Game, as an extension of the previous model, Axelrod proposed that if agent $`i`$ defects and agent $`j`$ does not punish $`i`$, then with probability $`v^{}`$ any other agent $`k`$ should punish $`j`$ *for failing to enforce the norm!* The meta-punishment $`P^{}`$ and cost associated with it, $`E^{}`$ were assumed to be equal to $`P`$ and $`E`$.

More importantly, the meta-vengefulness $`v^{}`$ was assumed to be equal to $`v`$. This decision was based on rather vague arguments. The link definitely simplified calculations, but, as noted by Axelrod himself, prohibited situations when certain agents would use different values of $`v`$ and $`v^{}`$ to their advantage. Results of simulations with $`v=v^{}`$ were uniform: in all cases the society converged to almost pure high vengefulness, low boldness, norm-observing one.

From his simulations Axelrod has concluded that, for pure Norms Game,

> …at first, boldness levels fell dramatically due to the vengefulness in the population. Then, gradually, the amount of vengefulness also fell because there was no direct incentive to pay the enforcement cost of punishing the defection. Once vengefulness became rare, the average level of boldness rose again, and then the norm completely collapsed. Moreover, the collapse was a stable outcome.

In contrast, the inclusion of meta-norm enforcement had provided the agents with a strong incentive to increase their vengefulness and this to promote the norm (decrease boldness). This was indeed confirmed by the simulation results.

From one point of view this is obvious consequence of the model, in which some characteristics are penalised (boldness in case of NG, boldness and lack of vengefulness in the case of MNG). Lack meta-vengefulness is not penalised at all, just the opposite — it carries its own cost, and is thus disadvantageous evolutionarily. This leads to minimization of the meta-vengefulness in the evolution process – unless, as in Axelrod’s model, it is forcibly fixed, for example by coupling it to vengefulness. By setting $`v=v^{}`$, $`P=P^{}`$ and $`E=E^{}`$ Axelrod has not only simplified the simulations, but also made a sort of self-policing model. He has thus excluded analysis treating meta-enforcement and enforcement as independent yet interacting phenomena.

Additionally, because in the original form, the Norms Game was played with a very small number of agents (20) no scaling effects were visible. One of such effects is the cumulation of penalties due to the model in which anyone is capable of punishing everyone, and individual punishments add up. This results in penalties easily "overcoming" benefits of trespassing. Despite the mentioned oversimplifications of the simulations themselves, the MNG model seems to be interesting enough to use it as a basis for more detailed studies. A significant and open question remains whether the results would be changed by larger and longer simulations and by more flexible model.

### 1.2 Modified form of Norms Game with single enforcers

To extend the validity of the model and to check if the assumptions made by Axelrod could be derived from more general principles we propose a modified form of MNG.

We assume that $`N_T`$ agents ($`N_T1`$) are connected via a network of connections, each interacting with a number of neighbors. In different simulations we have used four basic network topologies, typical for biological or social networks. These types included: fully connected network, in which every agent is connected with all others; random network, where agents are connected by randomly distributed links; nearest neighbour network, where the links are highly clustered and finally the small world networks. The three latter types of networks may be called local, as the number of agents that a given agent interacts with (actively or passively) is much smaller than the total number of agents. This limits the influence and sets stage for possible *i*ndirect influences. General description of such networks may be found in Albert and Barabási , Dorogovtsev and Mendes , Dorogovtsev et al. , Dorogovtsev and Mendes , Newman , Newman and Park . All these networks are characterized by filling ratio, which relates to the average number of links per node (agent). The larger the filling ratio, the more agents are directly connected. In our simulations we have used ratios from 0.005 to 0.02 (which corresponds to circa 5 to 20 links per agent. The difference between the various types of local networks are in the way the agents are connected.

In random (RAND) networks the links are distributed randomly, thus we have a meshed network of links, with agents differing in number of the neighbours, and no general structure.

The nearest neighbor (NN) nets are formed by linking certain number of closest neighbors together. The easiest way to visualize such network is to place the agents on an imaginary circle and connecting each agent to $`n`$ neighbors. The number of connections per agent is constant. Interesting property of NN networks is that for small filling factors, agents on the opposite points at the circle to communicate must go through many intermediaries. For 2000 agents and number of neighbors set at 10, the longest ‘distance’ is 100 ‘hops’. This suggests influence on the changes in norm adherence and enforcement: any change of behavior of the agent $`i`$ is seen immediately by his closest neighbors (who can see him through direct links) but only after considerable delay and filtering by members of society located far from $`i`$.

The Small World (SW) networks, introduced and popularized in recent years, reproduce a curious fact observed in many natural and human-produced networks, namely that the distance between any two nodes of the network, measured in links needed to connect them, is usually much smaller than that in nearest neighbor or even random networks of the same filling ratio. The name of the network category comes exactly from such observation. One of the models for SW networks is a simple reworking of the NN model: one takes a small (even very small!) fraction of the links from between nearest neighbors and applies them instead between random agents. Keeping in mind the visualization of NN networks as nodes along a circle, this corresponds to adding connections that criss-cross the circle at random. Due to such shortcuts, even if their number is very small, the average distance between any two nodes drops dramatically. Thus we have a network that for each agent, locally is very similar to NN model (as most of the neighbors are, in fact, the same), but globally the communication through the network is much faster.

Within our MNG scenario each agent has a choice between acting in accordance with a norm, in which case he gets no extra bonus, or breaking the norm (trespassing), in which case he gets extra bonus (treasure) $`T`$. If the agent $`i`$ misbehaves, every agent in the population is harmed by this action. We can imagine this corresponds, for example, to the agent stealing some ‘communal property’. The treasure $`T`$ is the benefit to the trespasser, but the harm $`H`$ to the community may be greater than $`T`$. It is then assumed that the harm $`H`$ is divided evenly among population, each agent ‘losing’ $`H/N_T`$.

Let’s denote the probability of misbehavior (boldness) of agent $`i`$ by $`b_i`$. The payoff $`X(i)`$ of agent $`i`$ is (without any punishments) given by

$$X(i)=b_iTN_Tb_TH/N_T=b_iTHb_T$$
(1)
where $`b_T`$ denotes average of $`b_i`$ over the total population. With $`H>T`$ the overall average payoff becomes negative $`b_T(TH)`$ and is disadvantageous to the community, but as any agent with boldness higher than the average gets more benefit than the rest, such behaviour gives evolutionary advantage, and thus multiplies. Natural end of such a process is a free-for-all, $`b_i=1`$ society.

As in Axelrod’s work we propose that to curb boldness, society enforces the norm through penalization of trespassing. Our model differs slightly in the way the process of detection and punishment is encoded.

Each agent watches his neighbours for acts of misbehavior and is may be willing to punish for such acts. For an agent $`j`$ with vengefulness $`v_j`$ the probability of detection that another agent $`i`$ has broken the norm is just $`v_j`$. The watchfulness and vengefulness result in enforcement effort $`E`$, the overall cost of being vigilant to an agent $`j`$ being $`v_jE`$. The first agent to notice that someone has broken the norm punishes the transgressor with a punishment of $`P`$. Thus only one punishment act take place and there is no cumulation effects, but the value of $`P`$ may be set freely. To calculate the expected penalty value we start with the probability of no-one detecting the misbehaviour of $`i`$ is

$$\underset{(ji)}{}(1v_j),$$
(2)
where the product is taken over all agents $`j`$ with links to $`i`$, that is over the neighbours of $`i`$. As a result the penalty for agent $`i`$ would be given by

$$\pi _P(i)=\left(1\underset{(ji)}{}(1v_j)\right)b_iP,$$
(3)
that is proportional to probability is at least one agent detecting the act of breaking the norm times probability of such act times the penalty value. Without meta-norms the payoff for player $`i`$ is thus

$$b_iTHb_Tv_iE\left(1\underset{(ji)}{}(1v_j)\right)b_iP$$
(4)
So far the difference between the modified model and Axelrod’s one was in the fact that the punishment takes place only once. The real change comes in the way the meta-norm is introduced. In the original model every failure to punish a norm breaker was meta-punished. Here, we have only one punishing agent (the first one to detect the trespasser). The problem is how should the meta-enforcer decide whether non-punishment by certain agent was intentional (due to lack of vengefulness or it’s small value) or just a result of coming in second, after the enforcement has already taken place? One way is to set the probability of meta-punishment as being proportional to the lack of vengefulness of the observed agent:

$$X_{MP}(j,i)=v_j^{}(1v_i),$$
(5)
with $`v_j^{}`$ being the meta-vengefulness of the observing agent $`j`$ and $`1v_i`$ being the lack of vengefulness of the observed agent $`i`$. Thus the expected value of meta-punishment for an agent $`i`$ is

$$\pi _{MP}(i)=\left(1\underset{(ji)}{}\left(1X_{MP}(j,i)\right)\right)P^{}$$
(6)
As our simulations were conceived as extension of the work of Axelrod, we have used similar convention of allowing discretized values of $`b_i`$, $`v_i`$ and $`v_i^{}`$ (from 0 to 1 with steps of 0.1). To see if such discretization does not produce any artificial effects we have also performed calculations in which boldness, vengefulness and meta-vengefulness were allowed to take arbitrary values between 0 and 1. To our slight surprise the results did not differ significantly. The specific results presented in this paper refer to the discrete model.

#### 1.2.1 Results of NG and MNG for large fully connected networks

Lets turn now to the NG and MNG played in relatively large, fully connected networks. Here every agent interacts with all agents as in the case of the small model of Axelrod, and the evolution seems rather simple. The most important fact here is that even a single enforcer or meta-enforcer (say an accidental mutant) would be able to influence the payoffs of all agents, and direct the course of evolution of population.

Lets consider first the situation where there is no meta-norm enforcement at all ($`P^{}=0`$). If one starts from a random population in all simulations the final average value of vengefulness was $`v=0`$. This is due to additional costs borne by the enforcers, and the pressure to minimize this cost. However, as initially there were some norm enforcers (and remember that in a fully connected network just one mutant would suffice to enforce the norm on the whole population!) the enforcement put a clear divide between two basic situations. When the treasure $`T`$ is greater than $`P`$ it pays to be bold, and $`b1`$. On the contrary, for $`T<P`$ we have $`b0`$. This corresponds to common sense observation that the norm would be obeyed if the penalty for breaking it would be sure (which is made probable by large number of observers) and greater than the benefits.

When the meta-norm enforcement is present the situation shifts to meta-level. In all simulations the final average value of meta-vengefulness goes to zero ($`v^{}0`$) — due to the same arguments as vengefulness dropped in previously described simulation. However, as already noted, in a fully connected network even a single meta-enforcer is able to meta-punish everyone in society! Thus, although nominally the number of these meta-enforcers goes is zero in the stable situation, their initial influence and occasional mutant presence are sufficient to distinguish the situation from the previously described case of $`P^{}=0`$.

The results of simulations differ in small details (such as the number of iterations it takes to achieve stability) but in principle are governed by two simple rules.

* The boldness is decided by relation of $`T`$ vs. $`P`$, as described above
* The vengefulness is decided by relation of the cost of vengefulness $`E`$ and the meta-punishment $`P^{}`$

Again, the above results are very much common sense. Decoupling $`v`$ from $`v^{}`$ results in $`v^{}0`$ for MNG, which contrasts the assumed $`v=v^{}`$ in Axelrod approach. Large size of the population leads to elimination of ‘undecided’ runs — all simulations ended with uniform populations.

#### 1.2.2 Local networks

In contrast to the fully connected case, the local networks present much more interesting norms dependence on the payoff parameters.

As in the original formulation of MNG, in addition to a breeding mechanism, in which certain percentage of agents with the best payoffs were allowed to breed (at the cost of the worst faring agents), we have used additional mutation process. This was simulated as allowing certain fraction of agents (given by mutation ratio) to change their characteristic $`b_i`$, $`v_i`$ and $`v_i^{}`$ values between interactions. Although the number of mutants per iteration was usually kept very small in relation with the population size (e.g. 2 mutants out of 2000 per iteration), its effects were quite important. Without mutation populations sometimes froze in a given configuration, ‘far’ from results obtained with even minimal mutation rates.

Another effect we set out to investigate was the dependence of the behaviour on the starting conditions. Two general types of starting conditions were used: in the first the initial populations had randomly distributed values of $`b_i,v_i`$ and $`v_i^{}`$. In the second, we have started from pure society of, for example, $`b_i=0`$, $`v_i=1`$ and $`v_i^{}=0`$, or $`b_i=1`$, $`v_i=1`$ and $`v_i^{}=0`$ — to observe if such pure societies were evolutionary stable with respect to small mutation rates. It turned out that regardless of the starting condition, if the mutation was present, the final state was very similar.

We present general results of simulation runs through a few examples, and discuss the role of the starting conditions. The simulations were run for 2000 agents using fixed values of $`P=20`$, $`P^{}=3`$, $`E^{}=1`$, $`H=12`$, with top 1% of the population being allowed to reproduce each iteration (at the expense of the bottom 1%). Other parameters ($`T`$, $`E`$) were used to control the simulations and the dependence on them is presented and discussed below.

#### 1.2.3 Similarities and differences between random, NN and SW networks

In our simulations we compare results for all three types of local networks. Despite the fact that the three types of locally connected networks have many different properties (see, for example ) the results of our simulations for random, small world and nearest neighbour networks are strikingly similar (see Figure 1). It turns out that the key factor is the common characteristic of all models, given in the first approximation by the average number of neighbours. The number of neighbors is related network filling factor $`\varphi `$, defined as ratio of the actual links in the network to the total possible number of links:

$$\varphi =\frac{N_{LINKS}}{N_T(N_T1)/2}.$$
(7)
For regular networks, such as NN network, the number of neighbors is the same for all agents, $`N_N=N_{LINKS}=\varphi N_{TOT}/2`$. For small world networks, constructed as in this work, the $`N_N`$ defined above is a very close measure of the number of neighbors, with very few exceptions. For random networks, the dispersion in number of links per agent is greater, but still $`N_N`$ gives the average value.

The differences between local network topologies for the same filling factor values are relatively small. Thus the use of the word ‘local’ is justified: behaviour of an agent is then determined by his immediate surrounding, mainly by the number of agents who observe it and vice-versa, the number of agents it observes.

The two vertical columns of plots in Figure 1 correspond to two cases of $`E/P^{}`$ ratio. The left column has $`E/P^{}=2/3<1`$, which in the case of fully connected network would correspond to large vengefulness $`v`$ – so, to strong pressure to decrease $`b`$. The right column has $`E/P^{}=5/3`$, with meta-penalty smaller than the cost of being vigilant there is little incentive to be a norm enforcer. Indeed, the results of simulations for the two cases are strikingly different. For $`E/P^{}<1`$ $`v`$ is relatively large and $`b`$ resembles the step-like function found in fully connected network. $`v^{}`$ remains relatively small. On the other hand, for $`E/P^{}>1`$ $`v`$ is smaller, $`v^{}`$ larger, and the step like character of $`b`$ is replaced by gradual growth. This corresponds to the presence of trespassers who, in local environments, find that locally the effective punishment is smaller than the treasure value. Thus even at $`T/P`$ much smaller than 1 we have relatively high average $`b`$.

#### 1.2.4 Discussion of the dependence on the filling factor

Figure 2 presents another set of comparisons of the average final $`b`$, $`v`$ and $`v^{}`$ values, this time for a single type of the network topology (Small World) but differing by the filling factor $`\varphi `$. The top row corresponds to $`\varphi =0.005`$, middle row to $`\varphi =0.01`$, bottom row to $`\varphi =0.02`$. Corresponding average number of neighbours changes from 5 to 20. For the case of $`E/P^{}=5/3`$ (right column) the difference is visible mostly in lowering values of $`v`$ and $`v^{}`$ below $`T/P=1`$, and in more pronounced kink in $`b`$ at $`T/P=1`$ with increasing $`\varphi `$.

For $`E/P^{}=2/3`$ one can observe two effects:

* a shift of the step-like increase of $`b`$ as function of $`T/P`$ to values lower than 1;
* an increase in vengefulness $`v`$ in the region of $`T`$ far below $`P`$ for $`E=2`$, where simple explanations suggest no dependence of $`v`$ over $`T`$.

The first effect can find a qualitative explanation through analogy with simplified, small size fully connected model. Then the finite size of neighbourhood changes the condition for the misbehaviour to be advantageous from $`T>P`$ to:

$$T_i>P\left(1\underset{ji}{}(1v_j)\right)P\left(1(1v)^{N_N}\right)$$
(8)
with $`v`$ being the average vengefulness. The expression $`P\left(1(1v)^{N_N}\right)`$ may be dubbed the effective deterrent, as it weights the penalty with probability of being punished. For large number of neighbors $`N_N`$ the term $`(1v)^{N_N}`$ goes to zero, and we recover the $`T>P`$ condition. For small $`N_N`$ the difference is visible as the shift of the step in $`b(T)`$ function.

### 1.3 Simulations with group enforcement

The previous sections have discussed the extreme opposite to Axelrod’s assumption of additive punishment — only a single agent was executing the punishment, and it was assumed that this single agent had enough power to effect the punishment on anyone it had links with. Within this model a single enforcer or meta-enforcer could ‘police’ the whole neighbourhood.

It is interesting to see if the results obtained in previous simulations would be changed when one would require some greater number of agents, acting in unison, to effect a punishment or meta-punishment. Such situation of ‘group action’, or coalition forming against a trespasser seems more similar to ‘real life’.

Lets assume that the punishment would take place if at least $`k_P`$ neighbors of trespassing agent $`i`$ would be willing to enforce the norm. Similarly, meta-punishment would require $`k_{MP}`$ neighbors to act together. In our calculations we would limit ourselves to low values of $`k_P`$ and $`k_{MP}`$ of 2 and 3.

To obtain the expressions for group punishment model one starts with $`\mathrm{\Phi }_m(i)`$ and $`\mathrm{\Psi }_m(i)`$ defined as: probability that exactly $`m`$ neighbors if agent $`i`$ would be willing to enforce the norm (or meta-norm) on $`i`$. We have:

$`\mathrm{\Phi }_0(i)`$ $`=`$ $`{\displaystyle \underset{ji}{}}(1v_j)`$ (9)
$`\mathrm{\Phi }_1(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}{\displaystyle \underset{ji,k}{}}v_k(1v_j)`$ (10)
$`\mathrm{\Phi }_2(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}{\displaystyle \underset{l>k,li}{}}v_kv_l{\displaystyle \underset{ji,k,l}{}}(1v_j)`$ (11)
$`\mathrm{\Phi }_3(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}{\displaystyle \underset{l>k,li}{}}{\displaystyle \underset{m>l,mi}{}}v_kv_lv_m{\displaystyle \underset{ji,k,l,m}{}}(1v_j)`$ (12)
$`\mathrm{\Psi }_0(i)`$ $`=`$ $`{\displaystyle \underset{ji}{}}\left(1v_j^{}(1v_i)\right)`$ (13)
$`\mathrm{\Psi }_1(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}v_k^{}(1v_i){\displaystyle \underset{ji,k}{}}\left(1v_j^{}(1v_i)\right)`$ (14)
$`\mathrm{\Psi }_2(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}{\displaystyle \underset{l>k,li}{}}v_k^{}v_l^{}(1v_i)^2{\displaystyle \underset{ji,k,l}{}}\left(1v_j^{}(1v_i)\right)`$ (15)
$`\mathrm{\Psi }_3(i)`$ $`=`$ $`{\displaystyle \underset{ki}{}}{\displaystyle \underset{l>k,li}{}}{\displaystyle \underset{m>l,mi}{}}v_k^{}v_l^{}v_m^{}(1v_i)^3{\displaystyle \underset{ji,k,l,m}{}}\left(1v_j^{}(1v_i)\right)`$ (16)
Thus, requiring that more than $`k_P`$ agents (who are neighbors) decide group together and punish transgressor $`i`$ (or respectively $`k_{MP}`$ agents group to meta-punish $`i`$) we have the expression for punishment (meta-punishment)

$`\pi _P^{k_P}(i)`$ $`=`$ $`\left(1{\displaystyle \underset{k=0}{\overset{k_P}{}}}\mathrm{\Phi }_k(i)\right)b_iP`$ (17)
$`\pi _{MP}^{k_P}(i)`$ $`=`$ $`\left(1{\displaystyle \underset{k=0}{\overset{k_{MP}}{}}}\mathrm{\Psi }_k(i)\right)P^{}`$ (18)
Obviously, for $`k_P=0`$ and $`k_{MP}=0`$ (which corresponds to just one agent being able to effect punishment) the above expressions reduce to Equations 3, 6.

Figure 3 compares results for SW network at low filling factor $`\varphi =0.005`$, where the problems related to the need to find enough enforcers for the group action should be most visible, for three values of $`k_P=1,2,3`$ (top, middle and bottom row respectively). The shift of high $`b`$ regime to regions of $`T/P<1`$ is even more pronounced than in Figure 2, the reason being clear: it is increasingly more difficult to find two and three agents acting in unison in neighborhoods of average size of 5. However, there is no radical change in behaviour of the system as the whole.

## 2 To enforce or to break the norm?

In Axelrod’s study as well as in the results presented above, boldness and vengefulness were treated as separate characteristics of an agent. Details of our simulations show that in most cases the best payoffs are achieved by agents who at the same time take advantage from breaking the norm and actively pursue enforcement on other trespassers, that is agents with simultaneously high $`b`$ and $`v`$. while perfectly acceptable from mathematical point of view, this is somewhat against psychological expectations: one can either be a policeman or a thief. We dub this case the exclusive model.

To see if such contradiction of choices has influence on the simulated societies we have performed series of runs in which $`b_i`$ and $`v_i`$ were coupled through $`v_i=1b_i`$, that is, the more the agent acted as trespasser (higher $`b_i`$) the smaller was his vengefulness. This simple relation does not reflect the complexity of the issue of the exclusivity of choices, but allows an insight into possible outcomes.

The simulations returned quite interesting results, as presented in Figure 4 for SW network at two filling factors ($`\varphi =0.005`$ and $`\varphi =0.02`$), with single agent enforcement. For small $`T/P`$ stable result was $`b=0`$ and $`v=1`$. For large $`T/P`$ the opposite result was true. The width of the transition region depends on the filling factor, for larger number of neighbors the region was much narrower. The change in the value of the of $`E`$ resulted in almost direct shift of the of the step-like increase of $`b`$ as function of $`T`$ toward lower $`T`$. The onset of the norm-breaking behaviour was given by simple relation $`T_0=PE`$. This corresponds to condition when the benefit of norm breaking ($`T`$) exceeds the optimum benefits of the enforcers ($`PE`$).

The exclusive $`b`$$`v`$ model could be enhanced in a way that instead of using relation $`b=1v`$ (which allows an agent to act ‘a bit as a policeman and a bit as a thief’) one would use more complicated rule. This may be subject of separate investigation.

## 3 Norm enforcement through bounty hunting

Within Axelrod’s model vengefulness and meta-vengefulness were directly coupled. In our extension of the approach described above we have decoupled them, and studied influence of the meta-enforcement of the norm persistence. There is, however another interesting possibility, which is based on simplified observations of ‘real life’ situations, where norm enforcement, although in itself costly (through constant watchfulness cost $`E`$), may be positively reinforced, without resorting to meta-norms.

The mechanism is quite simple. Suppose the trespasser is caught. Then the agent or agents who have caught him would get a bounty $`B`$. The bounty may be set arbitrarily, but in out simulations we use the following scenario. The effective penalty $`\pi _P^{k_P}(j)`$ may be treated as a fine imposed on the trespasser $`j`$. Part of it accounts for general costs, but the rest, described through fraction $`f`$, is divided among all agents that have participated in the act of enforcement on the trespasser (or, to formulate it within the model used here, among all agents that might have participated in the act of enforcement). This allows for a strategy using bounty hunting as means of increasing the payoff of an agent, and thus increases $`v`$, which might help in keeping $`b`$ down.

Suppose we want to establish the payoff for agent $`i`$ for catching a trespassing agent $`j`$. The whole bounty $`f\pi _P^{k_P}(j)`$ has to be divided among all agents $`k`$ that are neighbors of $`j`$, in a way that relates to their ability and willingness to enforce the norm. We introduce a simple (though arbitrary) formula that preserves the relation of the individual bounties $`B(i,j)`$ to the vengefulness effort. First we introduce the summed vengefulness around agent $`j`$:

$$v_j^S=\underset{k}{\overset{(\text{neighbors of }j)}{}}v_k,$$
(19)
where summation goes over all neighbors of $`j`$. Using this it is easy to divide the bounty: the bounty for agent $`i`$ for catching trespasser $`j`$ (who’s willingness to break the norm is given by boldness $`b_j`$) is

$$B(i,j)=f\pi _P^{k_P}(j)v_i/v_j^S.$$
(20)
Thus the total bounty payoff for agent $`i`$ is given by:

$$B(i)=\underset{j}{}B(i,j)=fPv_i\underset{j}{\overset{(\text{neighbors of }i)}{}}\frac{b_j_{n=0}^{k_P}\left(1\mathrm{\Phi }_n(j)\right)}{_k^{(\text{neighbors of }j)}v_k}$$
(21)
For large number of enforcers with high vengefulness, the bounty would be divided among many agents, and thus relatively small. But suppose there is just a few enforcers — they would get significant advantage for catching the trespassers. The model has, in place of $`P^{}`$ and $`E^{}`$ only a single parameter $`f`$, describing how much of the penalty $`P`$ is divided among the vengeful agents.

The new formula for payoffs replaces the meta-enforcement terms with the bounty term. The simulation space is now simpler (two-dimensional, with variables $`b`$ and $`v`$). Example of simulation results is presented in Figure 5. Left column shows average $`b`$ and $`v`$ values as functions of $`T/P`$, right column shows results for the exclusive model discussed in previous section.

Rows correspond to increasing fractions of the penalty divided as bounty among enforcers, from $`f=0.01`$, through $`f=0.05`$ and $`f=0.4`$ to $`f=1`$. It is worth noting that as each agent has on average 5 neighbors, these fractions allow, in principle to get quite significant payoff from bounty. With our simulation we have set the enforcement cost at $`E=2`$ and penalty $`P=20`$. Thus the bounty from 5 trespassers at $`f=0.05`$ might reach 5, overcoming the cost of being vengeful.

As a result in the model with independent $`b`$ and $`v`$ only for very low $`f`$ we observe behaviour resembling the one found in MNG simulations. For larger $`f`$ both $`b`$ and $`v`$ increase simultaneously, in a smooth way. These results form the best illustration of the statement that it pays best to be a thief and policemen at the same time.

For the exclusive case the growth of $`b`$ is quite similar, with $`v`$ mirroring it according to condition $`v=1b`$. The onset of the increase of $`b`$ as function of $`T/P`$ takes place at almost the same values for both the normal and exclusive models

## 4 Conclusions

The conclusions of our models and simulations may be summarized as follows:

* decoupling norm enforcement and meta-norm enforcement ($`v`$ and $`v^{}`$), set as equal in the original model of Axelrod, has shown that the quantities evolve in very different ways, thus the original assumption was probably introducing artificial effects into results presented in ;
* for large (more than 20 agents), fully connected societies, boldness and vengefulness followed simple rules, directly relating them to net payoffs for trespassers and enforcers ($`T>P`$, $`P^{}>E`$);
* locally networked societies show complicated interplay of various parameters, with the leading role of the $`T/P`$ and $`E/P^{}`$ ratios and strong dependence on the filling factor (number of neighbours). No discernible differences were observed for different network topologies;
* in all cases we have observed local ‘pockets’ of norm breakers (high $`b`$ agents) well below the limit of $`T/P=1`$;
* a variant of the model, aimed at replicating the dichotomy of choice norm breaking/norm enforcement, in the simplest form of setting $`v=1b`$ leads to a smoothed step like behaviour of $`b`$, with the position of the step decided by $`T`$ vs $`PE`$, and its width by the number of neighbours;
* a new model, Bounty Norms Game (BNG), alternative to Meta Norms Game, based on the division of part of the penalty among norm enforcers has been presented. Despite totally different nature of the process governing the evolutionary benefits of enforcement the results resemble in general way results from MNG simulations. Simplicity and ‘real life’ background of BNG make it a good candidate for at least some of the norm enforcement studies.

## 5 Figures


text:
# Indication on the universal hadron substructure — constituent quarks.

## Acknowledgements

We are grateful to N. E. Tyurin for the interest to the work and fruitful discussions. This work was partially supported by Russian Foundation for Basic Research grant 03-02-16919.


text:
# A Millennium Scale Sunspot Number Reconstruction: Evidence For an Unusually Active Sun Since the 1940’s

## I Introduction

The sunspot number (SN) series represents the longest running direct record of solar activity, with reliable observations starting in 1610, soon after the invention of the telescope. The behaviour of solar activity in the past, before the era of direct measurements, is of importance for a variety of reasons. For example, it allows an improved knowledge of the statistical behaviour of the solar dynamo process which generates the cyclically varying solar magnetic field. It also should help to produce superior estimates of the fraction of time the Sun spends in states of very low activity, the so-called Great Minima, such as the Maunder Minimum in the second half of the 17$`^{\text{th}}`$ century. This is of particular interest when comparing the behaviour of the Sun with that of other Sun-like stars bali90 . The level of solar activity also affects the Sun’s radiative output will91 , which in turn may influence the Earth’s climate sola98 . However, any such influence takes place on time scales longer than the solar cycle, so that a statistically significant comparison with paleoclimatic records requires a long time series of solar activity data.

We are specifically interested in the past evolution of sunspot activity. Sunspots lie at the heart of solar active regions and trace the emergence of large-scale magnetic flux, which is responsible for the various phenomena of solar activity. Consequently, sunspots are a good tracer for solar magnetic activity, particularly so during times of medium to high activity.

The sunspot number record shows intriguing contrasts between the extremes reached during the Maunder Minimum when practically no sunspots were seen on the face of the Sun eddy76 ; ribe93 , and in the period since the 1940s when SN reached the average value of about 75. Extensions to earlier times have been attempted in the past by extrapolating this record, based on mathematical modelling using statistical properties of the observed SN record nago97 ; rigo01  or adjusting them to fragmentary data on naked-eye sunspot and auroral observations shov55 . Such extrapolations suffer from rapidly increasing uncertainty for earlier times. Alternatively, the SN prior to 1610 has been estimated from archival proxies, such as the concentration of cosmogenic 14C isotope in tree rings or 10Be isotope in ice cores drilled in Greenland and Antarctica beer90 ; damo91 ; beer00 . For want of a physical relationship, a simple linear regression between the SN and the isotope concentration has generally been assumed.

Recently, detailed physical models have been developed for each individual link in the chain connecting the SN with the cosmogenic isotopes. This includes a physical model relating the heliospheric magnetic flux (the Sun’s open magnetic flux) to the SN sola00 ; sola02 , a model for the transport and modulation of galactic cosmic rays within the heliosphere usos02a , and a model describing the 10Be isotope production in the terrestrial atmosphere webb03 ; masa99 . We have combined these models, such that the output of one model becomes the input for the next step. It has thus become possible to model the complete sequence of processes and to calculate the expected 10Be concentration from 1610 onwards on the basis of the SN record usos02b ; usos03 . The inversion of this chain has been successfully demonstrated as well usos03 . For artificial, noise-free 10Be data the yearly SN could be reconstructed with an error of $`\pm 10`$ compared with the typical solar maximum SN of over 100 during recent cycles. For real 10Be data, the noise makes the reconstruction of yearly sunspot numbers impractical, but robust reconstructions of 11-year (solar cycle) averages of this quantity are shown to be possible.

Follow the recent approach usos03  we present here a reconstruction of the SN since the year 850, based upon the measured 10Be concentrations in ice cores at the Dye-3 site in Greenland (annual data for 1424–1985) beer90  and at the South Pole (roughly 8-year sampled data for 850–1900) bard97 .

## II Sunspot number reconstruction

The reconstruction comprises five steps, which are described in detail by Usoskin et al. usos02b ; usos03 :

$${}_{}{}^{10}\mathrm{Be}\stackrel{(1)}{}CR\stackrel{(2)}{}\mathrm{\Phi }\stackrel{(3)}{}F_o\stackrel{(4)}{}S\stackrel{(5)}{}SN$$
(1)
First the flux of cosmic rays (CR) impinging on the Earth’s atmosphere is derived from the measured 10Be concentration. This yields the modulation strength $`\mathrm{\Phi }`$ of cosmic rays through a model of heliospheric transport of cosmic rays usos02a , which in turn is used to determine the Sun’s open magnetic flux $`F_o`$. The model sola00 ; sola02  of Solanki et al. is employed to obtain the source term $`S`$ for the open magnetic flux and, finally, the sunspot number SN. We use a realistic 10Be yield function webb03  and include also heavier species of cosmic rays, in particular $`\alpha `$-particles, whose contribution to the 10Be production is about 30%. The advantage of our physical reconstruction is that it takes into account the non-linear nature of the relation between 10Be concentration and SN and thus returns more reliable values of SN once the parameters of the physical models are fixed on the basis of actual measurements. For the period of time from 1700 to 1940 when both SN and 10Be data sets are reliable, Fig. 1 illustrates the improvement achieved with our reconstruction method (left panel) compared to a linear regression between the SN and the 10Be production rate (right panel). The distinct excursions over limited time intervals (open circles: 1722–1743, 1761–1789, 1828–1845, 1882–1894) are probably due to local climatic effects in Greenland, as indicated by comparison with the 14C data (see Fig. 2), which are much less affected by local climatic variability. Most of the time, however, the reconstructed SN is closely related to the actual SN observations note1 , with the full dots showing a considerably smaller scatter (with a linear correlation coefficient of $`r=0.96`$) than the corresponding result using a linear regression ($`r=0.81`$). While our reconstruction reproduces the SN at intermediate to high levels of solar activity relatively well, it tends to overestimate the SN during periods of low activity. This partly reflects the existence of a residual level of magnetic activity which modulates the cosmic ray flux even in times of almost vanishing sunspots beer98 ; usos01 . Such activity could, e.g., be due to the emergence of small, spotless ephemeral active regions harvey92 ; sola02 . The SN based on our physical reconstruction can therefore be considered as upper limits during periods of low SN. The assumption that 10Be is deposited locally usos03  may also contribute to the offset at small SN seen in Fig. 1. Mixing in the terrestrial atmosphere would cause the 10Be concentration to be affected by geomagnetic variations. Since the Earth’s magnetic field was stronger during the Maunder minimum than it is today, neglecting this change could lead to an overestimate of the reconstructed SN. We have used the (nearly linear) relationship between the actual and the reconstructed values on the left panel of Fig. 1 to apply a simple correction to the reconstruction, which mainly affects the low SN values.

## III Discussion

Fig. 2 shows the (1-2-1 averaged) SN reconstructed from the 8-year sampled Antarctic 10Be record for the years 850–1900 bard97  and from the Greenland Dye-3 record for the period 1424–1985 beer90 . Also given is the similarly averaged group sunspot number hoyt98  based on observations after 1610 and the (scaled) 14C concentration in tree rings, corrected for the change of the geomagnetic field stui80 ; stui89 . For easier comparison, the latter curve has been scaled to match the mean and the range of the reconstructed SN. The reconstructed SN profiles depicted by the (red and green) coloured areas are bounded from above by the actual reconstruction results and from below by the low-SN corrected results described above.

The reconstructed SN series confirms the various Great Minima and also the Medieval Maximum (roughly between the years 1100 and 1250) identified in previous, statistical studies of the 10Be and 14C records stui89 ; bard97 . The two reconstructed and the measured SN series generally are in good agreement after the end of the Maunder minimum around 1700. The differences between the results from the Antarctica and the Greenland 10Be records are greater in 1450–1700, during the so-called ‘Little Ice Age’ eddy76 , and they can possibly be ascribed to local climatic effects in 10Be deposition. This interpretation is supported by the good correlation between the reconstruction from the Antarctica data and the 14C record during this period, once the phase shift of about 20 years due to the long attenuation time for 14C bard97  has been taken into account.

The most striking feature of the complete SN profile is the uniqueness of the steep rise of sunspot activity during the first half of the 20th century. Never during the eleven centuries prior to that was the Sun nearly as active. While the average value of the reconstructed SN between 850 and 1900 is about 30, it reaches values of 60 since 1900 and 76 since 1944. For the observed group SN series since 1610 these values are 25, 61, and 75, respectively. The largest 100-year average of the reconstructed SN prior to 1900 is 44, which occurs in 1140–1240, i.e., during the Medieval Maximum, but even this is significantly less than the level reached in the last century. The Medieval Maximum is remarkable, however, in the length of time that the Sun has consistently remained at the average SN level of about 40–50. Only during the recent period of high activity since about 1830, i.e., after the Dalton minimum, has the SN remained consistently above 30 for a similar length of time. We conclude that the high level of solar activity since the 1940s is unique since the year 850. This can be considered a robust conclusion since we have shown that our reconstruction is particularly reliable in phases of high and intermediate sunspot activity, while during periods of low activity the SN may be overestimated.

The good overall agreement of the reconstructed SN with the 14C data further supports the reliability of our reconstruction: the cross-correlation coefficient (taking into account the overall 20-year delay in 14C concentration) is $`0.83\pm 0.07`$. Since the globally mixed 14C is not affected by the vagaries of the local climate, the good correspondence between the 14C curve and the SN reconstructed from the Antarctic 10Be data indicates that long-term climatic variability does not strongly affect our results. This conclusion is reinforced by the good agreement between the measured 14C concentrations and the corresponding values derived from the 10Be data on the basis of a 14C redistribution model bard97 . The fact that the reconstruction based on the Greenland Dye-3 core shows stronger deviations from 14C during the Little Ice Age suggests that the Greenland 10Be record is more strongly affected by local climate fluctuations bard97 .

It is known that the geomagnetic field has decreased by about 30% during the last 1000 years (see, e.g., baum98 ). The stronger geomagnetic field in earlier times has led to a more effective shielding of cosmic rays and may, depending on the amount of atmospheric mixing of 10Be before precipitation, have caused a reduced 10Be production. Our calculations use the present geomagnetic field and neglect the possible effect of the changing geomagnetic field. This is probably well justified at least for the Antarctica record usos03 , as indicated by the good correspondence with the 14C record, which has been corrected for changes of the geomagnetic field stui80 . In any case, without such correction our reconstruction model ascribes any effect of a stronger geomagnetic field to a higher SN in the past. Consequently, our reconstructed SN values in the pre-telescopic era are to be considered as upper bounds, emphasizing even more the exceptional nature of the high solar activity during the last 60 years.

Although our SN reconstruction still covers a rather limited length of time (but nonetheless about 3 times longer than the telescopic sunspot record), the unusually high number of sunspots during the past century suggests that we currently may be seeing a state of the solar dynamo that is uncharacteristic of the Sun at middle age. Also, the higher activity level implies more coronal mass ejections and more solar energetic particles hitting the Earth. Thus we expect that the late 20th century has been particularly rich in phenomena like geomagnetic storms and aurorae. The flux of energetic galactic cosmic rays (in the neutron monitor energy range above several GeV) reaching the Earth is presently about 10% lower than it was around 1900 usos02b . The suppression of lower-energy cosmic rays (about 2 GeV), which are mainly responsible for the production of cosmogenic isotopes, is even stronger, reaching up to 40%.

The current high level of solar activity may also have an impact on the terrestrial climate. We note a general similarity between our long-term SN reconstruction and different reconstructions of temperature mann99 ; jone01 : (1) both SN and temperature show a slow decreasing trend just prior to 1900, followed by a steep rise that is unprecedented during the last millenium; (2) Great Minima in the SN data are accompanied by cool periods while the generally higher levels of solar activity between about 1100 and 1300 correspond to a relatively higher temperature (the Medieval Warm Period) brad01 . To clarify whether this similarity reflects a real physical connection requires a more detailed study of the various proposed mechanisms for a solar influence on climate issi00 .

We thank the anonymous referee for useful comments on improving this paper.


text:
# Fermi-Bose Correspondence and Bose-Einstein Condensation in The Two-Dimensional Ideal Gas

## I Introduction

Since their inception into the theoretical framework of physics, much work has been carried out on both the Fermi-Dirac as well as the Bose-Einstein distributions Huang (1987); Pathria (1996). Much of this focuses on the behaviour of fermions and bosons at low temperatures, where the deviations in their behaviour arising from the distributions are significant and quantum effects dominate. The advent of comprehensive experimental techniques for lowering the temperature of systems to previously theoretical values over the last few decades has contributed significantly to this and revealed a wealth of unusual behaviour at extremely low temperatures, in particular the unusual phenomenon of Bose-Einstein condensation - a macroscopic occupation of the ground state of a gas of bosons at finite temperature.

For most of this time, however, there has been only a cursory and passing interest in low-dimensional Fermi and Bose systems. This has been mainly because they have been viewed as little more than abstract theoretical and mathematical constructs. Moreover, it was felt that there was little behaviour of interest to be found in these lower dimensions. It is only over the last couple of decades that there has been a growing interest in the properties of gases obeying the Fermi and Bose distributions in lower dimensions. This has been driven by firstly, the growing practical importance of low-dimensional systems, such as for example in computers, as well as new experimental techniques (such as, for example, 2D adsorbed systems, in particular, 2D adsorbed helium films Lusher et al. (1991); Siqueira et al. (1993); Casey et al. (1998); Ray et al. (1998)) that raise the possibility of being able to trap and study low-dimensional quantum systems. Typical systems being studied include superfluid and superconducting films, quantum Hall and related two-dimensional electron gases and low-dimensional trapped Bose gases. Further, the renewed theoretical interest in low-dimensional physics over the last couple of decades has led to a growing realisation that these systems possess equally fascinating and unique quantum behaviour of their own (see, for example, Ref. Schakel, 1998). In particular, recent observations of exotic quantum behavior, such as Bose-Einstein condensation (BEC) Anderson et al. (1995); Bradley et al. (1997); Davis et al. (1995) and Fermi degeneracy DeMarco and Jin (1999), in three-dimensional systems have led to a growing interest in lower-dimensional systems.

One special theoretical interest in the 2D ideal gas is due to the fact that we can obtain completely analytic expressions for the chemical potential $`\mu `$ and other thermodynamic properties. This arises, in part, from the energy-independence of the 2D density of states within the semi-classical approximation. Furthermore, in two dimensions, there is an unusual correspondence between the Bose and Fermi cases, first noted by May May (1964) and subsequently, explored further in later papers Mullin and Fernandez (2003); Toda et al. (1983); Aldrovandi (1992); Viefers et al. (1995); Lee (1997).

In his seminal paper May (1964), May noted the unexpected result that the thermal capacities of the ideal 2D Bose and Fermi gases are identical in the thermodynamic limit, i.e. their internal energies differ only by a constant temperature-independent amount. This only holds in the thermodynamic limit, or equivalently, the semi-classical approximation that the energy level spacing is negligible compared with the temperature T, allowing us to replace the discrete (quantum) structure by a continuum. Real-world systems, however, are always finite, typically of the order of $`10^310^6`$ particles for trapped gases, and more for adsorbed systems. Moreover, as we shall see later, the chemical potential $`\mu `$ has a logarithmic dependence on the total number of particles $`N`$ in the system. Thus, the thermodynamic limit is never truly reached in two dimensions and is an abstract theoretical construct, making it essential to consider finite-size effects. These may be experimentally significant, and lead to deviations from May’s Theorem.

In the following sections, we shall consider the behaviour of the uniform 2D ideal Bose and Fermi gases. In the next section, we shall briefly outline the gases in the thermodynamic limit and derive May’s Theorem before examining the Bose and Fermi gases for a finite number of particles in Section III. Though there is no BEC in the thermodynamic limit, we will demonstrate that there is a pseudo-condensation at very low temperatures for the finite 2D Bose gas. This is then used to show in Section IV that the finite-size effects lead to differing heat capacities for the Bose and Fermi gases, and a breakdown in May’s Theorem for gases with a finite number of particles. Further, the inadequacy of heat capacities as a means of determining phase transitions is discussed.

## II The Thermodynamic Limit

As the 2D ideal gas and May’s Theorem have already been examined at length before in the thermodynamic limitMay (1964); Mullin and Fernandez (2003); Toda et al. (1983); Aldrovandi (1992); Viefers et al. (1995); Lee (1997), we shall simply draw on earlier approaches to provide a brief overview. Consider $`N`$ free and identical particles of mass $`m`$ in a 2D box of ‘volume’$`A`$, with rigid-box boundary conditions. The energy values for an arbitrary particle then are

$$\epsilon _i=\frac{\mathrm{}^2k^2}{2m},k^2=\frac{\pi ^2}{A}\left(n_x^2+n_y^2\right),$$
(1)
where $`k`$ is the momentum characterising the energy states $`\epsilon _i`$, $`i,n_x,n_y=1,2,3\mathrm{}`$ and $`(n_x,n_y)`$ are points in the positive quadrant of a unit square lattice. In the thermodynamic limit, the ground state energy tends to 0. For a finite system, however, due to the Heisenberg uncertainty principle, this is always greater than 0, as $`(n_x,n_y)=(1,1)`$. So, for simplicity, we shift the energy scale, taking the ground state energy as 0,

$$\epsilon _i=\frac{\mathrm{}^2k^2}{2m}\epsilon _0=\frac{\pi ^2\mathrm{}^2}{2mA}\left(n_x^2+n_y^22\right).$$
(2)
Assuming that the system is in thermal equilibrium at temperature $`T`$ and chemical potential $`\mu `$, then the mean number of particles $`n\left(\epsilon _i\right)`$ at a given energy $`\epsilon _i`$ is given by

$`n\left(\epsilon _i\right)={\displaystyle \frac{1}{e^{\left(\epsilon _i\mu \right)/k_BT}\pm 1}}.`$ (3)
Here and later, the upper (lower) signs are for fermions (bosons). The total number of particles N and the internal energy E of the system are now given by

$`N={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}n\left(\epsilon _i\right),`$ (4)
$`E={\displaystyle \underset{i=1}{\overset{\mathrm{}}{}}}\epsilon _in\left(\epsilon _i\right).`$ (5)
where the sum is taken over the single particle energy states. In the semiclassical approximation, we assume that the spacing between the energy levels of the particle states is negligible. This allows us to replace the sum with an integral over energy states using the density of states $`g(\epsilon )`$ and the spin degeneracy $`\alpha `$,

$`f=\alpha {\displaystyle _0^{\mathrm{}}}f\left(\epsilon \right)g\left(\epsilon \right)n\left(\epsilon \right)𝑑\epsilon ,`$ (6)
where $`f\left(\epsilon \right)`$ is a general single-particle function of energy. However, this approximation is only exact in the thermodynamic limit when the system size is infinite, and at very high temperatures $`T`$. As we shall see in Section III, for finite systems at low T, these assumptions are invalid and the discreteness of the energy levels cannot be ignored. In other words, we will be at a temperature $`T`$ below which the semiclassical approximation breaks down and modifications need to be made. However, it may be possible to compensate by modifying the semi-classical approximation to allow for finite-size effects.

We now evaluate the various thermodynamic quantities for arbitrary spin, i.e. arbitrary $`\alpha `$, for both the 2D Fermi and Bose gases in the thermodynamic limit. In 2D, $`g(\epsilon )`$ and thereby, the Fermi energy $`\epsilon _F`$ are given by,

$`g(\epsilon )={\displaystyle \frac{mA}{2\pi \mathrm{}^2}},\epsilon _F={\displaystyle \frac{2\pi \mathrm{}^2N}{\alpha mA}}.`$ (7)
For comparison, we adopt $`\epsilon _F`$ as a characteristic energy $`\epsilon _C`$ for the Bose gas also, allowing us to express later results in a unified form. In 2D, as $`g(\epsilon )`$ is a constant, we have a simple integral for N,

$$N=\frac{\alpha mA}{2\pi \mathrm{}^2}_0^{\mathrm{}}\frac{d\epsilon }{e^{\left(\epsilon \mu \right)/k_BT}\pm 1}.$$
(8)
This may be expressed using the Fermi-Dirac and Bose-Einstein functions Huang (1987); Pathria (1996), $`f_\nu (z)`$ and $`g_\nu (z)`$ respectively,

$`f_\nu \left(z\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\nu \right)}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{x^{\nu 1}dx}{z^1e^x+1}}`$ (9)
$`g_\nu \left(z\right)={\displaystyle \frac{1}{\mathrm{\Gamma }\left(\nu \right)}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{x^{\nu 1}dx}{z^1e^x1}},`$ (10)
where $`\mathrm{\Gamma }\left(\nu \right)`$ is the Euler gamma function. These can also be expressed using the polylogarithmic functions $`\text{Li}_\nu (z)`$ Abramowitz and Stegun (1964); Lewin (1958),

$`f_\nu \left(z\right)=\text{Li}_\nu \left(z\right),`$ (11)
$`g_\nu \left(z\right)=\text{Li}_\nu \left(z\right).`$ (12)
For N above, $`\nu =1`$ and $`z`$ is the fugacity ($`z=e^{\mu /k_BT}`$). Solving (8) in terms of $`\epsilon _C`$ and rearranging, we can now express the fugacity $`z`$ and thereby, the chemical potential $`\mu `$ as a completely analytic expression

$$z_{F,B}=e^{\mu /k_BT}=\left(1e^{\pm \epsilon _C/k_BT}\right),$$
(13)
where $`z_F`$ ($`z_B`$) is the fugacity for the ideal Fermi (Bose) gas. As $`T0`$, $`z_B`$ smoothly approaches 1. Denoting the fugacity $`z_B`$ by $`z_B=e^\alpha `$, as $`z_B1`$,

$$g_1\left(e^\alpha \right)=\mathrm{ln}\left(1e^\alpha \right)\underset{\alpha 0}{}\mathrm{ln}\alpha .$$
(14)
N diverges logarithmically as $`\alpha 0`$. Thus, it has no upper bound, and there is no temperature below which the ground state can be said to be macroscopically occupied in comparison to the excited states. Thus, there is no BEC in the thermodynamic limit. This is a qualification of the oft-paraded statement that there is no BEC for $`D2`$, where $`D`$ is the number of dimensions. However, as we shall see in Section III, the logarithmic divergence implies that at sufficiently low T, bosons still crowd into low-lying states to give a pseudo-condensate.

The internal energy may be similarly solved,

$`E_F=Nk_BT\left({\displaystyle \frac{k_BT}{\epsilon _C}}\right)f_2\left(e^{\epsilon _C/k_BT}1\right),`$ (15)
$`E_B=Nk_BT\left({\displaystyle \frac{k_BT}{\epsilon _C}}\right)g_2\left(1e^{\epsilon _C/k_BT}\right).`$ (16)
Now, using the following property of the dilogarithmic functionMay (1964); Lewin (1958) $`\text{Li}_2\left(x\right)`$,

$$f_2\left(\frac{1}{x}1\right)=g_2\left(1x\right)+\frac{1}{2}\left(\mathrm{ln}x\right)^2,$$
(17)
we can easily show

$`E_F=E_B+{\displaystyle \frac{1}{2}}N\epsilon _C.`$ (18)
i.e. for all $`T`$, $`E_B`$ differs from $`E_F`$ by a constant temperature independent energy, which does not contribute to the heat capacity. This gives May’s Theorem May (1964), which states that the heat capacities of the ideal 2D Fermi and Bose gases are identical at all T and N,

$$C_{V,Bose}=\frac{E_B}{T}|__V=\frac{E_F}{T}|__V=C_{V,Fermi}.$$
(19)
The above result indicates that the 2D ideal Bose gas does not undergo a phase transition, as any experimental measurement of the heat capacity would be the same as for a 2D Fermi gas, and would lack the characteristic bump that indicates the presence of BEC in the 3D Bose gas. However, as we shall see, this is because we have failed to consider the finite nature of any real systems.

## III The Finite Gas

The case of the finite ideal Bose gas in two dimensions has already been considered to varying degrees in earlier papers Osborne (1949); Ziman (1953); Mills (1964); Goble and Trainor (1967); Chester (1968); Krueger (1968); Irnry (1969); Barber (1983); Grossmann and Holthaus (1995), and more recently, by Pathria Pathria (1998). Most have commented that it may be possible to have a pseudo-condensation in the Bose gas in two dimensions, with the transition temperature typically going as $`T_C1/\mathrm{ln}N`$. Here, we present a simpler and more comprehensive analysis of the problem that not only offers far more insight into how the finite size of the system leads to the breakdown of May’s Theorem and a pseudo-condensation in the Bose gas, but further, in Section IV, discusses how the heat capacities for the Fermi and Bose gases differ for finite systems, and how such a pseudo-condensation may be experimentally measured.

At $`T=0`$, the 2D Bose gas is in its ground state, and we have the trivial case of ground state macroscopic occupation. Extrapolating, however, a range of finite temperatures must exist near $`T=0`$, which may depend on $`N`$, where the number of particles in the ground state is also of $`O\left(N\right)`$ \- a macroscopic occupation. At very low $`T`$, this can be a significant fraction of $`N`$, and the semi-classical approach is now flawed as $`N`$ is no longer continuous. Taking macroscopic occupation of the ground state as the defining characteristic of BEC, we consider the ground state separately analogous to standard treatments of the 3D Bose gas Huang (1987); Pathria (1996). For finite systems, the discrete level structure requires us to introduce a cutoff energy at the first excited level $`\epsilon _1`$. The integral is now bounded and we can remove the logarithmic divergence. Thus, we can have a pseudo-BEC in a finite two-dimensional system. This must be associated with a breakdown of May’s theorem at sufficiently low $`T`$ for finite systems. We thus write (8) for the Bose gas as

$$N=N_0+\frac{\alpha mA}{2\pi \mathrm{}^2}_{\epsilon _1}^{\mathrm{}}\frac{d\epsilon }{e^{\left(\epsilon \mu \right)/k_BT}1},$$
(20)
It should be noted that the above approximation is quite a simple and crude one, as it is only the leading term of the Euler-Maclaurin summation formula. However, it suits the purposes of our paper as our intention here is merely to highlight how even the simplest approximations show that there is a pseudo-condensation in the finite Bose gas. A more detailed approximation will only lead to minor modifications to the core results presented in the rest of this paper. Further, the leading correction term to the integral above in the Euler-Maclaurin expansion is of the order of $`N_1=n\left(\epsilon _1\right)=\left(e^{\left(\epsilon _1\mu \right)/k_BT}\right)^1`$. Taking the ratio $`N_0/N_1`$ of the number of particles in the ground state to that in the first excited state, we can see

$$\frac{N_0}{N_1}\frac{\epsilon _1\mu }{\mu }$$
(21)
As the chemical potential $`\mu `$ for the 2D Bose gas is always negative and less than $`0`$, the above ration is always greater than 1. Further, as $`T0`$, $`\mu 0`$ and the above expression clearly goes to $`\mathrm{}`$. Further, a more detailed evaluation of (21) gives us

$$\frac{N_0}{N_1}\frac{3\pi ^2\mathrm{}^2}{2mAk_BT}e^{\frac{2\pi \mathrm{}^2N}{\alpha mAk_BT}}$$
(22)
This clearly rapidly goes to infinity as both the temperature $`T`$ as well as the number of particles $`N`$ increase. Moreover, for the regime in question, we are interested in very low temperatures where from the above equation (22), the ratio $`N_0/N_1`$ is clearly significantly large.

Thus, we are justified in making our crude approximation and denoting the number of particles $`N`$ by (20) above. Solving as before, we now have

$$NN_0=\frac{\alpha mAk_BT}{2\pi \mathrm{}^2}\mathrm{ln}\left(1z_Be^{3\pi ^2\mathrm{}^2/2mAk_BT}\right)$$
(23)
and $`\epsilon _1`$ prevents the logarithmic divergence. This is because the excited states cannot accommodate all the particles at low $`T`$, and the excess are forced into the ground state, leading to a pseudo-BEC. Equivalently, changing variables $`\epsilon \epsilon +\epsilon _1`$ and resetting the integral from zero, we may also express (20) in terms of the Bose-Einstein functions $`g_\nu \left(z\right)`$ as

$$NN_0=\frac{\alpha mAk_BT}{2\pi \mathrm{}^2}g_1\left(z_Be^{3\pi ^2\mathrm{}^2/2mAk_BT}\right).$$
(24)
Thus, we now have an effective fugacity

$`z_Bz_B^{}=z_Be^{3\pi ^2\mathrm{}^2/2mAk_BT}`$
$`=e^{3\alpha \pi \epsilon _C/4Nk_BT}\left(1e^{\epsilon _C/k_BT}\right).`$ (25)
where $`z_B`$ is as before. This summarizes the finite-size effects on the Bose gas. The extra term is dependent on the system size and tends to 1 as $`N\mathrm{}`$. For finite $`N`$ and low $`T`$, it is less than 1, and so, the fugacity is modified. Thus, there is a pseudo-BEC here but none in the thermodynamic limit. This macroscopic occupation is a purely boson phenomenon, and so, the behaviour will be markedly different from that of the two-dimensional Fermi gas.

For the Fermi gas, Pauli’s principle forbids macroscopic occupation of any state. Thus, the first term is not divergent. $`z_F`$ (Eq. (13)) still holds, and finite-size effects only appear for the Bose gas, where the excited states take as many particles as possible at low $`T`$, with the excess forced into the ground state. This occurs at the maximal value of the fugacity $`z_B=1`$. Thus, we can define a characteristic temperature $`T_C`$ that marks the onset of condensation where $`N_00`$ and $`z1`$,

$`N{\displaystyle \frac{\alpha mAk_BT_C}{2\pi \mathrm{}^2}}\mathrm{ln}\left(1e^{3\pi ^2\mathrm{}^2/2mAk_BT_C}\right)`$
$`{\displaystyle \frac{Nk_BT_C}{\epsilon _C}}\mathrm{ln}\left({\displaystyle \frac{4Nk_BT_C}{3\pi \alpha \epsilon _C}}\right)`$ (26)
for very low $`T`$. We may solve (26) to give

$$T_C=\frac{\epsilon _C}{k_B}\frac{1}{W\left(4N/3\pi \alpha \right)},$$
(27)
where, $`W\left(x\right)`$ is the Lambert function Corless et al. (1996) and is the principal solution for $`w`$ in the equation $`x=we^w`$. This is plotted in Fig. 1. It may be viewed as a generalisation to a logarithm. Thus, our result (27) is a refinement of the earlier noted result that the transition temperature behaves like $`T_C1/\mathrm{ln}N`$. As $`N\mathrm{}`$, the reduced characteristic temperature $`k_BT_C/\epsilon _C0`$. Thus, in the thermodynamic limit, we recover the conventional result that there is no BEC in two dimensions. The condensate fraction is

$`{\displaystyle \frac{N_0}{N}}1{\displaystyle \frac{k_BT}{\epsilon _C}}\mathrm{ln}\left({\displaystyle \frac{4Nk_BT}{3\pi \alpha \epsilon _C}}\right)`$
$`1{\displaystyle \frac{T}{T_C}}\left[1+{\displaystyle \frac{k_BT_C}{\epsilon _C}}\mathrm{ln}\left({\displaystyle \frac{T}{T_C}}\right)\right].`$ (28)
This only holds for very low temperatures $`T<T_C`$. As can be seen in Fig. 2, the condensate occurs at lower and lower $`T`$ as the system size $`N`$ increases. In the thermodynamic limit, there is no BEC at all.

Having shown the pseudo-condensation, we now move on to consider how the Fermi-Bose correspondence predicted by May is affected when the gases are assumed finite.

## IV The Finite Gas: Internal Energy and Heat Capacity

The internal energies may be similarly found. The first term in (5) for the Bose gas is divergent as it is the internal energy $`E_{gs}`$ associated with the anomalously large ground state. Separating this out,

$$E_B=E_{gs}+\underset{i=1}{\overset{\mathrm{}}{}}\frac{\epsilon _i}{e^{\left(\epsilon _i\mu \right)/k_BT}1},$$
(29)
As we have set $`E_{gs}`$ as the zero-point energy, the first term is 0. Changing variables $`\epsilon \epsilon +\epsilon _1`$, and converting to an integral,

$`E_B={\displaystyle \frac{N}{\epsilon _C}}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\epsilon +\epsilon _1}{e^{\left(\epsilon +\epsilon _1\mu \right)/k_BT}1}}𝑑\epsilon `$
$`={\displaystyle \frac{N}{\epsilon _C}}\left(k_BT\right)^2g_2\left(z_B^{}\right)+{\displaystyle \frac{\epsilon _1}{\epsilon _C}}Nk_BTg_1\left(z_B^{}\right).`$ (30)
where $`z_B^{}`$ is defined by (25). We again have an effective fugacity $`z_B^{}`$ and additional terms, which are negligible in the thermodynamic limit. By contrast, for the finite Fermi gas, there is no divergence in the ground state term, and we have the same expression for the fugacity (13) as before. The internal energies of the finite two-dimensional Fermi and Bose gases are then given by

$`E_F`$ $`=`$ $`Nk_BT\left({\displaystyle \frac{k_BT}{\epsilon _C}}\right)f_2\left(e^{\epsilon _C/k_BT}1\right)`$ (31)
$`E_B`$ $`=`$ $`\left({\displaystyle \frac{Nk_B^2T^2}{\epsilon _C}}\right)g_2\left(z_B^{}\right){\displaystyle \frac{3\pi \alpha }{4}}k_BT\mathrm{ln}\left(1z_B^{}\right)`$ (32)
Eq. (18) is no longer true, as we cannot simply express $`E_F`$ in terms of $`E_B`$ as before. $`E_B`$ now has additional $`T`$-dependent terms that lead to a differing heat capacity.

May’s Theorem is thus a special limiting case for the 2D Bose and Fermi internal energies in the thermodynamic limit. Here, the additional factor $`e^{3\pi \alpha \epsilon _C/4Nk_BT}`$ tends to unity in the thermodynamic limit, recovering May’s Theorem. The same occurs in the high $`T`$ limit for a finite gas as above $`T_C`$, there is no macroscopic ground state occupation. Thus, above $`T_C`$, we may approximate the internal energies by (15) and (16) with differences becoming marked only at very low $`T`$ and/or very small system size $`N`$.

We now calculate the Bose and Fermi heat capacities by using the following polylogarithmic identity Pathria (1996); Abramowitz and Stegun (1964),

$$\frac{d}{dz}Li_{n+1}\left(z\right)=\frac{1}{z}Li_n\left(z\right)$$
(33)
The heat capacity for the finite gas is then given by

$`C_{V,Fermi}={\displaystyle \frac{Nk_B^2}{\epsilon _C}}\left[2Tg_2\left(1e^{\epsilon _C/k_BT}\right){\displaystyle \frac{\epsilon _C^2}{k_B^2T}}{\displaystyle \frac{1}{e^{\epsilon _C/k_BT}1}}\right]`$ (34)
$`C_{V,Bose}={\displaystyle \frac{Nk_B^2}{\epsilon _C}}\left[2Tg_2\left(z_B^{}\right)+{\displaystyle \frac{\epsilon _C}{k_B}}{\displaystyle \frac{\mathrm{ln}\left(1z_B^{}\right)}{e^{\epsilon _C/k_BT}1}}\right]`$
$`+{\displaystyle \frac{3\pi \alpha k_B}{16N}}\left[{\displaystyle \frac{\epsilon _C\left(3\pi \alpha z_B^{}4Ne^{\epsilon _C/k_BT3\alpha \pi \epsilon _C/4Nk_BT}\right)}{k_BT\left(1z_B^{}\right)}}8N\mathrm{ln}\left(1z_B^{}\right)\right]`$ (35)
The Bose heat capacity $`C_{V,Bose}`$ now contains several terms dependent on both $`N`$ and $`T`$. It should also be noted that the last two terms in the heat capacity are very much smaller than the first two terms and become increasingly negligible as the system size increases. At high $`T`$ and/or very large $`N`$, the additional terms tend to zero, i.e. as $`N\mathrm{}`$ and/or $`T\mathrm{}`$, $`\pi \alpha \epsilon _C/Nk_BT0`$, and (35) reduces to the standard expression in the thermodynamic limit,

$`C_{V,Bose}{\displaystyle \frac{2Nk_B^2T}{\epsilon _C}}g_2\left(1e^{\epsilon _C/k_BT}\right){\displaystyle \frac{N\epsilon _C}{T}}{\displaystyle \frac{1}{e^{\epsilon _C/k_BT}1}}=C_{V,Fermi}.`$ (36)
This is one of the most important results of our paper, and shows clearly that May’s theorem arises as the thermodynamic limit of (34) and (35), when the heat capacities of the Bose and Fermi gases become equal. Nevertheless, for most real-world 2D systems, above the characteristic temperature $`T_C`$, which marks the onset of Bose-Einstein condensation in the two-dimensional ideal Bose gas, $`\pi \alpha \epsilon _C/Nk_BT`$ is very small and so, even for finite gases, the thermal capacities of the Bose and Fermi systems are approximately equal above the condensation temperature.

In Fig. 3, we plot the reduced heat capacities ($`C_V/Nk_BT^{}`$) for the Fermi and Bose gases against the reduced temperature $`T^{}`$(where $`T^{}=k_BT/\epsilon _C`$) for varying system size $`N`$, using the full expression for both $`C_{V,Bose}`$ and $`C_{V,Fermi}`$. As the heat capacity is approximately linear close to $`T=0`$ in the thermodynamic limit, we would expect the plots for the 2D Fermi gas to tend to a horizontal line as $`T0`$. This limiting value will be at the critical value when $`T^{}=0`$, i.e.

$`{\displaystyle \frac{C_{V,Fermi}}{Nk_BT^{}}}=2g_2\left(1\right)=2\zeta \left(2\right)={\displaystyle \frac{\pi ^2}{3}}.`$ (37)
In contrast, for the finite 2D Bose gas, we expect to see a deviation from the above mentioned low-temperature behaviour. In particular, we expect that the reduced heat capacity ($`C_V/Nk_BT^{}`$) will fall to 0 following the formation of the condensate. Our plots indicate distinct differences between the reduced heat capacities for the Bose and Fermi gases, particularly at very low $`T`$ in the vicinity of and below the condensation temperature $`T_C`$. Then, $`C_{V,Bose}`$ is less than $`C_{V,Fermi}`$ and starts to deviate significantly around $`T_C`$, before rapidly going to zero as we approach $`T=0`$. This is in marked contrast to the 3D Bose gas, where the conventional heat capacity $`C_V/Nk_B`$ peaks at some finite $`T`$. In the two-dimensional case, however, the heat capacity $`C_V/Nk_B`$ is a smooth function of $`T`$, in keeping with the size-dependent nature of the condensation. Here, instead, we see a “bump”in the plot for the reduced heat capacity $`C_V/Nk_BT^{}`$ of the 2D Bose gas, indicating the presence of a condensation. Further, as N increases, the behaviour tends to that of the 2D Fermi gas as the condensation forms at lower temperatures, and consequently, the fall-off of the reduced heat capacity occurs at lower temperatures. Indeed, with $`N=500000`$, the two plots of the Fermi and Bose gases are virtually indistinguishable. This is significant as conventionally, the heat capacity $`C_V/Nk_B`$ is used experimentally to determine if a Bose-Einstein condensation has occurred. However, as can be seen from above, the heat capacity does not necessarily peak for a weak transition which is dependent on the system size, and we have to examine alternate expressions to determine the presence of a condensation.

## V Discussion and Conclusion

In this paper, we set out to discuss the two-dimensional ideal Bose and Fermi gases in the absence of any trapping potentials, and both in the thermodynamic limit and the finite case. In particular, we were interested in May’s Theorem, and its breakdown for finite systems.

As can be seen above in Sections III and IV, the 2D finite Fermi system displays the same behaviour as in the thermodynamic limit but the Bose gas deviates significantly for a finite system and undergoes a pseudo-Bose-Einstein condensation. It is not a “true”condensation as for the 3D Bose gas, where there is a distinct phase transition in the thermodynamic limit and a peak in the heat capacity. Rather, here we have a characteristic condensation temperature $`T_C`$ that is dependent on the system size and tends to zero in the thermodynamic limit. Moreover, though these is no heat capacity peak here, it still deviates significantly from the Fermi heat capacity below $`T_C`$ with the deviations tending to zero as the system size and/or the temperature increase. In other words, we have a peak in the reduced heat capacity, as shown in Fig. 3. This raises the question of what exactly is defined as a phase transition. The two-dimensional Bose gas undergoes a significant change at low temperatures but only in the finite case. Thus, it is not what may be termed a “conventional”phase transition as for the three-dimensional Bose gas, since it ceases to exist in the thermodynamic limit. However, from an experimental perspective, a system is never truly in the thermodynamic limit, and so, what we have termed a “pseudo-condensation”is expected. This also, therefore, casts doubt on the use of the heat capacity as a means of identifying phase transitions, as discussed above.

As a final note, we would like to consider some of the outstanding problems. The above analysis is a simple one and deals with only the ideal two-dimensional gas. There is still considerable debate whether such behaviour is to be found for an interacting two-dimensional gas. Some authorsPetrov et al. (2000) have claimed that there is a true condensation for the finite interacting 2D system at sufficiently low temperatures as the coherence length becomes larger than the condensate size. Others, however, have claimed that any phase transition would disappear in the presence of particle interactionsMullin (1997), as predicted by the Hohenberg theorem. Clearly, the presence of interactions will modify our results greatly. As yet, however, there is no indication if there is any correspondence between the two-dimensional Fermi and Bose gases, when they are interacting, and further, if a pseudo-condensation is to be found for the finite Bose gas. Further, similar analogues exist in other dimensions for ideal gases in a trapping potentialKetterle and van Druten (1996); Mullin (1997). It would be interesting to debate whether there are deeper implications due to the Fermi-Bose correspondence, and where such an situation occurs, it is a special unique case. This is an area of active work, with new experimental innovations making the resolution of the above problems and the elucidation of two-dimensional gases only a matter of time now.