format: text

text:
# Zero Temperature Properties of RNA Secondary Structures

## I Introduction

RNA plays a fundamental role in the biochemistry of all living systems BRANDEN\_TOOZE , and it is commonly believed to be at the origin of the pre-Darwinian epoch of life RNA\_WORLD . Much like for DNA, the RNA primary structure can be described in terms of strings of the four letter alphabet composed by Adenosine, Citosine, Guanine, Uracile (ACGU). Since RNA is usually found in the single stranded pattern, formation of double-helix regions is accomplished by the molecule folding back onto itself to form Watson-Crick (WC) base pairs G$``$C and A$`=`$U, or the slightly less stable G$``$U pair. One of the most intriguing features of RNA folded secondary structures is that in most cases the connectivity graph is planar: This property greatly reduces the computational efforts needed for calculating the ground state structure.

It might be asked whether secondary structures provide an adequate level of description for RNA real molecules TINOCO\_BUSTAMANTE . It is believed that secondary structure description is biologically relevant for a number of reasons: Base pairing and base pair stacking provide the major part of the free energy of folding VIENNA\_PACKAGE , secondary structures have been used successfully by biologists in the interpretation of RNA function and activity RNA\_WORLD , and because structures are conserved in evolutionary phylogeny. At the same time computer scientists find this level of description rather appealing since secondary structures are discrete and therefore easy to compare. Moreover, thanks to the planarity condition, efficient recursive algorithms for the computation of the native (ground state) structure are easily implemented NUSS\_JACOB .

Beside the genuine biological interest of RNA models, recently this subject has raised considerable attention as an intriguing problem in statistical mechanics of disordered systems. The focus is now set on the presence and the nature of a low temperature phase in ensembles of random sequences. In a series of recent papers HIGGS\_PRL ; NOSTRO\_PRL ; COMM\_HARTMANN ; NOSTRO\_REPLY ; BUND\_HWA\_1 ; BUND\_HWA\_2 ; KMM\_RNA  different authors have presented a number of evidences (mostly numerical) supporting the existence of a transition to a glassy phase (for a review see Ref. HIGGS\_REVIEW ). While a careful study of the equilibrium thermodynamics of the model suggests a smoother than second order phase transition NOSTRO\_PRL , there is still much debate about the nature of the low temperature phase, since finite size scaling corrections are very hard to keep under control. It has been shown in NOSTRO\_PRL ; NOSTRO\_REPLY  that, at least for systems up to 1000 bases, a broad overlap distribution characterizes the low temperature phase, but a safe extrapolation to the thermodynamic limit from the available data is still out of control COMM\_HARTMANN ; NOSTRO\_REPLY .

Bundschuh and Hwa BUND\_HWA\_1 ; BUND\_HWA\_2  have presented an extensive study on a variety of similar RNA models supporting the existence of a low temperature glassy phase. They were able to show analytically, via a two replica calculation, that weak quenched sequence disorder is equivalent to a high temperature phase in which all replicas are independent (molten phase), and that there must be a finite temperature below which replicas start feeling themselves as in a strong coupling regime. Numerically the authors established this glassy transition measuring the free energy cost of imposing a pinch between two bases. They observed that the energy of the pinching excitation (with respect to the ground state) increases with the sequence length following a logarithmic law (even if a power law with small exponent was not excluded).

In this paper we study the scaling regime of the lowest energy excitations in different models of random RNA secondary structures. Following an idea put forward in NOSTRO\_REPLY , we use a perturbing method CPPS ; MEZEMEZ  which has been very valuable in the study of low temperature properties of disordered systems EPS\_METHOD : In the following we will call it the $`\epsilon `$-coupling method. Very recently the same procedure has been followed by Krzakala, Mézard and Müller KMM\_RNA . We will comment on their results in the concluding section after having presented our data.

The goal of the $`\epsilon `$-coupling method is to calculate the energy cost of typical excitations above the ground state involving a finite fraction of the system. As in the droplet model DROPLET  these energy excitations are assumed to scale as $`\mathrm{\Delta }E(L)L^\theta `$, $`L`$ being the length of the molecule, and $`\theta `$ being a relevant exponent, that we would like to determine. It is well possible that also the local pinching of coupled replicas imposed in BUND\_HWA\_1 ; BUND\_HWA\_2  would generate “typical” configurations, but we believe that our method based on a bulk perturbation will surely do so (in this sense we find very illuminating its application). We apply a perturbation that is simply a repulsion term from the ground state structure, which forces the system to find new low-energy structures far away from the one of the original ground state, without any other constraint.

The paper is organized as follows: in section II we introduce the formalism for describing the RNA secondary structure and the different ways we have used to remove the ground state degeneracy intrinsic to the original model NOSTRO\_PRL . In section III we sketch the method we have applied for calculating the low-energy spectrum of the model. We also discuss the measurable observables. In section IV we present our results, focusing on the differences and similarities among the models we have introduced. Finally, in section V, we summarize our findings, we compare with previous work and we comment on further developments.

## II Models

The secondary structure of RNA is the set of base pairs that occur in its three-dimensional structure. Let us define a sequence of basis as $`\{r_1,r_2,\mathrm{},r_n\}`$, $`r_i`$ being the $`i^{th}`$ base of the chain and $`r_i\{A,C,G,U\}`$. A secondary structure on $``$ is now defined as a set $`𝒮`$ of $`(i,j)`$ pairs (with the convention that $`1i<jL`$) satisfying the following rules:

* $`ji4`$ : this restriction permits flexibility of the chain in its three-dimensional arrangement;
* two different base pairs $`(i,j),(i^{},j^{})𝒮`$ if and only if (without loss of generality we can assume that $`i<i^{}`$) $`i<j<i^{}<j^{}`$, i.e. the pair $`(i,j)`$ precedes $`(i^{},j^{})`$, or $`i<i^{}<j^{}<j`$, i.e. the pair $`(i,j)`$ includes $`(i^{},j^{})`$. This rule, called planarity condition, excludes the occurrence of the so-called pseudo-knots, which are very unlikely in real RNA.

We consider a simplified model for RNA folding, very similar to the one studied in NOSTRO\_PRL . The model is described in terms of the Hamiltonian:

$$=\underset{(i,j)𝒮}{}e_{ij}=\underset{(i,j)}{}e_{ij}\mathrm{}_{ij},$$
(1)
where $`e_{ij}`$ is the pairing energy between bases $`i`$ and $`j`$ and the variable $`\mathrm{}_{ij}`$ takes value 1 if $`(i,j)𝒮`$ and 0 otherwise. On a first approximation one can assume that the pairing energies depend only on the paired bases, $`e_{ij}=e(r_i,r_j)`$. Reasonable values for the energies $`e(r_i,r_j)`$ of the allowed base pairs (C-G, A-U and G-U) at room temperature are of $`𝒪(1)`$ kcal/mole 111In the models studied here we do no try a quantitative, detailed comparison to experimental results, and so we only select pairing energies of the correct experimental order of magnitude.. One could consider other phenomenological parameters in order to take into account the whole complexity of a realistic energy function VIENNA\_PACKAGE ; REALISTIC\_PARAM .

We have assumed a drastic approximation in order to get a tractable model both from a numerical and analytical point of view. We consider sequences made of $`4`$ symbols (A, C, G and U) and we assume that only Watson-Crick base pairs may occur: we use a strong C-G coupling of energy $`2`$ (in arbitrary units) and a weak A-U coupling of energy $`1`$. All the other possible couplings increase the energy, so that the system avoids these links. One of the advantages of this model is that the role of the disorder (encoded in the random sequence $``$) is clearly separated from that of the frustration (induced by the planarity condition on the structure $`𝒮`$).

This $`4`$ letters model has an exponentially large ground state degeneracy, which gives a finite $`T=0`$ entropy (as already found in NOSTRO\_PRL  for the $`2`$ letters model). For this reason we refer to it as the degenerate model (the D model).

The large ground state degeneracy that occurs in this D model is a pathology of a frustrated models with simple discrete interactions: Since the couplings can take only the two negative values $`2`$ and $`1`$ the same exact energetic situation can be realized in many ways. This accidental degeneracy will probably not play a relevant role in the physical RNA: Since real RNA energy function is far more complex than that, ground state degeneracy is unlikely to occur. Because of that we define two new models with modified pairing energies, in order to remove the degeneracy. In both models this aim is accomplished by adding a small random perturbing term $`\eta _{ij}`$ to the pairing energies: $`e_{ij}e_{ij}+\eta _{ij}`$.

In the quasi-degenerate model (the QD model) the $`\eta _{ij}`$ are i.i.d. variables extracted from a Gaussian distribution of zero mean $`\eta =0`$ and variance $`\eta ^2=\eta _0^2/L`$, with $`\eta _0`$ a small and finite constant of the order of $`0.1`$: When $`L\mathrm{}`$ the pairing energies are modified of an infinitesimal amount. The variance is chosen such that the energies of the ground states (which are degenerate for $`\eta _0=0`$) are split over an $`𝒪(\eta _0)`$ range. In this way we preserve somehow the structure of the original energy spectrum and the sequence still plays a key role, but the unphysical degeneracy is lifted and the ground state is now unique.

In the non-degenerate model (the ND model) the variance of the $`\eta _{ij}`$ variables is finite, $`\eta ^2=\eta _0^2`$ (we use $`\eta _0=0.1`$). This variance induces an $`𝒪(\sqrt{L})`$ splitting of the degenerated ground states, which has to be considered as a strong reshuffling of the original energy spectrum, since the energy gaps among levels in the original model were of $`𝒪(1)`$. The resulting energy landscape now depends very little on the sequence. Because of that the ND model is very similar to the “Gaussian disorder model” (the GD model), already discussed in BUND\_HWA\_2 , where the $`e_{ij}`$ are i.i.d. Gaussian variables of zero mean and unitary variance. In this model the sequence plays no role.

## III Methods

The $`\epsilon `$-coupling method we use to calculate low-energy excitations is the one already used in EPS\_METHOD  and in KMM\_RNA . It works as follows: First of all one calculates the ground state structure $`\mathbf{}_0=\{\mathrm{}_{ij}^{(0)}\}`$ which minimizes $``$. Then one adds a perturbation to the Hamiltonian, $`^{}=\epsilon (1q)`$, where $`q\frac{1}{L}_{ij}\mathrm{}_{ij}\mathrm{}_{ij}^{(0)}=\frac{1}{L}\mathbf{}\mathbf{}_0`$ is the overlap with the ground state structure (note that, with this definition, the overlap is always positive). The perturbation term penalizes the structures which are close to the ground state $`\mathbf{}_0`$ and thus acts as a repulsive term in the space of structures. Finally one calculates the ground state structure of $`^{}`$ for many values of $`\epsilon `$. Let us call these new structures $`\mathbf{}_\epsilon `$.

By definition, for any disorder realization $`𝒥=\{,𝜼\}`$, both the distance $`d_𝒥(\epsilon ,L)=1\frac{1}{L}\mathbf{}_\epsilon \mathbf{}_0`$ as well as the energy difference $`\mathrm{\Delta }E_𝒥(\epsilon ,L)=(\mathbf{}_\epsilon )(\mathbf{}_0)`$ between $`\mathbf{}_\epsilon `$ and $`\mathbf{}_0`$ are non-decreasing functions of $`\epsilon `$. Moreover $`\mathrm{\Delta }E_𝒥(\epsilon ,L)<\epsilon `$, since the Hamiltonian has been perturbed by a term whose absolute value is less than $`\epsilon `$, and the structures $`\mathbf{}_\epsilon `$ are then low energy excited states of the original Hamiltonian. We will indicate without the $`𝒥`$ subscript the observables averaged over the quenched disorder $`𝒥=\{,𝜼\}`$: $`d(\epsilon ,L)=\overline{d_𝒥(\epsilon ,L)}`$ and $`\mathrm{\Delta }E(\epsilon ,L)=\overline{\mathrm{\Delta }E_𝒥(\epsilon ,L)}`$.

The algorithm for finding the new ground states of $`\epsilon `$-coupled system is exactly the same one used for the original Hamiltonian: The repulsion from the first ground state is included by modifying the values of the original pairing energies $`e_{ij}`$.

In the thermodynamical limit structures differing by a finite $`\mathrm{\Delta }E`$ have the same intensive energy, and one could try to understand how they are organized in the configurational space. An interesting question is whether, in the large $`L`$ limit, these structures are extremely close together or spread over finite distances. The answer to this question can be given in terms of the asymptotic quantity

$$d_{\mathrm{}}(\epsilon )=\underset{L\mathrm{}}{lim}d(\epsilon ,L),$$
(2)
which is again a non-decreasing function of $`\epsilon `$. If $`d_{\mathrm{}}(\epsilon )=0`$ for any finite $`\epsilon `$ then structures with the same energy are close together, while if $`d_{\mathrm{}}(\epsilon )>0`$ for $`\epsilon >\epsilon ^{}𝒪(1)`$ then structures with the same intensive energy may have a broad probability distribution functions of their distances and overlaps.

In the case where $`d_{\mathrm{}}(\epsilon )=0`$ we can derive a relation describing the way $`d(\epsilon ,L)`$ vanishes. We assume, as in the droplet model DROPLET , that the energy cost of a typical excitation involving a finite fraction of the system (i.e. having finite $`d`$) scales with the system size as

$$\mathrm{\Delta }E_{\text{typ}}L^\theta .$$
(3)
We call $`\mathrm{\Pi }(\mathrm{\Delta }E,d,L)`$ the probability distribution (over the disorder) of excitations with energy $`\mathrm{\Delta }E`$ and size $`dL`$ in systems of size $`L`$. For any fixed and finite $`d(0,1]`$, we assume that $`\mathrm{\Pi }(\mathrm{\Delta }E,d,L)`$ has a finite weight in $`\mathrm{\Delta }E=0`$, and so, for normalization reasons, we must have $`\mathrm{\Pi }(0,d,L)=c(d)L^\theta `$ for large $`L`$ (unless there is a delta function in $`\mathrm{\Delta }E=0`$ as in the D model), where $`c(d)`$ is a smooth function in the scaling region.

Once we add the perturbing term $`\epsilon d`$ to the Hamiltonian, an excitation of size $`d`$ will be activated only if its energy satisfies $`\mathrm{\Delta }E<\epsilon d`$. Thus the average distance of the new ground state structure is given by

$$d(\epsilon ,L)=_0^1y\text{d}y_0^{\epsilon y}\mathrm{\Pi }(x,y,L)\text{d}x=\epsilon L^\theta _0^1y^2c(y)\text{d}y,$$
(4)
for small $`\epsilon `$ and large $`L`$ 222The D model has $`\theta =0`$, since, for a finite range of $`d`$ values, its $`\mathrm{\Pi }(\mathrm{\Delta }E,d,L)`$ has a delta function in $`\mathrm{\Delta }E=0`$ that implies the presence of a gap (energy levels are discretized), so that the integral in Eq. (4) is finite for any $`L`$ and $`\epsilon `$ small enough..

Then we can evaluate the $`\theta `$ exponent by two independent ways:

* from Eq. (4), $`d(\epsilon ,L)\epsilon L^\theta `$, by measuring the average distance $`d(\epsilon ,L)`$ for a fixed small $`\epsilon `$ as a function of the system size $`L`$;
* from Eq. (3), which can be equivalently rewritten for the average energy difference as

$$\mathrm{\Delta }E(d,L)L^\theta ,$$
(5)
by measuring the average energy difference for a fixed distance (not fixed $`\epsilon `$ !) as a function of the system size $`L`$.

## IV Results

We study zero-temperature properties of the models described above, i.e. we analyze ground states structures (GSS) of the original and of the perturbed Hamiltonians. We start showing the data for the $`T=0`$ overlap distribution in the D model (the only one with many different ground states). After that we present the results obtained with the $`\epsilon `$-coupling method for all the models defined in section II.

### IV.1 The D Model

The D model possesses an exponentially large number of GSS, which form a set that we call $`𝒢`$. In order to understand how they are distributed in the space of structures one can calculate the probability distribution function of the overlap, which is defined, for any pair of structures, as $`q\frac{1}{L}_{ij}\mathrm{}_{ij}^{(1)}\mathrm{}_{ij}^{(2)}`$.

Unfortunately the zero-temperature entropy of the D model is too large in order to list all the ground state structures for values of $`L`$ large enough to be interesting. Because of that we have added to the D model a further constraint, suggested by observations on biological RNA, which strongly reduces the entropy while keeping the $`P(q)`$ very similar in shape to the one of the unconstrained model. We avoid structures where a single base pairs is surrounded by non-paired bases, that is a structure with $`\mathrm{}_{i1,j+1}=0`$, $`\mathrm{}_{i,j}=1`$ and $`\mathrm{}_{i+1,j1}=0`$ is forbidden.

The overlap distribution, averaged over $`1000`$ samples, is shown in figure 1. Is worth noticing that the tail for small values of $`q`$ is disappearing very slowly with increasing system size. The peak location as well as the mean overlap converge somewhere around $`q0.87`$, while the variance decreases with the system size approximately as $`\sigma ^2L^{0.45}`$. So the $`P(q)`$ seems to converge, in the thermodynamic limit, to a delta function centered on a value of $`q`$ smaller than $`1`$. Such a value is compatible with the observation that in a typical ground state structure the paired bases are a finite fraction (smaller than $`1`$) of all bases. Nevertheless, as already explained in NOSTRO\_REPLY ; HRT , the triviality of the $`P(q)`$ at zero temperature does not imply a trivial behavior of the whole low-temperature phase, and so we resort to the study of low-energy density of states.

We have calculated the ground state of the Hamiltonian $`^{}`$ for $`18`$ values of $`\epsilon [0.001,131.072]`$ (equally spaced on a logarithmic scale), and many $`L`$ values. We have analyzed a minimum of $`500`$ disorder realizations for the largest chain ($`L=4096`$), and a maximum of circa $`510^4`$ samples for the smallest one ($`L=128`$).

The first GSS $`\mathbf{}_0`$ is chosen with uniform probability in $`𝒢`$ (the set of all the degenerate ground states of the D model). When we switch on the perturbation the new GSS will be the one in $`𝒢`$ having the smallest overlap with $`\mathbf{}_0`$. Nothing else will change as long as $`\epsilon 1`$, that is as long as $`\epsilon `$ is not large enough to make an excited state with $`\mathrm{\Delta }E=1`$ to become the new GSS. This explains the plateau for $`\epsilon 1`$ in figure 2, where we show the average distance between $`\mathbf{}_0`$ and $`\mathbf{}_\epsilon `$ as a function of $`\epsilon `$. The main information we get from figure 2 is the value of the plateau distance $`d0.38`$, corresponding to an overlap $`q=1d0.62`$. This distance can be viewed as the radius of a sphere containing the set $`𝒢`$ of all the GSS. Note that the GSS are not uniformly distributed in this sphere \[otherwise the $`P(q)`$ would be peaked on a much smaller overlap value\], but they are very dense in the central region and very sparse on the boundaries. This means that if one chooses two GSS at random they will typically be very close in the dense region, giving a value of $`q0.87`$, but if one forces the two GSS to be as far as possible the resulting minimum overlap $`q_{\text{min}}`$ will be much smaller, and will depend strongly on the specific disorder realization \[see figure 3, where we plot its probability distribution $`P(q_{\text{min}})`$\].

As it is made clear by the results shown in figure 2, for models with high degeneracy the perturbing method does not work properly in the interesting region of small $`\epsilon `$, giving only information about the minimal overlap among GSS. Now we remove the degeneracy and analyze the other models.

### IV.2 The QD Model

We have defined the QD model such as to keep as much as possible of the degenerate D model, even after removal of the accidental degeneracy. Here we can still distinguish two different regimes (see figure 4, where we plot the average distance $`d(\epsilon ,L)`$ as a function of $`\epsilon `$ for different chain lengths $`L`$): For $`\epsilon >1`$ the data coincide with those for the (unconstrained) D model, while for $`\epsilon 1`$ they have now a non trivial behavior.

In the interesting region of small $`\epsilon `$ the scaling of the data is very subtle and good results can be obtained either with $`\theta =0`$ or with $`\theta >0`$. Our finite size scaling analysis does not allow us to reach a quantitative estimate, and we cannot distinguish in a statistically significant way among a power law (droplet like) scaling and a logarithmic scaling. Further and longer studies are needed to understand better this model.

Despite the difficulties in the data analysis, we believe that the QD model has a large interest and relevance. Indeed it has the great advantage of a single non-degenerate ground state, but still the perturbation added in order to remove the original degeneracy modifies the energies of the structures by a quantity of order $`\eta _00.1`$, thus keeping a large amount of information about the original energy landscape of the degenerate D model of RNA.

### IV.3 The ND Model

In the ND model the ground state degeneracy has been removed by a random term which strongly reshuffles the energy levels.

In figure 5 we show the average distance $`d(\epsilon ,L)`$ for some values of $`\epsilon `$ and $`L`$ (the error is of the order of the symbol size). Data are now smooth functions of $`\epsilon `$, with no singular point at $`\epsilon 1`$, and a finite size scaling analysis can be performed in an easier way.

Following Eq. (4) we have rescaled the data plotting them versus $`\epsilon L^\theta `$. The results are shown in figure 6, where we have included all data points. The best collapse is achieved when using $`\theta 0.33`$. The dotted line has unitary slope and clearly shows that $`d(\epsilon ,L)\epsilon `$ for small $`\epsilon `$ and any fixed $`L`$. We notice that we are looking for a finite size scaling that works well only up to a given value of $`\epsilon L^\theta `$: The method we are using is based on the idea of having a “small” perturbation that acts as a probe. Very large values of $`\epsilon `$ drastically change the Hamiltonian and the energy landscape. It is interesting to note that the scaling seems to work well up to $`d0.4`$, that is the radius of the ball of degenerate ground states in the original D model.

The physical interpretation of this result is the following. For any fixed $`\epsilon `$ the scaling variable $`\epsilon L^\theta `$ vanishes in the thermodynamical limit, which implies that the unperturbed ground state is stable against this kind of perturbation. In order to have a different and more complex behavior one should perturb the original Hamiltonian by a term whose amplitude increases with the size at least as fast as $`L^\theta `$. The energy scale $`L^\theta `$ can thus be interpreted as the energy cost for reaching the first excited state.

A still clearer picture of this phenomenon is given in figure 7, where we plot the average energy difference $`\mathrm{\Delta }E(\epsilon ,L)`$ as a function of the average distance $`d(\epsilon ,L)`$ between the unperturbed and the perturbed ground states. It is evident that, for any fixed distance, the energy difference is growing with the system size, according to the argument given above. We have rescaled the data following equation (5) and the results are shown in figure 8. Again the best collapse is achieved for $`\theta 0.33`$, in perfect agreement with the previous analysis. Also here the scaling region extends over distances up to $`d0.5`$. For small distances the average energy increases as $`\mathrm{\Delta }Ed^2`$ (see the dotted line in figure 8). The same behavior has been observed also in the QD model and we will present a simple explanation in the next subsection, dedicated to the GD model.

Since the value of the $`\theta `$ exponent is small we also tried to fit the data under the assumption $`\mathrm{\Delta }E_{\text{typ}}(\mathrm{log}L)^a`$, that is $`\theta =0`$. This behavior is suggested by mean-field solutions of disordered models MPV  and more particularly from previous findings on similar models NOSTRO\_PRL ; BUND\_HWA\_1 . The conclusion is that a logarithmic fit with $`a=1.85\pm 0.15`$ still works rather well, but definitely the power law fit with $`\theta =0.33\pm 0.01`$ is more accurate and always has a smaller $`\chi ^2`$ value.

### IV.4 The GD Model

Under the application of the $`\epsilon `$-coupling perturbation, the GD model behaves very similarly to the ND model in that:

* All the data perfectly collapse with $`\theta 1/3`$;
* For small $`\epsilon `$, $`d\epsilon `$;
* For small $`d`$, $`\mathrm{\Delta }Ed^2`$.

Moreover the results listed here do not depend on the presence of the constraint $`\mathrm{}_{ij}=0`$ for $`|ij|<4`$ (see section II). These findings imply that we have at hand a very simplified model which share the same phenomenology with more realistic models and which is more amenable to an analytical treatment. For sake of clarity we recall its definition: We have an Hamiltonian of the form given in equation (1), where the pairing energies $`e_{ij}=e_{ji}`$ are $`L(L1)/2`$ independent Gaussian variables with zero mean and unitary variance, and the $`\mathrm{}_{ij}`$ satisfy the planarity condition.

Within this model it is easier, for example, to understand the behavior $`\mathrm{\Delta }Ed^2`$ for small $`d`$. First of all we observe from numerical simulations that in a typical GSS $`\mathbf{}_0`$ the fraction of paired bases is $`f<1`$ and the distribution of the pairing energies $`e_{ij}`$ of the active links, the ones with $`\mathrm{}_{ij}^{(0)}=1`$, can be very well approximated by a Gaussian of negative mean and finite width (the distribution is truncated since positive pairing energies are forbidden in the GSS). Let us call the distribution of the pairing energies absolute values $`P_e(e)`$. The only property we need for the proof is a finite weight in zero, $`P_e(0)>0`$, and this is the case for the GD model (and also for the ND model).

Now we construct a sequence of structures $`\mathbf{}_k`$ such that $`\mathrm{\Delta }Ed^2`$ for small $`d`$. $`\mathbf{}_k`$ is obtained from $`\mathbf{}_0`$ removing the $`k`$ weakest links, i.e. those with the smallest (in absolute value) pairing energies. So the distance between $`\mathbf{}_0`$ and $`\mathbf{}_k`$ is $`d=k/L`$ and the energy difference $`\mathrm{\Delta }E`$ is the sum of the smallest (in absolute value) $`k`$ pairing energies. For large $`L`$ we can write

$$\mathrm{\Delta }E=_0^uP_e(e^{})e^{}\text{d}e^{}=P_e(0)\frac{u^2}{2},$$
(6)
where the last equality only holds for small $`u`$. The upper integration limit $`u`$ is chosen such that $`k`$ pairing energies, or equivalently a fraction $`2k/(fL)`$ of pairing energies, are summed, that is

$$\frac{2k}{fL}=\frac{2d}{f}=_0^uP_e(e^{})\text{d}e^{}=P_e(0)u,$$
(7)
where the last equation is valid for small $`u`$. Combining the above equations we obtain $`\mathrm{\Delta }Ed^2`$.

Since we have chosen a sequence of structures which are not guaranteed to have the lowest possible energies, we can only argue that $`\mathrm{\Delta }Ed^\alpha `$ with $`\alpha 2`$. Nevertheless from numerical simulations the exponent turns out to be exactly $`2`$ (see figure 8).

In this very simplified GD model we can make one more analytical prediction, regarding the fraction of paired bases in the GSS. The number of planar structures with a fraction $`f`$ of paired bases can be easily calculated with the help of generating functions and turns out to be given by $`\mathrm{exp}[Ls(f)]`$, with an intensive entropy

$$s(f)=f\mathrm{log}f(1f)\mathrm{log}(1f)+f\mathrm{log}2.$$
(8)
$`s(f)`$ has a maximum for $`f=\frac{2}{3}`$, with $`s(\frac{2}{3})=\mathrm{log}3`$.

Let us now fix $`f`$ and see how the energies of the $`\mathrm{exp}[Ls(f)]`$ structures are distributed. They look random, but actually, since the independent Gaussian random variables are only $`L(L1)/2`$, there must be many correlations among them. Since any of these energies is the sum of $`fL/2`$ random Gaussian pairing energies, we make the approximation that the distribution of structure energies is also Gaussian with a variance proportional to $`fL`$, i.e. $`𝒫(E)\mathrm{exp}[E^2/(bfL)]`$. The evaluation of the coefficient $`b`$ is out of our present scopes. Given, for any fixed $`f`$, the number of structures and the distribution of the energies, we can estimate the most probable lowest energy $`E_{\text{min}}(f)`$ through

$$e^{Ls(f)}=_{\mathrm{}}^{E_{\text{min}}}𝒫(E^{})\text{d}E^{}e^{E_{\text{min}}^2/(bfL)},$$
(9)
where the last equality holds because $`E_{\text{min}}`$ is negative and large. The above equation implies $`E_{\text{min}}(f)=L\sqrt{bfs(f)}`$. In order to find the fraction of paired bases corresponding to the GSS one has to minimize $`E_{\text{min}}(f)`$, or equivalently maximize $`fs(f)`$, over $`f`$. Such an extremum is achieved for $`f=0.86`$ to be compared with the fraction of paired bases found numerically $`f=0.856`$. The rather small discrepancy tells us that the approximation made on the form of the structure energies distribution $`𝒫(E)`$ is not so bad.

The apparent correctness of such a simple approximation could suggest that the GD model has a trivial energy landscape. We have checked for this possibility with the following method: In a trivial energy landscape any reasonably smart greedy algorithm should be able to reach, or at least to closely approach, the ground state energy. We have used a greedy algorithm which builds up the structure in the following way: It starts with a structure with no links, at each step it chooses the lowest negative pairing energy (largest in absolute value) among the set of those allowed by the planarity condition, and adds the corresponding link to the growing structure. Using this greedy algorithm we can reach the energies shown in figure 9 which are more than 10% higher than the corresponding ground state energies. So, it seems that finding a structure with low energy in linear time is not an easy task. This suggests a complex energy landscape. A deeper analysis is obviously needed in order to say how much complex the energy landscape of the GD model is.

## V Conclusions

Our results allow to describe a clear and simple physical picture for the RNA-inspired models studied here. All of them possess a glassy phase at low enough temperatures (since we have analyzed very low energy density of states we cannot make predictions on the location of the critical temperature). We can claim that our study clearly selects a positive $`\theta `$ exponent value for all cases but for the quasi degenerate QD model, where our numerical results are not precise enough to allow us quantitative statements.

At variance with the results of BUND\_HWA\_1 ; BUND\_HWA\_2  we find that in the non-degenerate ND model the broken phase does not look marginal, but a standard droplet glassy phase with $`\theta 1/3>0`$. Our way of analyzing the data allows us to exclude (with good confidence) a simple logarithmic divergence of the energy difference between ground state and excited states. On this issue we agree with the results of KMM\_RNA : The difference in the estimate we give for $`\theta `$, as compared to the $`\theta 0.23`$ of KMM\_RNA , is probably due to the different fitting procedure, and to the number of free parameters used in the fit.

The $`\theta `$ exponent we find is perfectly compatible with that for directed polymers in random media in 1+1 dimensions MEZEMEZ , $`\theta _{DPRM}=\frac{1}{3}`$. Since the two models have some similarities this relation could indeed hide a deep connection.

The degenerate D model and the quasi degenerate QD model we have defined above are maybe the less trivial and the most intriguing from the theoretical point of view. Unfortunately we were not able to determine accurately the asymptotic scaling behavior in the latter.

It is probable, on the contrary, that the most part of the analytic developments will be obtained for the Gaussian disorder GD model, that is by far the simplest among all the models with a non-trivial behavior.

###### Acknowledgements.

We thank T. Hwa, F. Krzakala, M. Mézard, M. Müller and G. Parisi for a number of interesting conversations.


text:
# Fraction of uninfected walkers in the one-dimensional Potts model

## I Introduction

When a system is quenched from a homogeneous high-temperature disordered state into a low-temperature regime, well defined domains of equilibrium ordered phases form randomly and then grow with time (“coarsen”) in a self-similar way, until the domain size becomes comparable to the size of the system . A dynamic scaling hypothesis suggests that, at late times, the evolution of the system is characterized by a single length scale $`L(t)`$ that represents a typical linear size of the domains. It is well established, at least for systems with a scalar order parameter, that $`L(t)t^n`$ with $`n=1/2`$ for non-conserved dynamics and $`n=1/3`$ for conserved dynamics . The Ising spin model evolving with Glauber spin flip dynamics  demonstrates behavior characteristic of the former while evolution according to Kawasaki spin dynamics  exemplifies the latter. A generalization of the Glauber-Ising model is the $`q`$-state Potts model , which has $`q`$ distinct, but equivalent, ordered phases. Experimental realizations are known for $`q=2`$ i.e. the Ising model and for $`q=3,4,\mathrm{}`$ . The $`q=\mathrm{}`$ case describes several cellular structures , e.g. polycrystals , soap froths , magnetic bubbles  and foams .

Such coarsening systems are amongst those which have received considerable attention in recent years with regard to their persistence properties . The persistence probability $`P(t)`$ is simply the probability that a given stochastic variable $`\varphi (x,t)`$, with zero mean, does not change sign in the time interval $`[0,t]`$. Theoretical and computational studies of persistence include spin systems in one  and higher  dimensions, diffusion fields , fluctuating interfaces , phase-ordering dynamics , Lotka-Volterra models  and reaction-diffusion systems . Experimental studies include the coarsening dynamics of breath figures , foams , soap froths , twisted nematic liquid crystals , and one-dimensional gas diffusion . Persistence in non-equilibrium critical phenomena has also been studied in the context of the global order parameter $`M(t)`$ (e.g. the total magnetization of a ferromagnet), regarded as a stochastic process . In many systems of physical interest, the persistence decays algebraically according to $`P(t)t^\theta `$, where $`\theta `$ is the persistence exponent and is, in general, nontrivial. The non-triviality of $`\theta `$ emerges as a consequence of the coupling of the field $`\varphi (x,t)`$ to its neighbors, since such coupling implies that the stochastic process at a fixed point in space and time is non-Markovian. For the $`q`$-state Potts model at $`T=0`$, the fraction of spins which have never flipped up to time $`t`$ (i.e that fraction that have persisted in their original phase), has been observed  to obey a power law decay $`P(t)t^{\theta (q)}`$, where $`\theta (q)`$ has been obtained exactly, in one-dimension (1D), by Derrida et al.  for arbitrary $`q`$,

$$\theta (q)=\frac{1}{8}+\frac{2}{\pi ^2}\left[\mathrm{cos}^1\left(\frac{2q}{\sqrt{2}q}\right)\right]^2.$$
(1)
There is a very direct way  of relating the Glauber dynamics of the 1D $`q`$-state Potts model at zero temperature to reaction-diffusion models . Consider uncorrelated initial conditions where each of the $`q`$ phases is present with equal density $`1/q`$. During any time interval $`dt`$, each spin $`S_i(t)`$ has a probability $`dt/2`$ of becoming equal to its right neighbor, $`dt/2`$ of becoming equal to its left neighbor and a probability $`1dt`$ of retaining its own value. The domain walls therefore perform random walks. Upon contact, two domain walls react instantaneously, either annihilating $`A+A\mathrm{}`$ with probability $`1/(q1)`$ or coalescing $`A+AA`$ with probability $`(q2)/(q1)`$, these numbers being the probabilities that the states on the far sides of the walkers are the same (annihilation) or different (coalescence). This reaction-diffusion model, together with the fact that in the initial condition each bond is occupied by a domain wall with probability $`(q1)/q`$, is completely equivalent to the spin problem. The fraction of spins which have never flipped up to time $`t`$ is then simply the fraction of sites which have never been visited by a random walker. The probabilistic algorithm for implementing the Potts model through the annihilation or coalescence of random walkers allows $`q`$ to be a real number but is restricted to $`q2`$. However, an equivalent Ising spin representation of the Potts persistence problem  permits a study  of $`\theta (q)`$ in the regime $`1<q<2`$.

In the present work we study the fraction of random walkers which have never encountered another walker, using the reaction-diffusion representation of the 1D Glauber-Potts model at zero temperature. To facilitate our discussion, we consider all the random walkers at time $`t=0`$ to be uninfected. All the walkers are, however, considered to be infectious, so that for $`t>0`$ any contact between walkers (assuming they survive the encounter) leads to mutual contamination. Walkers that avoid all contact remain uninfected. Our goal, therefore, is to address the fraction of uninfected walkers $`U(t)`$ up to time $`t`$.

This problem has been addressed in some detail by a number of different authors , although under the guise of a different interpretation. The general approach used involves monitoring the motion of a tagged test particle, released in the system at time $`t=0`$. If the tagged particle diffuses with diffusion constant $`D^{}`$ and the other particles in the system diffuse with diffusion constant $`D`$, the tagged particle is viewed as an external impurity diffusing through a homogeneous background. One is then interested in addressing the survival probability of the tagged particle i.e., the probability is has not been absorbed by another particle in the system. Of course in our terminology, this corresponds to the probability that the tagged particle remains uninfected. It is well established  that $`U(t)`$ decays according to a power law $`U(t)t^{\varphi (q,c)}`$, where $`c=D^{}/D`$ is the ratio of the two diffusion constants in the problem. When $`D^{}=0`$, corresponding to a static impurity, the model clearly maps onto the problem of persistent spins and in this case, $`\varphi (q,0)=\theta (q)`$. The exponent $`\varphi (q,c)`$ is therefore considered a generalization of the standard persistence exponent $`\theta (q)`$, characterizing the survival probability of a mobile particle under the coarsening dynamics of the Potts model. A recent study of the survival probability of a mobile particle, moving according to either deterministic or stochastic rules, in given in  , where the external fluctuating field takes the form of either the solution to the diffusion equation, the coarsening dynamics of the $`q`$-state Potts model or spatially uncorrelated Brownian signals.

In the limit $`q\mathrm{}`$ the Potts model is equivalent to the $`A+AA`$ single species coalescence process. The dynamics of this system are particularly simple, since each particle sees only the “cage” formed by its two nearest neighbors. It is immaterial that these neighbors may coalesce with more distant particles. Consequently, the kinetics of each particle involves only itself and its two nearest neighbors and one is therefore able to calculate $`\varphi (\mathrm{},c)`$ exactly :

$$\varphi (\mathrm{},c)=\frac{\pi }{2\mathrm{cos}^1\left[\frac{c}{1+c}\right]}.$$
(2)
(Note that for $`c=0`$, this of course reduces to the static persistence exponent, $`\varphi (\mathrm{},0)=\theta (\mathrm{})=1`$.) In the limit of the Ising model ($`q=2`$), there is, however, no corresponding exact result for general $`c`$. An exact result is only known for $`c=0`$ i.e., $`\varphi (2,0)=\theta (2)=3/8`$. A mean-field Smoluchowski  approach predicts for $`q=2`$ and general $`c`$ that, $`\varphi (2,c)=\sqrt{(1+c)/8}`$, in good agreement with numerical simulations . However, when extended to the large $`q`$ limit, this Smoluchowski approach differs substantially from the exact $`q\mathrm{}`$ result Eq. (2). More recently, however, a perturbation theory developed by Monthus  has evaluated the exponent $`\varphi (q,c)`$ to first order in $`q1`$ for arbitrary $`c`$ and at first order in $`c`$ for arbitrary $`q`$. In this paper, we will not study the exponent $`\varphi (q,c)`$ in complete generality since we are primarily interested in the special case where all of the particles are equally mobile $`(D^{}=D)`$ and henceforth we therefore write $`\varphi (q,1)=\varphi (q)`$.

Within the context of the current work, where we are interested in the case $`c=1`$, the random walkers represent the motion of domain walls undergoing zero temperature coarsening. The fraction of uncollided domain walls (i.e., uninfected walkers) has previously been interpreted as “domain wall persistence” or equally, the probability that two adjacent domains survive. Under this interpretation, Majumdar et al.  have studied the quantity $`U(t)`$ numerically for various $`q`$. They observed a power law decay, $`U(t)t^{\varphi (q)}`$ in agreement with the results of Monthus  and obtained a spectrum of values for $`\varphi (q)`$: $`1/2\varphi (q)3/2`$, where $`2q\mathrm{}`$. In this paper, we extend and improve these numerical results by presenting (Sec. II) more accurate data than previously achieved  and exploring a larger range of $`q`$. We also relate the value of $`\varphi (q)`$ to the static persistent exponent $`\theta (q)`$ via $`\varphi (q)=\gamma (q)\theta (q)`$, where $`\gamma `$ is a new exponent that we recently introduced in .

We additionally consider (Sec. III) the coupled diffusion-limited reaction $`A+AB`$, $`B+BA`$ in one-dimension with equal initial densities, $`\rho _A(0)=\rho _B(0)`$. We find that the density of particles decays in this model according to $`\rho t^{1/2}`$, independent of the initial walker density, characteristic of “Potts like” behavior . The fraction of sites unvisited by either an $`A`$ or a $`B`$ particle decays as $`P(t)t^\theta `$ with $`\theta 1.33`$. We discuss these exponents within the context of the $`q`$-state Potts model and observe that the fraction of walkers which remain uninfected by their own species decays according to $`U(t)t^\varphi `$, where $`\varphi 1.13`$. We also study the fraction of walkers which remain uninfected by both species and find that this too obey a power law, with $`\varphi 1.93`$. In analogy to Sec. II we express the exponent $`\varphi `$ is terms of the static persistence exponent $`\theta `$ and similarly obtain $`\varphi =\gamma \theta `$. We conclude in Sec. IV with a summary of the results and a discussion of some open questions.

## II THE POTTS MODEL

We investigate numerically the fraction of uninfected walkers $`U(t)`$ in the 1D $`q`$-state Potts model, at zero temperature, as a function of $`q`$. Our simulations are performed on a 1D lattice of size $`10^7`$ with periodic boundary conditions. At $`t=0`$ a random walker is placed at each lattice site, so that $`\rho (0)=1`$. We choose this initial high density of walkers to accelerate the system’s evolution into the asymptotic regime. Our model is updated using the direct method  i.e. at each computational step, a particle is picked at random and moved with probability $`D=1/2`$ to a neighboring site, where $`D`$ is the diffusion constant. If the destination site is occupied, the two particles either annihilate, with probability $`1/(q1)`$, or coalesce, to become a single particle, with probability $`(q2)/(q1)`$, in accordance with the reaction-diffusion dynamics of the Potts model detailed in Sec. I. For each move of a particle, time $`t`$ is increased by $`dt=1/N`$, where $`N`$ is the current number of particles in the system. Such sequential dynamics can be chosen without loss of generality as parallel dynamics exhibit similar asymptotic behavior . Our simulations are performed for $`5000`$ time steps and our results are averaged over $`100`$ runs.

Fig. 1 clearly shows that the fraction of uninfected walkers decays according to a power-law, $`U(t)t^{\varphi (q)}`$. In Table I we present our numerical values for $`\varphi (q)`$ in the range $`2q\mathrm{}`$. These values were obtained by performing a linear regression on log-log curves such as those presented in Fig. 1. The regression was taken in the range $`10t1000`$. This range was chosen to avoid initial transients and to avoid statistical fluctuations between different runs which become prominent as $`t`$ increases, especially for large $`q`$. In Fig. 2 we plot $`\varphi `$ as a function of $`q`$, where we choose to represent the data on a linear-log plot merely for clarity of presentation.

| $`q`$ | $`\varphi `$ |
| --- | --- |
| 2 | 0.5006(7) |
| 4 | 0.867(5) |
| 8 | 1.121(5) |
| 16 | 1.283(6) |
| 25 | 1.352(6) |
| 32 | 1.380(7) |
| 50 | 1.419(7) |
| 64 | 1.435(8) |
| 100 | 1.461(8) |
| 128 | 1.464(8) |
| 256 | 1.481(9) |
| $`\mathrm{}`$ | 1.495(9) |

TABLE I. Numerical values of the exponent $`\varphi (q)`$ in the 1D, $`T=0`$, Glauber-Potts model where $`U(t)t^{\varphi (q)}`$, for various values of q.

The case $`q=2`$, which corresponds to the Ising model, reduces to the $`A+A\mathrm{}`$ reaction-diffusion process, where the density of uninfected walkers is clearly equal to the particle density. It is well known  that in the $`q`$-state Potts model the particle density obeys a power-law decay, $`\rho (t)t^{1/2}`$, which is consistent with our value of $`\varphi (2)=0.5006(7)`$. For the case $`q=\mathrm{}`$, one can obtain the exact value of $`\varphi (\mathrm{})`$ from Eq. (2) i.e. $`\varphi (\mathrm{},1)=3/2`$. This is again consistent with our numerical result, $`\varphi (\mathrm{})=1.495(9)`$. These limiting values, along with our numerics, suggest that the fraction of uninfected walkers decays according to $`U(t)t^{\varphi (q)}`$ where $`1/2\varphi (q)3/2`$, for $`2q\mathrm{}`$, consistent with the results of Monthus  and Majumdar et al. . Although the values of $`\varphi (q)=1/2`$ and $`\varphi (\mathrm{})=3/2`$ at the limit of the $`q`$ spectrum are understood, it was pointed out by Monthus  that the formalism that led to the exact determination of the static persistence exponent $`\theta (q)`$  cannot be easily applied to compute the exact value of $`\varphi (q,c)`$ for general $`c`$, including the limit $`c=1`$. The exact calculation of $`\varphi (q)`$ therefore remains an open problem.

Recently, we have explored the connection between uninfected walkers and unvisited sites in one-dimension within the context of a system of non-interacting randomly diffusing particles $`A\mathrm{}\mathrm{}A`$, and the $`A+B\mathrm{}`$ diffusion-limited reaction . We found that in both of these models,

$$U(t)[P(t)]^\gamma ,\mathrm{with}\gamma 1.39.$$
(3)
While this relationship is obeyed reasonably accurately over the range of times accessible to our simulations, there is evidence, for both processes, that $`\gamma `$ approaches a smaller value at late times . Clearly, we can write an analogous relationship for the $`q`$-state Potts model, where $`U(t)t^{\varphi (q)}`$, $`P(t)t^{\theta (q)}`$ and $`\gamma `$ is now some function of $`q`$, i.e. $`t^{\varphi (q)}=t^{\gamma (q)\theta (q)}`$, so that

$$\varphi (q)=\gamma (q)\theta (q).$$
(4)
We plot $`\varphi (q)`$ against $`\theta (q)`$ in Fig. 3, where our values for $`\varphi (q)`$ have been taken from Table I and the corresponding values of $`\theta (q)`$ have been calculated exactly from Eq. (1). The point at the origin corresponds to $`q=1`$, where $`\varphi `$ and $`\theta `$ both vanish. Fig. 3 indicates that $`\varphi (q)`$ increases monotonically as a function of $`\theta (q)`$. Therefore, given that we know the exact values of $`\theta (q)`$ and $`\varphi (q)`$ for $`q=2,\mathrm{}`$, and given the convex nature of the plot in Fig. 3, we can use Eq. (4) to determine precise bounds on $`\gamma (q)`$ in this range:

$$4/3\gamma (q)3/2,q[2,\mathrm{}].$$
(5)
Thus while $`\varphi (q)`$ and $`\theta (q)`$ are separately quite strongly dependent on $`q`$, their ratio $`\gamma (q)`$ depends only weakly on $`q`$.

## III THE EXCHANGE MODEL

The $`A+AB`$, $`B+BA`$ diffusion-limited reaction offers a curious combination of birth and death processes. For brevity, we dub this model the “exchange model”. The particles in the exchange model diffuse randomly. Whenever a particle encounters a walker belonging to its own species, the two particles annihilate instantaneously, to create a particle belonging to the other species at the point of annihilation. For the purposes of this study, we allow lattice sites to contain multiple particles, enabling two particles of differing species to cohabit a single lattice site. This feature of our model ensures that, in one-dimension, the lattice structure $`\mathrm{}ABABABABABAB\mathrm{}`$ continues to evolve. Confining our system to a maximum of one particle per site would otherwise mean that the preceding configuration would be stable for all time $`t`$. We focus our attention on the 1D exchange model with equal initial densities of $`A`$ and $`B`$ particles.

Our numerical simulations are performed on a one-dimensional lattice with periodic boundary conditions. At $`t=0`$ equal numbers of $`A`$ and $`B`$ particles are randomly distributed on the lattice with a maximum of one particle per site. Our system then evolves according to the dynamics described above. The system is updated in the same manner as for the $`q`$-state Potts model described in Sec. II.

### A Particle Density

The first question we address is the particle decay. The unimolecular reactions defining the exchange model are characteristic of the $`q`$-state Potts model and therefore we can expect the particle density to exhibit similar asymptotic behavior. Indeed, for every pair of reactions between two $`A`$ and two $`B`$ particles, a single $`A`$ and a single $`B`$ are produced, which is the same result as for two $`q=\mathrm{}`$ Potts models, $`A+AA`$, $`B+BB`$ operating on a single lattice. We therefore expect the particle density to decay in the exchange model according to $`\rho _A(t)Ct^\alpha `$, with $`\alpha =1/2`$, independent of the initial particle density . This is indeed confirmed by our simulations. In Fig. 4 we present our results for the particle density for three different initial values of $`\rho _A(0)`$ (recall $`\rho _B(0)=\rho _A(0)`$). The simulations were performed on a lattice of size $`10^6`$ for $`10^5`$ time steps and the results averaged over 100 runs.

In Table II we present our numerical values for the particle density decay exponent $`\alpha `$ and amplitude $`C`$, where $`\rho _A(t)Ct^\alpha `$. These values were obtained by performing a linear regression on the the plots in Fig. 4 in the asymptotic regime, $`10^4t10^5`$.

| $`\rho _A(0)`$ | $`\alpha `$ | $`C`$ |
| --- | --- | --- |
| 0.5 | 0.5007(9) | 0.419(5) |
| 0.25 | 0.5007(9) | 0.419(5) |
| 0.1 | 0.5000(5) | 0.416(9) |

TABLE II. Numerical values for the particle density decay exponent $`\alpha `$ and the amplitude $`C`$ in the 1D exchange model, where $`\rho _A(t)Ct^\alpha `$, for different initial densities $`\rho _A(0)`$.

For the 1D $`q`$-state Potts model, the amplitude $`C`$ is $`q`$-dependent ,

$$C(q)=\frac{q1}{q}\frac{1}{\sqrt{2\pi D}}.$$
(6)
Is there perhaps a value of $`q`$ which corresponds to the exchange model? A value of $`q=4`$ in Eq. (6), with $`D=1/2`$ returns a value $`C(4)=0.423`$, in good agreement with the values of $`C`$ recorded in Table II. Notice that this value lies exactly midway between the minimum value, $`C(2)=0.282`$ and the maximum value, $`C(\mathrm{})=0.564`$, set by the $`q`$-state Potts model. Our numerical results therefore suggest, that in terms of the particle density, the two reactions which define the exchange model ($`A+AB,B+BA`$) are equivalent, in some loose sense, to two independent $`q=4`$ Potts models, though we would not wish to stretch this comparison too far.

### B Fraction of Unvisited Sites

We now turn our attention to the persistence properties of the exchange model. Here, we are interested in the fraction of sites, $`P(t)`$, that have never been visited by either an $`A`$ or a $`B`$ particle up to time $`t`$. An obvious starting point is to see if $`P(t)`$ obeys a power-law decay with an exponent equal to that returned in the case of two $`q=4`$ Potts Models i.e. $`P(t)t^{2\times \theta (q=4)}`$ where, using Eq. (1), $`2\times \theta (q=4)=1.263`$. Our numerical simulations do indeed indicate a power-law decay, $`P(t)t^{\theta _{ex}}`$. However, the decay exponent $`\theta _{ex}1.33`$, in disagreement with what one might have naively expected from a study of the particle decay.

Our simulations are performed on a lattice of size $`10^6`$ for 5000 time steps and the results averaged over 200 runs. In Fig. 5 we present a log-log plot of $`P(t)`$. In Table III we show the results of a linear regression performed on the curves in Fig. 5, to obtain values of $`\theta _{ex}`$. The regression was performed in the regime $`500t5000`$, to avoid initial transients.

| $`\rho _A(0)`$ | $`\theta _{ex}`$ |
| --- | --- |
| 0.2 | 1.330(7) |
| 0.3 | 1.333(5) |

TABLE III. Numerical values for the persistence decay exponent $`\theta _{ex}`$ in the 1D exchange model, where $`P(t)t^{\theta _{ex}}`$, for initial densities $`\rho _A(0)`$.

### C Fraction of Uninfected Walkers

Finally, we consider the fraction of uninfected walkers in the exchange model. We focus, without loss of generality, on just the fraction of uninfected $`A`$ particles. In this model, we are only interested in original $`A`$ particles. Those $`A`$ particles which have been created through the reaction $`B+BA`$ after $`t=0`$, are also therefore considered to be infected. There are two distinct cases to consider: (i) the fraction of $`A`$ particles which remain uninfected by only their own species and (ii) the fraction of $`A`$ particles which have avoided infection by both types of species. Our approach is wholly numerical. We address the former case first.

Fig. 6 indicates that the fraction of uninfected $`A`$ walkers obeys a power-law decay according to $`U(t)t^{\varphi _A}`$, where the subscript $`A`$ on the exponent denotes infection by only $`A`$ particles (case (i)). The results for $`\rho _A(0)=1`$ were generated from 500 runs on a lattice of size $`5\times 10^5`$ over 5000 time steps. The results for $`\rho _A(0)=0.5`$ were generated using a lattice of size $`10^6`$, run for $`10^4`$ time steps and also averaged over 500 runs. In Table IV, we present our numerical values of $`\varphi _A`$ for two different initial starting densities. The values of the exponent were obtained by performing a linear regression on the curves in Fig. 6 in the regime $`100t5000`$.

| $`\rho _A(0)`$ | $`\varphi _A`$ |
| --- | --- |
| 1.0 | 1.133(5) |
| 0.5 | 1.138(8) |

TABLE IV. Numerical values for the exponent $`\varphi _A`$ in the 1D exchange model, where $`U(t)t^{\varphi _A}`$, for initial densities $`\rho _A(0)`$.

We now address case (ii), the fraction of $`A`$ particles which have avoided infection by both $`A`$ and $`B`$ walkers. Our simulations are performed on a lattice of size $`10^6`$ for $`2000`$ time steps and averaged over $`500`$ runs. In line with our previous data, we present our results for $`U(t)`$ on a log-log plot in Fig. 7.

Unfortunately, we do not obtain straight line graphs. However, an inspection of individual runs of the data (Fig. 8) shows that, in this model, statistical fluctuations between different runs become prominent very quickly, whereas earlier times are governed by transient effects. The dominance of these two regimes renders any determination of an exponent particularly difficult. Using a lower initial density, with the intention of delaying the onset of noisy data only extends the transient regime, whereas using a higher initial density to minimize transient effects accelerates the onset of large statistical fluctuations. In view of the power-law behavior we identified in case (i) and the general trend of our data, we suggest that, with sufficiently good statistics, the fraction of $`A`$ walkers uninfected by either $`A`$ or $`B`$ particles would also exhibit a power-law decay $`U(t)t^{\varphi _{AB}}`$, where the subscript on the exponent denotes infection is permitted by both species. There is, however, a cleaner method of extracting the exponent $`\varphi _{AB}`$. We suggest, in analogy with the other models we have studied, that $`U(t)[P(t)]^{\gamma _{ex}}`$ so that,

$$\varphi _{AB}=\gamma _{ex}\theta _{ex}.$$
(7)
If we can establish the validity of the relation $`U(t)[P(t)]^{\gamma _{ex}}`$, and determine the value of $`\gamma _{ex}`$, we can then read off the exponent $`\varphi _{AB}`$ from Eq. (7). Our initial condition specifies a maximum of one particle per site, so $`U(0)=1`$. While, in the continuum description of the model, $`P(0)=1`$ also, on a lattice this becomes $`P(0)=1\rho _0`$. Therefore, to place these two quantities on an equal footing in our numerical study, we define, $`p(t)=P(t)/P(0)`$, so that both $`U(0)`$ and $`p(0)`$ are equal to unity. Note, of course, that in the continuum limit, $`\rho _00`$, we recover $`p(t)P(t)`$. We present in Fig. 9 a log-log plot of $`U(t)\text{ vs. }p(t)`$.

The plot in Fig. 9 is quite linear, but does not go through the origin, suggesting that in fact $`U(t)Kp(t)^{\gamma _{ex}}`$ with $`K1`$. We note that in the numerical analysis of the relationship between noninfectedness and persistence in the $`A\mathrm{}\mathrm{}A`$ and $`A+B\mathrm{}`$ models  the equivalent prefactor $`K`$ was close enough to unity ($`K0.99`$) to tempt the speculation that it might become unity in the continuum (low density) limit, though we stress that in the latter models there was good evidence that some of the simulations had not reached the asymptotic regime. In order to evaluate $`\gamma _{ex}`$ and $`K`$ for the present model, we perform a linear regression on the data in Fig. 9 in the region indicated by the arrows, thereby avoiding initial transients associated with lattice effects, and the onset of statistical fluctuations between different runs of the data. We summarize our results for $`\gamma _{ex}`$ and $`K`$ in Table V.

| $`\rho _A(0)`$ | $`\gamma _{ex}`$ | $`K`$ |
| --- | --- | --- |
| 0.2 | 1.453(7) | 0.863(9) |
| 0.3 | 1.456(7) | 0.754(9) |

TABLE V. Numerical values of $`\gamma _{ex}`$ and $`K`$ in the 1D exchange model, where $`U(t)K[p(t)]^{\gamma _{ex}}`$, for two initial densities $`\rho _A(0)`$.

Note that the value of $`K`$ for $`\rho _A(0)=0.2`$ is larger than that for $`\rho _A(0)=0.3`$, suggesting the possibility that $`K1`$ as $`\rho _A(0)0`$, consistent with our results for other models , though we have no concrete arguments to support this idea. To test our results, we attempt to collapse the data for the uninfected fraction and the persistence on to a single curve by plotting $`\mathrm{ln}U(t)`$ and $`\gamma _{ex}\mathrm{ln}p(t)+\mathrm{ln}K`$ as functions of time $`t`$, in Fig. 10.

The data collapse is excellent over both transient and asymptotic regimes, suggesting that the relation $`U(t)K[p(t)]^{\gamma _{ex}}`$ holds to a high degree of accuracy for all times $`t`$. We attribute any deviations between the two curves at late times to the onset of statistical fluctuations between different runs of the data. Clearly we can improve the collapse at either early or late times by performing the linear regression in Fig. 9 on the appropriate part of the curve. We find that over the whole regime presented in Fig. 10, $`\gamma _{ex}1.45(1)`$ works well. For brevity and clarity of presentation, we have presented our results for $`\rho _A(0)=0.2`$ only, although we have achieved equally good results for $`\rho _A(0)=0.3`$.

The values $`\gamma _{ex}1.45`$ and $`\theta _{ex}1.33`$ yield, using Eq. (7), $`\varphi _{AB}1.45\times 1.331.93`$. We therefore argue that, when infection is caused by both species, $`U(t)t^{\varphi _{AB}}`$ with $`\varphi _{AB}1.93`$.

## IV Discussion and Summary

In this paper, we have investigated the fraction of uninfected walkers $`U(t)`$ in the 1D $`q`$-state Potts model evolving at zero temperature. Our numerical results are consistent with previous work , suggesting that $`U(t)t^{\varphi (q)}`$ with a nontrivial exponent $`\varphi (q)`$, where $`1/2\varphi (q)3/2`$ for $`2q\mathrm{}`$. Although the values of $`\varphi (q)`$ for $`q=2`$ and $`q=\mathrm{}`$ are understood, the exact calculation of $`\varphi (q)`$ for all $`q`$, appears to be a hard problem. In analogy with the other models we have studied , we have reduced the familiar study of persistent sites $`P(t)`$, to a limiting case in the study of uninfected walkers and noted that $`U(t)[P(t)]^{\gamma (q)}`$, with $`4/3\gamma 3/2`$ for $`2q\mathrm{}`$. We discussed in Sec. I that, with regard to the persistence properties of the Potts model, it is possible to probe the region $`1q2`$. Although we have included a point for the state $`q=1`$, where $`\theta =\varphi =0`$, in Fig. 2 and Fig. 3, we leave a study of the fraction of uninfected walkers in this latter regime as an open challenge.

We have also studied the 1D exchange model, defined by the reactions $`A+AB`$, $`B+BA`$. Our numerical simulations indicate similarities between this diffusion-limited system and the Potts class of models, in that the particle density decays according to $`\rho (t)t^{1/2}`$, independent of the initial particle density. The amplitude of the particle decay lies approximately halfway between the minimum and maximum values set by the Potts model suggesting that, within the context of the particle decay, the exchange model roughly mimics the behavior of two independent $`q=4`$ Potts models. We have also shown that the fraction of sites unvisited by either an $`A`$ or a $`B`$ particle decays according to $`P(t)t^{\theta _{ex}}`$ where $`\theta _{ex}1.33`$. Our study of the fraction of uninfected walkers in the exchange model focused on two cases. In the first instance, we addressed the fraction of $`A`$ walkers infected by only their own species and established that, in this case, $`U(t)t^{\varphi _A}`$ where $`\varphi _A1.13`$. When cross species contamination was also included in our simulation, we found $`U(t)t^{\varphi _{AB}}`$ with $`\varphi _{AB}1.93`$. We presented numerical evidence that $`U(t)K[p(t)]^{\gamma _{ex}}`$ with $`\gamma _{ex}1.45`$. An obvious goal for the future is to establish a theoretical framework for the exponents that occur in the exchange model and, in particular, for the observed simple relation between $`U(t)`$ and $`P(t)`$. Across the range of models we have studied, in both this and our previous paper , the fraction of uninfected walkers, $`U(t)`$, and the fraction of unvisited sites, $`P(t)`$, are seemingly related by $`U(t)[P(t)]^\gamma `$. This remarkable relationship, which seems to hold reasonably accurately across the entire time regime accessible to simulation, merits further attention. In Table VI, we summarize the values of $`\gamma `$ for the various models we have studied (note, however, that there is some evidence  that, for the first two models, the asymptotic value of $`\gamma `$ is smaller 1.39).

| Model | $`\gamma `$ |
| --- | --- |
| $`A\mathrm{}\mathrm{}A`$ | 1.39 |
| $`A+B\mathrm{}`$ | 1.39 |
| $`q`$-state Potts Model | $`[4/3,3/2]`$, $`q[2,\mathrm{}]`$ |
| exchange Model | 1.45 |

TABLE VI. Numerical results for the exponent $`\gamma `$ in various 1D models for which $`U(t)[P(t)]^\gamma `$ (results for the first two models from ref. ).

Finally, we note that all values of $`\gamma `$ obtained so far lie in the range $`4/3q3/2`$ spanned by the $`q`$-state Potts model with $`2q\mathrm{}`$. Up to now our analysis of $`\gamma `$ has been purely numerical, and naturally a more fundamental understanding of this exponent is required.

## V ACKNOWLEDGMENT

We thank Satya Majumdar for drawing our attention to some related work. This work was supported by EPSRC (UK).