text:
# Reconciling Conductance Fluctuations and the Scaling Theory of Localization

## Abstract

We reconcile the phenomenon of mesoscopic conductance fluctuations with the single parameter scaling theory of the Anderson transition. We calculate three averages of the conductance distribution: $`\mathrm{exp}\left(\mathrm{ln}g\right)`$, $`g`$ and $`1/R`$ where $`g`$ is the conductance in units of $`e^2/h`$ and $`R=1/g`$ is the resistance and demonstrate that these quantities obey single parameter scaling laws. We obtain consistent estimates of the critical exponent from the scaling of all these quantities.

In the original proposal of the scaling theory of localization the conductance $`g`$ (in units of $`e^2/h`$) is the relevant parameter . The scaling of this parameter is deduced by looking at the limiting cases of a good metal, where weak localization theory is applicable, and the strongly localized regime where an exponential dependence on system size is expected. A smooth monotonic interpolation is supposed between the two limits. The central equation of the theory is the $`\beta `$ function

$$\beta (g)\frac{\mathrm{d}\mathrm{ln}g}{\mathrm{d}\mathrm{ln}L}$$
(1)
which describes how the conductance of an $`L\times L\times L`$ cube is renormalized with system size $`L`$.

The scaling theory has been very influential. It is the basis for the predictions that the Anderson transition is continuous and that the lower critical dimension for the Anderson transition is two.

Shortly after its proposal, the pioneering numerical calculations of Pichard , and MacKinnon and Kramer  provided indirect support for the scaling theory. In these calculations the localization length of electrons on quasi- one dimensional bars was calculated. The theory of finite size scaling was then applied to deduce the critical parameters from the dependence of the localization length on the transverse dimension of the bars. This has proved useful for making quantitative estimates of critical parameters .

It became clear later, in particular after the investigation of universal conductance fluctuations in mesoscopic systems, that the conductance of a disordered system is a random variable . In the vicinity of the Anderson transition the fluctuations in $`g`$ are of the same order of magnitude as it’s mean value . The smooth scaling behavior predicted by (1) is clearly inconsistent with the fluctuating behavior of the conductance which occurs in practice and casts strong doubts over the soundness of the scaling theory . There seem to be two principal remedies which we discuss in turn .

First, we could attempt to establish that the distribution of conductance $`p_L(g)`$ (in the limit that the system size $`L`$ and the correlation length $`\xi `$ are much longer than any microscopic lengths) obeys a single parameter scaling law. The precise meaning of this statement is that it should be possible to parameterize the bulk of the conductance distribution with a single parameter. If we denote this parameter by $`X`$ then the bulk of the conductance distribution must be of the form

$$p_L(g)F(g;X(L))$$
(2)
and the parameter $`X`$ must obey a single parameter scaling law

$$\frac{\mathrm{d}\mathrm{ln}X(L)}{\mathrm{d}\mathrm{ln}L}=\beta (X)$$
(3)
Note that $`X`$ need not necessarily be one of the moments of the distribution. Indeed, single parameter scaling of a distribution does not necessarily imply scaling of its moments. The moments may be dominated by non- universal tails of the distribution or might not even exist. We shall refer to this first possibility as strong single parameter scaling.

Second, we could attempt to establish that some typical or average conductance obeys a single parameter scaling law. This is a somewhat weaker statement since one parameter scaling of some average or typical quantity does not imply single parameter scaling of the corresponding distribution. For example, it might happen that several independent parameters are needed to describe the distribution. For this reason we shall refer to this as weak single parameter scaling.

Our purpose in this paper is to firmly establish single parameter scaling in the weaker sense given above, for systems near the critical point of the Anderson transition in three dimensions. A secondary objective is to identify which averages or typical values obey single parameter scaling laws. By simulating the conductance distribution for large ensembles of disordered $`L\times L\times L`$ cubes we have achieved both objectives.

We supposed perfect leads were attached to a pair of opposite sides of each cube and used the Landauer formula to relate the (Landauer) conductance $`g_L`$ to the transmission matrix $`t`$ which describes the transmission of electrons from one lead to the other .

$$g_L=2\mathrm{t}rtt^{}$$
(4)
We use the notation $`g_L`$ to emphasize that this is the conductance that would be measured in a two probe measuring geometry. In the original work on the scaling theory the Thouless conductance $`g=E_C/\mathrm{\Delta }`$ was considered where $`E_c`$ is the Thouless energy and $`\mathrm{\Delta }`$ the mean energy level spacing. The Landauer conductance and the Thouless conductance are not completely equivalent. Their relation has been considered in detail by Braun et al  who state that the contact resistance, which is always present in a two terminal measurement, should be subtracted from $`g_L`$ . Hence we study the statistics of

$$\frac{1}{g}=\frac{1}{g_L}\frac{1}{2N}$$
(5)
Here $`NN(E_F)`$ is the number of propagating channels in the contacts at Fermi energy $`E_F`$ and $`1/2N`$ is the contact resistance appropriate for the situation we have simulated.

The motion of the electrons in the system is described by the Anderson model

$$H=V\underset{<i,j>}{}C_i^{}C_j+\underset{i}{}W_iC_i^{}C_i,$$
(6)
where $`C_i^{}(C_i)`$ denotes the creation (annihilation) operator of an electron at the site $`i`$ of a 3D cubic lattice. The amplitude of the random potential at site $`i`$ is $`W_i`$. Hopping is restricted to nearest neighbors and its amplitude is taken as the unit of energy, $`V=1`$. We supposed a box distribution with each $`W_i`$ uniformly distributed on the interval $`[W/2,W/2]`$. Previous work has verified the universality of the critical behavior in this model with respect to the choice of distribution of the random potential . We imposed fixed boundary conditions in the transverse directions since we have found in previous work that corrections to scaling vanish more quickly with system size in this case . The Hamiltonian (6) has both time reversal and spin rotation symmetries so that the observed critical behavior should be that of the orthogonal universality class. We used the method of Pendry et al  to calculate the transmission matrix which appears in (4). We set the Fermi energy $`E_F=0.5V`$ and for each combination of disorder $`W`$ and system size $`L`$ we generated an ensemble of $`1,000,000`$ samples (except for $`L=16`$ where 500,000 were generated). This allowed us to estimate the various averages of the conductance to roughly an accuracy of $`0.1\%`$.

We examined the behavior of three different averages of the conductance distribution: $`\mathrm{exp}\left(\mathrm{ln}g\right)`$, $`g`$ and $`1/R`$ where $`R=1/g`$ is the resistance. Note that for each average it is possible to define a different $`\beta `$ function. For example, for the typical conductance $`\mathrm{exp}\left(\mathrm{ln}g\right)`$

$$\beta (\mathrm{exp}\left(\mathrm{ln}g\right))\frac{\mathrm{d}\mathrm{ln}g}{\mathrm{d}\mathrm{ln}L}.$$
(7)
In the critical region the conductance fluctuations are of the same order of magnitude as the mean conductance so that these averages are not at all equivalent.

We fitted the disorder and system size dependence of the averages to the standard scaling form. Taking the typical conductance as an example we supposed that

$$\mathrm{ln}g=F(\psi L^{1/\nu },\varphi L^y),$$
(8)
where $`\psi `$ is a relevant scaling variable and $`\varphi `$ is an irrelevant scaling variable which allows us to take account of corrections to scaling. Such deviations from perfect scaling are always present in a simulation of a finite system and it is necessary to have some means of accounting for them. We approximated this scaling function by its first order expansion in the irrelevant scaling variable and fitted the numerical data to the form

$$\mathrm{ln}g=F_0(\psi L^{1/\nu })+\varphi L^yF_1(\psi L^{1/\nu }).$$
(9)
The scaling variables were approximated by expansions in terms of the dimensionless disorder $`w=(W_cW)/W_c`$ where $`W_c`$ is the critical disorder separating the insulating and metallic phases.

$$\psi =\psi _1w+\psi _2w^2,\varphi =\varphi _0.$$
(10)
The critical exponent $`\nu `$ describes the divergence of the localization (correlation) near the critical point

$$\xi =\xi _0\left|\psi \right|^\nu .$$
(11)
The absolute scale of the localization length $`\xi _0`$ cannot be determined from this fit. The system size dependence of the irrelevant scaling variable is described by an exponent $`y<0`$. The functions $`F_0`$ and $`F_1`$ were expanded to third order in $`w`$. We found that this fitting scheme was the simplest which still allowed for goodness of fit probabilities in excess of $`0.1`$.

The results of the analysis are displayed in Table I and Figures 1-4. A number of points can be noticed. First, acceptable fits to the single parameter scaling law are obtained for all three averages. Second, mutually consistent values for the critical exponent are obtained and these estimates are also consistent with estimates based on the scaling of the localization length of electrons on bars . Third, mutually consistent values of the critical disorder are obtained. Finally, the critical values of the quantities considered vary widely, indicative of the broadness of the distribution of the conductance fluctuations in the critical regime. This means that the precise form of the $`\beta `$ function is different for each average, though all should have the same slope $`(=1/\nu )`$ at their respective zeros . Note that the $`\beta `$ function will also depend on the boundary conditions since it is known that critical conductance distribution depends strongly on the boundary conditions even in the limit that $`L\mathrm{}`$ .

These results firmly establish that the typical conductance, the mean conductance and the mean resistance all obey single parameter scaling in the critical regime.

To what extent do these results also support scaling in the strong sense? Let us consider the vicinity of the metallic, insulating and critical fixed points in turn.

The metallic fixed point can be reached by taking the limit $`L\mathrm{}`$ with $`\xi `$ fixed from any starting point on the metallic side of the transition. In this limit $`p_L(g)`$ approaches a normal distribution with a size independent variance . Therefore a single parameter, $`g`$, is sufficient to parameterize the distribution. Further it can be established that $`g`$ obeys a single parameter scaling law using weak localization theory near the metallic fixed point .

The insulating fixed point can be reached by taking the limit $`L\mathrm{}`$ with $`\xi `$ fixed from any starting point on the insulating side of the transition. In the insulating regime exact results are available only for one dimensional systems . It is found that distribution of conductance is log- normal with a width which is related to the mean of the distribution (provided certain conditions are satisfied .) Hence a single parameter $`\mathrm{ln}g`$ is sufficient to parameterize the distribution. Also $`\mathrm{ln}g`$ obeys a single parameter scaling law in the limit $`L\xi `$. Calculations of the average conductance and the average resistance, in this limit, in strictly one dimension, give

$$g=\frac{1}{2}\left(\frac{\pi \xi }{2L}\right)^{3/2}\mathrm{exp}\left(L/2\xi \right)$$
(12)
and

$$R=\frac{1}{2}\mathrm{exp}\left(4L/\xi \right)$$
(13)
showing that both these quantities also obey single parameter scaling laws. However, the average conductance (and also its higher moments) are not consistent with the log- normal form of the bulk of the conductance distribution. They turn out to be determined by the probability of the occurrence of rare resonance states with anomalously large conductances which are not described by the the log- normal distribution. It is possible, but not certain, that these results are qualitatively correct for the insulating fixed point in three dimensions.

At the critical fixed point $`\xi `$ diverges and the single parameter scaling hypothesis predicts a scale independent universal critical conductance distribution $`p_c(g)`$. This has been confirmed in numerical simulations .

Thus, it seems likely that the conductance distribution also obeys single parameter scaling in the strong, as well as the weak, senses. However, the available results are not completely conclusive. It seems to us that a demonstration that $`p_L(g)`$ obeys $`(\text{2})`$ and $`(\text{3})`$ in the limit $`L\mathrm{}`$ and $`\xi \mathrm{}`$ with $`L/\xi `$ fixed for any value of the ratio $`L/\xi `$ is required.

We would like to thank the ISSP and the Slovak Academy of Sciences for the use of their computer facilities. PM would like to thank the Japan Society for the Promotion of Science and Sophia University for their hospitality and VEGA for financial support.


format: text

text:
# Spinor Bosonic Atoms in Optical Lattices: Symmetry Breaking and Fractionalization

\[

## Abstract

We study superfluid and Mott insulator phases of cold spin-1 Bose atoms with antiferromagnetic interactions in an optical lattice, including a usual polar condensate phase, a condensate of singlet pairs, a crystal spin nematic phase, and a spin singlet crystal phase. We suggest a possibility of exotic fractionalized phases of spinor BEC and discuss them in the language of topological defect condensation and $`Z_2`$ lattice gauge theory.

\]

Recent experiments on the condensation of Bose gases with internal degrees of freedom initiated a considerable amount of work on understanding the nature of spinor BEC. In particular Bose condensation in alkali atoms with nuclear spin $`I=3/2`$ that have three low energy hyperfine states and therefore behave as $`F=1`$ bosons, have been a subject of active experimental  and theoretical research. It was pointed out that the ground state of atoms with antiferromagnetic interactions in an optical trap is not a usual condensate of a macroscopic number of atoms in a single quantum state, but a complicated many body state of atoms arranged into a total singlet. In this article we examine the problem of Bose $`F=1`$ atoms with antiferromagnetic interaction (spin symmetric interactions are assumed throughout the paper) in optical lattices - arrays of microscopic potentials created by interfering laser beams . The dynamics of spinless bosonic atoms in arrays of optical wells may be described by a Bose-Hubbard model with the possibility of quantum phase transitions between insulating and superfluid phases induced by varying the properties of the laser light . Here, we propose several new phases that will appear for $`F=1`$ bosons due to an interplay between spin and charge degrees of freedom: singlet pair condensate that only breaks charge symmetry, a spin nematic crystal phase that only breaks spin symmetry, and a novel “strong coupling pairing phase” that breaks both spin and charge symmetry but is distinct from a simple polar BEC. These phases may be considered as “fragmented condensates” , and have fractionalized topological excitations: half-vortices and $`\pi `$-disclinations. Josephson type experiments may be used to distinguish among insulating and various superfluid phases. We also conjecture the possibility of fractionalized phases in spinor Bose gases. Fractionalization is characterized by a topological order and may coexist with any broken continuous symmetry .

In what follows we consider the case of antiferromagnetic interaction between the atoms $`a_0>a_2`$, where $`a_F`$ is an $`s`$-wave scattering length in the $`F`$ channel. In the individual wells we take wavefunctions that are superpositions of the mean-field solutions of the interacting problem in the absence of tunneling between the wells  $`|\psi _N`$=$`_𝐧\psi (𝐧)|N,𝐧`$, with $`|N,𝐧=\frac{1}{𝒩}(n_xa_x^{}+n_ya_y^{}+n_za_z^{})^N|0`$. Here $`𝐧`$ is a real unit vector, $`a_\alpha ^{}`$ are boson creation operators in the lowest vibrational state, $`|0`$ is the vacuum state, and $`𝒩`$ is a normalization factor. Two properties of these wavefunctions are important: (i) they satisfy a symmetry condition $`|N,𝐧=()^N|N,𝐧`$ ; (ii) in the limit of large $`N`$ states that correspond to different $`𝐧`$’s are orthogonal to each other. Hamiltonian of an array of identical optical wells with spin symmetric tunneling is given by

$``$ $`=`$ $`{\displaystyle \underset{i}{}}_i+{\displaystyle \underset{ij}{}}_{ij}`$ (1)
$`_i`$ $`=`$ $`{\displaystyle \frac{u}{2}}N_i^2\mu N_i+{\displaystyle \frac{g}{2}}𝐋_i^2`$ (2)
$`_{ij}`$ $`=`$ $`2t𝐧_i𝐧_j(b_i^{}b_j+b_j^{}b_i)`$ (3)
where we defined the charge creation and annihilation operators that change the number of particles $`N_i`$, but not the direction of $`𝐧_i`$: $`b_i^{}|N_i,𝐧_i=(N_i+1)^{1/2}|N_i+1,𝐧_i`$; the number of particles in each well $`N_i=b_i^{}b_i`$; and the angular momentum operators $`𝐋_i=i𝐧_i\times \frac{}{𝐧_i}`$ that describe a collective spin in well $`i`$ . Parity condition on the wavefunctions implies constraint $`N_i+L_i=`$even for all $`i`$, and for simplicity in (3) we chose to work in the grand canonical ensemble using a chemical potential $`\mu `$. For a single well Hamiltonian (3) correctly gives a spin singlet ground state if the number of particles is even . When optical lattices are produced by lasers with wavelength $`\lambda `$, they create an optical potential $`V(𝐱)=V_0_i\mathrm{sin}^2(kx_i)`$. We can estimate parameters of the Hamiltonian (3) as $`g=\frac{2\pi ^2}{3}E_R\frac{(a_0a_2)}{\lambda }(\frac{V_0}{E_R})^{3/4}`$, $`U=\frac{2\pi ^2}{3}E_R\frac{(a_0+2a_2)}{\lambda }(\frac{V_0}{E_R})^{3/4}`$, $`t=(4V_0E_R)^{1/2}Exp\{2\pi (\frac{V_0}{E_R})^{1/2}(12(\frac{E_R}{V_0})^{1/2})\}`$, where $`E_R=\mathrm{}^2k^2/2M`$ is the recoil energy. The average number of particles in each well will be related to the atomic gas density $`n`$ as $`\overline{N}=n(\frac{\lambda }{2\pi })^3(\frac{E_R}{V_0})^{3/4}`$. The obvious parameters to tune in experiments is $`\overline{N}`$ (e.g. by varying boson density) and the strength of periodic potential $`V_0`$ (by tuning the laser intensity). When the latter is varied it affects tunneling exponentially and the interaction coefficients only as a power law. In what follows, we therefore take $`U`$ and $`g`$ as constant and consider a phase diagram in the coordinates $`t`$ and $`\mu `$.

There are two limits when construction of the phase diagram of (3) is simple. When $`g=0`$ there is no energetic penalty for having high angular momentum (non zero spin) states of quantum rotors in (3) and spin symmetry is broken. For integer filling factors we have spin nematic insulator phases (NI) for small $`t`$ and a superfluid polar condensate phase (PC) for larger values of $`t`$. For fractional filling factors the system is always superfluid. The other simple case is when $`g`$ is large and states with odd number of atoms in individual wells are energetically costly: they are not allowed to have $`L=0`$ state of a rotor due to parity constraint and have a higher rotational energy. It is energetically favorable to have even numbers of atoms in each well arranged into singlet combinations, which can be interpreted as binding of atoms into singlet pairs with a pair binding energy $`E_p(Q)=2E(Q+1)E(Q+2)E(Q)=g`$. In the limit when $`g\mathrm{}`$ the crystal phases are possible for even numbers of atoms per well only and correspond to spin singlet insulators (SSI). The superfluid phase is a condensate of singlet pairs (SSC), in which tunneling of individual atoms between the wells is suppressed and only singlet pairs are delocalized. The origin of pairing in this case is not the attraction between individual bosons, but a singlet formation on the scale of individual wells. This is reminiscent of the “attraction from repulsion” mechanism of electron pairing proposed by Chakravarty and Kivelson for high Tc cuprates, $`C_{60}`$, and polyacenes . In the case of finite $`g`$ we expect Mott crystal phases for all integer filling factors. For even $`N`$’s and small $`t`$ the rotors kinetic energy dominates and we have singlet crystal phases (SSI). When $`t`$ is increased the system goes to a spin nematic insulator (NI) phase that has admixture of $`L0`$ in the ground state. For odd $`N`$’s when all charge fluctuations are frozen there is always at least $`L=1`$ on each site, so we expect the system to be analogous to spin 1 lattice model with a broken spin symmetry. The superfluid phase will be spin singlet fluid of pairs (SSC) for small $`t`$, and a polar condensate (PC) for larger $`t`$’s.

Simple Josephson type experiment may be used to distinguish superfluid phases suggested above. We imagine spinor atoms in an optical lattice in the presense of gravity. Gravity acts as voltage in a conventional Josephson junction arrays and gives rise to Josephson oscillations . Josephson oscillations may be measured selectively in spin (1), or measuring all the particles regardless of their spin (2). In the case of single atom condensate (PC) we expect oscillations for both kinds of measurements with frequent $`\mathrm{}\omega _J=mg\lambda /2`$. A singlet pair condensate will only give oscillations in channel (2) with frequency $`2\mathrm{}\omega _J=mg\lambda `$. Insulating phases will not show any Josephson oscillations, and nematic vs spin singlet may be distinguished further through different correlations for particles of different spin.

Order parameters for various phases suggested above and shown on Figure 1 may be conveniently discussed using a Euclidean space-time action. We implement a constraint $`N_i+L_i=`$even using a projection operator $`P_i=\frac{1}{2}_{\sigma _i=\pm 1}e^{i\frac{\pi }{2}(1\sigma _i)(N_i+L_i)}`$ . To calculate the partition function $`Z=Tr[e^\beta P]`$ we divide imaginary time $`\beta `$ into $`M`$ slices of length $`ϵ=\beta /M`$ and after a few standard manipulations that include Hubbard-Stratanovich transformation of the last term in (3) we find

$`S={\displaystyle \underset{rr^{}}{}}J_{rr^{}}^c\sigma _{rr^{}}cos\varphi _{rr^{}}{\displaystyle \underset{rr^{}}{}}J_{rr^{}}^{2c}cos(2\varphi _{rr^{}})`$ (4)
$`{\displaystyle \underset{rr^{}}{}}J_{rr^{}}^s\sigma _{rr^{}}𝐧_r𝐧_r^{}{\displaystyle \underset{rr^{}}{}}J_{rr^{}}^{2s}Q_r^{ab}Q_r^{}^{ab}`$ (5)
Here summation goes over sites $`r=(i,\tau )`$ in the space-time lattice, $`\varphi _r`$ is the phase variable conjugate to $`N_r`$, $`\varphi _{rr^{}}=\varphi _r\varphi _r^{}`$, $`Q_r^{ab}=n_r^an_r^b\delta ^{ab}/3`$ is a nematic order parameter, and $`\sigma _{rr^{}}=\pm 1`$ is an Ising field that lives on the links rather than sites. The coupling constants are $`J_{r,r\pm \widehat{\tau }}^c=(ϵu)^1`$, $`J_{r,r\pm \widehat{\tau }}^s=(ϵg)^1`$, $`J_{r,r\pm \{\widehat{x},\widehat{y},\widehat{z}\}}^{c,s}=ϵJ|\chi |`$, $`J_{r,r\pm \widehat{\tau }}^{2c}=J_{r,r\pm \widehat{\tau }}^{2s}=0`$, $`J_{r,r\pm \{\widehat{x},\widehat{y},\widehat{z}\}}^{2c,2s}=ϵJ/4`$, with $`\chi `$ being the saddle point value of the Hubbard-Stratanovich field. The usual periodic boundary conditions for $`\varphi `$, $`𝐧`$, and $`\sigma `$ are assumed at $`\tau =0`$ and $`\tau =\beta `$. In writing (5) we omitted the Berry phase terms that are important for quantitative calculations of the phase diagram but not for the symmetry arguments discussed below.

Action (5) has a $`Z_2`$ gauge symmetry $`\varphi _r\varphi _r+\frac{\pi }{2}(1ϵ_r)`$, $`𝐧_rϵ_r𝐧_r`$, $`\sigma _{r,r+\widehat{a}}ϵ_r\sigma _{r,r+\widehat{a}}ϵ_{r+\widehat{a}}`$, where $`ϵ_r=\pm 1`$. Such a symmetry has been pointed out earlier in , where it was observed that the physical order parameter is a complex vector $`𝐝=e^{i\varphi }𝐧`$ (not $`\varphi `$ and $`𝐧`$ separately) that has a symmetry of $`e^{i\varphi }e^{i\varphi }`$, $`𝐧𝐧`$. Equation (5) is the simplest action consistent with the charge $`U(1)`$, spin $`SO(3)`$, and gauge $`Z_2`$ symmetries of the model. Another term allowed by the $`Z_2`$ symmetry is the analogue of the Maxwell terms for the lattice gauge models $`S_\sigma =K_{\mathrm{}}_{\mathrm{}}\sigma _{ij}`$, where summation goes over plaquettes in $`d+1`$ dimensional lattice. This term may be generated by integrating out the high energy degrees of freedom.

The existence of the local symmetry imposes important constraints on the possible order parameters and symmetry breaking states of the system. Only order parameters that are gauge invariant may acquire expectation values. For example, the expectation values $`𝐧0`$ or $`e^{i\varphi }0`$ are not allowed in the models described by (5), since $`𝐧_i`$ and $`e^{i\varphi _i}`$ are not invariant under $`Z_2`$ gauge transformations. Physically this restriction arises from the fact that the wavefunction $`\mathrm{\Psi }(𝐧)`$ should be an even or odd function of $`𝐧`$ for $`N`$ even or odd respectively. Below we consider several examples of the order parameters that are gauge invariant and review broken symmetry phases that they lead to (see Figure 1). Polar BEC phase (PC) has $`𝐝=e^{i\varphi }𝐧0`$. It breaks both charge and spin symmetries and is the phase where atoms are condensed directly. Nematic Insulator phase (NI) has $`Q_{ab}=n_an_b\delta _{ab}/30`$ and breaks only the spin SO(3) symmetry. This is a phase where not the atoms but certain particle-hole composites are condensed. Spin singlet condensate (SSC) $`e^{i2\varphi }0`$ breaks the $`U(1)`$ symmetry but not $`SO(3)`$. It corresponds to a condensate of singlet pairs of atoms. Strong coupling pairing phase (SCP) is a phase where both $`e^{i2\varphi }0`$ and $`Q_{ab}0`$. This phase breaks both $`U(1)`$ and $`SO(3)`$ symmetries, however, it is fundamentally different from the polar condensate phase. Transition between phases PC and SCP should be Ising type transition. Similar phase has been recently suggested in the context of bilayer quantum Hall systems and called “strong coupling pairing phase” . Finally, spin singlet insulator (SSI) is a phase where no field has expectation value, so neither U(1) nor SO(3) symmetries are broken. We note that for a finite number of traps, when spontaneous symmetry breaking is not possible, one can observe signatures of the transitions by measuring fluctuations in the number of atoms, atom pairs, and spins .

It is interesting to point out that phases NI, SSC, and SCP correspond to the fragmented BEC discussed by P. Nozieres and D. Saint James. The part emphasized in  is that fragmented state has no macroscopic population of k=0 state for the bosons. Another interesting aspect of these states is that they achieve fractionalization of topological objects. The unfragmented PC phase has topological excitations that are composites of 1/2 vortex in the charge sector and $`\pi `$-disclination in the spin sector. Here 1/2 charge vortex is a topological defect around which the phase $`\varphi `$ winds by $`\pi `$ and $`e^{i\varphi }`$ changes sign. $`\pi `$-disclination in the spin sector is introduced as a topological defect around which $`𝐧`$ changes sign (vortex like objects in $`𝐧`$ are called merons, so $`\pi `$-disclination is 1/2 of a meron). Both 1/2 vortex and $`\pi `$-disclination are pointlike objects in two dimensions and lines in three dimensions. No individual 1/2 charge vortices or $`\pi `$-disclinations are allowed in the PC phase. The spin and charge parts of the vortex are always glued together. A condensate of singlet pairs only, the SSC phase, has 1/2 charge vortices as separate excitations, but no spin vortices; and the nematic phase NI has $`\pi `$ disclinations and no charge vortices. A condensate of S=0 pairs and S=2 pairs, a SCP phase, will have 1/2 charge vortices and $`\pi `$ disclinations as separate excitations.

Establishing the nature of the broken symmetry does not fully characterize models with gauge symmetry. Another important aspect of the system described by (5) is the possibility of confining and deconfining phases of the $`Z_2`$ gauge theory in dimension $`d+13`$ . Difference between the confining and deconfining phases relies on the concept of topological order  and has important implications on the nature of excitations in the system. Topological order may coexist with any true long range order, so any of the phases reviewed earlier may be confining or deconfining. For the pure $`Z_2`$ gauge models the confining and the deconfining phases may be distinguished using Wilson loops, which are characterized by an area law for the confining phase and by a perimeter law for the deconfining phase. With matter field present Wilson loops may no longer be used to discriminate between the two phases , however distinction between them survives, and there is a phase transition between the two. A simple way to understand the transition is to think of it as condensation of visons, topological excitations of the $`Z_2`$ gauge theory that describe frustrated plaquettes $`_{\mathrm{}}\sigma _{rr^{}}=1`$. Particles that carry a $`Z_2`$ charge are frustrated when traveling around a vison, so when visons are condensed such particles may not propagate coherently and we have a confining phase. When visons are not condensed, they are finite energy topological excitations and we have a deconfined phase. Presense of finite energy visons leads to an additional degeneracy of the gauge models in a deconfined state on manifolds with non-trivial degeneracy. Systems with $`Z_2`$ gauge symmetries in condensed matter have already been discussed by Lammert et.al.  in connection with liquid crystals and Senthil and Fisher in connection with high temperature superconductors . For our system the discussion above implies that all the broken symmetry states , as well as a state SSI, where neither U(1) nor SO(3) symmetries are broken, come in two varieties: confining and deconfining. Detailed comparison between confining and deconfining versions of various will be given in . Here we only note that some of the most striking implications of deconfinement appear for the SSI phase. The deconfining SSI phase allows charge $`b=e^{i\varphi }`$ and spin $`𝐧`$ carrying excitations propagate independently .

Another perspective on the nature of deconfining phases comes from considering phase transitions as condensation of topological defects. Discussion will be given in the case of two spatial dimensions, where condensation of point-like topological objects (vortices, $`\pi `$-disclinations, merons, and etc.) drive transitions into the phases of restored symmetry. We expect, however, that most of the conclusions reached for $`d=2`$ will also hold in $`d=3`$. Let us consider a transition between the PC phase that breaks charge and spin symmetry and NI phase that breaks only the spin symmetry. In the phase PC we have topological objects that are 1/2 vortex $`\pi `$-disclination composites. Half-vortex and half-meron parts may come with different signs, so we have 4 types of these composites ($`\pm `$1/2, $`\pm `$1/2), where the first and the second numbers stand for the charge and the spin parts respectively. To go to a phase that has broken spin symmetry and restored charge symmetry one expects that we need to condense topological objects in the charge sector, with topological objects in the spin sector remaining at finite energy. This can be achieved by condensing composites $`(1/2,1/2)+(1/2,1/2)=(1,0)`$ and $`(1/2,1/2)+(1/2,1/2)=(1,0)`$, but having individual topological defects $`(\pm 1/2,\pm 1/2)`$ gapped. Such transition would correspond to condensation of integer vortices and would be equivalent to order-disorder transition for the charge order parameter $`e^{i\varphi }`$ in the usual Bose-Hubbard model. An important observation is that the insulating phase obtained this way is a fractionalized phase NI! It does not have a condensate of minimal topological object in the charge sector (1/2 vortex), but only a condensate of integer vortices. So, only the $`PCNI^{}`$ transition may be a continuous second order transition, but $`PCNI`$ must be a first order transition. The existence of fractionalized phases may, therefore, be established through study of phase transitions . We note, however, that starting from a strong coupling pairing phase SCP (that has unbound 1/2- vortices and $`\pi `$-disclinations) we will have $`SCPNI`$ as a second order transition.

To summarize, we reviewed several novel phenomena that will take place for $`F=1`$ Bose gases with antiferromagnetic interactions. We showed that in optical lattices they will have several distinct superfluid and Mott insulator phases depending on the density and parameters of the system. Another possibility suggested in this paper is a class of fractionalized phases of spinor BEC, characterized by the topological order that may coexist with any true long range order in the system. Quantitative analysis of experimental systems will be given in subsequent publications . Discussion in this paper may be generalized to systems with higher hyperfine spins. We also expect some of the results to be applicable to triplet superconductors.

We thank M. Fisher, H.Y. Kee, Y.B. Kim, M. Lukin, C. Nayak, S. Sachdev, J. Sethna, T. Senthil, H. T. Stoof, P. Wiegnman and particularly F.D. Haldane for useful discussions. ED is supported by the Harvard Society of Fellows.


text:
# Nonequilibrium phase transitions in a model for the origin of life

## 1 Introduction

The most fundamental event in the history of life was probably the appearance via spontaneous creation of a molecule capable of replicating itself (replicator). Given a possible mechanism of replication, which in this case is some form of template activity, the evolution of such replicators has been extensively investigated through the chemical kinetics formalism put forward by Eigen and co-workers in the seventies . Those studies have raised a series of objections to the simplistic view of the emergence of a complex organism from a collection of competing species of replicators. For instance, the finding that the length of a molecule (polynucleotide) is limited due to the finite replication accuracy per nucleotide has prompted the proposal of models that incorporated some sort of cooperation between the replicators, such as the molecular catalytic feedback networks termed hypercycles . These models, however, have attracted their own criticisms since, as pointed out by Maynard Smith, giving catalytic support in such molecular networks is in fact an altruistic behavior and so they are extremely vulnerable to the presence of parasites, i.e., molecules that do not reciprocate the catalytic support they receive . A possible solution to this problem is provided by the structured deme formulation of group selection , where it is assumed that the replicators are spatially localized, say, in rock crevices or water droplets, so that the benefits accrued from cooperation are directed mostly to the members of the catalytic network . Yet another successful approach to the problem of resistance against parasites is based on a reaction-diffusion system where replication and diffusion taking place on an adsorbing surface generate self-organized spiral structures . Interestingly, as these spatial structures, which greatly increase the stability of the hypercycles against parasites, can be viewed as super-organisms that approach is also related to the group selection theory .

Since in the prebiotic or chemical evolution context, natural selection is essentially the dynamics of replicators, it is not surprising that most of the studies in this subject have focused almost exclusively on the competition between replicators, among which the so-called malthusian and hypercyclic replicators are the most important . The former corresponds to the simplest reproduction process, namely, the binary fission of a parent replicator and is modelled by the chemical reaction

$$A+E\stackrel{s}{}2A$$
(1)
where $`A`$ is the replicator and $`E`$ is the source materials (mononucleotide resources). It is well-known that the concentration of $`A`$ grows exponentially with the rate constant $`s`$, provided that the concentration of $`E`$ is kept constant, hence the name malthusian replicator. To avoid this explosive growth, one usually imposes a constraint on the total concentration of replicators which can be implemented in practice by a dilution flux . Alternatively, one can allow the replicators to be degraded by hydrolysis into its mononucleotide components $`E`$ according to the reaction

$$A\stackrel{\gamma }{}E$$
(2)
which seems a more natural approach to limit the growth of $`A`$.

As best exemplified by sexual reproduction, there are situations that cannot be described by (1) since two replicators are necessary to produce a third one. In this case the corresponding chemical reaction is

$$2A+E\stackrel{c}{}3A$$
(3)
which leads to a hyperbolic growth of the concentration of replicator $`A`$ . A hypercyclic replicator is defined as one that can replicate itself using both (1) and (3) reaction schemes. Actually, the term hypercycle derives from the superimposition of the catalytic replication cycle (3) on the self-replication cycle (1). Of course, the limit $`c=0`$ corresponds to the malthusian replicator while $`s=0`$ can be associated to an obligatory sexual replicator.

In contrast to previous works that have concentrated on the competition between replicators either of the same kind but with different production rates  or of different kinds , in this paper we address a more fundamental problem that has received comparatively little attention, namely, the stability of the different kinds of replicators, viewed here as an active (ordered) phase of the molecular system, against the empty (disordered) phase composed of the resource materials only. This lack of interest was probably due to the fact that the usual kinetics formalism used to study the dynamics of replicators does not represent the mononucleotide resource dynamics explicitly (see however ), thus precluding the study of the issues addressed in the present contribution. More pointedly, we consider the dynamics of a population of identical hypercyclic replicators on a lattice space both in the deterministic infinite diffusion (mean-field) limit and in the stochastic position-fixed (contact process) limit where each replicator on a lattice cell never moves. The last limit is particularly interesting because it allows the connection between the replicator models and some standard models of nonequilibrium phase transition in a lattice (e.g. directed percolation) . As a result, the powerful analytical tools of statistical mechanics can be readily used to advance our understanding of the evolution of replicators. Of particular relevance is the so-called dynamic Monte Carlo method whose idea is to set the system initially in the empty state with a seed of replicators in the center of the lattice and then study the subsequent spreading of activity . More importantly, the thorough analysis of both limits exposes the limitations of the widely used deterministic chemical kinetics or mean-field formalism to study the problem of the emergence of life.

The remainder of the paper is organized as follows. In Sec. 2 we present the set of rules that govern the evolution of a population of hypercyclic replicators in a two-dimensional square lattice. The mean-field or infinite diffusion limit, which models an ideally mixed medium, is studied analytically in Sec. 3. The results are summarized in a phase diagram showing the regions of stability of the empty and active regimes in the space of the control parameters of the model. Those regions are delimitated by continuous as well as discontinuous transition lines that join at a tricritical point. In Sec. 4 we study the position-fixed or contact process limit using mainly the dynamic Monte Carlo method which allows the computation of the critical dynamic exponents that describe quantitatively the spreading of a vanishingly small population of replicators. Finally, some concluding remarks are presented in Sec. 5. In particular, we compare the hypercyclic replicator model with Schlögl’s models of nonequilibrium phase transition in reaction-diffusion systems .

## 2 The model

It has been suggested that chemical evolution started with a surface-bonded autocatalytic chemical network as there are enormous thermodynamics and kinetics advantages of surface binding reactions, especially in the case of reactions that require unlikely collisions of many reactants . The binding must be strong enough to keep the reactants on the surface but also flexible enough to allow their slow migration on it. Interestingly, if this proposal proves correct it will probably lead to the replacement of the popular notion of a primitive soup by that of a primitive pizza instead . Accordingly, we define our replicator model in a two-dimensional space consisting of $`L\times L`$ cells in a square toroidal lattice. Each cell is either empty or occupied by a replicator and it is assumed that an empty cell contains all source material required to assemble a new replicator. The evolution of the population of replicators is governed by the following local rules:

* A replicator has a probability $`\gamma `$ of decaying; after decay the cell becomes empty. This rule is motivated by the hydrolysis reaction (2).
* A replicator in one of the four first neighbor cells (von Neumann neighborhood) of an empty cell can replicate into that cell with probability $`s`$. This process is referred to as non-catalised self-replication and is motivated by the reaction (1).
* Regardless of the previous rule, a replicator in the von Neumann neighborhood of an empty cell can replicate into that cell if there are other replicators in the intersection of the Moore neighborhoods of both cells. The probability of this type of replication, which is motivated by reaction (3), is $`c`$ for each pair of replicators. We recall that the Moore neighborhod of a given cell consists of its first and second nearest neighbors, adding up to eight cells.

Hence in the extreme situation where an empty cell is surrounded by eight replicators, it can become occupied with probability $`4s+16c`$. To carry out the simulations we choose the parameters $`c`$ and $`s`$ such that $`4s+16c1`$. These rules are applied simultaneously to all cells in the lattice so our model can be viewed as a two-dimensional stochastic cellular automaton. Actually the model is essentially an adaptation to one-membered hypercycles of the spatial cellular automaton model of multi-membered hypercycles proposed by Boerlijst and Hogeweg . The dynamics defined by the rules given above is manifestly irreversible and, in particular, the state characterized by empty cells only is an absorbing state, i.e., a configuration from which the system cannot escape. In this sense, the principle of detailed balance is broken and the active stationary state is in fact in nonequilibrium. Although the more realistic situation is a diffusion-controled reaction where each reactant can move randomly on the lattice, in this paper we choose to study in detail the simpler extreme cases of infinite diffusion (mean field) and no diffusion (contact process). Of course, we hope that features common to both limits will be present in the finite diffusion situation as well.

## 3 The mean-field limit

The mean-field limit describes exactly an infinite population of reactants in an ideally mixed medium and so it is equivalent to the usual chemical kinetics formulation. Neglecting spatial correlations among cells, i.e., assuming that at any time the molecules are distributed randomly over the lattice cells, it is straightforward to write the evolution equation for the density of replicators or occupied cells at time $`t`$, namely,

$$\rho _{t+1}=\rho _t\left(1\gamma \right)+4\rho _t\left(1\rho _t\right)\left(s+4c\rho _t\right).$$
(4)
We will consider only the stationary solutions $`\rho _{t+1}=\rho _t=\rho `$ of this equation. The absorbing (empty) state $`\rho =0`$ is always a solution, while the non-zero solutions are given by the roots of the quadratic equation

$$\stackrel{~}{c}\rho ^2\rho \left(\stackrel{~}{c}\stackrel{~}{s}\right)+1\stackrel{~}{s}=0$$
(5)
where we have introduced the dimensionless parameters

$$\stackrel{~}{c}=\frac{16c}{\gamma }\text{and}\stackrel{~}{s}=\frac{4s}{\gamma }.$$
(6)
This equation has real roots provided that the condition $`\left(\stackrel{~}{c}+\stackrel{~}{s}\right)^24\stackrel{~}{c}`$ is satisfied. In addition we can easily show that: (i) for $`\stackrel{~}{s}<1`$ and $`\stackrel{~}{c}>\stackrel{~}{s}`$ both roots are positive; (ii) for $`\stackrel{~}{s}<1`$ and $`\stackrel{~}{c}<\stackrel{~}{s}`$ both roots are negative; and (iii) for $`\stackrel{~}{s}>1`$ only one of the roots is negative. Furthermore, though the solution $`\rho =0`$ exists in the entire plane $`(\stackrel{~}{c},\stackrel{~}{s})`$, it is stable only if the condition

$$\frac{\rho _{t+1}}{\rho _t}|_{\rho =0}=1\gamma +4s<1,$$
(7)
which reduces to $`\stackrel{~}{s}<1`$, is satisfied. In the region where the non-zero roots are physical (i.e., real and positive), the stable root is always the largest one. In Fig. 1 we show the steady-state density of replicators for two different choices of initial density. We identify three distinct phases: the absorbing or empty phase $`\left(E\right)`$ associated to the solution $`\rho =0`$; the replicating or active phase $`\left(A\right)`$ associated to the solution $`\rho >0`$; and the phase labeled $`\left(EA\right)`$ where both solutions $`\rho =0`$ and $`\rho >0`$ are stable. In this phase, the outcome of the dynamics is not determined by the control parameters only but also by the initial abundance of replicators. ¿From this figure it is clear that the system undergoes a continuous non-equilibrium phase transition from phase $`\left(E\right)`$ to phase $`\left(A\right)`$ at $`\stackrel{~}{s}=\stackrel{~}{s}_c=1`$ and $`\stackrel{~}{c}<1`$. Explicitly, near this transition the density of replicators can be written as

$$\rho \frac{\stackrel{~}{s}1}{1\stackrel{~}{c}}.$$
(8)
The continuous transition ends at the tricritical point $`\stackrel{~}{s}_t=\stackrel{~}{c}_t=1`$ so that for $`\stackrel{~}{c}>\stackrel{~}{c}_t`$ the transition between phases $`\left(E\right)`$ and $`\left(EA\right)`$, that takes place at $`\stackrel{~}{s}=2\sqrt{\stackrel{~}{c}}\stackrel{~}{c}`$, as well as the transition between phases $`\left(EA\right)`$ and $`\left(A\right)`$, that occurs at $`\stackrel{~}{s}=1`$, are discontinuous. In particular, the jumps of the densities of replicators are

$$\mathrm{\Delta }\rho =11/\stackrel{~}{c}^{1/2}\stackrel{~}{c}>1$$
(9)
at the transition $`\left(E\right)\left(EA\right)`$, and

$$\mathrm{\Delta }\rho =11/\stackrel{~}{c}\stackrel{~}{c}>1$$
(10)
at the transition $`\left(EA\right)\left(A\right)`$. These results are conveniently summarized in the phase diagram shown in Fig. 2. It is interesting to note that setting $`\stackrel{~}{s}=\stackrel{~}{c}`$ we find, close to the tricritical point,

$$\rho \left(\stackrel{~}{s}1\right)^{1/2}$$
(11)
so that the tricritical point is not in the same universality class as the transition observed in the absence of catalytically assisted replication.

The interpretation of our results within the prebiotic evolution context leads to the conclusion that for finite values of $`\stackrel{~}{c}`$ an obligatory sexual replicator cannot emerge spontaneously (i.e., appear at vanishingly small concentrations). For instance, for $`\stackrel{~}{s}=0`$, the minimal initial density of replicators necessary to engender a prosperous population is $`\rho =1/2`$ at the transition point $`\stackrel{~}{c}=4`$, vanishing as $`1/\stackrel{~}{c}`$ for large $`\stackrel{~}{c}`$. Actually, an initial colony of replicators is certain to grow from vanishingly small concentrations provided that $`\stackrel{~}{s}>1`$. These conclusions, however, must be taken with caution since in the deterministic limit a vanishingly small concentration means an infinite population of replicators, while one would expect the first replicators to show up as a single or a few copies at most. Of course, a proper understanding of this emergence phenomenon calls for a stochastic approach, which is the subject of the next section.

## 4 The position-fixed limit

The primary aim of this section is to determine what features of the rich phase-diagram obtained in the mean-field limit show up also in the opposite limit, where the replicators are fixed on the lattice cells. This is a rather challenging enterprise as, at least in the case of equilibrium phase-transitions, there is no totally unambiguous way by which one can detect the order of the transition through the analysis of finite systems alone . Nevertheless, we tackle this problem using both a steady-state approach for finite lattices and the dynamical Monte Carlo method for lattices of effectively infinite size.

First we measure de density $`\rho `$ of replicators in the steady state. Our results for two system sizes ($`L=100`$ and $`L=200`$) are shown in Figs. 3 and 4. For each set of the control parameters we made runs of $`10^5`$ generations, neglecting the first $`210^4`$ generations and recording $`\rho `$ at steps of $`200`$ generations. A generation corresponds to the simultaneous update of all lattice cells. Each data point is the arithmetic mean of these recorded data. Provided that the population is not extinct, the results are independent of the choice of the initial configuration. In particular, we have made runs starting from all cells occupied or from a seed of only four clustered occupied cells. Furthermore, as in the mean-field limit we have verified that our results depend only on the ratios $`s/\gamma `$ and $`c/\gamma `$ so throughout the remaining of this paper the value of the decay constant is held fixed at $`\gamma =0.05`$. These results clearly indicate the existence of a phase transition separating the active ($`\rho >0`$) and empty ($`\rho =0`$) phases of the model. Assuming that in the neighborhood of the transition points the density of replicators goes to zero as $`\rho \left(\stackrel{~}{s}\stackrel{~}{s}_c\right)^\beta `$ and using the least-square method we can estimate both the critical replication rate $`\stackrel{~}{s}_c`$ and the critical exponent $`\beta >0`$. These estimates are presented in table 1 and the quality of the fitting can be appreciated from Fig. 4. The statistical errors are of order of $`10^3`$ but the systematic errors, which are due mainly to the difficulty to carry out long runs very close to the transition point, are probably much larger. We note that our estimates of $`\beta `$ for small values of $`\stackrel{~}{c}`$ indicate that in this regime the replicator model belongs to the so-called $`\left(2+1\right)`$ directed-percolation universality class for which $`\beta =0.59\pm 0.02`$ . It should be pointed out that in finite systems the active regime is a meta-stable one as there is always a finite probability that the colony becomes extinct due to fluctuations in the stochastic dynamics. Since this extinction probability increases towards $`1`$ as the critical point is approached, it is very difficult to obtain reliable estimates of $`\stackrel{~}{s}_c`$ and $`\beta `$ by means of numerical simulations in the steady-state regime. Of course we are aware that if the transition happens to become discontinuous at some value of $`\stackrel{~}{c}>0`$ then the assumption of a power-law singularity at the transition point breaks down (actually $`\beta =0`$ in this case). In fact, the (anomalous) continuous decrease of $`\beta `$ as $`\stackrel{~}{c}`$ increases (see table 1) is an indication that this might be the case. Our careless use of the power-law assumption is intentional and aims to illustrate the difficulty of detecting the order of a phase transition using results of steady-state simulations: the abrupt variation of $`\rho `$ at criticality observed in Fig. 3 for certain values of $`\stackrel{~}{c}`$, which might indicate the occurrence of a first-order irreversible transition, can be explained as a continuous transition with a small exponent $`\beta `$ as well.

We turn now to the analysis of the spreading behavior of a small colony of replicators settled initially in the center of an otherwise empty lattice of infinite size. More pointedly, the initial colony is composed of four replicators located in the von Neumann neighborhod of the central empty cell. Finite size effects are absent because the lattice size is taken large enough so that during the time we follow the evolution of the colony the replicators can never reach the lattice boundaries. This of course sets an upper limit to the number of generations we can follow the colony and so, in particular, we let the population evolve up to typically $`t=10^4`$. As usual, we concentrate on the time dependence of the following key quantities : (i) the average number of replicators $`n\left(t\right)`$; (ii) the survival probability of the colony $`P\left(t\right)`$; and (iii) the average mean-square distance over which the replicators have spread $`R^2\left(t\right)`$. For each time $`t`$ we carry out $`M=210^5`$ independent runs, all starting with the same initial colony. Hence $`P\left(t\right)`$ is simply the fraction of runs for which there is at least one replicator in the lattice at time $`t`$. Since $`n\left(t\right)`$ is an average taken over all runs including those which have already been extinct at generation $`t`$, the average number of replicators per surviving run is given by the ratio $`N\left(t\right)=n\left(t\right)/P\left(t\right)`$. Furthermore, noting that $`R^2\left(t\right)`$ is averaged only over the surviving runs, we can define the fractal dimension $`d_f`$ of the surviving colonies of replicators at a given time $`t`$ as $`NR^{d_f}`$.

At the transition points we expect that the measured quantities obey the following scaling laws

$`n\left(t\right)`$ $``$ $`t^\eta `$ (12)
$`P\left(t\right)`$ $``$ $`t^\delta `$ (13)
$`R^2\left(t\right)`$ $``$ $`t^z`$ (14)
where $`\delta `$, $`\eta `$ and $`z`$ are dynamic exponents which are related with the fractal dimension of the clusters of replicators through the equation

$$d_f=2\frac{\eta +\delta }{z}.$$
(15)
In principle, these scaling laws are valid for continuous as well as discontinuous phase transitions, though the scaling relations between the exponents, such as the ‘hyperscaling’ relation

$$\frac{1}{2}dz\eta =2\delta $$
(16)
where $`d`$ is the lattice dimension, hold only in the case of continuous transitions.

In Figs. 5 and 6 we present log-log plots of $`n\left(t\right)`$ and $`P\left(t\right)`$, respectively, as functions of $`t`$ in the vicinity of the critical point for $`\stackrel{~}{c}=0`$. The dependence of $`R^2\left(t\right)`$ on $`t`$ is not shown since, near criticality, the curves for different values of $`\stackrel{~}{s}`$ are clustered together and do not reveal any qualitatively relevant information on the colony evolution. The asymptotic straight lines observed in these figures are the signature of critical behavior while upward and downward deviations indicate supercritical and subcritical behaviors, respectively. A precise estimate for the critical exponents is obtained by considering the local slopes of the curves shown in the previous figures. For instance, the local slope $`\delta \left(t\right)`$ is defined by

$$\delta \left(t\right)=\frac{\mathrm{ln}\left[P\left(t\right)/P\left(t/5\right)\right]}{\mathrm{ln}5},$$
(17)
which for large $`t`$ behaves as

$$\delta \left(t\right)\delta +\frac{a}{t}$$
(18)
where $`a`$ is a constant. Analogous expressions hold for $`\eta \left(t\right)`$ and $`z\left(t\right)`$. Hence plots of the local slopes as functions of $`1/t`$ allow the calculation of the critical exponents. Using this procedure we find $`\stackrel{~}{s}_c=1.628\pm 0.001`$ which yields the exponents $`\eta =0.23\pm 0.01`$, $`\delta =0.45\pm 0.01`$ and $`z=1.13\pm 0.01`$. The error in $`\stackrel{~}{s}_c`$ is estimated by determining two values of $`\stackrel{~}{s}`$ as close as possible to the critical point for which upward and downward deviations can be observed, while the errors in the critical exponents are, as usual, the statistical errors obtained by fitting the local slopes by straigh lines in the large $`t`$ regime. Our exponents are in good agreement with those of the $`\left(2+1\right)`$ directed-percolation  and satisfy very well the hyperscaling relation (16) thus indicating that the transition in the limit $`\stackrel{~}{c}=0`$ is continuous, as expected.

The results of the spreading analysis for the other extreme case, $`\stackrel{~}{s}=0`$, which models a population of obligatory sexual replicators are shown in Figs. 7 and 8. Although the dependence of $`\mathrm{ln}P\left(t\right)`$ on $`t`$ is similar to that observed in the previous case, the behavior pattern of $`\mathrm{ln}n\left(t\right)`$ (see Fig. 7) is rather different: in the supercritical regime (i.e., $`\stackrel{~}{c}>\stackrel{~}{c}_c`$) $`n\left(t\right)`$ first increases reaching a maximum, then decreases reaching a minimum and finally starts to increase monotonically again. This pattern is illustrated better in Fig. 9, which shows the spreading results for $`\stackrel{~}{s}=0.2`$. Analysis of these figures, which exceptionally show the colony evolution up to $`210^4`$ generations, suggests that a flat line separates the supercritical and the subcritical regimes implying thus the vanishing of the exponent $`\eta `$ at the transition point. Furthermore, the qualitatively distinct behavior patterns of $`\mathrm{ln}n\left(t\right)`$ illustrated in Figs. 5, 7 and 9 can be used to identify unambiguously the order of the nonequilibrium phase transition and hence to estimate the location of the tricritical point. To appreciate how the time dependence of $`\mathrm{ln}n\left(t\right)`$ in the supercritical regime changes continuously from the simple monotonic increase for $`\stackrel{~}{s}=0`$ to the complicated behavior described above for $`\stackrel{~}{c}=0`$ we present in Figs. 10 and 11 log-log plots of $`n\left(t\right)`$ as a function of $`t`$ for $`\stackrel{~}{s}=0.6`$ and $`\stackrel{~}{s}=1.0`$, respectively. In fact, analysis of Figs. 7 to 11 suggests that the turning point between those distinct behavior patterns occurs when the maximum and the minimum of $`n\left(t\right)`$ coincide, i.e., the critical curve $`\mathrm{ln}n\left(t\right)`$ vs. $`\mathrm{ln}t`$ has an inflection point (see Fig. 10). The value of $`\stackrel{~}{s}`$ and the corresponding $`\stackrel{~}{c}_c`$ at which this behavior occurs are then identified as the coordinates of the tricritical point. Applying this procedure we find $`\stackrel{~}{s}_t=0.60\pm 0.04`$ and $`\stackrel{~}{c}_t=4.00\pm 0.20`$.

The transition points and the dynamic exponents obtained via the scaling laws (12)-(14) and via the analysis of the local slopes are summarized in table 2. The errors in the estimates of the transition points are calculated as described before for the case $`\stackrel{~}{c}=0`$. Except for that case, we refrain from presenting the (statistical) errors in the exponents since the systematic errors are unusually large, due probably to the crossover behavior among (at least) three different universality classes. For instance, in the vicinity of the transition point for $`\stackrel{~}{s}=0.6`$, analysis of the local slope $`\delta \left(t\right)`$ up to $`t=500`$ indicates a clear tendency to the asymptotic value $`\delta 1`$ while for $`t>500`$ the tendency suddenly changes towards the asymptotic value $`\delta 0.6`$. A similar phenomenon occurs for $`\stackrel{~}{s}=1.0`$ also: the initial tendency is towards $`\delta 0.6`$ and then changes towards $`\delta 0.45`$ for larger times. As a result the estimate of the exponents becomes strongly dependent on the precise location of the transition points, which requires even better statistics as well as much longer runs. We leave this interesting research line which includes, for instance, the identification of the universality class of the tricritical point to a future, more technical contribution.

The evidences in support of our claim that, similarly to the mean-field limit, in the position-fixed limit the nonequilibrium phase transition between the empty and active phases is discontinuous for small $`\stackrel{~}{s}`$ and so there is a tricritical point in the phase diagram of the model are threefold: First, the vanishing of the exponent $`\eta `$ in this range signalizes a distinct asymptotic behavior of the average number of replicators $`n\left(t\right)`$. Second, the hyperscaling relation (16) is clearly violated for small $`\stackrel{~}{s}`$ while it is satisfied in the regime where we expect the transition to be a second-order one. Third, using Eq. (15) and the data of table 2 we find $`d_fd=2`$ in the region of small $`\stackrel{~}{s}`$ indicating that the clusters of replicators are not fractal objects, in contrast to the clusters observed in the vicinity of a second-order transition (we find, for instance, $`d_f1.21`$ for $`\stackrel{~}{c}=0`$). This point is illustrated in Figs. 12 and 13 which show snapshots of typical colonies at $`t=10^4`$ for the two extreme cases. In both figures the relative distances to the critical points are the same. The reason for the colonies of obligatory replicators to be much denser than those of malthusian replicators is that the number of replicators per surviving colony ($`N\left(t\right)t^{\eta +\delta }`$) increases roughly as $`t`$ for the former and as $`t^{0.7}`$ for the latter. In addition, the average square distance over which the initial colony has spread from the center of the lattice at time $`t`$ is roughly of the same order in both extremes as indicated by the values of the exponent $`z`$.

Finally, in Fig. 14 we present the phase diagram for the position-fixed limit. In this case there are only two phases, namely, the empty phase $`\left(E\right)`$ characterized by a vanishing probability of survival $`P_{\mathrm{}}lim_t\mathrm{}P\left(t\right)=0`$ and the phase $`\left(EA\right)`$ where the active and empty states can occur with probabilities $`P_{\mathrm{}}`$ and $`1P_{\mathrm{}}`$ (see Figs. 6 and 8), respectively. We note that the population size is effectively unlimited so that extinction is not certain to occur as in the steady-state analysis of finite systems. It is worth emphasizing that for finite production rates one has $`P_{\mathrm{}}<1`$ and so there is always a nonvanishing probability of extinction. (Actually, $`P_{\mathrm{}}`$ depends on the number and location of the replicators in the initial colony, but the conditions we have chosen are the most relevant for the emergence of life problem.) This result contrasts to our findings in the mean-field limit that the fixed point associated to the empty phase becomes unstable for $`\stackrel{~}{s}>1`$ so that even starting with a vanishingly small concentration of replicators the population never dies out. As mentioned before, the reason for this discrepancy is not the difference in mobility of the replicators but the fact that the mean-field analysis actually considers an infinite population and hence it fails to take into account the stochastic fluctuations that could drive a small population to extinction.

For sake of completeness, we should mention that we have also carried out a similar analysis for one-dimensional lattices (chains). While the results for the mean-field limit are of course the same (provided we properly redefine the dimensionless parameters $`\stackrel{~}{c}`$ and $`\stackrel{~}{s}`$), the fixed-point limit has some distinct features that are worth-mentioning. In particular, we find no evidence for a first-order transition; instead, we find that in both extremes $`\stackrel{~}{c}=0`$ and $`\stackrel{~}{s}=0`$ the empty and active phases are separated by a second-order phase transition that belongs to the $`\left(1+1\right)`$ directed-percolation universality class . Interestingly, in one-dimension the steady-state analysis of finite chains is rendered practically useless by the very pronounced finite-size effects, which are probably due to the proximity to the lower critical dimension of the model.

## 5 Conclusion

In this paper we have focussed on the prior step in the evolution of life: What are the necessary conditions for small colonies of molecules capable of making copies of themselves via some template mechanism to persist? This step must be passed before one can consider issues such as the outcome of the competition between the replicators and their defective copies  or between different kinds of replicators . As a model of replicator we have considered the well-established hypercyclic replicator (one-membered hypercycle) which incorporates two independent mechanisms of replication, namely, the direct template replication reaction (1) and the catalytically-assisted template replication reaction (3), whose rates are proportional to the parameters $`s`$ and $`c`$, respectively. Furthermore, motivated by the modern theories for the evolution of life that suggest a scenario of diffusion-controlled chemical reactions taking place on adsorbing surfaces (probably pyrite) where each reactant can move randomly on the surface , we have considered a two-dimensional lattice model where each reactant can occupy one of the lattice cells.

Since the diffusion process of reactants complicates considerably the analysis, we have focussed on the two extreme situations: the infinite diffusion or mean-field limit and the position-fixed or contact process limit. The expectation is that features common to both limits should also be present in the more realistic, intermediate situations. In both cases we found rich phase diagrams showing the regions in the plane $`(c,s)`$ where the replicators persist (active phase) and die out (empty phase): these regions are separated by second-order nonequilibrium phase transitions which turn into first-order ones at tricritical points. The dynamic Monte Carlo method has proven very well suited to our investigation of the position-fixed limit not only because the method is based on the analysis of the spreading behavior of a small colony of active cells, which is exactly the problem we are interested in, but because, rather surprisingly, it allows an unambiguous identification of the order of the nonequilibrium phase transition. In addition, we show that the continuous transition is in the universality class of the $`\left(2+1\right)`$ directed percolation.

Some remarks on the apparent similarity between Schlögl’s first and second models and the replicator models studied in this paper are in order . In fact, irreversible versions of Schlögl’s first and second models are recovered when the reactant $`E`$ is eliminated from reactions (1) and (3), respectively, so that the existence of an empty cell containing source materials is not required for replication. Interestingly, this difference is not important in the case of malthusian replicators ($`c=0`$) since this model has the same critical behavior as Schlögl’s first model, namely, a second-order phase transition which is in the same universality class of the $`\left(d+1\right)`$ directed percolation (the mean-field limit is obtained for $`d4`$) . The comparison between Schlögl’s second model and obligatory sexual replicators ($`s=0`$) is more involved. On the one hand, the mean-field analysis of Schlögl’s second model predicts a first-order phase transition  but Monte Carlo calculations indicate that the transition for $`d<4`$ is a second-order one which is actually in the same universality class as the transition in Schlögl’s first model . On the other hand, our results for the case $`s=0`$ show that a directed-percolation-like, second-order phase transition takes place for $`d=1`$ only, the transition being discontinuous for higher dimensions ($`d2`$). We should mention, however, that the Monte Carlo implementation of Schlögl’s second model actually allows the diffusion of reactants to neighboring cells and, in addition, allows a cell to shelter more than one reactant  so that a comparison with the position-fixed limit may not be appropriate. In any event, it is our opinion that the hypercyclic replicator model should not be viewed as a mere variant of Schlögl’s models; rather it is a well-established model of chemical evolution  which, as far as we know, has not been studied beyond the mean-field limit.

Although the simple replicator model considered in this paper turned out to be a quite exciting model of nonequilibrium phase transitions, we should not lose sight of the original purpose of this work and so, at this stage, it is important to highlight the relevance of our results to the origin of life issue. In fact, the mere existence of a phase transition between the empty and the active regimes poses a difficulty to our scenario of the emergence of life since the scaled production rate $`\stackrel{~}{s}`$ of the spontaneously created self-replicating molecule must be larger than some threshold value already at the outset. Though increasing the mobility of the reactants decreases this threshold somewhat, that scenario would be more plausible if replicators with vanishingly small production rates could also thrive. The situation becomes even worse in the case of first-order transitions: in the deterministic mean-field limit the initial abundance of the spontaneously created replicators should be large as well, while in the stochastic position-fixed limit the probability of survival in the vicinity of the transition point is some orders of magnitude smaller than in the case of a second-order transition (see Figs. 6 and 8). Furthermore, our results indicate that some important conclusions, such as the certainty of survival for $`\stackrel{~}{s}>1`$ or the role played by the initial concentration of replicators near a discontinuous transition, are actually artifacts of the deterministic formalism commonly used to study chemical evolution.

In summary, our results show the necessity of adding some new elements to the standard scenario for the emergence of life whose effect would be to avoid the phase-transitions, allowing thus inefficient replicators to thrive at this first stage of life. Only then one can invoke natural selection and imperfect replication to boost the replication rates. In addition, our results point to the inadequacy of the deterministic mean-field or chemical kinetics formulation to address the origin of life issue and suggest as an alternative stochastic formulation the dynamic Monte Carlo method that has been extensively used in the physics literature to study nonequilibrium phase transitions.

The work of J.F.F. is supported in part by Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP), Proj. No. 99/09644-9. C.P.F. is supported by FAPESP.


text:
# Acknowledgement.

## Acknowledgement.

During the preparation of this paper we became aware from G. Dvali that he is also working out a similar proposal for the solution of the strong CP problem. We thank him for correspondence. This work was supported by the Academy of Finland under the Project No. 163394.


text:
# Information preserved guided scan pixel difference coding for medical images

## I Introduction

This research attempts at developing an information preserving image coding suitable for handling medical images which usually do not allow any degree of image degradation in comparison with the source images saghri ; takaya . JPEG, one of the well accepted compression standards, exhibits an excellent compression gain, but does not insure that the decoded result is the exact replica of the original. A variation of DPCM or pixel difference coding, namely ”Guided Scan Pixel Difference Coding was studied. A major modification to the classical DPCM is the inclusion of scan information, in terms of the direction used to calculate the pixel difference and step size. This additional information does not significantly increase the volume of coded data. It actually opens the extensibility for the DPCM to be able to handle multi-dimensional image data, in a similar manner as JPEG is extended to MPEG to handle movie presentation. Another important asp ect of medical image coding, which is the ability to access the full details of a local image (often referred to as ROI, Region of Interest), is also satisfied. Neither the location of a ROI nor the scanning method is restricted so that a local area of ROI can be selected anywhere in the original image. The flexibility in scanning makes it possible to use the same code for real time movie presentation in the same way MPEG uses JPEG. This research considers two stage redundancy removal, one by taking pixel differences and the other by an entropic source coding similar to the Huffman coding to achieve a redundancy removal of approximately 4-5 bpp (bits per pixel).

The sizes of medical images have increased with the advancement in various medical imaging modalities, typically in MRI and X-ray CAT scan. MRI can scan a whole body within a reasonable time of 30 minutes and X-ray CAT has adopted helical scanning for a higher spatial resolution. In the case of 3D scanning, the image size could be as large as 100 million pixels. Images of this size are usually stored on writable laser disks (CD-WORM), but this imposes serious problems when the retrieval of such data is attempted through a communication channel in the wide area network (WAN) environment using 64 kbits/sec bit-rate telephone channels.

The problems associated with the compression of medical images are two fold. The first requirement is that compression must be absolutely non-degradable. Regardless of whether a picture has been affected by measurement noise that occurs in the process of physical measurement and in image reconstruction, the source data acquired at an imager must be retained without any kind of losses. The second condition is the integrity of the objects contained in a picture data. Physicians examine very closely the image in a ROI (Region of Interest) or a number of ROI’s. A selected ROI needs to be immediately accessible and it must be translated into a serial stream of compressed data for transmission via a communication channel.

## II Technical Background

The method is best represented by the name Guided Scan Pixel Difference Coding takaya , a variation of DPCM widely used in compressing voice signals. It is well known that a PCM voice signal of 6-8 bits can be compressed down to a 3 bit DPCM. Even discounting the fact that DPCM uses a nonlinear quantization scale having a greater step size for a larger signal level, it is easy to remove 3 bits per code word, reducing the total volume down to half of the original. In comparison, JPEG performs a compression rate of 10 and is better than DPCM, because it permits the encoded image to degrade somewhat from the original. DPCM is a non-degradable source coding if the quantization step size is fixed to the same linear scale as already done in the source images. Furthermore, the encoding/decoding method is basically the simple operation of subtraction/summation. Further compression to remove the redundancy leftover from DPCM is accomplished by an entropy coding such as Huffman code.

In order to satisfy the other aspect of medical image compression, accessibility to a ROI, additional information describing how an individual code representing a pixel difference is derived, must be attached at the expense of compression gain. In the case of 2D images, there are four (up, down, left and right) or eight (if diagonal differences are considered) possible directions by which a pixel difference can be calculated. If an ROI is specified by a person viewing an image specifying several guide marks surrounding the ROI, the length of each raster along a specified scan direction must be registered. Also, the step size used to calculate the pixel difference must be included, if coarsely sampled images need to be transmitted for viewing an overall image at a reduced image quality. These additional pieces of information could be 4 bits for directions and about 3 bits to indicate the step size per picture. Information regarding a length of scanning depends on how irregular a shape is set to be an ROI. Nevertheless, no significant increase in the total volume of coded words is expected. The term ”Guided Scan is meant to include the information as to how an image or a ROI is scanned. The way of scanning an image should be left for the image coding system to determine, so that the total volume of information can be minimized for individual selected ROI’s. Another interesting aspect of the Guided Scan DPCM is the structural similarity observed in the forward or backward difference interpolation algorithm. Since the tree of high order differences, necessary for the interpolation, can be calculated easily from the first order differences coded by the Guided Scan DPCM, there is potential to artificially increase the spatial resolution of the source image without distorting the image.

## III Information Entropy

The information entropy is defined by:

$$H=\underset{i=1}{\overset{I}{}}p_i\text{ log}_2p_i$$
(1)
where $`p_i`$ is the probability associated with the occurrence of a value i when the total number of values used is $`I`$. If eight bits are used to code an image and the probability for any number between 0 and 255 to occur is the same, the information entropy is 8. If only one value, say 128, occurs all times, the entropy is zero. Adjacent pixels in a digitized image are highly correlated. If two adjacent pixels are considered, the probability for the first pixel $`x`$ to take a value i, that for the second pixel $`y`$ to take j, and the joint probability for the pixel $`x`$ to take a value i and the pixel $`y`$ to take j are given respectively as follows:

$$p_i=\frac{N_i}{N},p_j=\frac{N_j}{N}\text{ and }p_{ij}=\frac{N_{ij}}{N}$$
(2)
The entropies based on the joint probability and the conditional probabilities are then obtained as follows:

$`H(x,y)`$ $`=`$ $`{\displaystyle \underset{i,j}{}}p_{ij}\text{ log}_2p_{ij}`$
$`H(x|y)`$ $`=`$ $`{\displaystyle \underset{i,j}{}}p_{ij}\text{ log}_2{\displaystyle \frac{p_{ij}}{p_j}}`$
$`H(y|x)`$ $`=`$ $`{\displaystyle \underset{i,j}{}}p_{ij}\text{ log}_2{\displaystyle \frac{p_{ij}}{p_i}}`$ (3)
We also know that $`H(x,y)=H(x)+H(y|x)=H(y)+H(x|y)`$. Since theoretically:

$$H(y|x)H(y)0$$
(4)
The entropy calculated for the second pixel $`y`$ knowing the occurrence of the first pixel $`x`$ is smaller than the entropy calculated from $`y`$ alone. The uncertainty coefficient of $`y`$

$$U(y|x)=\frac{H(y)H(y|x)}{H(y)}$$
(5)
is indicative of the dependency of $`y`$ on $`x`$. The symmetrical uncertainty:

$$U(x,y)=\frac{H(x)U(x|y)+H(y)U(y|x)}{H(x)+H(y)}$$
(6)
yields 0 if $`x`$ and $`y`$ are completely independent and 1 if they are completely dependent.

## IV Entropy Calculated from Stretched Exponential PDF

When the probability density function (PDF), obtained from a histogram of pixel intensities, is well approximated by the stretched exponential probability density function mallat91 ; mallat92 ; mallat89 :

$$p(x)=Ke^{(\frac{|x|}{\alpha })^\beta }\text{ with }K=\frac{\beta }{2\alpha \mathrm{\Gamma }(\frac{1}{\beta })}$$
(7)
the theoretical entropy can be calculated from:

$$H=2K\frac{\alpha }{\beta }\frac{(lnK)\mathrm{\Gamma }(\frac{1}{\beta })+\mathrm{\Gamma }(1+\frac{1}{\beta })}{ln2}$$
(8)
After evaluating the absolute mean $`\mu `$ and standard deviation $`\sigma `$, the parameters $`\alpha `$ and $`\beta `$ are extracted from the equation press :

$$F(\beta )=\frac{\mu ^2}{\sigma ^2+\mu ^2}$$
(9)
where, the function $`F`$ is defined by:

$$F(x)=\frac{\mathrm{\Gamma }(\frac{2}{x})^2}{\mathrm{\Gamma }(\frac{1}{x})\mathrm{\Gamma }(\frac{3}{x})}$$
(10)
and $`\alpha `$ is defined by:

$$\alpha =\frac{(\sigma ^2+\mu ^2)\mathrm{\Gamma }(\frac{1}{\beta })}{\mathrm{\Gamma }(\frac{3}{\beta })}$$
(11)
Fig. 1 depicts the variation of the entropy $`H`$ as a function of $`\alpha `$ in the interval \[0.1-5\] for a given $`\beta `$ picked in the interval \[0.8-1.5\].

## V Preliminary Studies on Information Entropy

A 3-D MRI brain section image which consists of 64 slices was analyzed with an attempt to find the lower bound of the information entropy for medical images takaya . The lower bound is a definite measure that tells exactly how many bits can be removed per pixel as redundancy. One slice of the 3-D image is shown in Fig. 2. Two different approaches are taken to calculate the entropies for the image shown in Fig. 2.

1. Use Eq. 1 which does not consider the dependency between adjacent pixels for n-th order difference images (n = 0,…,10).
2. Use Eq. 3 which takes the dependency of a pixel on its immediate neighbours into account (equivalent to a first order Markov model).

The first approach finds a crude information entropy, if it is applied to the original picture since the redundancy due to pixel-to-pixel correlation is still intact. In order to find a more accurate information entropy, it is necessary to remove the redundancy by some means that allows the recovery of the original image. A difference image is introduced for this purpose. When a 2-D image is denoted by $`A`$ and $`z^1`$ represents the one-step shift operator towards the positive side of the horizontal axis (or vertical or diagonal axis), the first difference is:

$$D_1=Az^1A=(1z^1)A$$
(12)
The second difference is given by

$$D_2=(1z^1)D_1=(1z^1)^2A$$
(13)
Thus the n-th difference is given by

$$D_n=(1z^1)^nA$$
(14)
In order to make reconstruction possible, the first column of the first difference must retains the first column of the original if horizontal shift is used. The second difference must retain the first column of the original and the second column of the first difference in its second column. Additional DC restoration columns must be progressively added when a new difference image is created. Successive applications of this difference operator remove the pixel-to-pixel correlation and resulting images become gradually more random. Table 1 shows how the entropy associated with such a difference image varies when n is increased.

2-D auto-correlation functions for the difference images of n= 0, 1 and 2 are shown in Fig. 3. The top figure is the 2-D auto-correlation of the original image. Comparing the middle n =1 and the bottom n =2, it is observed that the central peak of n=2 is sharper (less correlated) than that of n=1. As for the entropies calculated, the entropy drastically drops to the minimum at the very first difference operation then it starts increasing as n increases. Successive difference operations seem to decorrelate the image and make it more random, but the entropy monotonously increases after the first difference. This phenomenon can be explained from the frequency response of the n-th difference operation described by the transfer function of a high pass filter,

$$G(z)=(1z^1)^n$$
(15)
Since the magnitude response is given by:

$$|G(e^{j\omega })|=2^nsin^n(\frac{\omega }{2})\text{ where }0\omega <\pi $$
(16)
and the bandwidth becomes wider as n increases. The power of the n-th difference image is greater than that of (n-1)th difference image, and so is the variance.

The second approach using Eq. 3 considers the dependency between adjacent pixels in calculating information entropies so that it is no longer necessary to decorrelate images. Since the first difference achieves the minimum entropy value, Eq. 3 is applied to the image of the first difference. Table 2 summarizes the results. The entropy $`H(x,y)`$ based on the joint probability $`p_{ij}`$ is divided by 2 to translate it into the entropy per pixel. This value is significantly smaller than $`H(x)\text{ or }H(y)`$.

Another interesting result is found by fitting the stretched exponential PDF daub  to the histogram of the original image shown in Fig. 4. The entropy for the large peak located in the lower gray scale range of $`<50`$ is calculated to be 2.749977 without including the range $`>50`$. This is nearly one half of the total entropy found for the original image. A whole picture of a medical image usually contains a significant portion of dark background.

## VI Discussion

Reviewing the experimental results presented in the previous sections, some conclusive remarks can be made to determine the strategy for information preserving source coding. Information preserving coding or lossless coding generally means that the picture received is the picture sent in terms of the bit structure of the image, not in terms of the visual impression of the image before and after image compression/transmission. In this narrow sense of lossless coding, altering pixels is prohibited. No alterations in quantization are permitted. With these strict constraints, the only possible source of redundancy that could be removed is limited to the pixel-to-pixel correlation. Pixel difference coding (DPCM) alone brings the information entropy down to its near minimum. As seen in the entropy $`H(x,y)`$, the result of the pixel difference coding can be further trimmed but not to a large extent. The MR images analyzed are all 8 bit images. According to the sample calculations shown in this paper and other tested results, a bare-bone information entropy per pixel is slightly greater than 4 bits/pixel. An optimistic compression ratio CR is therefore CR=2:1, as long as a near optimum compression algorithm, typically Huffman coding, is used. Blending the Huffman code and the codes which consider state transitions that frequently occur within a near zero range of the pixel difference scale, a small improvement in compression ratio will be made. For example, a run of successive zeros up to a length of 10 can be treated as a code word if the frequencies of such occurrences are sufficiently high.

Further improvement of CR requires removing a constraint on quantization. As used in DPCM in voice signal coding, there are several methods to set up a quantization scale which minimizes the information entropy. If the methods of transform coding are allowed, CR can be improved drastically. Approaches to control the quality of medical images, for example texture, appearance of speckles, maintaining repeated basic patterns as fractal images do, etc… If one coding method can assure the fidelity of a certain image attribute, it might be considered as a better coding scheme. The medical community may be prepared to accept it as a better replacement for the totally reconstructible information preserving coding.

If some loss is permitted in medical image coding, the n-th order difference images discussed in this paper has another potential application for image coding. Recalling that the power spectrum $`S(\omega _x,\omega _y)`$ of an n-th order difference image $`D_n(x,y)`$ relates to its auto-correlation $`R(x,y)`$ with the 2-D Fourier transform,

$`S(\omega _x,\omega _y)=R(x,y)`$
$`=\{D_n(x,y)\}\{D_n(x,y)\}^{}`$ (17)
it is apparent that the magnitude information of $`D_n(x,y)`$ is contained in the highly concentrated peak of $`R(x,y)`$. Since the phase information is lost in calculating $`R(x,y)`$, it is necessary to preserve the phase of $`D_n(x,y)`$. The auto-correlation of $`\mathrm{}D_n(x,y)`$ shown in Fig. 5 indicates that the magnitude of phase is also highly concentrated at zero of the 2-D phase auto-correlation. The n-th order difference image produces a highly concentrated auto-correlation both in magnitude and in phase. By preserving the profile of the peak and discarding the rest, it seems possible to achieve a high compression gain. It is however known that image degradation is usually enhanced by inaccurate phase estimation.

Acknowledgement Li Yuan acknowledges support through a TRlabs fellowship. C. Tannous acknowledges support from NSERC president’s fund grant.

Figure Captions

* Theoretical entropy calculated for stretched exponential probability density functions. The graphs show the entropy versus $`\alpha `$ with $`\beta `$ as a parameter.
* A Slice of a 3-D MRI Head Image (256x256)
* 2D Auto-correlation, original, 1st difference, 2nd difference image and phase image of the original (from top to bottom). Top vertical scale is in units of $`10^7`$ , middle is in units of $`10^5`$ whereas bottom part is in $`10^6`$ units.
* Histograms for the original, 1st difference and 2nd difference images.
* Auto-correlation of the phase image. Vertical scale is in units of $`10^4`$.